static Vertices conformation;
static double transform[D3][D3]= {{1,0,0}, {0,1,0}, {0,0,1}};
-STATIC_GSL_MATRIX(transform);
+static GSL_MATRIX(transform);
const char *input_filename;
}
static void transform_coordinates(void) {
- /*
- static double result[D3];
- static const gsl_vector result_gsl= { D3,1,&result[0]; }
+ double result[D3];
+ GSL_VECTOR(result);
+ gsl_vector input_gsl= { D3,1 };
- int v;
+ int v, k;
-
FOR_VERTEX(v) {
- GA( gsl_blas_dgemv(CblasNoTrans,
- }*/
+ input_gsl.data= &conformation[v][0];
+ GA( gsl_blas_dgemv(CblasNoTrans, 1.0,&transform_gsl, &input_gsl,
+ 0.0, &result_gsl) );
+ K conformation[v][k]= result[k];
+ }
}
static void addtriangle(int va, int vb, int vc) {
static void drag_none_abandon(void) { }
DRAG(none);
+static void pvectorcore(const char *n, double v[D3]) {
+ int k;
+ printf("%10s [ ",n);
+ K printf("%# 10.10f ",v[k]);
+ printf("]\n");
+}
+static void pvector(const char *n, double v[D3]) {
+ pvectorcore(n,v);
+ putchar('\n');
+}
+static void pmatrix(const char *n, double m[D3][D3]) {
+ int j;
+ 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);
- static double qr[D3][D3];
- STATIC_GSL_MATRIX(qr);
+ double d= hypot(dx,dy);
+ if (d < 1e-6) return;
- static double tau[D3];
- STATIC_GSL_VECTOR(tau);
+ double ang= d*2.0;
- int k;
+ 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);
- K rotateby[k][k]= 1;
- rotateby[0][2]= dx;
- rotateby[1][2]= dy;
+ 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,
- &rotateby_gsl,&transform_gsl, 0.0,&qr_gsl) );
- GA( gsl_linalg_QR_decomp(&qr_gsl, &tau_gsl) );
- GA( gsl_linalg_QR_unpack(&qr_gsl, &tau_gsl,
- &transform_gsl, &rotateby_gsl /*dummy*/) );
+ &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();
~ian / moebius2.git / blobdiff