From bb0e48bffe4373fabc17a954bbae651bf41def2d Mon Sep 17 00:00:00 2001
From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Mon, 31 Dec 2007 00:57:54 +0000
Subject: [PATCH] can rotate but does not normalise

---
 project.c | 92 +++++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 69 insertions(+), 23 deletions(-)

diff --git a/project.c b/project.c
index 764a4dd..07d5fa7 100644
--- a/project.c
+++ b/project.c
@@ -232,30 +232,90 @@ static void pvectorcore(const char *n, double v[D3]) {
   printf("]\n");
 }
 static void pvector(const char *n, double v[D3]) {
-  pvectorcore(m,v);
+  pvectorcore(n,v);
   putchar('\n');
 }
 static void pmatrix(const char *n, double m[D3][D3]) {
   int j;
-  for (j=0; j<D3; j++) { pvector(n,m[j]); n=""; }
+  for (j=0; j<D3; j++) { pvectorcore(n,m[j]); n=""; }
   putchar('\n');
 }
 #define PMATRIX(x) pmatrix(#x,x);
 
 static void drag_rotate_delta(double dx, double dy) {
-  /* We multiple our transformation matrix by this matrix:
-   *     [  1   0   0 ]
-   *     [  0   1   0 ]
-   *     [ dx  dy   1 ]
-   * and then renormalise.
+  /* We multiple our transformation matrix by a matrix:
+   *
+   * If we just had y movement, we would rotate about x axis:
+   *  rotation X = [  1    0   0 ]
+   *               [  0   cy  sy ]
+   *               [  0  -sy  cy ]
+   *  where cy and sy are sin and cos of y rotation
+   *
+   * But we should pre-rotate this by a rotation about the z axis
+   * to get it to the right angle (to include x rotation).  So
+   * we make cy and sy be cos() and sin(hypot(x,y)) and use
+   * with cr,sr as cos() and sin(atan2(y,y)):
+   *
+   * Ie we would do  T' = R^T X R T   where
+   *             or  T' =    C    T   where  C = R^T X R  and
+   *
+   *  adjustment R = [  cr  sr  0 ]
+   *                 [ -sr  cr  0 ]
+   *                 [  0    0  1 ]
    */
 
-  static double rotateby[D3][D3]= {{1,0,0},{0,1,0},{-20,-20,1}};
-  static GSL_MATRIX(rotateby);
+  double rotx[D3][D3], adjr[D3][D3];
+  GSL_MATRIX(rotx);
+  GSL_MATRIX(adjr);
+
+  static double temp[D3][D3], change[D3][D3];
+  static GSL_MATRIX(temp);
+  static GSL_MATRIX(change);
+
+  double d= hypot(dx,dy);
+  if (d < 1e-6) return;
+
+  double ang= d*2.0;
+  
+  double cy= cos(ang);
+  double sy= sin(ang);
+  double cr= -dy / d;
+  double sr=  dx / d;
+  printf("\n d=%g cy,sy=%g,%g cr,sr=%g,%g\n\n", d,cy,sy,cr,sr);
+  
+  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);
 
@@ -263,18 +323,6 @@ static void drag_rotate_delta(double dx, double dy) {
   GSL_VECTOR(sigma);
   GSL_VECTOR(tau);
   
-  rotateby[2][0]= dx;
-  rotateby[2][1]= dy;
-  PMATRIX(rotateby);
-  
-  GA( gsl_blas_dgemm(CblasNoTrans,CblasNoTrans, 1.0,
-		     &rotateby_gsl,&transform_gsl,
-		     0.0, &skew_gsl) );
-
-  /* Now we want to normalise skew, the result becomes new transform */
-
-  pmatrix("skew", skew);
-
   /* We use notation from Wikipedia Polar_decomposition
    *     Wikipedia's  W      is GSL's U
    *     Wikipedia's  Sigma  is GSL's S
@@ -405,11 +453,9 @@ int main(int argc, const char *const *argv) {
     case ConfigureNotify:   event_config(&event.xconfigure);  break;
 
     case MotionNotify:
-#if 0
       motion_x= event.xmotion.x;
       motion_y= event.xmotion.y;
       motion_deferred= 1;
-#endif
       continue;
       
     default:
-- 
2.30.2