+ rotx[0][0]= 1; rotx[0][1]= 0; rotx[0][2]= 0;
+ rotx[1][0]= 0; rotx[1][1]= cy; rotx[1][2]= sy;
+ rotx[2][0]= 0; rotx[2][1]= -sy; rotx[2][2]= cy;
+ PMATRIX(rotx);
+
+ adjr[0][0]= cr; adjr[0][1]= sr; adjr[0][2]= 0;
+ adjr[1][0]= -sr; adjr[1][1]= cr; adjr[1][2]= 0;
+ adjr[2][0]= 0; adjr[2][1]= 0; adjr[2][2]= 1;
+ PMATRIX(adjr);
+
+ GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
+ &rotx_gsl,&adjr_gsl,
+ 0.0, &temp_gsl) );
+ PMATRIX(temp);
+
+ GA( gsl_blas_dgemm(CblasTrans,CblasNoTrans, 1.0,
+ &adjr_gsl,&temp_gsl,
+ 0.0, &change_gsl) );
+ PMATRIX(change);
+
+ static double skew[D3][D3];
+ static GSL_MATRIX(skew);
+
+ GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
+ &change_gsl,&transform_gsl,
+ 0.0, &skew_gsl) );
+ PMATRIX(skew);
+
+ memcpy(&transform,&skew,sizeof(transform));
+ show();
+ return;
+
+ /* Now we want to normalise skew, the result becomes new transform */
+ double svd_v[D3][D3];
+ GSL_MATRIX(svd_v);
+
+ double sigma[D3], tau[D3];
+ GSL_VECTOR(sigma);
+ GSL_VECTOR(tau);
+
+ /* We use notation from Wikipedia Polar_decomposition
+ * Wikipedia's W is GSL's U
+ * Wikipedia's Sigma is GSL's S
+ * Wikipedia's V is GSL's V
+ * Wikipedia's U is our desired result
+ * Wikipedia which says if the SVD is A = W Sigma V*
+ * then the polar decomposition is A = U P
+ * where P = V Sigma V*
+ * and U = W V*
+ */
+
+ GA( gsl_linalg_SV_decomp(&skew_gsl, &svd_v_gsl, &sigma_gsl, &tau_gsl) );
+ pmatrix("W", skew);
+ pvector("Sigma",sigma);
+ pmatrix("V", svd_v);
+
+ /* We only need U, not P. */
+ GA( gsl_blas_dgemm(CblasNoTrans,CblasTrans, 1.0,
+ &skew_gsl,&svd_v_gsl,
+ 0.0,&transform_gsl) );
+
+ pmatrix("U", transform);
+
+ printf("drag_rotate_delta...\n");
+ show();
+}
+DRAG_SAVING(rotate, transform);
+
+static void drag_sizeadj_delta(double dx, double dy) {
+ sizeadj_scale *= pow(3.0, -dy);
+ show();
+}
+DRAG_SAVING(sizeadj, sizeadj_scale);
+
+static void drag_3d_delta(double dx, double dy) {
+ const double min_eyes_apart= -0.02;
+ eyes_apart += dx * 0.1;
+ if (eyes_apart < min_eyes_apart) eyes_apart= min_eyes_apart;
+ printf("sizeadj eyes_apart %g\n", eyes_apart);
+ show();
+}
+DRAG_SAVING(3d, eyes_apart);
+
+static const Drag *drag= &drag_none;
+
+static int drag_last_x, drag_last_y;
+
+static void drag_position(int x, int y) {
+ drag->delta((x - drag_last_x) * 1.0 / wmaxdim,
+ (y - drag_last_y) * 1.0 / wmaxdim);
+ drag_last_x= x;
+ drag_last_y= y;
+}
+
+static void event_button(XButtonEvent *e) {
+ if (e->window != window || !e->same_screen) return;
+ if (e->type == ButtonPress) {
+ if (e->state || drag != &drag_none) {
+ printf("drag=%s press state=0x%lx abandon\n",
+ drag->name, (unsigned long)e->state);
+ drag->abandon();
+ drag= &drag_none;
+ return;
+ }
+ switch (e->button) {
+ case Button1: drag= &drag_rotate; break;
+ case Button2: drag= &drag_sizeadj; break;
+ case Button3: drag= &drag_3d; break;
+ default: printf("unknown drag start %d\n", e->button);
+ }
+ printf("drag=%s press button=%lu start %d,%d\n",
+ drag->name, (unsigned long)e->button, e->x, e->y);
+ drag_last_x= e->x;
+ drag_last_y= e->y;
+ drag->start();
+ }
+ if (e->type == ButtonRelease) {
+ printf("drag=%s release %d,%d\n", drag->name, e->x, e->y);
+ drag_position(e->x, e->y);
+ drag->conclude();
+ drag= &drag_none;
+ }
+}
+
+static void event_motion(int x, int y) {
+ printf("drag=%s motion %d,%d\n", drag->name, x, y);
+ drag_position(x,y);
+}
+
+static void event_config(XConfigureEvent *e) {
+ if (e->width == wwidth && e->height == wheight)
+ return;
+
+ wwidth= e->width; wheight= e->height;
+ wmaxdim= wwidth > wheight ? wwidth : wheight;
+ wmindim= wwidth < wheight ? wwidth : wheight;
+
+ XA( XSetWindowBackground(display,window,BlackPixel(display,x11screen)) );
+ for (currentbuffer=0; currentbuffer<2; currentbuffer++)
+ XA( XFreePixmap(display, doublebuffers[currentbuffer]) );
+
+ mkpixmaps();
+ show();
+}
+