- 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);