Astro-PAL
view release on metacpan or search on metacpan
palsrc/palDmoon.c view on Meta::CPAN
/*
*+
* Name:
* palDmoon
* Purpose:
* Approximate geocentric position and velocity of the Moon
* Language:
* Starlink ANSI C
* Type of Module:
* Library routine
* Invocation:
* void palDmoon( double date, double pv[6] );
* Arguments:
* date = double (Given)
* TDB as a Modified Julian Date (JD-2400000.5)
* pv = double [6] (Returned)
* Moon x,y,z,xdot,ydot,zdot, mean equator and
* equinox of date (AU, AU/s)
* Description:
* Calculate the approximate geocentric position of the Moon
* using a full implementation of the algorithm published by
* Meeus (l'Astronomie, June 1984, p348).
* Authors:
* TIMJ: Tim Jenness (JAC, Hawaii)
* PTW: Patrick T. Wallace
* {enter_new_authors_here}
* Notes:
* - Meeus quotes accuracies of 10 arcsec in longitude, 3 arcsec in
* latitude and 0.2 arcsec in HP (equivalent to about 20 km in
* distance). Comparison with JPL DE200 over the interval
* 1960-2025 gives RMS errors of 3.7 arcsec and 83 mas/hour in
* longitude, 2.3 arcsec and 48 mas/hour in latitude, 11 km
* and 81 mm/s in distance. The maximum errors over the same
* interval are 18 arcsec and 0.50 arcsec/hour in longitude,
* 11 arcsec and 0.24 arcsec/hour in latitude, 40 km and 0.29 m/s
* in distance.
* - The original algorithm is expressed in terms of the obsolete
* timescale Ephemeris Time. Either TDB or TT can be used, but
* not UT without incurring significant errors (30 arcsec at
* the present time) due to the Moon's 0.5 arcsec/sec movement.
* - The algorithm is based on pre IAU 1976 standards. However,
* the result has been moved onto the new (FK5) equinox, an
* adjustment which is in any case much smaller than the
* intrinsic accuracy of the procedure.
* - Velocity is obtained by a complete analytical differentiation
* of the Meeus model.
* History:
* 2012-03-07 (TIMJ):
* Initial version based on a direct port of the SLA/F code.
* Adapted with permission from the Fortran SLALIB library.
* {enter_further_changes_here}
* Copyright:
* Copyright (C) 1998 Rutherford Appleton Laboratory
* Copyright (C) 2012 Science and Technology Facilities Council.
* All Rights Reserved.
* Licence:
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation; either version 3 of
* the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be
* useful, but WITHOUT ANY WARRANTY; without even the implied
* warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
* PURPOSE. See the GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
* MA 02110-1301, USA.
* Bugs:
* {note_any_bugs_here}
*-
*/
#include "pal.h"
#include "pal1sofa.h"
#include "palmac.h"
/* Autoconf can give us -DPIC */
#undef PIC
void palDmoon( double date, double pv[6] ) {
/* Seconds per Julian century (86400*36525) */
const double CJ = 3155760000.0;
/* Julian epoch of B1950 */
const double B1950 = 1949.9997904423;
/* Earth equatorial radius in AU ( = 6378.137 / 149597870 ) */
const double ERADAU=4.2635212653763e-5;
double T,THETA,SINOM,COSOM,DOMCOM,WA,DWA,WB,DWB,WOM,
DWOM,SINWOM,COSWOM,V,DV,COEFF,EMN,EMPN,DN,FN,EN,
DEN,DTHETA,FTHETA,EL,DEL,B,DB,BF,DBF,P,DP,SP,R,
DR,X,Y,Z,XD,YD,ZD,SEL,CEL,SB,CB,RCB,RBD,W,EPJ,
EQCOR,EPS,SINEPS,COSEPS,ES,EC;
double ELP, DELP;
double EM, DEM, EMP, DEMP, D, DD, F, DF, OM, DOM, E, DESQ, ESQ, DE;
int N,I;
/*
* Coefficients for fundamental arguments
*
* at J1900: T**0, T**1, T**2, T**3
* at epoch: T**0, T**1
*
* Units are degrees for position and Julian centuries for time
*
*/
/* Moon's mean longitude */
const double ELP0=270.434164;
const double ELP1=481267.8831;
const double ELP2=-0.001133;
const double ELP3=0.0000019;
/* Sun's mean anomaly */
const double EM0=358.475833;
const double EM1=35999.0498;
const double EM2=-0.000150;
const double EM3=-0.0000033;
/* Moon's mean anomaly */
const double EMP0=296.104608;
const double EMP1=477198.8491;
const double EMP2=0.009192;
const double EMP3=0.0000144;
/* Moon's mean elongation */
const double D0=350.737486;
const double D1=445267.1142;
const double D2=-0.001436;
const double D3=0.0000019;
/* Mean distance of the Moon from its ascending node */
const double F0=11.250889;
const double F1=483202.0251;
const double F2=-0.003211;
const double F3=-0.0000003;
/* Longitude of the Moon's ascending node */
const double OM0=259.183275;
const double OM1=-1934.1420;
const double OM2=0.002078;
const double OM3=0.0000022;
/* Coefficients for (dimensionless) E factor */
const double E1=-0.002495;
const double E2=-0.00000752;
/* Coefficients for periodic variations etc */
const double PAC=0.000233;
const double PA0=51.2;
const double PA1=20.2;
const double PBC=-0.001778;
const double PCC=0.000817;
const double PDC=0.002011;
const double PEC=0.003964;
const double PE0=346.560;
const double PE1=132.870;
const double PE2=-0.0091731;
const double PFC=0.001964;
const double PGC=0.002541;
const double PHC=0.001964;
const double PIC=-0.024691;
const double PJC=-0.004328;
const double PJ0=275.05;
const double PJ1=-2.30;
const double CW1=0.0004664;
const double CW2=0.0000754;
/*
* Coefficients for Moon position
*
* Tx(N) = coefficient of L, B or P term (deg)
* ITx(N,1-5) = coefficients of M, M', D, F, E**n in argument
*/
#define NL 50
#define NB 45
#define NP 31
/*
* Longitude
*/
const double TL[NL] = {
6.28875,1.274018,.658309,.213616,-.185596,
-.114336,.058793,.057212,.05332,.045874,.041024,-.034718,-.030465,
.015326,-.012528,-.01098,.010674,.010034,.008548,-.00791,-.006783,
.005162,.005,.004049,.003996,.003862,.003665,.002695,.002602,
.002396,-.002349,.002249,-.002125,-.002079,.002059,-.001773,
-.001595,.00122,-.00111,8.92e-4,-8.11e-4,7.61e-4,7.17e-4,7.04e-4,
6.93e-4,5.98e-4,5.5e-4,5.38e-4,5.21e-4,4.86e-4
};
const int ITL[NL][5] = {
/* M M' D F n */
{ +0, +1, +0, +0, 0 },
{ +0, -1, +2, +0, 0 },
{ +0, +0, +2, +0, 0 },
{ +0, +2, +0, +0, 0 },
{ +1, +0, +0, +0, 1 },
{ +0, +0, +0, +2, 0 },
{ +0, -2, +2, +0, 0 },
{ -1, -1, +2, +0, 1 },
{ +0, +1, +2, +0, 0 },
{ -1, +0, +2, +0, 1 },
{ -1, +1, +0, +0, 1 },
{ +0, +0, +1, +0, 0 },
{ +1, +1, +0, +0, 1 },
{ +0, +0, +2, -2, 0 },
{ +0, +1, +0, +2, 0 },
{ +0, -1, +0, +2, 0 },
{ +0, -1, +4, +0, 0 },
{ +0, +3, +0, +0, 0 },
{ +0, -2, +4, +0, 0 },
{ +1, -1, +2, +0, 1 },
{ +1, +0, +2, +0, 1 },
{ +0, +1, -1, +0, 0 },
{ +1, +0, +1, +0, 1 },
{ -1, +1, +2, +0, 1 },
{ +0, +2, +2, +0, 0 },
{ +0, +0, +4, +0, 0 },
{ +0, -3, +2, +0, 0 },
{ -1, +2, +0, +0, 1 },
{ +0, +1, -2, -2, 0 },
{ -1, -2, +2, +0, 1 },
{ +0, +1, +1, +0, 0 },
{ -2, +0, +2, +0, 2 },
{ +1, +2, +0, +0, 1 },
{ +2, +0, +0, +0, 2 },
{ -2, -1, +2, +0, 2 },
{ +0, +1, +2, -2, 0 },
{ +0, +0, +2, +2, 0 },
{ -1, -1, +4, +0, 1 },
{ +0, +2, +0, +2, 0 },
{ +0, +1, -3, +0, 0 },
{ +1, +1, +2, +0, 1 },
{ -1, -2, +4, +0, 1 },
{ -2, +1, +0, +0, 2 },
{ -2, +1, -2, +0, 2 },
{ +1, -2, +2, +0, 1 },
{ -1, +0, +2, -2, 1 },
{ +0, +1, +4, +0, 0 },
{ +0, +4, +0, +0, 0 },
{ -1, +0, +4, +0, 1 },
{ +0, +2, -1, +0, 0 }
};
/*
* Latitude
*/
const double TB[NB] = {
5.128189,.280606,.277693,.173238,.055413,
.046272,.032573,.017198,.009267,.008823,.008247,.004323,.0042,
.003372,.002472,.002222,.002072,.001877,.001828,-.001803,-.00175,
.00157,-.001487,-.001481,.001417,.00135,.00133,.001106,.00102,
8.33e-4,7.81e-4,6.7e-4,6.06e-4,5.97e-4,4.92e-4,4.5e-4,4.39e-4,
4.23e-4,4.22e-4,-3.67e-4,-3.53e-4,3.31e-4,3.17e-4,3.06e-4,
-2.83e-4
};
const int ITB[NB][5] = {
/* M M' D F n */
{ +0, +0, +0, +1, 0 },
{ +0, +1, +0, +1, 0 },
{ +0, +1, +0, -1, 0 },
{ +0, +0, +2, -1, 0 },
{ +0, -1, +2, +1, 0 },
{ +0, -1, +2, -1, 0 },
{ +0, +0, +2, +1, 0 },
{ +0, +2, +0, +1, 0 },
{ +0, +1, +2, -1, 0 },
{ +0, +2, +0, -1, 0 },
{ -1, +0, +2, -1, 1 },
{ +0, -2, +2, -1, 0 },
{ +0, +1, +2, +1, 0 },
{ -1, +0, -2, +1, 1 },
{ -1, -1, +2, +1, 1 },
{ -1, +0, +2, +1, 1 },
{ -1, -1, +2, -1, 1 },
{ -1, +1, +0, +1, 1 },
{ +0, -1, +4, -1, 0 },
{ +1, +0, +0, +1, 1 },
{ +0, +0, +0, +3, 0 },
{ -1, +1, +0, -1, 1 },
{ +0, +0, +1, +1, 0 },
{ +1, +1, +0, +1, 1 },
{ -1, -1, +0, +1, 1 },
{ -1, +0, +0, +1, 1 },
{ +0, +0, -1, +1, 0 },
{ +0, +3, +0, +1, 0 },
{ +0, +0, +4, -1, 0 },
{ +0, -1, +4, +1, 0 },
{ +0, +1, +0, -3, 0 },
{ +0, -2, +4, +1, 0 },
{ +0, +0, +2, -3, 0 },
{ +0, +2, +2, -1, 0 },
{ -1, +1, +2, -1, 1 },
{ +0, +2, -2, -1, 0 },
{ +0, +3, +0, -1, 0 },
{ +0, +2, +2, +1, 0 },
{ +0, -3, +2, -1, 0 },
{ +1, -1, +2, +1, 1 },
{ +1, +0, +2, +1, 1 },
{ +0, +0, +4, +1, 0 },
{ -1, +1, +2, +1, 1 },
{ -2, +0, +2, -1, 2 },
{ +0, +1, +0, +3, 0 }
};
/*
* Parallax
*/
const double TP[NP] = {
.950724,.051818,.009531,.007843,.002824,
8.57e-4,5.33e-4,4.01e-4,3.2e-4,-2.71e-4,-2.64e-4,-1.98e-4,1.73e-4,
1.67e-4,-1.11e-4,1.03e-4,-8.4e-5,-8.3e-5,7.9e-5,7.2e-5,6.4e-5,
-6.3e-5,4.1e-5,3.5e-5,-3.3e-5,-3e-5,-2.9e-5,-2.9e-5,2.6e-5,
-2.3e-5,1.9e-5
};
const int ITP[NP][5] = {
/* M M' D F n */
{ +0, +0, +0, +0, 0 },
{ +0, +1, +0, +0, 0 },
{ +0, -1, +2, +0, 0 },
{ +0, +0, +2, +0, 0 },
{ +0, +2, +0, +0, 0 },
{ +0, +1, +2, +0, 0 },
{ -1, +0, +2, +0, 1 },
{ -1, -1, +2, +0, 1 },
{ -1, +1, +0, +0, 1 },
{ +0, +0, +1, +0, 0 },
{ +1, +1, +0, +0, 1 },
{ +0, -1, +0, +2, 0 },
{ +0, +3, +0, +0, 0 },
{ +0, -1, +4, +0, 0 },
{ +1, +0, +0, +0, 1 },
{ +0, -2, +4, +0, 0 },
{ +0, +2, -2, +0, 0 },
{ +1, +0, +2, +0, 1 },
{ +0, +2, +2, +0, 0 },
{ +0, +0, +4, +0, 0 },
{ -1, +1, +2, +0, 1 },
{ +1, -1, +2, +0, 1 },
{ +1, +0, +1, +0, 1 },
{ -1, +2, +0, +0, 1 },
{ +0, +3, -2, +0, 0 },
{ +0, +1, +1, +0, 0 },
{ +0, +0, -2, +2, 0 },
{ +1, +2, +0, +0, 1 },
{ -2, +0, +2, +0, 2 },
{ +0, +1, -2, +2, 0 },
{ -1, -1, +4, +0, 1 }
};
/* Centuries since J1900 */
T=(date-15019.5)/36525.;
/*
* Fundamental arguments (radians) and derivatives (radians per
* Julian century) for the current epoch
*/
/* Moon's mean longitude */
ELP=PAL__DD2R*fmod(ELP0+(ELP1+(ELP2+ELP3*T)*T)*T,360.);
DELP=PAL__DD2R*(ELP1+(2.*ELP2+3*ELP3*T)*T);
/* Sun's mean anomaly */
EM=PAL__DD2R*fmod(EM0+(EM1+(EM2+EM3*T)*T)*T,360.);
DEM=PAL__DD2R*(EM1+(2.*EM2+3*EM3*T)*T);
/* Moon's mean anomaly */
EMP=PAL__DD2R*fmod(EMP0+(EMP1+(EMP2+EMP3*T)*T)*T,360.);
DEMP=PAL__DD2R*(EMP1+(2.*EMP2+3*EMP3*T)*T);
( run in 1.619 second using v1.01-cache-2.11-cpan-8f98c5d2c55 )