#include "util.h" #include "NRVects.h" NRMatrix33P Rot_ZtoU_YtoZxV(NRVectP u1,NRVectP u2, NRVectP v1,NRVectP v2) /* U = u2-u1 V = v2-v1 */ { NRMatrix33P M1 = dmatrix(1,3,1,3); AllocateNRVect3 yhat, xhat ; NRVectP zhat, vhat; AllocateNRVect3 o = {0,0,0,0}; if (u2 == NULL) { zhat = direction_from_2vects(o,u1); } else { /* (u2-u1)/|u2-u1| */ zhat = direction_from_2vects(u1,u2); } if (v2 == NULL) { vhat = direction_from_2vects(o,v1); } else { vhat = direction_from_2vects(v1,v2); } NR_perp_dir_1eq2x3(yhat,zhat,vhat); NR_perp_dir_1eq2x3(xhat,yhat,zhat); NR_MakeMat1WithVect234(M1,xhat,yhat,zhat); Transpose(M1); free_dvector(zhat,1,3); free_dvector(vhat,1,3); return M1 ; } void NR_SimpPrintMat33 (NRMatrix33P M) { NRVect_Print(M[1],3); NRVect_Print(M[2],3); NRVect_Print(M[3],3); } double ** NR_mkrotmat(NRMatrix33P M, NRVectP v, double t) /* From Altman (1986), "Rotations, Quaternions, & Double Groups" : 1 - 2(mm+nn)S -ns + 2lmS ms + 2nlS ns + 2lmS 1 - 2(nn+ll)S -ls + 2mnS -ms + 2nlS ls + 2mnS 1 - 2(ll+mm)S This gives a matrix with the opposite sense of rotation to the XPLOR convention. Consequently, I change the sign of the rotation angle t in the routine. */ { double T = TORAD(-t), l = v[1], m = v[2], n = v[3], s = sin(T), S = sin(T/2)*sin(T/2); M[1][1]= 1-2*(m*m+n*n)*S ;M[1][2]= -n*s+2*l*m*S ;M[1][3]= m*s+2*n*l*S ; M[2][1]= n*s+2*l*m*S ;M[2][2]= 1-2*(n*n+l*l)*S ;M[2][3]= -l*s+2*m*n*S ; M[3][1]= -m*s+2*n*l*S ;M[3][2]= l*s+2*m*n*S ;M[3][3]= 1-2*(l*l+m*m)*S ; return (M); } double * VectRotate_1is2times3 (double * v, double ** R, double * u) { v[1] = R[1][1] * u[1] + R[1][2] * u[2] + R[1][3] * u[3] ; v[2] = R[2][1] * u[1] + R[2][2] * u[2] + R[2][3] * u[3] ; v[3] = R[3][1] * u[1] + R[3][2] * u[2] + R[3][3] * u[3] ; return(v); } double * VectRotate_1by2 (double * u, double ** R) { AllocateNRVect3 v; v[1] = R[1][1] * u[1] + R[1][2] * u[2] + R[1][3] * u[3] ; v[2] = R[2][1] * u[1] + R[2][2] * u[2] + R[2][3] * u[3] ; v[3] = R[3][1] * u[1] + R[3][2] * u[2] + R[3][3] * u[3] ; u[1]=v[1]; u[2]=v[2]; u[3]=v[3]; return(u); } double * VectTranslate_1by2 (double *v, double * u) { v[1] += u[1]; v[2] += u[2]; v[3] += u[3]; return(v); } double * VectSub_1is2minus3 (double *w, double * u, double * v) { w[1] = u[1] - v[1] ; w[2] = u[2] - v[2] ; w[3] = u[3] - v[3] ; return(w); } NRMatrix33P IdentityMatrix(NRMatrix33P R) { R[1][1]= 1 ;R[1][2]= 0 ;R[1][3]= 0 ; R[2][1]= 0 ;R[2][2]= 1 ;R[2][3]= 0 ; R[3][1]= 0 ;R[3][2]= 0 ;R[3][3]= 1 ; return R ; } double ** Transpose(double ** M) { double save = M[1][2] ; M[1][2] = M[2][1]; M[2][1] = save ; save = M[1][3] ; M[1][3] = M[3][1]; M[3][1] = save ; save = M[2][3] ; M[2][3] = M[3][2]; M[3][2] = save ; return M ; } double NR_dotprod(NRVectP u, NRVectP v) { return (u[1]*v[1] + u[2]*v[2] + u[3]*v[3]); } double NR_length(NRVectP v) { return sqrt(NR_dotprod(v,v)); } double NRVect_Mean(NRVectP v, int n) { int i; double s = 0; for (i=1; i<=n ; i++) { s += v[i]; } return s / n ; } double NRVect_SumSq(NRVectP v, int n) { int i; double s = 0; for (i=1; i<=n ; i++) { s += v[i]*v[i]; } return s ; } /* Useful for Chris */ double NRVect_Covar(NRVectP v, NRVectP u, int n) { double v_mean = NRVect_Mean(v,n); double u_mean = NRVect_Mean(u,n); int i; double s = 0; for (i=1; i<=n ; i++) { double dv = v[i] - v_mean ; double du = u[i] - u_mean ; s += dv*du; } return s / n ; } void NRVect_Print(NRVectP v, int n) { int i; for (i=1; i<=n ; i++) printf("%12.6f ",v[i]); printf("\n"); } /* For all nr routines I will use ROT[row][col] */ NRMatrix33P NR_MakeMat1WithVect234(NRMatrix33P M, NRVectP v1, NRVectP v2, NRVectP v3) { M[1][1] = v1[1] ; M[2][1] = v1[2] ; M[3][1] = v1[3] ; M[1][2] = v2[1] ; M[2][2] = v2[2] ; M[3][2] = v2[3] ; M[1][3] = v3[1] ; M[2][3] = v3[2] ; M[3][3] = v3[3] ; return M ; } NRVectP NR_Normalize (NRVectP v) { double len = NR_length(v); v[1] = v[1] / len ; v[2] = v[2] / len ; v[3] = v[3] / len ; } NRVectP CopyNRVect3_1to2 (NRVectP from, NRVectP to) { to[1]=from[1];to[2]=from[2];to[3]=from[3]; return(to); } NRVectP NR_cross_1eq2x3(NRVectP prod, NRVectP v1, NRVectP v2) { prod[1] = v1[2]*v2[3] - v2[2]*v1[3] ; prod[2] = v1[3]*v2[1] - v2[3]*v1[1] ; prod[3] = v1[1]*v2[2] - v2[1]*v1[2] ; return prod ; } NRVectP NR_perp_dir_1eq2x3(NRVectP prod,NRVectP v1, NRVectP v2) { NR_cross_1eq2x3(prod,v1,v2); NR_Normalize(prod); return prod ; } NRVectP NR_GetXYZVector(NRVectP M, FILE * stream, char * string, char * TheDefault) { # define WHOLE_ARRAY(M) \ M[1],M[2],M[3] char * fmt = "%lg%lg%lg" ; char * from ; if (stream != NULL) { from = "a stream" ; if (fscanf(stream,fmt,WHOLE_ARRAY(&M)) != 3) return NR_GetXYZVector(M,NULL,string,TheDefault); } else if (string != NULL) { from = "a string" ; if (sscanf(string,fmt,WHOLE_ARRAY(&M)) != 3) return NR_GetXYZVector(M,stream,NULL,TheDefault); } else { from = "the default" ; sscanf(TheDefault,fmt,WHOLE_ARRAY(&M)); } fmt = "[Got 3x1 from %12s] " "( %7.4f %7.4f %7.4f )\n" ; fprintf(stderr,fmt,from,WHOLE_ARRAY(M)); return (double *) M ; # undef WHOLE_ARRAY } NRMatrix33P NR_CalcI33(NRMatrix33P C, NRVectP x, NRVectP y, NRVectP z, int n) { C[1][1] = n * (NRVect_Covar(y,y,n) + NRVect_Covar(z,z,n)); C[2][2] = n * (NRVect_Covar(x,x,n) + NRVect_Covar(z,z,n)); C[3][3] = n * (NRVect_Covar(y,y,n) + NRVect_Covar(x,x,n)); C[1][2] = - n * NRVect_Covar(x,y,n); C[1][3] = - n * NRVect_Covar(x,z,n); C[2][3] = - n * NRVect_Covar(y,z,n); C[2][1] = C[1][2] ; C[3][1] = C[1][3] ; C[3][2] = C[2][3] ; return C ; } NRMatrix33P NR_CalcCovarMat33(NRMatrix33P C, NRVectP x, NRVectP y, NRVectP z, int n) { C[1][1] = NRVect_Covar(x,x,n); C[1][2] = NRVect_Covar(x,y,n); C[1][3] = NRVect_Covar(x,z,n); C[2][2] = NRVect_Covar(y,y,n); C[2][3] = NRVect_Covar(y,z,n); C[3][3] = NRVect_Covar(z,z,n); C[2][1] = C[1][2] ; C[3][1] = C[1][3] ; C[3][2] = C[2][3] ; return C ; } void RotateXYZby1(NRMatrix33P R, NRVectP x, NRVectP y, NRVectP z, int n) { AllocateNRVect3 v; int i; for (i=1; i<=n; i++) { v[1] = x[i]; v[2] = y[i]; v[3] = z[i]; VectRotate_1by2(v,R); x[i] = v[1]; y[i] = v[2]; z[i] = v[3]; } } void TranslateXYZby1(NRVectP v, NRVectP x, NRVectP y, NRVectP z, int n) { int i; for (i=1; i<=n; i++) { x[i] += v[1]; y[i] += v[2]; z[i] += v[3]; } } void jacpark(float **a, int n, float d[], float **v, int *nrot); void DiagMatrix33 (NRMatrix33P Cd, NRVectP Eigenvalues, NRMatrix33P Eigenvectors) { float ** C = matrix(1,3,1,3) ; float * val = vector(1,3) ; float ** vec = matrix(1,3,1,3) ; int i,j; int nrot ; C[1][1] = Cd[1][1] ; C[1][2] = Cd[1][2] ; C[1][3] = Cd[1][3] ; C[2][2] = Cd[2][2] ; C[2][3] = Cd[2][3] ; C[3][3] = Cd[3][3] ; C[2][1] = Cd[2][1] ; C[3][1] = Cd[3][1] ; C[3][2] = Cd[3][2] ; jacpark(C,3,val,vec,&nrot); for (i=1; i<4; i++) { Eigenvalues[i] = val[i] ; for (j=1; j<4; j++) { Eigenvectors[i][j] = vec[i][j]; } } free_matrix(C,1,3,1,3); free_matrix(vec,1,3,1,3); free_vector(val,1,3); } void CyclicPermuteVecMat (NRVectP EigenVals, NRMatrix33P EigenVecs) { int i; NRMatrix33P M = EigenVecs ; AllocateNRVect3 T ; T[1]=EigenVals[3];T[2]=EigenVals[1];T[3]=EigenVals[2]; CopyNRVect3_1to2(T,EigenVals); for(i=1;i<=3;i++){ T[1]=M[3][i]; T[2]=M[1][i]; T[3]=M[2][i]; M[1][i]=T[1]; M[2][i]=T[2]; M[3][i]=T[3]; } } void PermuteEigenValsVecsSoLargestLast (NRVectP EigenVals, NRMatrix33P EigenVecs) { if (EigenVals[3] >= EigenVals[2] && EigenVals[3] >= EigenVals[1]) { return ; } else if (EigenVals[2] >= EigenVals[3] && EigenVals[2] >= EigenVals[1]) { CyclicPermuteVecMat(EigenVals,EigenVecs); return ; } else if (EigenVals[1] >= EigenVals[3] && EigenVals[1] >= EigenVals[2]) { CyclicPermuteVecMat(EigenVals,EigenVecs); CyclicPermuteVecMat(EigenVals,EigenVecs); return ; } } int LargestElementIndex (NRVectP v, int N) { int maxi = 1; int i; for (i=2; i<=N; i++) { if (v[i] >= v[maxi]) { maxi = i; } } return maxi; }