home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / 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[