| 1 | /*
|
|---|
| 2 | *+
|
|---|
| 3 | * Name:
|
|---|
| 4 | * palDmoon
|
|---|
| 5 |
|
|---|
| 6 | * Purpose:
|
|---|
| 7 | * Approximate geocentric position and velocity of the Moon
|
|---|
| 8 |
|
|---|
| 9 | * Language:
|
|---|
| 10 | * Starlink ANSI C
|
|---|
| 11 |
|
|---|
| 12 | * Type of Module:
|
|---|
| 13 | * Library routine
|
|---|
| 14 |
|
|---|
| 15 | * Invocation:
|
|---|
| 16 | * void palDmoon( double date, double pv[6] );
|
|---|
| 17 |
|
|---|
| 18 | * Arguments:
|
|---|
| 19 | * date = double (Given)
|
|---|
| 20 | * TDB as a Modified Julian Date (JD-2400000.5)
|
|---|
| 21 | * pv = double [6] (Returned)
|
|---|
| 22 | * Moon x,y,z,xdot,ydot,zdot, mean equator and
|
|---|
| 23 | * equinox of date (AU, AU/s)
|
|---|
| 24 |
|
|---|
| 25 | * Description:
|
|---|
| 26 | * Calculate the approximate geocentric position of the Moon
|
|---|
| 27 | * using a full implementation of the algorithm published by
|
|---|
| 28 | * Meeus (l'Astronomie, June 1984, p348).
|
|---|
| 29 |
|
|---|
| 30 | * Authors:
|
|---|
| 31 | * TIMJ: Tim Jenness (JAC, Hawaii)
|
|---|
| 32 | * PTW: Patrick T. Wallace
|
|---|
| 33 | * {enter_new_authors_here}
|
|---|
| 34 |
|
|---|
| 35 | * Notes:
|
|---|
| 36 | * - Meeus quotes accuracies of 10 arcsec in longitude, 3 arcsec in
|
|---|
| 37 | * latitude and 0.2 arcsec in HP (equivalent to about 20 km in
|
|---|
| 38 | * distance). Comparison with JPL DE200 over the interval
|
|---|
| 39 | * 1960-2025 gives RMS errors of 3.7 arcsec and 83 mas/hour in
|
|---|
| 40 | * longitude, 2.3 arcsec and 48 mas/hour in latitude, 11 km
|
|---|
| 41 | * and 81 mm/s in distance. The maximum errors over the same
|
|---|
| 42 | * interval are 18 arcsec and 0.50 arcsec/hour in longitude,
|
|---|
| 43 | * 11 arcsec and 0.24 arcsec/hour in latitude, 40 km and 0.29 m/s
|
|---|
| 44 | * in distance.
|
|---|
| 45 | * - The original algorithm is expressed in terms of the obsolete
|
|---|
| 46 | * timescale Ephemeris Time. Either TDB or TT can be used, but
|
|---|
| 47 | * not UT without incurring significant errors (30 arcsec at
|
|---|
| 48 | * the present time) due to the Moon's 0.5 arcsec/sec movement.
|
|---|
| 49 | * - The algorithm is based on pre IAU 1976 standards. However,
|
|---|
| 50 | * the result has been moved onto the new (FK5) equinox, an
|
|---|
| 51 | * adjustment which is in any case much smaller than the
|
|---|
| 52 | * intrinsic accuracy of the procedure.
|
|---|
| 53 | * - Velocity is obtained by a complete analytical differentiation
|
|---|
| 54 | * of the Meeus model.
|
|---|
| 55 |
|
|---|
| 56 | * History:
|
|---|
| 57 | * 2012-03-07 (TIMJ):
|
|---|
| 58 | * Initial version based on a direct port of the SLA/F code.
|
|---|
| 59 | * Adapted with permission from the Fortran SLALIB library.
|
|---|
| 60 | * {enter_further_changes_here}
|
|---|
| 61 |
|
|---|
| 62 | * Copyright:
|
|---|
| 63 | * Copyright (C) 1998 Rutherford Appleton Laboratory
|
|---|
| 64 | * Copyright (C) 2012 Science and Technology Facilities Council.
|
|---|
| 65 | * All Rights Reserved.
|
|---|
| 66 |
|
|---|
| 67 | * Licence:
|
|---|
| 68 | * This program is free software; you can redistribute it and/or
|
|---|
| 69 | * modify it under the terms of the GNU General Public License as
|
|---|
| 70 | * published by the Free Software Foundation; either version 3 of
|
|---|
| 71 | * the License, or (at your option) any later version.
|
|---|
| 72 | *
|
|---|
| 73 | * This program is distributed in the hope that it will be
|
|---|
| 74 | * useful, but WITHOUT ANY WARRANTY; without even the implied
|
|---|
| 75 | * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|---|
| 76 | * PURPOSE. See the GNU General Public License for more details.
|
|---|
| 77 | *
|
|---|
| 78 | * You should have received a copy of the GNU General Public License
|
|---|
| 79 | * along with this program; if not, write to the Free Software
|
|---|
| 80 | * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
|
|---|
| 81 | * MA 02110-1301, USA.
|
|---|
| 82 |
|
|---|
| 83 | * Bugs:
|
|---|
| 84 | * {note_any_bugs_here}
|
|---|
| 85 | *-
|
|---|
| 86 | */
|
|---|
| 87 |
|
|---|
| 88 | #include "pal.h"
|
|---|
| 89 | #include "pal1sofa.h"
|
|---|
| 90 | #include "palmac.h"
|
|---|
| 91 |
|
|---|
| 92 | /* Autoconf can give us -DPIC */
|
|---|
| 93 | #undef PIC
|
|---|
| 94 |
|
|---|
| 95 | void palDmoon( double date, double pv[6] ) {
|
|---|
| 96 |
|
|---|
| 97 | /* Seconds per Julian century (86400*36525) */
|
|---|
| 98 | const double CJ = 3155760000.0;
|
|---|
| 99 |
|
|---|
| 100 | /* Julian epoch of B1950 */
|
|---|
| 101 | const double B1950 = 1949.9997904423;
|
|---|
| 102 |
|
|---|
| 103 | /* Earth equatorial radius in AU ( = 6378.137 / 149597870 ) */
|
|---|
| 104 | const double ERADAU=4.2635212653763e-5;
|
|---|
| 105 |
|
|---|
| 106 | double T,THETA,SINOM,COSOM,DOMCOM,WA,DWA,WB,DWB,WOM,
|
|---|
| 107 | DWOM,SINWOM,COSWOM,V,DV,COEFF,EMN,EMPN,DN,FN,EN,
|
|---|
| 108 | DEN,DTHETA,FTHETA,EL,DEL,B,DB,BF,DBF,P,DP,SP,R,
|
|---|
| 109 | DR,X,Y,Z,XD,YD,ZD,SEL,CEL,SB,CB,RCB,RBD,W,EPJ,
|
|---|
| 110 | EQCOR,EPS,SINEPS,COSEPS,ES,EC;
|
|---|
| 111 | double ELP, DELP;
|
|---|
| 112 | double EM, DEM, EMP, DEMP, D, DD, F, DF, OM, DOM, E, DESQ, ESQ, DE;
|
|---|
| 113 | int N,I;
|
|---|
| 114 |
|
|---|
| 115 | /*
|
|---|
| 116 | * Coefficients for fundamental arguments
|
|---|
| 117 | *
|
|---|
| 118 | * at J1900: T**0, T**1, T**2, T**3
|
|---|
| 119 | * at epoch: T**0, T**1
|
|---|
| 120 | *
|
|---|
| 121 | * Units are degrees for position and Julian centuries for time
|
|---|
| 122 | *
|
|---|
| 123 | */
|
|---|
| 124 |
|
|---|
| 125 | /* Moon's mean longitude */
|
|---|
| 126 | const double ELP0=270.434164;
|
|---|
| 127 | const double ELP1=481267.8831;
|
|---|
| 128 | const double ELP2=-0.001133;
|
|---|
| 129 | const double ELP3=0.0000019;
|
|---|
| 130 |
|
|---|
| 131 | /* Sun's mean anomaly */
|
|---|
| 132 | const double EM0=358.475833;
|
|---|
| 133 | const double EM1=35999.0498;
|
|---|
| 134 | const double EM2=-0.000150;
|
|---|
| 135 | const double EM3=-0.0000033;
|
|---|
| 136 |
|
|---|
| 137 | /* Moon's mean anomaly */
|
|---|
| 138 | const double EMP0=296.104608;
|
|---|
| 139 | const double EMP1=477198.8491;
|
|---|
| 140 | const double EMP2=0.009192;
|
|---|
| 141 | const double EMP3=0.0000144;
|
|---|
| 142 |
|
|---|
| 143 | /* Moon's mean elongation */
|
|---|
| 144 | const double D0=350.737486;
|
|---|
| 145 | const double D1=445267.1142;
|
|---|
| 146 | const double D2=-0.001436;
|
|---|
| 147 | const double D3=0.0000019;
|
|---|
| 148 |
|
|---|
| 149 | /* Mean distance of the Moon from its ascending node */
|
|---|
| 150 | const double F0=11.250889;
|
|---|
| 151 | const double F1=483202.0251;
|
|---|
| 152 | const double F2=-0.003211;
|
|---|
| 153 | const double F3=-0.0000003;
|
|---|
| 154 |
|
|---|
| 155 | /* Longitude of the Moon's ascending node */
|
|---|
| 156 | const double OM0=259.183275;
|
|---|
| 157 | const double OM1=-1934.1420;
|
|---|
| 158 | const double OM2=0.002078;
|
|---|
| 159 | const double OM3=0.0000022;
|
|---|
| 160 |
|
|---|
| 161 | /* Coefficients for (dimensionless) E factor */
|
|---|
| 162 | const double E1=-0.002495;
|
|---|
| 163 | const double E2=-0.00000752;
|
|---|
| 164 |
|
|---|
| 165 | /* Coefficients for periodic variations etc */
|
|---|
| 166 | const double PAC=0.000233;
|
|---|
| 167 | const double PA0=51.2;
|
|---|
| 168 | const double PA1=20.2;
|
|---|
| 169 | const double PBC=-0.001778;
|
|---|
| 170 | const double PCC=0.000817;
|
|---|
| 171 | const double PDC=0.002011;
|
|---|
| 172 | const double PEC=0.003964;
|
|---|
| 173 | const double PE0=346.560;
|
|---|
| 174 | const double PE1=132.870;
|
|---|
| 175 | const double PE2=-0.0091731;
|
|---|
| 176 | const double PFC=0.001964;
|
|---|
| 177 | const double PGC=0.002541;
|
|---|
| 178 | const double PHC=0.001964;
|
|---|
| 179 | const double PIC=-0.024691;
|
|---|
| 180 | const double PJC=-0.004328;
|
|---|
| 181 | const double PJ0=275.05;
|
|---|
| 182 | const double PJ1=-2.30;
|
|---|
| 183 | const double CW1=0.0004664;
|
|---|
| 184 | const double CW2=0.0000754;
|
|---|
| 185 |
|
|---|
| 186 | /*
|
|---|
| 187 | * Coefficients for Moon position
|
|---|
| 188 | *
|
|---|
| 189 | * Tx(N) = coefficient of L, B or P term (deg)
|
|---|
| 190 | * ITx(N,1-5) = coefficients of M, M', D, F, E**n in argument
|
|---|
| 191 | */
|
|---|
| 192 | #define NL 50
|
|---|
| 193 | #define NB 45
|
|---|
| 194 | #define NP 31
|
|---|
| 195 |
|
|---|
| 196 | /*
|
|---|
| 197 | * Longitude
|
|---|
| 198 | */
|
|---|
| 199 | const double TL[NL] = {
|
|---|
| 200 | 6.28875,1.274018,.658309,.213616,-.185596,
|
|---|
| 201 | -.114336,.058793,.057212,.05332,.045874,.041024,-.034718,-.030465,
|
|---|
| 202 | .015326,-.012528,-.01098,.010674,.010034,.008548,-.00791,-.006783,
|
|---|
| 203 | .005162,.005,.004049,.003996,.003862,.003665,.002695,.002602,
|
|---|
| 204 | .002396,-.002349,.002249,-.002125,-.002079,.002059,-.001773,
|
|---|
| 205 | -.001595,.00122,-.00111,8.92e-4,-8.11e-4,7.61e-4,7.17e-4,7.04e-4,
|
|---|
| 206 | 6.93e-4,5.98e-4,5.5e-4,5.38e-4,5.21e-4,4.86e-4
|
|---|
| 207 | };
|
|---|
| 208 |
|
|---|
| 209 | const int ITL[NL][5] = {
|
|---|
| 210 | /* M M' D F n */
|
|---|
| 211 | { +0, +1, +0, +0, 0 },
|
|---|
| 212 | { +0, -1, +2, +0, 0 },
|
|---|
| 213 | { +0, +0, +2, +0, 0 },
|
|---|
| 214 | { +0, +2, +0, +0, 0 },
|
|---|
| 215 | { +1, +0, +0, +0, 1 },
|
|---|
| 216 | { +0, +0, +0, +2, 0 },
|
|---|
| 217 | { +0, -2, +2, +0, 0 },
|
|---|
| 218 | { -1, -1, +2, +0, 1 },
|
|---|
| 219 | { +0, +1, +2, +0, 0 },
|
|---|
| 220 | { -1, +0, +2, +0, 1 },
|
|---|
| 221 | { -1, +1, +0, +0, 1 },
|
|---|
| 222 | { +0, +0, +1, +0, 0 },
|
|---|
| 223 | { +1, +1, +0, +0, 1 },
|
|---|
| 224 | { +0, +0, +2, -2, 0 },
|
|---|
| 225 | { +0, +1, +0, +2, 0 },
|
|---|
| 226 | { +0, -1, +0, +2, 0 },
|
|---|
| 227 | { +0, -1, +4, +0, 0 },
|
|---|
| 228 | { +0, +3, +0, +0, 0 },
|
|---|
| 229 | { +0, -2, +4, +0, 0 },
|
|---|
| 230 | { +1, -1, +2, +0, 1 },
|
|---|
| 231 | { +1, +0, +2, +0, 1 },
|
|---|
| 232 | { +0, +1, -1, +0, 0 },
|
|---|
| 233 | { +1, +0, +1, +0, 1 },
|
|---|
| 234 | { -1, +1, +2, +0, 1 },
|
|---|
| 235 | { +0, +2, +2, +0, 0 },
|
|---|
| 236 | { +0, +0, +4, +0, 0 },
|
|---|
| 237 | { +0, -3, +2, +0, 0 },
|
|---|
| 238 | { -1, +2, +0, +0, 1 },
|
|---|
| 239 | { +0, +1, -2, -2, 0 },
|
|---|
| 240 | { -1, -2, +2, +0, 1 },
|
|---|
| 241 | { +0, +1, +1, +0, 0 },
|
|---|
| 242 | { -2, +0, +2, +0, 2 },
|
|---|
| 243 | { +1, +2, +0, +0, 1 },
|
|---|
| 244 | { +2, +0, +0, +0, 2 },
|
|---|
| 245 | { -2, -1, +2, +0, 2 },
|
|---|
| 246 | { +0, +1, +2, -2, 0 },
|
|---|
| 247 | { +0, +0, +2, +2, 0 },
|
|---|
| 248 | { -1, -1, +4, +0, 1 },
|
|---|
| 249 | { +0, +2, +0, +2, 0 },
|
|---|
| 250 | { +0, +1, -3, +0, 0 },
|
|---|
| 251 | { +1, +1, +2, +0, 1 },
|
|---|
| 252 | { -1, -2, +4, +0, 1 },
|
|---|
| 253 | { -2, +1, +0, +0, 2 },
|
|---|
| 254 | { -2, +1, -2, +0, 2 },
|
|---|
| 255 | { +1, -2, +2, +0, 1 },
|
|---|
| 256 | { -1, +0, +2, -2, 1 },
|
|---|
| 257 | { +0, +1, +4, +0, 0 },
|
|---|
| 258 | { +0, +4, +0, +0, 0 },
|
|---|
| 259 | { -1, +0, +4, +0, 1 },
|
|---|
| 260 | { +0, +2, -1, +0, 0 }
|
|---|
| 261 | };
|
|---|
| 262 |
|
|---|
| 263 | /*
|
|---|
| 264 | * Latitude
|
|---|
| 265 | */
|
|---|
| 266 | const double TB[NB] = {
|
|---|
| 267 | 5.128189,.280606,.277693,.173238,.055413,
|
|---|
| 268 | .046272,.032573,.017198,.009267,.008823,.008247,.004323,.0042,
|
|---|
| 269 | .003372,.002472,.002222,.002072,.001877,.001828,-.001803,-.00175,
|
|---|
| 270 | .00157,-.001487,-.001481,.001417,.00135,.00133,.001106,.00102,
|
|---|
| 271 | 8.33e-4,7.81e-4,6.7e-4,6.06e-4,5.97e-4,4.92e-4,4.5e-4,4.39e-4,
|
|---|
| 272 | 4.23e-4,4.22e-4,-3.67e-4,-3.53e-4,3.31e-4,3.17e-4,3.06e-4,
|
|---|
| 273 | -2.83e-4
|
|---|
| 274 | };
|
|---|
| 275 |
|
|---|
| 276 | const int ITB[NB][5] = {
|
|---|
| 277 | /* M M' D F n */
|
|---|
| 278 | { +0, +0, +0, +1, 0 },
|
|---|
| 279 | { +0, +1, +0, +1, 0 },
|
|---|
| 280 | { +0, +1, +0, -1, 0 },
|
|---|
| 281 | { +0, +0, +2, -1, 0 },
|
|---|
| 282 | { +0, -1, +2, +1, 0 },
|
|---|
| 283 | { +0, -1, +2, -1, 0 },
|
|---|
| 284 | { +0, +0, +2, +1, 0 },
|
|---|
| 285 | { +0, +2, +0, +1, 0 },
|
|---|
| 286 | { +0, +1, +2, -1, 0 },
|
|---|
| 287 | { +0, +2, +0, -1, 0 },
|
|---|
| 288 | { -1, +0, +2, -1, 1 },
|
|---|
| 289 | { +0, -2, +2, -1, 0 },
|
|---|
| 290 | { +0, +1, +2, +1, 0 },
|
|---|
| 291 | { -1, +0, -2, +1, 1 },
|
|---|
| 292 | { -1, -1, +2, +1, 1 },
|
|---|
| 293 | { -1, +0, +2, +1, 1 },
|
|---|
| 294 | { -1, -1, +2, -1, 1 },
|
|---|
| 295 | { -1, +1, +0, +1, 1 },
|
|---|
| 296 | { +0, -1, +4, -1, 0 },
|
|---|
| 297 | { +1, +0, +0, +1, 1 },
|
|---|
| 298 | { +0, +0, +0, +3, 0 },
|
|---|
| 299 | { -1, +1, +0, -1, 1 },
|
|---|
| 300 | { +0, +0, +1, +1, 0 },
|
|---|
| 301 | { +1, +1, +0, +1, 1 },
|
|---|
| 302 | { -1, -1, +0, +1, 1 },
|
|---|
| 303 | { -1, +0, +0, +1, 1 },
|
|---|
| 304 | { +0, +0, -1, +1, 0 },
|
|---|
| 305 | { +0, +3, +0, +1, 0 },
|
|---|
| 306 | { +0, +0, +4, -1, 0 },
|
|---|
| 307 | { +0, -1, +4, +1, 0 },
|
|---|
| 308 | { +0, +1, +0, -3, 0 },
|
|---|
| 309 | { +0, -2, +4, +1, 0 },
|
|---|
| 310 | { +0, +0, +2, -3, 0 },
|
|---|
| 311 | { +0, +2, +2, -1, 0 },
|
|---|
| 312 | { -1, +1, +2, -1, 1 },
|
|---|
| 313 | { +0, +2, -2, -1, 0 },
|
|---|
| 314 | { +0, +3, +0, -1, 0 },
|
|---|
| 315 | { +0, +2, +2, +1, 0 },
|
|---|
| 316 | { +0, -3, +2, -1, 0 },
|
|---|
| 317 | { +1, -1, +2, +1, 1 },
|
|---|
| 318 | { +1, +0, +2, +1, 1 },
|
|---|
| 319 | { +0, +0, +4, +1, 0 },
|
|---|
| 320 | { -1, +1, +2, +1, 1 },
|
|---|
| 321 | { -2, +0, +2, -1, 2 },
|
|---|
| 322 | { +0, +1, +0, +3, 0 }
|
|---|
| 323 | };
|
|---|
| 324 |
|
|---|
| 325 | /*
|
|---|
| 326 | * Parallax
|
|---|
| 327 | */
|
|---|
| 328 | const double TP[NP] = {
|
|---|
| 329 | .950724,.051818,.009531,.007843,.002824,
|
|---|
| 330 | 8.57e-4,5.33e-4,4.01e-4,3.2e-4,-2.71e-4,-2.64e-4,-1.98e-4,1.73e-4,
|
|---|
| 331 | 1.67e-4,-1.11e-4,1.03e-4,-8.4e-5,-8.3e-5,7.9e-5,7.2e-5,6.4e-5,
|
|---|
| 332 | -6.3e-5,4.1e-5,3.5e-5,-3.3e-5,-3e-5,-2.9e-5,-2.9e-5,2.6e-5,
|
|---|
| 333 | -2.3e-5,1.9e-5
|
|---|
| 334 | };
|
|---|
| 335 |
|
|---|
| 336 | const int ITP[NP][5] = {
|
|---|
| 337 | /* M M' D F n */
|
|---|
| 338 | { +0, +0, +0, +0, 0 },
|
|---|
| 339 | { +0, +1, +0, +0, 0 },
|
|---|
| 340 | { +0, -1, +2, +0, 0 },
|
|---|
| 341 | { +0, +0, +2, +0, 0 },
|
|---|
| 342 | { +0, +2, +0, +0, 0 },
|
|---|
| 343 | { +0, +1, +2, +0, 0 },
|
|---|
| 344 | { -1, +0, +2, +0, 1 },
|
|---|
| 345 | { -1, -1, +2, +0, 1 },
|
|---|
| 346 | { -1, +1, +0, +0, 1 },
|
|---|
| 347 | { +0, +0, +1, +0, 0 },
|
|---|
| 348 | { +1, +1, +0, +0, 1 },
|
|---|
| 349 | { +0, -1, +0, +2, 0 },
|
|---|
| 350 | { +0, +3, +0, +0, 0 },
|
|---|
| 351 | { +0, -1, +4, +0, 0 },
|
|---|
| 352 | { +1, +0, +0, +0, 1 },
|
|---|
| 353 | { +0, -2, +4, +0, 0 },
|
|---|
| 354 | { +0, +2, -2, +0, 0 },
|
|---|
| 355 | { +1, +0, +2, +0, 1 },
|
|---|
| 356 | { +0, +2, +2, +0, 0 },
|
|---|
| 357 | { +0, +0, +4, +0, 0 },
|
|---|
| 358 | { -1, +1, +2, +0, 1 },
|
|---|
| 359 | { +1, -1, +2, +0, 1 },
|
|---|
| 360 | { +1, +0, +1, +0, 1 },
|
|---|
| 361 | { -1, +2, +0, +0, 1 },
|
|---|
| 362 | { +0, +3, -2, +0, 0 },
|
|---|
| 363 | { +0, +1, +1, +0, 0 },
|
|---|
| 364 | { +0, +0, -2, +2, 0 },
|
|---|
| 365 | { +1, +2, +0, +0, 1 },
|
|---|
| 366 | { -2, +0, +2, +0, 2 },
|
|---|
| 367 | { +0, +1, -2, +2, 0 },
|
|---|
| 368 | { -1, -1, +4, +0, 1 }
|
|---|
| 369 | };
|
|---|
| 370 |
|
|---|
| 371 | /* Centuries since J1900 */
|
|---|
| 372 | T=(date-15019.5)/36525.;
|
|---|
| 373 |
|
|---|
| 374 | /*
|
|---|
| 375 | * Fundamental arguments (radians) and derivatives (radians per
|
|---|
| 376 | * Julian century) for the current epoch
|
|---|
| 377 | */
|
|---|
| 378 |
|
|---|
| 379 | /* Moon's mean longitude */
|
|---|
| 380 | ELP=PAL__DD2R*fmod(ELP0+(ELP1+(ELP2+ELP3*T)*T)*T,360.);
|
|---|
| 381 | DELP=PAL__DD2R*(ELP1+(2.*ELP2+3*ELP3*T)*T);
|
|---|
| 382 |
|
|---|
| 383 | /* Sun's mean anomaly */
|
|---|
| 384 | EM=PAL__DD2R*fmod(EM0+(EM1+(EM2+EM3*T)*T)*T,360.);
|
|---|
| 385 | DEM=PAL__DD2R*(EM1+(2.*EM2+3*EM3*T)*T);
|
|---|
| 386 |
|
|---|
| 387 | /* Moon's mean anomaly */
|
|---|
| 388 | EMP=PAL__DD2R*fmod(EMP0+(EMP1+(EMP2+EMP3*T)*T)*T,360.);
|
|---|
| 389 | DEMP=PAL__DD2R*(EMP1+(2.*EMP2+3*EMP3*T)*T);
|
|---|
| 390 |
|
|---|
| 391 | /* Moon's mean elongation */
|
|---|
| 392 | D=PAL__DD2R*fmod(D0+(D1+(D2+D3*T)*T)*T,360.);
|
|---|
| 393 | DD=PAL__DD2R*(D1+(2.*D2+3.*D3*T)*T);
|
|---|
| 394 |
|
|---|
| 395 | /* Mean distance of the Moon from its ascending node */
|
|---|
| 396 | F=PAL__DD2R*fmod(F0+(F1+(F2+F3*T)*T)*T,360.);
|
|---|
| 397 | DF=PAL__DD2R*(F1+(2.*F2+3.*F3*T)*T);
|
|---|
| 398 |
|
|---|
| 399 | /* Longitude of the Moon's ascending node */
|
|---|
| 400 | OM=PAL__DD2R*fmod(OM0+(OM1+(OM2+OM3*T)*T)*T,360.);
|
|---|
| 401 | DOM=PAL__DD2R*(OM1+(2.*OM2+3.*OM3*T)*T);
|
|---|
| 402 | SINOM=sin(OM);
|
|---|
| 403 | COSOM=cos(OM);
|
|---|
| 404 | DOMCOM=DOM*COSOM;
|
|---|
| 405 |
|
|---|
| 406 | /* Add the periodic variations */
|
|---|
| 407 | THETA=PAL__DD2R*(PA0+PA1*T);
|
|---|
| 408 | WA=sin(THETA);
|
|---|
| 409 | DWA=PAL__DD2R*PA1*cos(THETA);
|
|---|
| 410 | THETA=PAL__DD2R*(PE0+(PE1+PE2*T)*T);
|
|---|
| 411 | WB=PEC*sin(THETA);
|
|---|
| 412 | DWB=PAL__DD2R*PEC*(PE1+2.*PE2*T)*cos(THETA);
|
|---|
| 413 | ELP=ELP+PAL__DD2R*(PAC*WA+WB+PFC*SINOM);
|
|---|
| 414 | DELP=DELP+PAL__DD2R*(PAC*DWA+DWB+PFC*DOMCOM);
|
|---|
| 415 | EM=EM+PAL__DD2R*PBC*WA;
|
|---|
| 416 | DEM=DEM+PAL__DD2R*PBC*DWA;
|
|---|
| 417 | EMP=EMP+PAL__DD2R*(PCC*WA+WB+PGC*SINOM);
|
|---|
| 418 | DEMP=DEMP+PAL__DD2R*(PCC*DWA+DWB+PGC*DOMCOM);
|
|---|
| 419 | D=D+PAL__DD2R*(PDC*WA+WB+PHC*SINOM);
|
|---|
| 420 | DD=DD+PAL__DD2R*(PDC*DWA+DWB+PHC*DOMCOM);
|
|---|
| 421 | WOM=OM+PAL__DD2R*(PJ0+PJ1*T);
|
|---|
| 422 | DWOM=DOM+PAL__DD2R*PJ1;
|
|---|
| 423 | SINWOM=sin(WOM);
|
|---|
| 424 | COSWOM=cos(WOM);
|
|---|
| 425 | F=F+PAL__DD2R*(WB+PIC*SINOM+PJC*SINWOM);
|
|---|
| 426 | DF=DF+PAL__DD2R*(DWB+PIC*DOMCOM+PJC*DWOM*COSWOM);
|
|---|
| 427 |
|
|---|
| 428 | /* E-factor, and square */
|
|---|
| 429 | E=1.+(E1+E2*T)*T;
|
|---|
| 430 | DE=E1+2.*E2*T;
|
|---|
| 431 | ESQ=E*E;
|
|---|
| 432 | DESQ=2.*E*DE;
|
|---|
| 433 |
|
|---|
| 434 | /*
|
|---|
| 435 | * Series expansions
|
|---|
| 436 | */
|
|---|
| 437 |
|
|---|
| 438 | /* Longitude */
|
|---|
| 439 | V=0.;
|
|---|
| 440 | DV=0.;
|
|---|
| 441 | for (N=NL-1; N>=0; N--) { /* DO N=NL, 1, -1 */
|
|---|
| 442 | COEFF=TL[N];
|
|---|
| 443 | EMN=(double)(ITL[N][0]);
|
|---|
| 444 | EMPN=(double)(ITL[N][1]);
|
|---|
| 445 | DN=(double)(ITL[N][2]);
|
|---|
| 446 | FN=(double)(ITL[N][3]);
|
|---|
| 447 | I=ITL[N][4];
|
|---|
| 448 | if (I == 0) {
|
|---|
| 449 | EN=1.;
|
|---|
| 450 | DEN=0.;
|
|---|
| 451 | } else if (I == 1) {
|
|---|
| 452 | EN=E;
|
|---|
| 453 | DEN=DE;
|
|---|
| 454 | } else {
|
|---|
| 455 | EN=ESQ;
|
|---|
| 456 | DEN=DESQ;
|
|---|
| 457 | }
|
|---|
| 458 | THETA=EMN*EM+EMPN*EMP+DN*D+FN*F;
|
|---|
| 459 | DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF;
|
|---|
| 460 | FTHETA=sin(THETA);
|
|---|
| 461 | V=V+COEFF*FTHETA*EN;
|
|---|
| 462 | DV=DV+COEFF*(cos(THETA)*DTHETA*EN+FTHETA*DEN);
|
|---|
| 463 | }
|
|---|
| 464 | EL=ELP+PAL__DD2R*V;
|
|---|
| 465 | DEL=(DELP+PAL__DD2R*DV)/CJ;
|
|---|
| 466 |
|
|---|
| 467 | /* Latitude */
|
|---|
| 468 | V=0.;
|
|---|
| 469 | DV=0.;
|
|---|
| 470 | for (N=NB-1; N>=0; N--) { /* DO N=NB,1,-1 */
|
|---|
| 471 | COEFF=TB[N];
|
|---|
| 472 | EMN=(double)(ITB[N][0]);
|
|---|
| 473 | EMPN=(double)(ITB[N][1]);
|
|---|
| 474 | DN=(double)(ITB[N][2]);
|
|---|
| 475 | FN=(double)(ITB[N][3]);
|
|---|
| 476 | I=ITB[N][4];
|
|---|
| 477 | if (I == 0 ) {
|
|---|
| 478 | EN=1.;
|
|---|
| 479 | DEN=0.;
|
|---|
| 480 | } else if (I == 1) {
|
|---|
| 481 | EN=E;
|
|---|
| 482 | DEN=DE;
|
|---|
| 483 | } else {
|
|---|
| 484 | EN=ESQ;
|
|---|
| 485 | DEN=DESQ;
|
|---|
| 486 | }
|
|---|
| 487 | THETA=EMN*EM+EMPN*EMP+DN*D+FN*F;
|
|---|
| 488 | DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF;
|
|---|
| 489 | FTHETA=sin(THETA);
|
|---|
| 490 | V=V+COEFF*FTHETA*EN;
|
|---|
| 491 | DV=DV+COEFF*(cos(THETA)*DTHETA*EN+FTHETA*DEN);
|
|---|
| 492 | }
|
|---|
| 493 | BF=1.-CW1*COSOM-CW2*COSWOM;
|
|---|
| 494 | DBF=CW1*DOM*SINOM+CW2*DWOM*SINWOM;
|
|---|
| 495 | B=PAL__DD2R*V*BF;
|
|---|
| 496 | DB=PAL__DD2R*(DV*BF+V*DBF)/CJ;
|
|---|
| 497 |
|
|---|
| 498 | /* Parallax */
|
|---|
| 499 | V=0.;
|
|---|
| 500 | DV=0.;
|
|---|
| 501 | for (N=NP-1; N>=0; N--) { /* DO N=NP,1,-1 */
|
|---|
| 502 | COEFF=TP[N];
|
|---|
| 503 | EMN=(double)(ITP[N][0]);
|
|---|
| 504 | EMPN=(double)(ITP[N][1]);
|
|---|
| 505 | DN=(double)(ITP[N][2]);
|
|---|
| 506 | FN=(double)(ITP[N][3]);
|
|---|
| 507 | I=ITP[N][4];
|
|---|
| 508 | if (I == 0) {
|
|---|
| 509 | EN=1.;
|
|---|
| 510 | DEN=0.;
|
|---|
| 511 | } else if (I == 1) {
|
|---|
| 512 | EN=E;
|
|---|
| 513 | DEN=DE;
|
|---|
| 514 | } else {
|
|---|
| 515 | EN=ESQ;
|
|---|
| 516 | DEN=DESQ;
|
|---|
| 517 | }
|
|---|
| 518 | THETA=EMN*EM+EMPN*EMP+DN*D+FN*F;
|
|---|
| 519 | DTHETA=EMN*DEM+EMPN*DEMP+DN*DD+FN*DF;
|
|---|
| 520 | FTHETA=cos(THETA);
|
|---|
| 521 | V=V+COEFF*FTHETA*EN;
|
|---|
| 522 | DV=DV+COEFF*(-sin(THETA)*DTHETA*EN+FTHETA*DEN);
|
|---|
| 523 | }
|
|---|
| 524 | P=PAL__DD2R*V;
|
|---|
| 525 | DP=PAL__DD2R*DV/CJ;
|
|---|
| 526 |
|
|---|
| 527 | /*
|
|---|
| 528 | * Transformation into final form
|
|---|
| 529 | */
|
|---|
| 530 |
|
|---|
| 531 | /* Parallax to distance (AU, AU/sec) */
|
|---|
| 532 | SP=sin(P);
|
|---|
| 533 | R=ERADAU/SP;
|
|---|
| 534 | DR=-R*DP*cos(P)/SP;
|
|---|
| 535 |
|
|---|
| 536 | /* Longitude, latitude to x,y,z (AU) */
|
|---|
| 537 | SEL=sin(EL);
|
|---|
| 538 | CEL=cos(EL);
|
|---|
| 539 | SB=sin(B);
|
|---|
| 540 | CB=cos(B);
|
|---|
| 541 | RCB=R*CB;
|
|---|
| 542 | RBD=R*DB;
|
|---|
| 543 | W=RBD*SB-CB*DR;
|
|---|
| 544 | X=RCB*CEL;
|
|---|
| 545 | Y=RCB*SEL;
|
|---|
| 546 | Z=R*SB;
|
|---|
| 547 | XD=-Y*DEL-W*CEL;
|
|---|
| 548 | YD=X*DEL-W*SEL;
|
|---|
| 549 | ZD=RBD*CB+SB*DR;
|
|---|
| 550 |
|
|---|
| 551 | /* Julian centuries since J2000 */
|
|---|
| 552 | T=(date-51544.5)/36525.;
|
|---|
| 553 |
|
|---|
| 554 | /* Fricke equinox correction */
|
|---|
| 555 | EPJ=2000.+T*100.;
|
|---|
| 556 | EQCOR=PAL__DS2R*(0.035+0.00085*(EPJ-B1950));
|
|---|
| 557 |
|
|---|
| 558 | /* Mean obliquity (IAU 1976) */
|
|---|
| 559 | EPS=PAL__DAS2R*(84381.448+(-46.8150+(-0.00059+0.001813*T)*T)*T);
|
|---|
| 560 |
|
|---|
| 561 | /* To the equatorial system, mean of date, FK5 system */
|
|---|
| 562 | SINEPS=sin(EPS);
|
|---|
| 563 | COSEPS=cos(EPS);
|
|---|
| 564 | ES=EQCOR*SINEPS;
|
|---|
| 565 | EC=EQCOR*COSEPS;
|
|---|
| 566 | pv[0]=X-EC*Y+ES*Z;
|
|---|
| 567 | pv[1]=EQCOR*X+Y*COSEPS-Z*SINEPS;
|
|---|
| 568 | pv[2]=Y*SINEPS+Z*COSEPS;
|
|---|
| 569 | pv[3]=XD-EC*YD+ES*ZD;
|
|---|
| 570 | pv[4]=EQCOR*XD+YD*COSEPS-ZD*SINEPS;
|
|---|
| 571 | pv[5]=YD*SINEPS+ZD*COSEPS;
|
|---|
| 572 |
|
|---|
| 573 | }
|
|---|
