/* +++++++++++ */ # define FILENAME "nrvect.c" /* +++++++++++ */ # include "util.h" # define DOT_PROD(u,v) (u[1]*v[1] + u[2]*v[2] + u[3]*v[3]) # define LENGTH(u) (sqrt(DOT_PROD(u,u))) /* ========== */ void FreeVector(double v[3]) /* ========== */ { free_dvector(v,1,3); } /* U = ai + bj + ck V = Ai + Bj + Ck U x V = (bC - cB) i + (cA - aC) j + (aB - bA) k */ /* Fortran convention for array indicies. WRITE(6,'(/A,3(/3F12.6))') & ' Rotation matrix =',((ROT(I,J),J=1,3),I=1,3) ROT(row,col) so for all nr routines I will use ROT[row][col] */ /* In fortran matrices are stored so that the first index is fast moving: A(1,1) A(2,1) A(3,1) A(1,2) A(2,2) A(3,2) A(1,3) A(2,3) A(3,3) This is opposite the order in C: A[1][1] A[1][2] A[1][3] A[2][1] A[2][2] A[2][3] A[3][1] A[3][2] A[3][3] */ /* =========== */ void dump_matrix (double ** m) /* =========== */ { int r,c; for (r=1; r<4; r++) for (c=1; c<4; c++) printf("%f\t",m[r][c]); } /* ========= */ double * argv2dvect (int *argc, char ***argv,char *s) /* ========= */ { double * v = dvector(1,3) ; if (*argc < 3) ITOD(*argc); v[1] = (double) atof((*argv)[1]) ; v[2] = (double) atof((*argv)[2]) ; v[3] = (double) atof((*argv)[3]) ; *argc -= 3; *argv += 3; printf(__FILE__" [%3d]: %s = ( %f %f %f )\n",__LINE__,s,v[1],v[2],v[3]); return v ; } /* ====================== */ double ** rotate_v_around_z_to_x (double * u) /* ====================== */ { double theta = atan2(u[2],u[1]); double * v = dvector(1,3); v[1]=0; v[2]=0; v[3]=1; return nrmkrotmat(v, theta); } /* ===================== */ double * direction_from_2vects (double * nterm,double * cterm) /* ===================== */ { double * nnterm = nrscale(-1,nterm); double * dir = nrsum(nnterm,cterm); double * norm = normalize(dir); return norm ; } /* ========= */ double * normalize (double * u) /* ========= */ { double dot = u[1]*u[1] + u[2]*u[2] + u[3]*u[3]; double * norm = nrscale(1/sqrt(dot),u); return norm ; } /* ====== */ double * nrcopy (double * u) /* ====== */ { double * v = dvector(1,3); v[1] = u[1] ; v[2] = u[2] ; v[3] = u[3] ; return(v); } /* ===== */ double * nrsum (double * v, double * u) /* ===== */ { double * w = dvector(1,3); w[1] = u[1] + v[1] ; w[2] = u[2] + v[2] ; w[3] = u[3] + v[3] ; return(w); } /* ======== */ double * midpoint (double * v, double * u) /* ======== */ { double * w = dvector(1,3); w[1] = (u[1] + v[1])/2 ; w[2] = (u[2] + v[2])/2 ; w[3] = (u[3] + v[3])/2 ; return(w); } /* ======= */ double * nrscale (double a, double * u) /* ======= */ { double * v = dvector(1,3); v[1] = a * u[1] ; v[2] = a * u[2] ; v[3] = a * u[3] ; return(v); } /* ===== */ double * nrrot (double ** R, double * u) /* ===== */ { double * v = dvector(1,3); 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 ** matmult33 (double ** R, double ** S) /* ========= */ { double ** T = dmatrix(1,3,1,3); int i,j,k; for (i=1;i<=3;i++) for (j=1;j<=3;j++) { T[j][i]=0; for (k=1;k<=3;k++) { T[j][i] += R[j][k]*S[k][i] ; } } return(T); } /* ELSE DO I=1,N DO J=1,M MAT3(J,I)=ZERO DO K=1,L MAT3(J,I)=MAT3(J,I)+MAT1(J,K)*MAT2(K,I) END DO END DO END DO */ /* T[1][1] = R[1][1] * S[1][1] + R[1][2] * S[2][1] + R[1][3] * S[3][1] ;*/ /* T[2][1] = R[2][1] * S[1][1] + R[2][2] * S[2][1] + R[2][3] * S[3][1] ;*/ /* ========== */ double ** nrmkrotmat (double * v0, double theta) /* ========== */ { # define AXIS(i) (v0[i]) # define RAD (PI/180.0) # define ZERO 0 # define ONE 1 # define SMALL 0.0000001 # define MAX(a,b) (a > b ? a : b) # define MIN(a,b) (a < b ? a : b) # define ROT(i,j) (R[i][j]) # define RSMALL 0.0000000001 double ** R = dmatrix(1,3,1,3) ; double NN=sqrt(AXIS(1)*AXIS(1)+AXIS(2)*AXIS(2)+AXIS(3)*AXIS(3)) ; double T1,T2,T3 = TODEG(theta); fprintf(stderr,"WARNING! The routine nrmkrotmat() has been entered. \n"); fprintf(stderr,"WARNING! This routine may be defective and has been \n"); fprintf(stderr,"WARNING! superceded as of 9 August 1994 by \n"); fprintf(stderr,"WARNING! NR_mkrotmat() in NRVects.c. \n"); if (NN < RSMALL) { STDERR("rotation vector has zero length"); T1=ZERO ; T2=ZERO ; T3=ZERO ; } else { AXIS(1)=AXIS(1)/NN ; AXIS(2)=AXIS(2)/NN ; AXIS(3)=AXIS(3)/NN ; T1=acos(AXIS(2))/RAD ; if (AXIS(1)*AXIS(1)+AXIS(3)*AXIS(3) < SMALL) { T2=ZERO ; } else { T2= acos(MAX(-ONE, MIN(ONE, AXIS(1)/sqrt(AXIS(1)*AXIS(1)+AXIS(3)*AXIS(3))))) /RAD ; if (AXIS(3) > ZERO) T2=-T2 ; } } { double S1=sin(RAD*T1) , S2=sin(RAD*T2) , S3=sin(RAD*T3) , C1=cos(RAD*T1) , C2=cos(RAD*T2) , C3=cos(RAD*T3) , S1SQ = S1*S1 , C3C =1.0 - C3 ; ROT(1,1) = C3 + S1SQ*C2*C2*C3C ; ROT(1,2) = S1*C1*C2*C3C - S1*S2*S3 ; ROT(1,3) =-S1SQ*C2*S2*C3C - C1*S3 ; ROT(2,1) = S1*C1*C2*C3C + S1*S2*S3 ; ROT(2,2) = C3 + C1*C1*C3C ; ROT(2,3) =-S1*C1*S2*C3C + S1*C2*S3 ; ROT(3,1) =-S1SQ*S2*C2*C3C + C1*S3 ; ROT(3,2) =-S1*C1*S2*C3C - S1*C2*S3 ; ROT(3,3) = C3 + S1SQ*S2*S2*C3C ; } return (R) ; # undef AXIS # undef RAD # undef ZERO # undef ONE # undef SMALL # undef MAX # undef MIN # undef ROT # undef RSMALL } /* =================== */ double ** nr_rotate_to_z_axis (double * N0) /* =================== */ { /* rotang = acos (N0 . z) */ /* rotaxis = N0 cross z = N0[0] x cross z + N[1] y cross z */ /* = -N0[0] y + N[1] x */ double rotang = - acos(N0[3]) ; double * rotaxis = dvector(1,3) ; rotaxis[1] = N0[2]; rotaxis[2] = -N0[1]; rotaxis[3] = 0; return nrmkrotmat(rotaxis, rotang); } /* =========== */ void dump_vector (double * v) /* =========== */ { int i; for (i=1; i<4; i++) printf("%f\t",v[i]); } /* ===================== */ void GraphicsVector3dCross (double * v, double L) /* ===================== */ { printf("%10.5f%10.5f%10.5f%5d\n",v[1]-L,v[2], v[3], 0); printf("%10.5f%10.5f%10.5f%5d\n",v[1]+L,v[2], v[3], 1); printf("%10.5f%10.5f%10.5f%5d\n",v[1] ,v[2]-L,v[3], 0); printf("%10.5f%10.5f%10.5f%5d\n",v[1] ,v[2]+L,v[3], 1); printf("%10.5f%10.5f%10.5f%5d\n",v[1] ,v[2], v[3]-L,0); printf("%10.5f%10.5f%10.5f%5d\n",v[1] ,v[2], v[3]+L,1); } /* ========================= */ void GraphicsVectorKludgeArrow (double * v, double *u) /* ========================= */ { double L = LENGTH(second_from_first(u,v)); printf("%10.5f%10.5f%10.5f%5d\n",v[1],v[2],v[3],0); printf("%10.5f%10.5f%10.5f%5d\n",u[1],u[2],u[3],1); GraphicsVector3dCross(u,L/10); } /* ================ */ void graphics_vectors (double * v, double * u) /* ================ */ { printf("%10.5f%10.5f%10.5f%5d\n",v[1],v[2],v[3],0); printf("%10.5f%10.5f%10.5f%5d\n",u[1],u[2],u[3],1); } /* =========== */ void atom_vector /* =========== */ (double * v, char * resnam, int iresn, double b) { static int iser = 1; printf(PDBFMT,iser++," CA ",resnam,' ', iresn,v[1],v[2],v[3],1.0,b," "); } /* ============== */ double angle_on_plane (double * norm, double * u, double * v) /* ============== */ { double * upar = nrscale(DOT_PROD(norm,u),norm); double * vpar = nrscale(DOT_PROD(norm,v),norm); double * uperp = second_from_first(u,upar); double * vperp = second_from_first(v,vpar); double dot = DOT_PROD(uperp,vperp) / ( LENGTH(uperp)*LENGTH(vperp) ); double angle = TODEG(acos(dot)); /* FTOD(dot); */ return angle ; } /* ============= */ double dist_on_plane (double * norm, double * u, double * v) /* ============= */ { double * w = second_from_first(u,v); double * wpar = nrscale(DOT_PROD(norm,w),norm); double * wperp = second_from_first(w,wpar); double dist = LENGTH(wperp); /* FTOD(dot); */ return dist ; } /* ================= */ double * second_from_first (double * u, double * v) /* ================= */ { return vect(u[1]-v[1], u[2]-v[2], u[3]-v[3]); } /* ==== */ double * vect (double a, double b, double c) /* ==== */ { double * w = dvector(1,3); w[1] = a ; w[2] = b ; w[3] = c ; return w; } /* ================== */ double * ProjectPtToSegment (double * Pt, double * Seg1,double * Seg2) /* ================== */ { STDERR("ProjectPtToSegment BEG"); { double * Seg12 = second_from_first(Seg2,Seg1); double LenSeg12 = LENGTH(Seg12); double * dirSeg12 = nrscale(1/LenSeg12,Seg12); double * dir_check = direction_from_2vects(Seg1,Seg2); double * Seg1Pt = second_from_first(Pt,Seg1); /* double LenSeg1Pt = LENGTH(Seg1Pt); double Proj = DOT_PROD(Seg12,Seg1Pt)/(LenSeg1Pt * LenSeg12) ; double * ProjSeg12 = nrscale(Proj,Seg12); */ double Proj = DOT_PROD(dirSeg12,Seg1Pt); double ratio = Proj / LenSeg12 ; double * ProjSeg12 = nrscale(ratio,Seg12); double * ProjPt = nrsum(Seg1,ProjSeg12); FTOE(Proj); NRTOE(dirSeg12); NRTOE(dir_check); # ifdef OFF assert(nreq(dirSeg12,dir_check)); # endif STDERR("ProjectPtToSegment END"); return ProjPt ; } } /* ==== */ bool nreq (double * v, double * w) /* ==== */ { # define TINY 0.000001 return ( (v[1]-w[1]