home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / de118i.zip / OBLATE.C < prev    next >
C/C++ Source or Header  |  1993-02-06  |  16KB  |  562 lines
  1. /* oblate.c
  2.  */
  3.  
  4. #ifndef NOINCS
  5. #include "mconf.h"
  6. #include "prec.h"
  7. #include "ssystem.h"
  8. #include "const.h"
  9. #endif
  10.  
  11. #ifdef DEBUG
  12. #undef DEBUG
  13. #endif
  14. #define DEBUG 0
  15.  
  16. extern DOUBLE GMs[];
  17. extern DOUBLE Rij[NTOTAL][NTOTAL];
  18. extern DOUBLE Rij3[NTOTAL][NTOTAL];
  19.  
  20. extern DOUBLE AM, AE, LOVENO, PHASE, EMRAT;
  21. extern DOUBLE PSLP1, PSLPI, PSLPIA, PSLPIB, JDEPOCH;
  22. extern DOUBLE LBET, LGAM, LGAMBET;
  23. extern DOUBLE Jm[], Je[], Cnm[], Snm[];
  24.  
  25. DOUBLE ofdate[3];
  26. DOUBLE phi, th, psi, phi1, th1, psi1, Nx, Ny, Nz;
  27. DOUBLE sinpsi, cospsi, sinth, costh, sinphi, cosphi;
  28. DOUBLE X00, X01, X02, X10, X11, X12, X20, X21, X22;
  29. extern DOUBLE B1950;
  30. DOUBLE SQRT(), SIN(), COS(), FLOOR();
  31.  
  32. #if LDOUBLE
  33. #if IBMPC
  34. short Cr375i[]  = {0x0000,0x0000,0x0000,0xc000,0x3ffd, XPD};
  35. short C2r5i[]   = {0x0000,0x0000,0x0000,0xa000,0x4000, XPD};
  36. short C3r75i[]  = {0x0000,0x0000,0x0000,0xf000,0x4000, XPD};
  37. short C4r375i[] = {0x0000,0x0000,0x0000,0x8c00,0x4001, XPD};
  38. short C7r5i[]   = {0x0000,0x0000,0x0000,0xf000,0x4001, XPD};
  39. short C15i[]    = {0x0000,0x0000,0x0000,0xf000,0x4002, XPD};
  40. short C17r5i[]  = {0x0000,0x0000,0x0000,0x8c00,0x4003, XPD};
  41. short C52r5i[]  = {0x0000,0x0000,0x0000,0xd200,0x4004, XPD};
  42. short C105i[]   = {0x0000,0x0000,0x0000,0xd200,0x4005, XPD};
  43. #endif
  44. #if MIEEE
  45. long Cr375i[]   = {0x3ffd0000,0xc0000000,0x00000000};
  46. long C2r5i[]    = {0x40000000,0xa0000000,0x00000000};
  47. long C3r75i[]   = {0x40000000,0xf0000000,0x00000000};
  48. long C4r375i[]  = {0x40010000,0x8c000000,0x00000000};
  49. long C7r5i[]    = {0x40010000,0xf0000000,0x00000000};
  50. long C15i[]     = {0x40020000,0xf0000000,0x00000000};
  51. long C17r5i[]   = {0x40030000,0x8c000000,0x00000000};
  52. long C52r5i[]   = {0x40040000,0xd2000000,0x00000000};
  53. long C105i[]    = {0x40050000,0xd2000000,0x00000000};
  54. #endif
  55. #define Cr375 (*(long double *)Cr375i)
  56. #define C2r5 (*(long double *)C2r5i)
  57. #define C3r75 (*(long double *)C3r75i)
  58. #define C4r375 (*(long double *)C4r375i)
  59. #define C7r5 (*(long double *)C7r5i)
  60. #define C15 (*(long double *)C15i)
  61. #define C17r5 (*(long double *)C17r5i)
  62. #define C52r5 (*(long double *)C52r5i)
  63. #define C105 (*(long double *)C105i)
  64. #else /* LDOUBLE */
  65. DOUBLE Cr375 =  3.75000000000000000000E-1;
  66. DOUBLE C2r5 =  2.50000000000000000000E0;
  67. DOUBLE C3r75 =  3.75000000000000000000E0;
  68. DOUBLE C4r375 = 4.37500000000000000000E0;
  69. DOUBLE C7r5 = 7.50000000000000000000E0;
  70. DOUBLE C15 = 1.50000000000000000000E1;
  71. DOUBLE C17r5 = 1.75000000000000000000E1;
  72. DOUBLE C52r5 = 5.25000000000000000000E1;
  73. DOUBLE C105 = 1.05000000000000000000E2;
  74. #endif /* not LDOUBLE */
  75.  
  76.  
  77. oblate( JD, yp, v, ymoon )
  78. DOUBLE JD;
  79. DOUBLE yp[], v[], ymoon[];
  80. {
  81. DOUBLE f, a, b, rem, res;
  82. DOUBLE xyz[3], c, s;
  83. int i;
  84.  
  85. /* x, y, z are the the (equatorial J2000)
  86.  * solar system barycentric Cartesian coordinates
  87.  * of the external point mass relative to the center of
  88.  * the extended body.
  89.  * Bx, By, Bz are the principal axes of inertia of the extended body.
  90.  * The Euler angles are
  91.  * theta, the angle from the z axis to the Bz axis
  92.  * phi, measured in the x-y plane from the x axis to the line of nodes
  93.  * psi, measured in the body's equatorial (Bx-By) plane
  94.  *      from the line of nodes to the Bx axis.
  95.  */
  96. phi1 = yp[0];
  97. phi = yp[1];
  98. th1 = yp[2];
  99. th = yp[3];
  100. psi1 = yp[4];
  101. v[1] = phi1; /* Copy the first derivative vector */
  102. v[3] = th1;
  103. v[5] = psi1;
  104. /* Add slope representing the Moon's mean rotation */
  105. psi1 += PSLP1;
  106. /* Calculate psi = yp[5] + PSLP1*(JD - JDEPOCH).
  107.  * Since only the sine and cosine of psi are needed, psi may
  108.  * here be reduced accurately modulo the mean rotation period
  109.  * of 2 pi / PSLP1 days.
  110.  */
  111. b = JD - JDEPOCH; /* elapsed days assumed arithmetically exact */
  112. /* Subtract nearest integer number of rotation periods from the time */
  113. f = FLOOR( b/PSLPI ); /* integer number of periods elapsed */
  114. a = b - f * PSLPIA; /* arithmetic is exact while  f < 2^26  periods */
  115. b = a - f * PSLPIB; /* note, PSLPIA + PSLPIB = 1 period */
  116. psi = yp[5] + PSLP1*b; /* yp[5] is the integrator output for psi */
  117. #if DEBUG
  118. printf( "phi %.5e, th %.5e, psi %.5e\n", (double) phi, (double) th, (double) psi );
  119. #endif
  120.  
  121.  
  122. /* Construct matrix elements that will be used by mfigure()
  123.  * to convert the extended-body-to-point-mass vector
  124.  * from space (xyz) coordinates to body (Bxyz) coordinates.
  125.  * This is a composite of three rotations:
  126.  *   first by an angle phi about the z axis,
  127.  *   second by an angle theta about the new x axis,
  128.  *   third by an angle psi about the final z axis.
  129.  * Each rotation is counter-clockwise, looking toward the
  130.  * origin along the indicated axis.  See H. Goldstein,
  131.  * _Classsical Mechanics_, Addison-Wesley, 1950, pp 107-109.
  132.  */
  133. sinpsi = SIN(psi);
  134. cospsi = COS(psi);
  135. sinth = SIN(th);
  136. costh = COS(th);
  137. sinphi = SIN(phi);
  138. cosphi = COS(phi);
  139. a = costh*sinphi;
  140. b = costh*cosphi;
  141. X00 = cospsi*cosphi - a*sinpsi;
  142. X01 = cospsi*sinphi + b*sinpsi;
  143. X02 = sinpsi*sinth;
  144. X10 = -sinpsi*cosphi - a*cospsi;
  145. X11 = -sinpsi*sinphi + b*cospsi;
  146. X12 = cospsi*sinth;
  147. X20 = sinth*sinphi;
  148. X21 = -sinth*cosphi;
  149. X22 = costh;
  150.  
  151. /* Initialize the torque accumulators.
  152.  */
  153. Nx = Zero;
  154. Ny = Zero;
  155. Nz = Zero;
  156.  
  157. /* Accelerations due to the figure of the Moon.
  158.  */
  159. /* effect of the point-mass Earth : 2e-12 au/d^2 */
  160. xyz[0] = -ymoon[1];
  161. xyz[1] = -ymoon[3];
  162. xyz[2] = -ymoon[5];
  163. mfigure( IEARTH, xyz, v );
  164.  
  165. /* effect of the point-mass Sun: 1e-17 au/d^2 */
  166. xyz[0] = yp[6*ISUN+1] - yp[6*IMOON+1];
  167. xyz[1] = yp[6*ISUN+3] - yp[6*IMOON+3];
  168. xyz[2] = yp[6*ISUN+5] - yp[6*IMOON+5];
  169. #if DEBUG
  170. printf( "Moon->Sun %.5e %.5e %.5e\n", (double) xyz[0], (double) xyz[1], (double) xyz[2] );
  171. #endif
  172. /* a = Rij[ISUN][IMOON]; */
  173. mfigure( ISUN, xyz, v );
  174.  
  175. /* Update accelerations of libration angles */
  176. librate( v );
  177.  
  178.  
  179. /* Oblateness of the earth
  180.  */
  181.  
  182. /* effect of the point-mass Moon : 7e-11 au/d^2
  183.  */
  184. for( i=0; i<3; i++ )
  185.     ofdate[i] = ymoon[2*i+1];
  186. /* Precess from basic epoch to date */
  187. precess( ofdate, JD, -1 );
  188. nutate( JD, ofdate );
  189. rem = Rij[IEARTH][IMOON];
  190. legendre( ofdate, AE, rem, Je, xyz, 0 );
  191. invnut( xyz );
  192. precess( xyz, JD, 1 );
  193. a = -GMs[IMOON];
  194. i = 6 * IEARTH;
  195. v[i] += a * xyz[0];
  196. v[i+2] += a * xyz[1];
  197. v[i+4] += a * xyz[2];
  198.  
  199. a = GMs[IEARTH];
  200. #if DEBUG
  201. printf( "E,M %.5e %.5e %.5e\n",
  202.    (double )(a*xyz[0]), (double )(a*xyz[1]), (double )(a*xyz[2]) );
  203. #endif
  204. i = 6 * IMOON;
  205. v[i] += a * xyz[0];
  206. v[i+2] += a * xyz[1];
  207. v[i+4] += a * xyz[2];
  208.  
  209. /* tides: 1e-15 au/d^2
  210.  */
  211. a = AE * rem;
  212. b = a*a;
  213. a = b*b*a*(One + One/EMRAT);
  214. a = -Three*LOVENO*a*GMs[IMOON] * Rij3[IEARTH][IMOON];
  215. xyz[0] = a * (ofdate[0] + PHASE * ofdate[1]);
  216. xyz[1] = a * (ofdate[1] - PHASE * ofdate[0]);
  217. xyz[2] = a * ofdate[2];
  218. invnut( xyz );
  219. precess( xyz, JD, 1 );
  220. c = One/(One+EMRAT);
  221. s = EMRAT*c;
  222. i = 6 * IEARTH;
  223. v[i] -= c * xyz[0];
  224. v[i+2] -= c * xyz[1];
  225. v[i+4] -= c * xyz[2];
  226. #if DEBUG
  227. printf( "T,M %.5e %.5e %.5e\n",
  228.     (double )(s*xyz[0]), (double )(s*xyz[1]), (double )(s*xyz[2]) );
  229. #endif
  230. i = 6 * IMOON;
  231. v[i] += s * xyz[0];
  232. v[i+2] += s * xyz[1];
  233. v[i+4] += s * xyz[2];
  234.  
  235.  
  236. /* effect of the Sun: 7e-16 au/d^2
  237.  */
  238. res = Rij[ISUN][IEARTH];
  239. ofdate[0] = yp[6*ISUN+1] - yp[6*IEARTH+1];
  240. ofdate[1] = yp[6*ISUN+3] - yp[6*IEARTH+3];
  241. ofdate[2] = yp[6*ISUN+5] - yp[6*IEARTH+5];
  242. precess( ofdate, JD, -1 );
  243. nutate( JD, ofdate );
  244. legendre( ofdate, AE, res, Je, xyz, 0 );
  245. invnut( xyz );
  246. precess( xyz, JD, 1 );
  247. a = -GMs[ISUN];
  248. #if DEBUG
  249. printf( "E,S %.5e %.5e %.5e\n",
  250.  (double )(a*xyz[0]), (double )(a*xyz[1]), (double )(a*xyz[2]) );
  251. #endif
  252. i = 6 * IEARTH;
  253. v[i] += a * xyz[0];
  254. v[i+2] += a * xyz[1];
  255. v[i+4] += a * xyz[2];
  256. a = GMs[IEARTH];
  257. i = 6 * ISUN;
  258. v[i] += a * xyz[0];
  259. v[i+2] += a * xyz[1];
  260. v[i+4] += a * xyz[2];
  261. }
  262.  
  263.  
  264.  
  265. mfigure( iobj, xyz, v )
  266. int iobj; /* Index of distance from center of extended body to point mass */
  267. DOUBLE xyz[]; /* space coordinates of point mass */
  268. DOUBLE v[]; /* Velocity and acceleration state vector */
  269. {
  270. DOUBLE a, rem;
  271. DOUBLE Bxyz[3], Fxyz[3];
  272. int i;
  273.  
  274. /* Apply the transformation from space to body coordinates.
  275.  */
  276. Bxyz[0] = X00*xyz[0] + X01*xyz[1] + X02*xyz[2];
  277. Bxyz[1] = X10*xyz[0] + X11*xyz[1] + X12*xyz[2];
  278. Bxyz[2] = X20*xyz[0] + X21*xyz[1] + X22*xyz[2];
  279.  
  280. rem = Rij[iobj][IMOON];
  281.  
  282. /* Compute acceleration due to Moon's oblateness
  283.  */
  284. legendre( Bxyz, AM, rem, Jm, Fxyz, 1 );
  285.  
  286. /* The torque on the extended body, in body coordinates is r X F.
  287.  * -GMiobj Fxyz is the acceleration of the extended body,
  288.  * so the torque = -GMiobj r X Fxyz M
  289.  * where M is the mass of the extended body.  M is not actually
  290.  * included, since it will cancel out later.
  291.  */
  292. a = -GMs[iobj];
  293. Nx += a*(Bxyz[1] * Fxyz[2] - Bxyz[2] * Fxyz[1]);
  294. Ny += a*(Bxyz[2] * Fxyz[0] - Bxyz[0] * Fxyz[2]);
  295. Nz += a*(Bxyz[0] * Fxyz[1] - Bxyz[1] * Fxyz[0]);
  296.  
  297. /* Transform from body coordinates to space coordinates
  298.  * by the inverse (i.e., the transpose) of the previous transformation.
  299.  */
  300. xyz[0] = X00*Fxyz[0] + X10*Fxyz[1] + X20*Fxyz[2];
  301. xyz[1] = X01*Fxyz[0] + X11*Fxyz[1] + X21*Fxyz[2];
  302. xyz[2] = X02*Fxyz[0] + X12*Fxyz[1] + X22*Fxyz[2];
  303.  
  304. #if DEBUG
  305. printf( "%d: %.5e %.5e %.5e\n",
  306.  iobj, (double )(a*xyz[0]), (double )(a*xyz[1]), (double )(a*xyz[2]) );
  307. #endif
  308. i = 6 * IMOON;
  309. v[i] += a * xyz[0];
  310. v[i+2] += a * xyz[1];
  311. v[i+4] += a * xyz[2];
  312. a = GMs[IMOON];
  313. i = 6 * iobj;
  314. v[i] += a * xyz[0];
  315. v[i+2] += a * xyz[1];
  316. v[i+4] += a * xyz[2];
  317. }
  318.  
  319.  
  320.  
  321.  
  322. /* Expansion in Legendre functions.
  323.  * The potential function is
  324.  *                inf.  n    n
  325.  *      GM  [      -    -   a      m           m           m        ]
  326.  *  U = --- [ 1 -  >    >   --  ( C  cos mL + S  sin mL ) P (sin d) ]
  327.  *       r  [      -    -    n     n           n           n        ]
  328.  *                n=2  m=0  r
  329.  *
  330.  * where
  331.  * a = mean radius
  332.  * r = distance
  333.  * L = longitude re principal axes
  334.  * d = latitude  re principal axes
  335.  * P = Legendre function
  336.  *
  337.  * The part GM/r has already been evaluated by the main program.
  338.  * The part under the summations is converted to an acceleration
  339.  * by taking its gradient in spherical coordinates and transforming
  340.  * to Cartesian coordinates xi, eta, zeta.  The program then transforms
  341.  * this into the principal axis body coordinate system.
  342.  */
  343. legendre( xyz, arad, r, Jn, Bxyz, assoc )
  344. DOUBLE xyz[];  /* body coordinates of the point mass */
  345. DOUBLE arad;   /* mean radius of extended body */
  346. DOUBLE r;      /* center-to-center distance = |x,y,z| */
  347. DOUBLE Jn[];   /* zonal harmonic coefficients */
  348. DOUBLE Bxyz[]; /* output acceleration in body coordinates */
  349. int assoc;     /* 1 = do tesseral harmonics */
  350. {
  351. DOUBLE P[3], dP[3], Pnm[8], dPnm[8], cosm[4], sinm[4];
  352. DOUBLE coslat, sinlat, coslon, sinlon;
  353. DOUBLE x, x2, x2p, sx2p, a, b, c, s, aor, aorn, np1;
  354. DOUBLE xez0, xez1, xez2, t0, t1, t2;
  355. int i, k, m;
  356.  
  357. /* Latitude of the point mass, in body coordinates:
  358.  * Bz = r sin(lat)
  359.  */
  360. sinlat = xyz[2]*r;
  361. x = sinlat;
  362. x2 = sinlat * sinlat;
  363. x2p = One - x2;
  364. sx2p = SQRT( x2p );
  365. coslat = sx2p;
  366.  
  367. /* Longitude of the point mass expressed in body coordinates:
  368.  * By = r_xy sin(lon)
  369.  * Bx = r_xy cos(lon)
  370.  */
  371. a = coslat/r; /* r is 1/distance */
  372. coslon = xyz[0]/a;
  373. sinlon = xyz[1]/a;
  374.  
  375. #if DEBUG
  376. printf( "sinlat %.5e, sinlon %.5e\n", (double )sinlat, (double )sinlon );
  377. #endif
  378. /* Legendre polynomials.
  379.  * P0(x) = 1
  380.  * P1(x) = x
  381.  * (n+1) Pn+1(x) = (2n+1) x Pn(x) - n Pn-1(x)
  382.  */
  383. P[0] = OneandaHalf*x2 - Half;          /* P2(x) */
  384. P[1] = (C2r5 * x2 - OneandaHalf)*x;    /* P3(x) */
  385. P[2] = (C4r375*x2 - C3r75)*x2 + Cr375; /* P4(x) */
  386.  
  387. /* Derivatives of Legendre polynomials
  388.  * P'0(x) = 0
  389.  * P'1(x) = 1
  390.  * (x^2 - 1) P'n(x) = n[ x Pn(x) - Pn-1(x) ]
  391.  */
  392. dP[0] = Three * x;             /* dP2(x)/dx */
  393. dP[1] = C7r5*x2 - OneandaHalf; /* dP3(x)/dx */
  394. dP[2] = (C17r5*x2 - C7r5)*x;   /* dP4(x)/dx */
  395.  
  396. if( assoc )
  397.     {
  398. /*
  399.  * Associated Legendre functions
  400.  * for integer u:
  401.  *
  402.  *   u                  u/2   u        u
  403.  *  P  (x)  =  (1 - x^2)     d Pv(x)/dx
  404.  *   v
  405.  *                                                    u
  406.  * (Caution: some math textbooks multiply this by (-1)  ).
  407.  * For v also an integer,
  408.  *   u
  409.  *  P  (x)  = 0   if u > v.
  410.  *   v
  411.  *
  412.  */
  413.  
  414.     Pnm[0] = Three * x2p;               /* P^2_2(x) */
  415.     Pnm[1] = sx2p * dP[1];              /* P^1_3(x) */
  416.     Pnm[2] = C15 * x * x2p;             /* P^2_3(x) */
  417.     b = x2p * sx2p;
  418.     Pnm[3] = C15 * b;                   /* P^3_3(x) */
  419.     Pnm[4] = sx2p * dP[2];              /* P^1_4(x) */
  420.     Pnm[5] = x2p * (C52r5 * x2 - C7r5); /* P^2_4(x) */
  421.     Pnm[6] = C105 * x * b;              /* P^3_4(x) */
  422.     Pnm[7] = C105 * x2p * x2p;          /* P^4_4(x) */
  423.  
  424. /* Derivatives of associated Legendre functions
  425.  */
  426.     a = -x/x2p;
  427.     b = One/sx2p;
  428.     dPnm[0] = -Six*x;                    /* dP^2_2(x)/dx */
  429.     dPnm[1] =     a*Pnm[1] + b*Pnm[2];   /* dP^1_3(x)/dx */
  430.     dPnm[2] = Two*a*Pnm[2] + b*Pnm[3];   /* dP^2_3(x)/dx */
  431.     dPnm[3] = Three*a*Pnm[3];            /* dP^3_3(x)/dx */
  432.     dPnm[4] =     a*Pnm[4] + b*Pnm[5];   /* dP^1_4(x)/dx */
  433.     dPnm[5] = Two*a*Pnm[5] + b*Pnm[6];   /* dP^2_4(x)/dx */
  434.     dPnm[6] = Three*a*Pnm[6] + b*Pnm[7]; /* dP^3_4(x)/dx */
  435.     dPnm[7] = Four*a*Pnm[7];             /* dP^4_4(x)/dx */
  436.  
  437.  
  438. /* sines and cosines of m * longitude, m = 1, 2, 3, 4
  439.  */
  440.     c = coslon;
  441.     s = sinlon;
  442.     cosm[0] = c;
  443.     sinm[0] = s;
  444.     for( i=1; i<4; i++ )
  445.         {
  446.         b = sinlon * c + coslon * s;
  447.         c = coslon * c - sinlon * s;
  448.         s = b;
  449.         cosm[i] = c;
  450.         sinm[i] = s;
  451.         }
  452.     }
  453.  
  454. xez0 = Zero; /* xi */
  455. xez1 = Zero; /* eta */
  456. xez2 = Zero; /* zeta */
  457. np1 = Two;
  458. aor = arad * r;
  459. aorn = aor;
  460. k = 0;
  461. for( i=0; i<3; i++ )
  462.     {
  463.     aorn *= aor;
  464.     np1 += One;
  465. /* zonal harmonics */
  466.     a = Jn[i];
  467.     t0 = a * np1 * P[i];
  468.     t1 = Zero;
  469.     t2 = -coslat * a * dP[i];
  470.  
  471.     if( assoc )
  472.         {
  473. /* tesseral harmonics */
  474.         for( m=0; m<=(i+1); m++ )
  475.             {
  476.             if( (i == 0) && (m == 0) )
  477.                 continue;
  478.             a = Pnm[k];
  479.             b = Cnm[k]*cosm[m] + Snm[k]*sinm[m];
  480.             t0 -= np1*a*b;
  481.             t1 += ((m+1)*a*(-Cnm[k]*sinm[m]+Snm[k]*cosm[m]))/coslat;
  482.             t2 += coslat*dPnm[k]*b;
  483.             k += 1;
  484.             }
  485.         }
  486.     xez0 += aorn * t0;
  487.     xez1 += aorn * t1;
  488.     xez2 += aorn * t2;
  489.     }
  490. r = r*r; /* 1/r^2 */
  491. xez0 *= r;
  492. if( assoc )
  493.     xez1 *= r;
  494. xez2 *= r;
  495. /* Rotate xi,eta,zeta into the Bx,By,Bz body coordinate system.
  496.  */
  497. /* First rotate counterclockwise about the eta axis by the latitude. */
  498. Bxyz[0] = coslat*xez0 - sinlat*xez2;
  499. Bxyz[1] = xez1;
  500. Bxyz[2] = sinlat*xez0 + coslat*xez2;
  501. /* Then rotate clockwise about the new z axis by the longitude. */
  502. b  = coslon*Bxyz[0] - sinlon*Bxyz[1];
  503. Bxyz[1] = sinlon*Bxyz[0] + coslon*Bxyz[1];
  504. Bxyz[0] = b;
  505. }
  506.  
  507.  
  508. /* Differential equation for Lunar librations
  509.  */
  510. extern DOUBLE AMR2, BMR2, CMR2;
  511.  
  512. librate( v )
  513. DOUBLE v[];
  514. {
  515. DOUBLE omBx, omBy, omBz, omB1x, omB1y, omB1z;
  516. DOUBLE phi2, th2, psi2, a;
  517.  
  518. /* The angular velocity of the extended body,
  519.  * expressed in body coordinates.  The suffix 1 indicates
  520.  * differentiation with respect to time.
  521.  */
  522. a = phi1 * sinth;
  523. omBx = a*sinpsi + th1*cospsi;
  524. omBy = a*cospsi - th1*sinpsi;
  525. omBz = phi1*costh + psi1;
  526.  
  527. /*
  528.  * The ratio of the torque N to the principal moment of
  529.  * inertia C is calculated by
  530.  *
  531.  *    N          r X F          r X (M a)         r X a
  532.  *   ----  =  -----------  =  ------------  =  ------------
  533.  *    C       C/MR^2 MR^2     C/MR^2 MR^2       C/MR^2 R^2
  534.  *
  535.  * where a is the acceleration due to oblateness,
  536.  * M is the mass of the extended body,
  537.  * R is the mean radius of the extended body,
  538.  * r is the vector from the center of the body to the point mass.
  539.  * r X a = (Nx,Ny,Nz), calculated earlier.
  540.  */
  541.  
  542. /* Euler's equations for the change in angular velocity.
  543.  */
  544. omB1x = LGAMBET*omBy*omBz + Nx/AMR2;
  545. omB1y = LBET*omBz*omBx + Ny/BMR2;
  546. omB1z = -LGAM*omBx*omBy + Nz/CMR2;
  547.  
  548. /* Differential equations for the Euler angles
  549.  */
  550. phi2 = (omB1x*sinpsi + omB1y*cospsi + th1*(psi1 - phi1*costh))/sinth;
  551. th2 = omB1x*cospsi - omB1y*sinpsi - a*psi1;
  552. psi2 = omB1z - phi2*costh + a*th1;
  553.  
  554. v[0] = phi2;
  555. v[2] = th2;
  556. v[4] = psi2;
  557. #if DEBUG
  558. printf( "phi2 %.12e, th2 %.12e, psi2 %.12e\n",
  559.    (double) phi2, (double) th2, (double) psi2 );
  560. #endif
  561. }
  562.