#include "bbltyp.h" /*----------------------------------------------------- Notice of blatant theft. This code was totally lifted from Roger Sayle's Internal2Cartesian routine in RasMol 2.6. Thanks Roger. Roger's code was a whole lot cleaner than the earlier version which was stolen from MOPAC 5. However, by looking at Roger's varialbe names it looks like his code was derived from MOPAC also. Oh well, what comes around goes around :-) WPW 4-1-96 April Fools !!!!! ------------------------------------------------------*/ int int_to_cart(ums_type *mol) { register double cosph,sinph,costh,sinth,coskh,sinkh; register double cosa,sina,cosd,sind; register double dist,angle,dihed; register double xpd,ypd,zpd,xqd,yqd,zqd; register double xa,ya,za,xb,yb,zb; register double rbc,xyb,yza,temp; register double xpa,ypa,zqa; register double xd,yd,zd; register int flag; int i, na, nb, nc; /* Atom #1 */ X(1) = 0.0; Y(1) = 0.0; Z(1) = 0.0; if (Atoms == 1) return( TRUE ); /* Atom #2 */ X(2) = R(2); Y(2) = 0.0; Z(2) = 0.0; if( Atoms == 2 ) return( TRUE ); /* Atom #3 */ dist = R(3); angle = W(3) * DEG_TO_RAD; cosa = cos(angle); sina = sin(angle); if( NA(3) == 1 ) { X(3) = X(1) + cosa*dist; } else { X(3) = X(2) - cosa*dist; } Y(3) = sina*dist; Z(3) = 0.0; for (i = 4; i <= Atoms; i++) { dist = R(i); angle = W(i) * DEG_TO_RAD; dihed = T(i) * DEG_TO_RAD; na = NA(i); nb = NB(i); nc = NC(i); xb = X(nb) - X(na); yb = Y(nb) - Y(na); zb = Z(nb) - Z(na); rbc = xb*xb + yb*yb + zb*zb; if( rbc < 0.0001 ) return( FALSE ); rbc = 1.0/sqrt(rbc); cosa = cos(angle); sina = sin(angle); if( fabs(cosa) >= 0.999999 ) { /* Colinear */ temp = dist*rbc*cosa; X(i) = X(na) + temp*xb; Y(i) = Y(na) + temp*yb; Z(i) = Z(na) + temp*zb; } else { xa = X(nc) - X(na); ya = Y(nc) - Y(na); za = Z(nc) - Z(na); sind = -sin(dihed); cosd = cos(dihed); xd = dist*cosa; yd = dist*sina*cosd; zd = dist*sina*sind; xyb = sqrt(xb*xb + yb*yb); if( xyb < 0.1 ) { /* Rotate about y-axis! */ temp = za; za = -xa; xa = temp; temp = zb; zb = -xb; xb = temp; xyb = sqrt(xb*xb + yb*yb); flag = TRUE; } else flag = FALSE; costh = xb/xyb; sinth = yb/xyb; xpa = costh*xa + sinth*ya; ypa = costh*ya - sinth*xa; sinph = zb*rbc; cosph = sqrt(1.0 - sinph*sinph); zqa = cosph*za - sinph*xpa; yza = sqrt(ypa*ypa + zqa*zqa); if( yza > 1.0E-10 ) { coskh = ypa/yza; sinkh = zqa/yza; ypd = coskh*yd - sinkh*zd; zpd = coskh*zd + sinkh*yd; } else { /* coskh = 1.0; */ /* sinkh = 0.0; */ ypd = yd; zpd = zd; } xpd = cosph*xd - sinph*zpd; zqd = cosph*zpd + sinph*xd; xqd = costh*xpd - sinth*ypd; yqd = costh*ypd + sinth*xpd; if( flag ) { /* Rotate about y-axis! */ X(i) = X(na) - zqd; Y(i) = Y(na) + yqd; Z(i) = Z(na) + xqd; } else { X(i) = X(na) + xqd; Y(i) = Y(na) + yqd; Z(i) = Z(na) + zqd; } } } return( TRUE ); }