GPS ToolKit Software Library: MoonPosition.cpp Source File

00001 #pragma ident "$Id: MoonPosition.cpp 1161 2008-03-27 17:16:22Z ckiesch $"
00002 
00009 //============================================================================
00010 //
00011 //  This file is part of GPSTk, the GPS Toolkit.
00012 //
00013 //  The GPSTk is free software; you can redistribute it and/or modify
00014 //  it under the terms of the GNU Lesser General Public License as published
00015 //  by the Free Software Foundation; either version 2.1 of the License, or
00016 //  any later version.
00017 //
00018 //  The GPSTk is distributed in the hope that it will be useful,
00019 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00020 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00021 //  GNU Lesser General Public License for more details.
00022 //
00023 //  You should have received a copy of the GNU Lesser General Public
00024 //  License along with GPSTk; if not, write to the Free Software Foundation,
00025 //  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00026 //  
00027 //  Dagoberto Salazar - gAGE ( http://www.gage.es ). 2007
00028 //
00029 //============================================================================
00030 
00031 
00032 #include "MoonPosition.hpp"
00033 
00034 
00035 namespace gpstk
00036 {
00037 
00038 
00039       // Time of the first valid time
00040    const DayTime MoonPosition::initialTime(1900, 3, 1, 0, 0, 0.0);
00041 
00042       // Time of the last valid time
00043    const DayTime MoonPosition::finalTime(2100, 2, 28, 0, 0, 0.0);
00044 
00045 
00046       // Coefficients for fundamental arguments
00047       // Units are degrees for position and Julian centuries for time
00048 
00049       // Moon's mean longitude
00050    const double MoonPosition::ELP0(270.434164),
00051                 MoonPosition::ELP1(481267.8831),
00052                 MoonPosition::ELP2(-0.001133),
00053                 MoonPosition::ELP3(0.0000019);
00054 
00055       // Sun's mean anomaly
00056    const double MoonPosition::EM0(358.475833),
00057                 MoonPosition::EM1(35999.0498),
00058                 MoonPosition::EM2(-0.000150),
00059                 MoonPosition::EM3(-0.0000033);
00060 
00061       // Moon's mean anomaly
00062    const double MoonPosition::EMP0(296.104608),
00063                 MoonPosition::EMP1(477198.8491),
00064                 MoonPosition::EMP2(0.009192),
00065                 MoonPosition::EMP3(0.0000144);
00066 
00067       // Moon's mean elongation
00068    const double MoonPosition::D0(350.737486),
00069                 MoonPosition::D1(445267.1142),
00070                 MoonPosition::D2(-0.001436),
00071                 MoonPosition::D3(0.0000019);
00072 
00073       // Mean distance of the Moon from its ascending node
00074    const double MoonPosition::F0(11.250889),
00075                 MoonPosition::F1(483202.0251),
00076                 MoonPosition::F2(-0.003211),
00077                 MoonPosition::F3(-0.0000003);
00078 
00079       // Longitude of the Moon's ascending node
00080    const double MoonPosition::OM0(259.183275),
00081                 MoonPosition::OM1(-1934.1420),
00082                 MoonPosition::OM2(0.002078),
00083                 MoonPosition::OM3(0.0000022);
00084 
00085 
00086       // Coefficients for (dimensionless) E factor
00087    const double MoonPosition::E1(-0.002495),
00088                 MoonPosition::E2(-0.00000752);
00089 
00090       // Coefficients for periodic variations, etc
00091    const double MoonPosition::PAC(0.000233),
00092                 MoonPosition::PA0(51.2),
00093                 MoonPosition::PA1(20.2);
00094    const double MoonPosition::PBC(-0.001778);
00095    const double MoonPosition::PCC(0.000817);
00096    const double MoonPosition::PDC(0.002011);
00097    const double MoonPosition::PEC(0.003964),
00098                 MoonPosition::PE0(346.560),
00099                 MoonPosition::PE1(132.870),
00100                 MoonPosition::PE2(-0.0091731);
00101    const double MoonPosition::PFC(0.001964);
00102    const double MoonPosition::PGC(0.002541);
00103    const double MoonPosition::PHC(0.001964);
00104    const double MoonPosition::cPIC(-0.024691);
00105    const double MoonPosition::PJC(-0.004328),
00106                 MoonPosition::PJ0(275.05),
00107                 MoonPosition::PJ1(-2.30);
00108    const double MoonPosition::CW1(0.0004664);
00109    const double MoonPosition::CW2(0.0000754);
00110 
00111       // Coefficients for Moon position
00112       //      Tx(N): coefficient of L, B or P term (deg)
00113       //      ITx(N,0-4): coefficients of M, M', D, F, E**n in argument
00114       //
00115    const size_t MoonPosition::NL(50),
00116                 MoonPosition::NB(45),
00117                 MoonPosition::NP(31);
00118    Vector<double> MoonPosition::TL(NL,0.0),
00119                   MoonPosition::TB(NB,0.0),
00120                   MoonPosition::TP(NP,0.0);
00121    Matrix<int> MoonPosition::ITL(5,NL,0),
00122                MoonPosition::ITB(5,NB,0),
00123                MoonPosition::ITP(5,NP,0);
00124 
00125 
00126 
00127       // Returns the position of Moon ECEF coordinates (meters) at 
00128       // the indicated time.
00129       //
00130       // @param[in] t   the time to look up
00131       //
00132       // @return the position of the Moon at time (as a Triple)
00133       //
00134       // @throw InvalidRequest If the request can not be completed for 
00135       // any reason, this is thrown. The text may have additional
00136       // information as to why the request failed.
00137    Triple MoonPosition::getPosition(const DayTime& t) const
00138       throw(InvalidRequest)
00139    {
00140 
00141          // Test if the time interval is correct
00142       if ( (t < MoonPosition::initialTime) ||
00143            (t > MoonPosition::finalTime) )
00144       {
00145          InvalidRequest ir("Provided epoch is out of bounds.");
00146          GPSTK_THROW(ir);
00147       }
00148 
00149 
00150          // We will store here the results
00151       Triple res;
00152 
00153       res = MoonPosition::getPositionCIS(t);
00154       res = CIS2CTS(res, t);
00155 
00156       return res;
00157 
00158    } // End MoonPosition::getPosition
00159 
00160 
00161 
00162       /* Function to compute Moon position in CIS system (coordinates 
00163        * in meters)
00164        *
00165        * @param t Epoch
00166        */
00167    Triple MoonPosition::getPositionCIS(const DayTime& t) const
00168       throw(InvalidRequest)
00169    {
00170 
00171          // Test if the time interval is correct
00172       if ( (t < MoonPosition::initialTime) ||
00173            (t > MoonPosition::finalTime) )
00174       {
00175          InvalidRequest ir("Provided epoch is out of bounds.");
00176          GPSTK_THROW(ir);
00177       }
00178 
00179 
00180 
00181          // Coefficients for Moon position
00182          //      Tx(N): coefficient of L, B or P term (deg)
00183          //      ITx(N,0-4): coefficients of M, M', D, F, E**n in argument
00184          //
00185 
00186          // Longitude
00187       TL( 0) = +6.288750;
00188       //          M              M'             D              F             n
00189 
00190       ITL(0, 0)= +0; ITL(1, 0)= +1; ITL(2, 0)= +0; ITL(3, 0)= +0, ITL(4, 0)= 0;
00191       TL( 1) = +1.274018;
00192       ITL(0, 1)= +0; ITL(1, 1)= -1; ITL(2, 1)= +2; ITL(3, 1)= +0, ITL(4, 1)= 0;
00193       TL( 2) = +0.658309;
00194       ITL(0, 2)= +0; ITL(1, 2)= +0; ITL(2, 2)= +2; ITL(3, 2)= +0, ITL(4, 2)= 0;
00195       TL( 3) = +0.213616;
00196       ITL(0, 3)= +0; ITL(1, 3)= +2; ITL(2, 3)= +0; ITL(3, 3)= +0, ITL(4, 3)= 0;
00197       TL( 4) = -0.185596;
00198       ITL(0, 4)= +1; ITL(1, 4)= +0; ITL(2, 4)= +0; ITL(3, 4)= +0, ITL(4, 4)= 1;
00199       TL( 5) = -0.114336;
00200       ITL(0, 5)= +0; ITL(1, 5)= +0; ITL(2, 5)= +0; ITL(3, 5)= +2, ITL(4, 5)= 0;
00201       TL( 6) = +0.058793;
00202       ITL(0, 6)= +0; ITL(1, 6)= -2; ITL(2, 6)= +2; ITL(3, 6)= +0, ITL(4, 6)= 0;
00203       TL( 7) = +0.057212;
00204       ITL(0, 7)= -1; ITL(1, 7)= -1; ITL(2, 7)= +2; ITL(3, 7)= +0, ITL(4, 7)= 1;
00205       TL( 8) = +0.053320;
00206       ITL(0, 8)= +0; ITL(1, 8)= +1; ITL(2, 8)= +2; ITL(3, 8)= +0, ITL(4, 8)= 0;
00207       TL( 9) = +0.045874;
00208       ITL(0, 9)= -1; ITL(1, 9)= +0; ITL(2, 9)= +2; ITL(3, 9)= +0, ITL(4, 9)= 1;
00209       TL(10) = +0.041024;
00210       ITL(0,10)= -1; ITL(1,10)= +1; ITL(2,10)= +0; ITL(3,10)= +0, ITL(4,10)= 1;
00211       TL(11) = -0.034718;
00212       ITL(0,11)= +0; ITL(1,11)= +0; ITL(2,11)= +1; ITL(3,11)= +0, ITL(4,11)= 0;
00213       TL(12) = -0.030465;
00214       ITL(0,12)= +1; ITL(1,12)= +1; ITL(2,12)= +0; ITL(3,12)= +0, ITL(4,12)= 1;
00215       TL(13) = +0.015326;
00216       ITL(0,13)= +0; ITL(1,13)= +0; ITL(2,13)= +2; ITL(3,13)= -2, ITL(4,13)= 0;
00217       TL(14) = -0.012528;
00218       ITL(0,14)= +0; ITL(1,14)= +1; ITL(2,14)= +0; ITL(3,14)= +2, ITL(4,14)= 0;
00219       TL(15) = -0.010980;
00220       ITL(0,15)= +0; ITL(1,15)= -1; ITL(2,15)= +0; ITL(3,15)= +2, ITL(4,15)= 0;
00221       TL(16) = +0.010674;
00222       ITL(0,16)= +0; ITL(1,16)= -1; ITL(2,16)= +4; ITL(3,16)= +0, ITL(4,16)= 0;
00223       TL(17) = +0.010034;
00224       ITL(0,17)= +0; ITL(1,17)= +3; ITL(2,17)= +0; ITL(3,17)= +0, ITL(4,17)= 0;
00225       TL(18) = +0.008548;
00226       ITL(0,18)= +0; ITL(1,18)= -2; ITL(2,18)= +4; ITL(3,18)= +0, ITL(4,18)= 0;
00227       TL(19) = -0.007910;
00228       ITL(0,19)= +1; ITL(1,19)= -1; ITL(2,19)= +2; ITL(3,19)= +0, ITL(4,19)= 1;
00229       TL(20) = -0.006783;
00230       ITL(0,20)= +1; ITL(1,20)= +0; ITL(2,20)= +2; ITL(3,20)= +0, ITL(4,20)= 1;
00231       TL(21) = +0.005162;
00232       ITL(0,21)= +0; ITL(1,21)= +1; ITL(2,21)= -1; ITL(3,21)= +0, ITL(4,21)= 0;
00233       TL(22) = +0.005000;
00234       ITL(0,22)= +1; ITL(1,22)= +0; ITL(2,22)= +1; ITL(3,22)= +0, ITL(4,22)= 1;
00235       TL(23) = +0.004049;
00236       ITL(0,23)= -1; ITL(1,23)= +1; ITL(2,23)= +2; ITL(3,23)= +0, ITL(4,23)= 1;
00237       TL(24) = +0.003996;
00238       ITL(0,24)= +0; ITL(1,24)= +2; ITL(2,24)= +2; ITL(3,24)= +0, ITL(4,24)= 0;
00239       TL(25) = +0.003862;
00240       ITL(0,25)= +0; ITL(1,25)= +0; ITL(2,25)= +4; ITL(3,25)= +0, ITL(4,25)= 0;
00241       TL(26) = +0.003665;
00242       ITL(0,26)= +0; ITL(1,26)= -3; ITL(2,26)= +2; ITL(3,26)= +0, ITL(4,26)= 0;
00243       TL(27) = +0.002695;
00244       ITL(0,27)= -1; ITL(1,27)= +2; ITL(2,27)= +0; ITL(3,27)= +0, ITL(4,27)= 1;
00245       TL(28) = +0.002602;
00246       ITL(0,28)= +0; ITL(1,28)= +1; ITL(2,28)= -2; ITL(3,28)= -2, ITL(4,28)= 0;
00247       TL(29) = +0.002396;
00248       ITL(0,29)= -1; ITL(1,29)= -2; ITL(2,29)= +2; ITL(3,29)= +0, ITL(4,29)= 1;
00249       TL(30) = -0.002349;
00250       ITL(0,30)= +0; ITL(1,30)= +1; ITL(2,30)= +1; ITL(3,30)= +0, ITL(4,30)= 0;
00251       TL(31) = +0.002249;
00252       ITL(0,31)= -2; ITL(1,31)= +0; ITL(2,31)= +2; ITL(3,31)= +0, ITL(4,31)= 2;
00253       TL(32) = -0.002125;
00254       ITL(0,32)= +1; ITL(1,32)= +2; ITL(2,32)= +0; ITL(3,32)= +0, ITL(4,32)= 1;
00255       TL(33) = -0.002079;
00256       ITL(0,33)= +2; ITL(1,33)= +0; ITL(2,33)= +0; ITL(3,33)= +0, ITL(4,33)= 2;
00257       TL(34) = +0.002059;
00258       ITL(0,34)= -2; ITL(1,34)= -1; ITL(2,34)= +2; ITL(3,34)= +0, ITL(4,34)= 2;
00259       TL(35) = -0.001773;
00260       ITL(0,35)= +0; ITL(1,35)= +1; ITL(2,35)= +2; ITL(3,35)= -2, ITL(4,35)= 0;
00261       TL(36) = -0.001595;
00262       ITL(0,36)= +0; ITL(1,36)= +0; ITL(2,36)= +2; ITL(3,36)= +2, ITL(4,36)= 0;
00263       TL(37) = +0.001220;
00264       ITL(0,37)= -1; ITL(1,37)= -1; ITL(2,37)= +4; ITL(3,37)= +0, ITL(4,37)= 1;
00265       TL(38) = -0.001110;
00266       ITL(0,38)= +0; ITL(1,38)= +2; ITL(2,38)= +0; ITL(3,38)= +2, ITL(4,38)= 0;
00267       TL(39) = +0.000892;
00268       ITL(0,39)= +0; ITL(1,39)= +1; ITL(2,39)= -3; ITL(3,39)= +0, ITL(4,39)= 0;
00269       TL(40) = -0.000811;
00270       ITL(0,40)= +1; ITL(1,40)= +1; ITL(2,40)= +2; ITL(3,40)= +0, ITL(4,40)= 1;
00271       TL(41) = +0.000761;
00272       ITL(0,41)= -1; ITL(1,41)= -2; ITL(2,41)= +4; ITL(3,41)= +0, ITL(4,41)= 1;
00273       TL(42) = +0.000717;
00274       ITL(0,42)= -2; ITL(1,42)= +1; ITL(2,42)= +0; ITL(3,42)= +0, ITL(4,42)= 2;
00275       TL(43) = +0.000704;
00276       ITL(0,43)= -2; ITL(1,43)= +1; ITL(2,43)= -2; ITL(3,43)= +0, ITL(4,43)= 2;
00277       TL(44) = +0.000693;
00278       ITL(0,44)= +1; ITL(1,44)= -2; ITL(2,44)= +2; ITL(3,44)= +0, ITL(4,44)= 1;
00279       TL(45) = +0.000598;
00280       ITL(0,45)= -1; ITL(1,45)= +0; ITL(2,45)= +2; ITL(3,45)= -2, ITL(4,45)= 1;
00281       TL(46) = +0.000550;
00282       ITL(0,46)= +0; ITL(1,46)= +1; ITL(2,46)= +4; ITL(3,46)= +0, ITL(4,46)= 0;
00283       TL(47) = +0.000538;
00284       ITL(0,47)= +0; ITL(1,47)= +4; ITL(2,47)= +0; ITL(3,47)= +0, ITL(4,47)= 0;
00285       TL(48) = +0.000521;
00286       ITL(0,48)= -1; ITL(1,48)= +0; ITL(2,48)= +4; ITL(3,48)= +0, ITL(4,48)= 1;
00287       TL(49) = +0.000486;
00288       ITL(0,49)= +0; ITL(1,49)= +2; ITL(2,49)= -1; ITL(3,49)= +0, ITL(4,49)= 0;
00289 
00290 
00291          // Latitude
00292       TB( 0) = +5.128189;
00293       //          M              M'             D              F             n
00294       ITB(0, 0)= +0; ITB(1, 0)= +0; ITB(2, 0)= +0; ITB(3, 0)= +1, ITB(4, 0)= 0;
00295       TB( 1) = +0.280606;
00296       ITB(0, 1)= +0; ITB(1, 1)= +1; ITB(2, 1)= +0; ITB(3, 1)= +1, ITB(4, 1)= 0;
00297       TB( 2) = +0.277693;
00298       ITB(0, 2)= +0; ITB(1, 2)= +1; ITB(2, 2)= +0; ITB(3, 2)= -1, ITB(4, 2)= 0;
00299       TB( 3) = +0.173238;
00300       ITB(0, 3)= +0; ITB(1, 3)= +0; ITB(2, 3)= +2; ITB(3, 3)= -1, ITB(4, 3)= 0;
00301       TB( 4) = +0.055413;
00302       ITB(0, 4)= +0; ITB(1, 4)= -1; ITB(2, 4)= +2; ITB(3, 4)= +1, ITB(4, 4)= 0;
00303       TB( 5) = +0.046272;
00304       ITB(0, 5)= +0; ITB(1, 5)= -1; ITB(2, 5)= +2; ITB(3, 5)= -1, ITB(4, 5)= 0;
00305       TB( 6) = +0.032573;
00306       ITB(0, 6)= +0; ITB(1, 6)= +0; ITB(2, 6)= +2; ITB(3, 6)= +1, ITB(4, 6)= 0;
00307       TB( 7) = +0.017198;
00308       ITB(0, 7)= +0; ITB(1, 7)= +2; ITB(2, 7)= +0; ITB(3, 7)= +1, ITB(4, 7)= 0;
00309       TB( 8) = +0.009267;
00310       ITB(0, 8)= +0; ITB(1, 8)= +1; ITB(2, 8)= +2; ITB(3, 8)= -1, ITB(4, 8)= 0;
00311       TB( 9) = +0.008823;
00312       ITB(0, 9)= +0; ITB(1, 9)= +2; ITB(2, 9)= +0; ITB(3, 9)= -1, ITB(4, 9)= 0;
00313       TB(10) = +0.008247;
00314       ITB(0,10)= -1; ITB(1,10)= +0; ITB(2,10)= +2; ITB(3,10)= -1, ITB(4,10)= 1;
00315       TB(11) = +0.004323;
00316       ITB(0,11)= +0; ITB(1,11)= -2; ITB(2,11)= +2; ITB(3,11)= -1, ITB(4,11)= 0;
00317       TB(12) = +0.004200;
00318       ITB(0,12)= +0; ITB(1,12)= +1; ITB(2,12)= +2; ITB(3,12)= +1, ITB(4,12)= 0;
00319       TB(13) = +0.003372;
00320       ITB(0,13)= -1; ITB(1,13)= +0; ITB(2,13)= -2; ITB(3,13)= +1, ITB(4,13)= 1;
00321       TB(14) = +0.002472;
00322       ITB(0,14)= -1; ITB(1,14)= -1; ITB(2,14)= +2; ITB(3,14)= +1, ITB(4,14)= 1;
00323       TB(15) = +0.002222;
00324       ITB(0,15)= -1; ITB(1,15)= +0; ITB(2,15)= +2; ITB(3,15)= +1, ITB(4,15)= 1;
00325       TB(16) = +0.002072;
00326       ITB(0,16)= -1; ITB(1,16)= -1; ITB(2,16)= +2; ITB(3,16)= -1, ITB(4,16)= 1;
00327       TB(17) = +0.001877;
00328       ITB(0,17)= -1; ITB(1,17)= +1; ITB(2,17)= +0; ITB(3,17)= +1, ITB(4,17)= 1;
00329       TB(18) = +0.001828;
00330       ITB(0,18)= +0; ITB(1,18)= -1; ITB(2,18)= +4; ITB(3,18)= -1, ITB(4,18)= 0;
00331       TB(19) = -0.001803;
00332       ITB(0,19)= +1; ITB(1,19)= +0; ITB(2,19)= +0; ITB(3,19)= +1, ITB(4,19)= 1;
00333       TB(20) = -0.001750;
00334       ITB(0,20)= +0; ITB(1,20)= +0; ITB(2,20)= +0; ITB(3,20)= +3, ITB(4,20)= 0;
00335       TB(21) = +0.001570;
00336       ITB(0,21)= -1; ITB(1,21)= +1; ITB(2,21)= +0; ITB(3,21)= -1, ITB(4,21)= 1;
00337       TB(22) = -0.001487;
00338       ITB(0,22)= +0; ITB(1,22)= +0; ITB(2,22)= +1; ITB(3,22)= +1, ITB(4,22)= 0;
00339       TB(23) = -0.001481;
00340       ITB(0,23)= +1; ITB(1,23)= +1; ITB(2,23)= +0; ITB(3,23)= +1, ITB(4,23)= 1;
00341       TB(24) = +0.001417;
00342       ITB(0,24)= -1; ITB(1,24)= -1; ITB(2,24)= +0; ITB(3,24)= +1, ITB(4,24)= 1;
00343       TB(25) = +0.001350;
00344       ITB(0,25)= -1; ITB(1,25)= +0; ITB(2,25)= +0; ITB(3,25)= +1, ITB(4,25)= 1;
00345       TB(26) = +0.001330;
00346       ITB(0,26)= +0; ITB(1,26)= +0; ITB(2,26)= -1; ITB(3,26)= +1, ITB(4,26)= 0;
00347       TB(27) = +0.001106;
00348       ITB(0,27)= +0; ITB(1,27)= +3; ITB(2,27)= +0; ITB(3,27)= +1, ITB(4,27)= 0;
00349       TB(28) = +0.001020;
00350       ITB(0,28)= +0; ITB(1,28)= +0; ITB(2,28)= +4; ITB(3,28)= -1, ITB(4,28)= 0;
00351       TB(29) = +0.000833;
00352       ITB(0,29)= +0; ITB(1,29)= -1; ITB(2,29)= +4; ITB(3,29)= +1, ITB(4,29)= 0;
00353       TB(30) = +0.000781;
00354       ITB(0,30)= +0; ITB(1,30)= +1; ITB(2,30)= +0; ITB(3,30)= -3, ITB(4,30)= 0;
00355       TB(31) = +0.000670;
00356       ITB(0,31)= +0; ITB(1,31)= -2; ITB(2,31)= +4; ITB(3,31)= +1, ITB(4,31)= 0;
00357       TB(32) = +0.000606;
00358       ITB(0,32)= +0; ITB(1,32)= +0; ITB(2,32)= +2; ITB(3,32)= -3, ITB(4,32)= 0;
00359       TB(33) = +0.000597;
00360       ITB(0,33)= +0; ITB(1,33)= +2; ITB(2,33)= +2; ITB(3,33)= -1, ITB(4,33)= 0;
00361       TB(34) = +0.000492;
00362       ITB(0,34)= -1; ITB(1,34)= +1; ITB(2,34)= +2; ITB(3,34)= -1, ITB(4,34)= 1;
00363       TB(35) = +0.000450;
00364       ITB(0,35)= +0; ITB(1,35)= +2; ITB(2,35)= -2; ITB(3,35)= -1, ITB(4,35)= 0;
00365       TB(36) = +0.000439;
00366       ITB(0,36)= +0; ITB(1,36)= +3; ITB(2,36)= +0; ITB(3,36)= -1, ITB(4,36)= 0;
00367       TB(37) = +0.000423;
00368       ITB(0,37)= +0; ITB(1,37)= +2; ITB(2,37)= +2; ITB(3,37)= +1, ITB(4,37)= 0;
00369       TB(38) = +0.000422;
00370       ITB(0,38)= +0; ITB(1,38)= -3; ITB(2,38)= +2; ITB(3,38)= -1, ITB(4,38)= 0;
00371       TB(39) = -0.000367;
00372       ITB(0,39)= +1; ITB(1,39)= -1; ITB(2,39)= +2; ITB(3,39)= +1, ITB(4,39)= 1;
00373       TB(40) = -0.000353;
00374       ITB(0,40)= +1; ITB(1,40)= +0; ITB(2,40)= +2; ITB(3,40)= +1, ITB(4,40)= 1;
00375       TB(41) = +0.000331;
00376       ITB(0,41)= +0; ITB(1,41)= +0; ITB(2,41)= +4; ITB(3,41)= +1, ITB(4,41)= 0;
00377       TB(42) = +0.000317;
00378       ITB(0,42)= -1; ITB(1,42)= +1; ITB(2,42)= +2; ITB(3,42)= +1, ITB(4,42)= 1;
00379       TB(43) = +0.000306;
00380       ITB(0,43)= -2; ITB(1,43)= +0; ITB(2,43)= +2; ITB(3,43)= -1, ITB(4,43)= 2;
00381       TB(44) = -0.000283;
00382       ITB(0,44)= +0; ITB(1,44)= +1; ITB(2,44)= +0; ITB(3,44)= +3, ITB(4,44)= 0;
00383 
00384 
00385          // Parallax
00386       TP( 0) = +0.950724;
00387       //          M              M'             D              F             n
00388       ITP(0, 0)= +0; ITP(1, 0)= +0; ITP(2, 0)= +0; ITP(3, 0)= +0, ITP(4, 0)= 0;
00389       TP( 1) = +0.051818;
00390       ITP(0, 1)= +0; ITP(1, 1)= +1; ITP(2, 1)= +0; ITP(3, 1)= +0, ITP(4, 1)= 0;
00391       TP( 2) = +0.009531;
00392       ITP(0, 2)= +0; ITP(1, 2)= -1; ITP(2, 2)= +2; ITP(3, 2)= +0, ITP(4, 2)= 0;
00393       TP( 3) = +0.007843;
00394       ITP(0, 3)= +0; ITP(1, 3)= +0; ITP(2, 3)= +2; ITP(3, 3)= +0, ITP(4, 3)= 0;
00395       TP( 4) = +0.002824;
00396       ITP(0, 4)= +0; ITP(1, 4)= +2; ITP(2, 4)= +0; ITP(3, 4)= +0, ITP(4, 4)= 0;
00397       TP( 5) = +0.000857;
00398       ITP(0, 5)= +0; ITP(1, 5)= +1; ITP(2, 5)= +2; ITP(3, 5)= +0, ITP(4, 5)= 0;
00399       TP( 6) = +0.000533;
00400       ITP(0, 6)= -1; ITP(1, 6)= +0; ITP(2, 6)= +2; ITP(3, 6)= +0, ITP(4, 6)= 1;
00401       TP( 7) = +0.000401;
00402       ITP(0, 7)= -1; ITP(1, 7)= -1; ITP(2, 7)= +2; ITP(3, 7)= +0, ITP(4, 7)= 1;
00403       TP( 8) = +0.000320;
00404       ITP(0, 8)= -1; ITP(1, 8)= +1; ITP(2, 8)= +0; ITP(3, 8)= +0, ITP(4, 8)= 1;
00405       TP( 9) = -0.000271;
00406       ITP(0, 9)= +0; ITP(1, 9)= +0; ITP(2, 9)= +1; ITP(3, 9)= +0, ITP(4, 9)= 0;
00407       TP(10) = -0.000264;
00408       ITP(0,10)= +1; ITP(1,10)= +1; ITP(2,10)= +0; ITP(3,10)= +0, ITP(4,10)= 1;
00409       TP(11) = -0.000198;
00410       ITP(0,11)= +0; ITP(1,11)= -1; ITP(2,11)= +0; ITP(3,11)= +2, ITP(4,11)= 0;
00411       TP(12) = +0.000173;
00412       ITP(0,12)= +0; ITP(1,12)= +3; ITP(2,12)= +0; ITP(3,12)= +0, ITP(4,12)= 0;
00413       TP(13) = +0.000167;
00414       ITP(0,13)= +0; ITP(1,13)= -1; ITP(2,13)= +4; ITP(3,13)= +0, ITP(4,13)= 0;
00415       TP(14) = -0.000111;
00416       ITP(0,14)= +1; ITP(1,14)= +0; ITP(2,14)= +0; ITP(3,14)= +0, ITP(4,14)= 1;
00417       TP(15) = +0.000103;
00418       ITP(0,15)= +0; ITP(1,15)= -2; ITP(2,15)= +4; ITP(3,15)= +0, ITP(4,15)= 0;
00419       TP(16) = -0.000084;
00420       ITP(0,16)= +0; ITP(1,16)= +2; ITP(2,16)= -2; ITP(3,16)= +0, ITP(4,16)= 0;
00421       TP(17) = -0.000083;
00422       ITP(0,17)= +1; ITP(1,17)= +0; ITP(2,17)= +2; ITP(3,17)= +0, ITP(4,17)= 1;
00423       TP(18) = +0.000079;
00424       ITP(0,18)= +0; ITP(1,18)= +2; ITP(2,18)= +2; ITP(3,18)= +0, ITP(4,18)= 0;
00425       TP(19) = +0.000072;
00426       ITP(0,19)= +0; ITP(1,19)= +0; ITP(2,19)= +4; ITP(3,19)= +0, ITP(4,19)= 0;
00427       TP(20) = +0.000064;
00428       ITP(0,20)= -1; ITP(1,20)= +1; ITP(2,20)= +2; ITP(3,20)= +0, ITP(4,20)= 1;
00429       TP(21) = -0.000063;
00430       ITP(0,21)= +1; ITP(1,21)= -1; ITP(2,21)= +2; ITP(3,21)= +0, ITP(4,21)= 1;
00431       TP(22) = +0.000041;
00432       ITP(0,22)= +1; ITP(1,22)= +0; ITP(2,22)= +1; ITP(3,22)= +0, ITP(4,22)= 1;
00433       TP(23) = +0.000035;
00434       ITP(0,23)= -1; ITP(1,23)= +2; ITP(2,23)= +0; ITP(3,23)= +0, ITP(4,23)= 1;
00435       TP(24) = -0.000033;
00436       ITP(0,24)= +0; ITP(1,24)= +3; ITP(2,24)= -2; ITP(3,24)= +0, ITP(4,24)= 0;
00437       TP(25) = -0.000030;
00438       ITP(0,25)= +0; ITP(1,25)= +1; ITP(2,25)= +1; ITP(3,25)= +0, ITP(4,25)= 0;
00439       TP(26) = -0.000029;
00440       ITP(0,26)= +0; ITP(1,26)= +0; ITP(2,26)= -2; ITP(3,26)= +2, ITP(4,26)= 0;
00441       TP(27) = -0.000029;
00442       ITP(0,27)= +1; ITP(1,27)= +2; ITP(2,27)= +0; ITP(3,27)= +0, ITP(4,27)= 1;
00443       TP(28) = +0.000026;
00444       ITP(0,28)= -2; ITP(1,28)= +0; ITP(2,28)= +2; ITP(3,28)= +0, ITP(4,28)= 2;
00445       TP(29) = -0.000023;
00446       ITP(0,29)= +0; ITP(1,29)= +1; ITP(2,29)= -2; ITP(3,29)= +2, ITP(4,29)= 0;
00447       TP(30) = +0.000019;
00448       ITP(0,30)= -1; ITP(1,30)= -1; ITP(2,30)= +4; ITP(3,30)= +0, ITP(4,30)= 1;
00449 
00450 
00451          // Centuries since J1900
00452       double tt((t.MJD()-15019.5)/36525.0);
00453 
00454 
00455          // Fundamental arguments (radians) and derivatives (radians per
00456          // Julian century) for the current epoch
00457 
00458          // Moon's mean longitude
00459       double ELP(D2R*fmod( (ELP0+(ELP1+(ELP2+ELP3*tt)*tt)*tt), 360.0));
00460       double DELP(D2R*(ELP1+(2.0*ELP2+3.0*ELP3*tt)*tt));
00461 
00462          // Sun's mean anomaly
00463       double EM(D2R*fmod( (EM0+(EM1+(EM2+EM3*tt)*tt)*tt), 360.0));
00464       double DEM(D2R*(EM1+(2.0*EM2+3.0*EM3*tt)*tt));
00465 
00466          // Moon's mean anomaly
00467       double EMP(D2R*fmod( (EMP0+(EMP1+(EMP2+EMP3*tt)*tt)*tt), 360.0));
00468       double DEMP(D2R*(EMP1+(2.0*EMP2+3.0*EMP3*tt)*tt));
00469 
00470          // Moon's mean elongation
00471       double D(D2R*fmod( (D0+(D1+(D2+D3*tt)*tt)*tt), 360.0));
00472       double DD(D2R*(D1+(2.0*D2+3.0*D3*tt)*tt));
00473 
00474          // Mean distance of the Moon from its ascending node
00475       double F(D2R*fmod( (F0+(F1+(F2+F3*tt)*tt)*tt), 360.0));
00476       double DF(D2R*(F1+(2.0*F2+3.0*F3*tt)*tt));
00477 
00478          // Longitude of the Moon's ascending node
00479       double OM(D2R*fmod( (OM0+(OM1+(OM2+OM3*tt)*tt)*tt), 360.0));
00480       double DOM(D2R*(OM1+(2.0*OM2+3.0*OM3*tt)*tt));
00481       double SINOM(std::sin(OM));
00482       double COSOM(std::cos(OM));
00483       double DOMCOM(DOM*COSOM);
00484 
00485 
00486          // Let's add the periodic variations
00487       double  THETA(D2R*(PA0+PA1*tt));
00488       double  WA(std::sin(THETA));
00489       double  DWA(D2R*PA1*std::cos(THETA));
00490       THETA = D2R*(PE0+(PE1+PE2*tt)*tt);
00491       double  WB(PEC*std::sin(THETA));
00492       double  DWB(D2R*PEC*(PE1+2.0*PE2*tt)*std::cos(THETA));
00493       ELP   = ELP+D2R*(PAC*WA+WB+PFC*SINOM);
00494       DELP  = DELP+D2R*(PAC*DWA+DWB+PFC*DOMCOM);
00495       EM    = EM+D2R*PBC*WA;
00496       DEM   = DEM+D2R*PBC*DWA;
00497       EMP   = EMP+D2R*(PCC*WA+WB+PGC*SINOM);
00498       DEMP  = DEMP+D2R*(PCC*DWA+DWB+PGC*DOMCOM);
00499       D     = D+D2R*(PDC*WA+WB+PHC*SINOM);
00500       DD    = DD+D2R*(PDC*DWA+DWB+PHC*DOMCOM);
00501       double  WOM(OM+D2R*(PJ0+PJ1*tt));
00502       double  DWOM(DOM+D2R*PJ1);
00503       double  SINWOM(std::sin(WOM));
00504       double  COSWOM(std::cos(WOM));
00505       F     = F+D2R*(WB+cPIC*SINOM+PJC*SINWOM);
00506       DF    = DF+D2R*(DWB+cPIC*DOMCOM+PJC*DWOM*COSWOM);
00507 
00508 
00509          // E-factor, and square
00510       double  E(1.0+(E1+E2*tt)*tt);
00511       double  DE(E1+2.0*E2*tt);
00512       double  ESQ(E*E);
00513       double  DESQ(2.0*E*DE);
00514 
00515 
00516          // Series expansions
00517 
00518          // Longitude
00519       double V(0.0);
00520       double DV(0.0);
00521       for (int n=(NL-1); n>=0; n--)
00522       {
00523          double EN, DEN;
00524          double COEFF(TL(n));
00525          double EMN(static_cast<double>(ITL(0,n)));
00526          double EMPN(static_cast<double>(ITL(1,n)));
00527          double DN(static_cast<double>(ITL(2,n)));
00528          double FN(static_cast<double>(ITL(3,n)));
00529          int I=ITL(4,n);
00530          if (I == 0)
00531          {
00532             EN  = 1.0;
00533             DEN = 0.0;
00534          }
00535          else
00536          {
00537             if (I == 1)
00538             {
00539                EN  = E;
00540                DEN = DE;
00541             }
00542             else
00543             {
00544                EN  = ESQ;
00545                DEN = DESQ;
00546             }
00547          }
00548          THETA = EMN*EM+EMPN*EMP+DN*D+FN*F;
00549          double DTHETA(EMN*DEM+EMPN*DEMP+DN*DD+FN*DF);
00550          double FTHETA(std::sin(THETA));
00551          V  = V+COEFF*FTHETA*EN;
00552          DV = DV+COEFF*(std::cos(THETA)*DTHETA*EN+FTHETA*DEN);
00553       }
00554       double EL(ELP+D2R*V);
00555 
00556 
00557          // Latitude
00558       V  = 0.0;
00559       DV = 0.0;
00560       for (int n=(NB-1); n>=0; n--)
00561       {
00562          double EN, DEN;
00563          double COEFF(TB(n));
00564          double EMN(static_cast<double>(ITB(0,n)));
00565          double EMPN(static_cast<double>(ITB(1,n)));
00566          double DN(static_cast<double>(ITB(2,n)));
00567          double FN(static_cast<double>(ITB(3,n)));
00568          int I=ITB(4,n);
00569          if (I == 0)
00570          {
00571             EN  = 1.0;
00572             DEN = 0.0;
00573          }
00574          else
00575          {
00576             if (I == 1)
00577             {
00578                EN  = E;
00579                DEN = DE;
00580             }
00581             else
00582             {
00583                EN  = ESQ;
00584                DEN = DESQ;
00585             }
00586          }
00587          THETA = EMN*EM+EMPN*EMP+DN*D+FN*F;
00588          double DTHETA(EMN*DEM+EMPN*DEMP+DN*DD+FN*DF);
00589          double FTHETA(std::sin(THETA));
00590          V  = V+COEFF*FTHETA*EN;
00591          DV = DV+COEFF*(std::cos(THETA)*DTHETA*EN+FTHETA*DEN);
00592       }
00593       double BF(1.0-CW1*COSOM-CW2*COSWOM);
00594       double B(D2R*V*BF);
00595 
00596 
00597          // Parallax
00598       V  = 0.0;
00599       DV = 0.0;
00600       for (int n=(NP-1); n>=0; n--)
00601       {
00602          double EN, DEN;
00603          double COEFF(TP(n));
00604          double EMN(static_cast<double>(ITP(0,n)));
00605          double EMPN(static_cast<double>(ITP(1,n)));
00606          double DN(static_cast<double>(ITP(2,n)));
00607          double FN(static_cast<double>(ITP(3,n)));
00608          int I=ITP(4,n);
00609          if (I == 0)
00610          {
00611             EN  = 1.0;
00612             DEN = 0.0;
00613          }
00614          else
00615          {
00616             if (I == 1)
00617             {
00618                EN  = E;
00619                DEN = DE;
00620             }
00621             else
00622             {
00623                EN  = ESQ;
00624                DEN = DESQ;
00625             }
00626          }
00627          THETA = EMN*EM+EMPN*EMP+DN*D+FN*F;
00628          double DTHETA(EMN*DEM+EMPN*DEMP+DN*DD+FN*DF);
00629          double FTHETA(std::cos(THETA));
00630          V  = V+COEFF*FTHETA*EN;
00631          DV = DV+COEFF*(-std::sin(THETA)*DTHETA*EN+FTHETA*DEN);
00632       }
00633       double P(D2R*V);
00634 
00635 
00636          // Transformation into final form
00637 
00638          // Parallax to distance (AU, AU/sec)
00639       double SP(std::sin(P));
00640       double R(ERADAU/SP);
00641 
00642          // Longitude, latitude to x,y,z (AU)
00643       double SEL(std::sin(EL));
00644       double CEL(std::cos(EL));
00645       double SB(std::sin(B));
00646       double CB(std::cos(B));
00647       double RCB(R*CB);
00648       double X(RCB*CEL);
00649       double Y(RCB*SEL);
00650       double Z(R*SB);
00651 
00652 
00653          // Julian centuries since J2000
00654       tt=(t.MJD()-51544.5)/36525.0;
00655 
00656          // Fricke equinox correction
00657       double EPJ(2000.0+tt*100.0);
00658       double EQCOR(DS2R*(0.035+0.00085*(EPJ-B1950)));
00659 
00660          // Mean obliquity (IAU 1976)
00661       double EPS(DAS2R*(84381.448 + (-46.8150 + 
00662                         (-0.00059+0.001813*tt)*tt)*tt));
00663 
00664          // Change to equatorial system, mean of date, FK5 system
00665       double SINEPS(std::sin(EPS));
00666       double COSEPS(std::cos(EPS));
00667       double ES(EQCOR*SINEPS);
00668       double EC(EQCOR*COSEPS);
00669 
00670       Triple res;
00671 
00672       res.theArray[0] = (X-EC*Y+ES*Z)*AU_CONST;
00673       res.theArray[1] = (EQCOR*X+Y*COSEPS-Z*SINEPS)*AU_CONST;
00674       res.theArray[2] = (Y*SINEPS+Z*COSEPS)*AU_CONST;
00675 
00676 
00677       return res;
00678 
00679    } // End MoonPosition::getPositionCIS()
00680 
00681 
00682 } // end namespace gpstk