home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / de118i.zip / NUT1T.C < prev    next >
C/C++ Source or Header  |  1993-02-06  |  7KB  |  270 lines
  1. /* One-term nutation program uses only the 18.6 year term.
  2.  *
  3.  *
  4.  * References:
  5.  * "Summary of 1980 IAU Theory of Nutation (Final Report of the
  6.  * IAU Working Group on Nutation)", P. K. Seidelmann et al., in
  7.  * Transactions of the IAU Vol. XVIII A, Reports on Astronomy,
  8.  * P. A. Wayman, ed.; D. Reidel Pub. Co., 1982.
  9.  *
  10.  * "Nutation and the Earth's Rotation",
  11.  * I.A.U. Symposium No. 78, May, 1977, page 256.
  12.  * I.A.U., 1980.
  13.  *
  14.  * Woolard, E.W., "A redevelopment of the theory of nutation",
  15.  * The Astronomical Journal, 58, 1-3 (1953).
  16.  */
  17.  
  18. #ifndef NOINCS
  19. #include "mconf.h"
  20. #include "prec.h"
  21. #include "ssystem.h"
  22. #include "const.h"
  23. #endif
  24.  
  25. #if LDOUBLE
  26. #if IBMPC
  27. /* 1980 series */
  28. static short nuAi[] = {
  29. 0xd17f,0x30cc,0x4f41,0x94fb,0x3fec, XPD
  30. 0xc469,0x5b0e,0x9761,0x8834,0x3ff6, XPD
  31. 0x41eb,0x9bf7,0x5c3f,0xf1c4,0xc009, XPD
  32. 0x5656,0x21e4,0xcb9d,0xfa16,0x4005, XPD
  33. 0x14a5,0x7b74,0x6349,0x8eb4,0xbff9, XPD
  34. 0xb780,0x8240,0xc7e2,0x8998,0xc003, XPD
  35. 0x73c1,0xe1ef,0xe392,0xe94e,0x3ff4, XPD
  36. 0x3d71,0xd70a,0x70a3,0x933d,0x4002, XPD
  37. };
  38. short C3600i[] = {0x0000,0x0000,0x0000,0xe100,0x400a, XPD};
  39. short C360i[] = {0x0000,0x0000,0x0000,0xb400,0x4007, XPD};
  40. #endif /* long double IBMPC */
  41. #if MIEEE
  42. /* 1980 series */
  43. long nuAi[24] = {
  44. 0x3fec0000,0x94fb4f41,0x30ccd17f,
  45. 0x3ff60000,0x88349761,0x5b0ec469,
  46. 0xc0090000,0xf1c45c3f,0x9bf741eb,
  47. 0x40050000,0xfa16cb9d,0x21e45656,
  48. 0xbff90000,0x8eb46349,0x7b7414a5,
  49. 0xc0030000,0x8998c7e2,0x8240b780,
  50. 0x3ff40000,0xe94ee392,0xe1ef73c1,
  51. 0x40020000,0x933d70a3,0xd70a3d71,
  52. };
  53. long C3600i[] = {0x400a0000,0xe1000000,0x00000000};
  54. long C360i[] = {0x40070000,0xb4000000,0x00000000};
  55. #endif /* long double MIEEE */
  56. #define nuA ((DOUBLE *)nuAi)
  57. #define C3600 (*(DOUBLE *)C3600i)
  58. #define C360 (*(DOUBLE *)C360i)
  59. #else /* not LDOUBLE: */
  60. /* 1977 series */
  61. /*
  62. DOUBLE nuA[8] = {
  63.  8.00000000000000000000E-3,
  64.  7.48000000000000000000E0,
  65. -6.96291123000000000000E6,
  66.  9.33059790000000000000E5,
  67. -1.73700000000000000000E-2,
  68. -1.71860000000000000000E1,
  69.  9.10000000000000000000E-4,
  70.  9.20159000000000000000E0,
  71. };
  72. */
  73. /* 1980 series */
  74. #if UNK
  75. static double nuA[8] = {
  76.  2.22000000000000000000E-6,
  77.  2.07833000000000000000E-3,
  78. -1.93413626080000000000E3,
  79.  1.25044522200000000000E2,
  80. -1.74200000000000000000E-2,
  81. -1.71996000000000000000E1,
  82.  8.90000000000000000000E-4,
  83.  9.20250000000000000000E0,
  84. };
  85. double C3600 = 3600.0;
  86. double C360 = 360.0;
  87. #endif /* double UNK */
  88. #if DEC
  89. static short nuAi[32] = {
  90. 0033424,0175517,0040460,0146321,
  91. 0036010,0032227,0060533,0007304,
  92. 0142761,0142134,0037633,0173502,
  93. 0041772,0013313,0116441,0162126,
  94. 0136616,0132143,0044573,0072025,
  95. 0141211,0114307,0161202,0040270,
  96. 0035551,0047343,0111341,0167564,
  97. 0041023,0036560,0121727,0005075,
  98. };
  99. short C3600i[] = {0043141,0000000,0000000,0000000};
  100. short C360i[] = {0042264,0000000,0000000,0000000};
  101. #define nuA ((DOUBLE *)nuAi)
  102. #define C3600 (*(DOUBLE *)C3600i)
  103. #define C360 (*(DOUBLE *)C360i)
  104. #endif /* double DEC */
  105. #if IBMPC
  106. static short nuAi[40] = {
  107. 0x199a,0xe826,0x9f69,0x3ec2,
  108. 0x61d9,0xec2b,0x0692,0x3f61,
  109. 0x7ee8,0x87f3,0x388b,0xc09e,
  110. 0x3c8b,0x73a4,0x42d9,0x405f,
  111. 0x6e83,0x692f,0xd68c,0xbf91,
  112. 0x4817,0xfc50,0x3318,0xc031,
  113. 0x3dee,0x725c,0x29dc,0x3f4d,
  114. 0xe148,0x147a,0x67ae,0x4022,
  115. };
  116. short C3600i[] = {0x0000,0x0000,0x2000,0x40ac};
  117. short C360i[] = {0x0000,0x0000,0x8000,0x4076};
  118. #define nuA ((DOUBLE *)nuAi)
  119. #define C3600 (*(DOUBLE *)C3600i)
  120. #define C360 (*(DOUBLE *)C360i)
  121. #endif /* double IBMPC */
  122. #if MIEEE
  123. static short nuAi[40] = {
  124. 0x3ec2,0x9f69,0xe826,0x199a,
  125. 0x3f61,0x0692,0xec2b,0x61d9,
  126. 0xc09e,0x388b,0x87f3,0x7ee8,
  127. 0x405f,0x42d9,0x73a4,0x3c8b,
  128. 0xbf91,0xd68c,0x692f,0x6e83,
  129. 0xc031,0x3318,0xfc50,0x4817,
  130. 0x3f4d,0x29dc,0x725c,0x3dee,
  131. 0x4022,0x67ae,0x147a,0xe148,
  132. };
  133. short C3600i[] = {0x40ac,0x2000,0x0000,0x0000};
  134. short C360i[] = {0x4076,0x8000,0x0000,0x0000};
  135. #define nuA ((DOUBLE *)nuAi)
  136. #define C3600 (*(DOUBLE *)C3600i)
  137. #define C360 (*(DOUBLE *)C360i)
  138. #endif /* double MIEEE */
  139. #endif /* not LDOUBLE */
  140.  
  141. /* The answers are posted here by nutlo():
  142.  */
  143. DOUBLE jdnut = 0.0;    /* time to which the nutation applies */
  144. DOUBLE nutl = 0.0;    /* nutation in longitude (radians) */
  145. DOUBLE nuto = 0.0;    /* nutation in obliquity (radians) */
  146. extern DOUBLE eps, coseps, sineps, DTR, STR, J1900, J2000, Jcentury;
  147. DOUBLE FLOOR();
  148. #define mod360(x) (x - C360*FLOOR(x/C360))
  149.  
  150. nutlo(J)
  151. DOUBLE J;
  152. {
  153. DOUBLE T, OM, C, D;
  154. DOUBLE SIN(), COS();
  155.  
  156. if( jdnut == J )
  157.     return(0);
  158. jdnut = J;
  159.  
  160. /* 1977 IAU series */
  161. /*
  162.  * T = (J-J1900)/Jcentury;
  163.  * OM = (((nuA[0]*T + nuA[1])*T + nuA[2])*T + nuA[3])/C3600;
  164.  */
  165.  
  166. /* 1980 IAU series */
  167. T = (J-J2000)/Jcentury;
  168. OM = ((nuA[0]*T + nuA[1])*T + nuA[2])*T + nuA[3];
  169. OM = DTR * mod360(OM);
  170.  
  171. /* Woolard series
  172.  * C = (-0.01737*T - 17.2327)*SIN(OM);
  173.  * D = ( 0.00091*T +  9.2100)*COS(OM);
  174.  */
  175.  
  176. /*
  177.  * C = (nuA[4]*T + nuA[5])*SIN(OM);
  178.  * D = ( nuA[6]*T + nuA[7])*COS(OM);
  179.  */
  180.  
  181. C = nuA[5]*SIN(OM);
  182. D = nuA[7]*COS(OM);
  183.  
  184. /* Save answers, expressed in radians */
  185. nutl = STR * C;
  186. nuto = STR * D;
  187. return(1);
  188. }
  189.  
  190.  
  191.  
  192. /* Nutation -- AA page B20
  193.  * using nutation in longitude and obliquity from nutlo()
  194.  * and obliquity of the ecliptic from epsiln()
  195.  * both calculated for Julian date J.
  196.  *
  197.  * p[] = equatorial rectangular position vector of object for
  198.  * mean ecliptic and equinox of date.
  199.  */
  200. DOUBLE nuN00, nuN01, nuN02, nuN10, nuN11, nuN12, nuN20, nuN21, nuN22;
  201.  
  202. nutate( J, p )
  203. DOUBLE J;
  204. DOUBLE p[];
  205. {
  206. DOUBLE ce, se, cl, sl, sino, f;
  207. DOUBLE p1[3];
  208. int i;
  209. DOUBLE SIN(), COS();
  210.  
  211. if( J == jdnut )
  212.     goto matmpy;
  213. nutlo(J); /* be sure we calculated nutl and nuto */
  214. epsiln(J); /* and also the obliquity of date */
  215.  
  216. f = eps + nuto;
  217. ce = COS( f );
  218. se = SIN( f );
  219. sino = SIN(nuto);
  220. cl = COS( nutl );
  221. sl = SIN( nutl );
  222.  
  223. /* Apply adjustment
  224.  * to equatorial rectangular coordinates of object.
  225.  *
  226.  * This is a composite of three rotations: rotate about x axis
  227.  * to ecliptic of date; rotate about new z axis by the nutation
  228.  * in longitude; rotate about new x axis back to equator of date
  229.  * plus nutation in obliquity.
  230.  */
  231. nuN00 =   cl;
  232. nuN01 =  -sl*coseps;
  233. nuN02 =  -sl*sineps;
  234.  
  235. nuN10 = sl*ce;
  236. nuN11 = cl*coseps*ce + sineps*se;
  237. nuN12 = -( sino + (One-cl)*sineps*se );
  238.  
  239. nuN20 = sl*se;
  240. nuN21 = sino + (cl-One)*se*coseps;
  241. nuN22 = cl*sineps*se + coseps*ce;
  242.  
  243. matmpy:
  244.  
  245. p1[0] = nuN00*p[0] + nuN01*p[1] + nuN02*p[2];
  246. p1[1] = nuN10*p[0] + nuN11*p[1] + nuN12*p[2];
  247. p1[2] = nuN20*p[0] + nuN21*p[1] + nuN22*p[2];
  248.  
  249. for( i=0; i<3; i++ )
  250.     p[i] = p1[i];
  251. }
  252.  
  253.  
  254. /* Inverse transformation using previously
  255.  * computed matrix elements.
  256.  */
  257. invnut( p1 )
  258. DOUBLE p1[];
  259. {
  260. DOUBLE p[3];
  261. int i;
  262.  
  263. p[0] = nuN00*p1[0] + nuN10*p1[1] + nuN20*p1[2];
  264. p[1] = nuN01*p1[0] + nuN11*p1[1] + nuN21*p1[2];
  265. p[2] = nuN02*p1[0] + nuN12*p1[1] + nuN22*p1[2];
  266.  
  267. for( i=0; i<3; i++ )
  268.     p1[i] = p[i];
  269. }
  270.