1 //============================================================================== 2 // 3 // This file is part of GPSTk, the GPS Toolkit. 4 // 5 // The GPSTk is free software; you can redistribute it and/or modify 6 // it under the terms of the GNU Lesser General Public License as published 7 // by the Free Software Foundation; either version 3.0 of the License, or 8 // any later version. 9 // 10 // The GPSTk is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU Lesser General Public License for more details. 14 // 15 // You should have received a copy of the GNU Lesser General Public 16 // License along with GPSTk; if not, write to the Free Software Foundation, 17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA 18 // 19 // This software was developed by Applied Research Laboratories at the 20 // University of Texas at Austin. 21 // Copyright 2004-2020, The Board of Regents of The University of Texas System 22 // 23 //============================================================================== 24 25 //============================================================================== 26 // 27 // This software was developed by Applied Research Laboratories at the 28 // University of Texas at Austin, under contract to an agency or agencies 29 // within the U.S. Department of Defense. The U.S. Government retains all 30 // rights to use, duplicate, distribute, disclose, or release this software. 31 // 32 // Pursuant to DoD Directive 523024 33 // 34 // DISTRIBUTION STATEMENT A: This software has been approved for public 35 // release, distribution is unlimited. 36 // 37 //============================================================================== 38 39 /** 40 * @file MoonPosition.cpp 41 * Returns the approximate position of the Moon at the given epoch in the 42 * ECEF system. 43 */ 44 45 #include "MoonPosition.hpp" 46 #include "MJD.hpp" 47 #include "CivilTime.hpp" 48 #include "EpochDataStore.hpp" 49 50 51 namespace gpstk 52 { 53 54 55 // Time of the first valid time 56 const CommonTime MoonPosition::initialTime = CivilTime(1900, 3, 1, 0, 0, 0.0,TimeSystem::Any); 57 58 // Time of the last valid time 59 const CommonTime MoonPosition::finalTime = CivilTime(2100, 2, 28, 0, 0, 0.0,TimeSystem::Any); 60 61 62 // Coefficients for fundamental arguments 63 // Units are degrees for position and Julian centuries for time 64 65 // Moon's mean longitude 66 const double MoonPosition::ELP0(270.434164), 67 MoonPosition::ELP1(481267.8831), 68 MoonPosition::ELP2(-0.001133), 69 MoonPosition::ELP3(0.0000019); 70 71 // Sun's mean anomaly 72 const double MoonPosition::EM0(358.475833), 73 MoonPosition::EM1(35999.0498), 74 MoonPosition::EM2(-0.000150), 75 MoonPosition::EM3(-0.0000033); 76 77 // Moon's mean anomaly 78 const double MoonPosition::EMP0(296.104608), 79 MoonPosition::EMP1(477198.8491), 80 MoonPosition::EMP2(0.009192), 81 MoonPosition::EMP3(0.0000144); 82 83 // Moon's mean elongation 84 const double MoonPosition::D0(350.737486), 85 MoonPosition::D1(445267.1142), 86 MoonPosition::D2(-0.001436), 87 MoonPosition::D3(0.0000019); 88 89 // Mean distance of the Moon from its ascending node 90 const double MoonPosition::F0(11.250889), 91 MoonPosition::F1(483202.0251), 92 MoonPosition::F2(-0.003211), 93 MoonPosition::F3(-0.0000003); 94 95 // Longitude of the Moon's ascending node 96 const double MoonPosition::OM0(259.183275), 97 MoonPosition::OM1(-1934.1420), 98 MoonPosition::OM2(0.002078), 99 MoonPosition::OM3(0.0000022); 100 101 102 // Coefficients for (dimensionless) E factor 103 const double MoonPosition::E1(-0.002495), 104 MoonPosition::E2(-0.00000752); 105 106 // Coefficients for periodic variations, etc 107 const double MoonPosition::PAC(0.000233), 108 MoonPosition::PA0(51.2), 109 MoonPosition::PA1(20.2); 110 const double MoonPosition::PBC(-0.001778); 111 const double MoonPosition::PCC(0.000817); 112 const double MoonPosition::PDC(0.002011); 113 const double MoonPosition::PEC(0.003964), 114 MoonPosition::PE0(346.560), 115 MoonPosition::PE1(132.870), 116 MoonPosition::PE2(-0.0091731); 117 const double MoonPosition::PFC(0.001964); 118 const double MoonPosition::PGC(0.002541); 119 const double MoonPosition::PHC(0.001964); 120 const double MoonPosition::cPIC(-0.024691); 121 const double MoonPosition::PJC(-0.004328), 122 MoonPosition::PJ0(275.05), 123 MoonPosition::PJ1(-2.30); 124 const double MoonPosition::CW1(0.0004664); 125 const double MoonPosition::CW2(0.0000754); 126 127 // Coefficients for Moon position 128 // Tx(N): coefficient of L, B or P term (deg) 129 // ITx(N,0-4): coefficients of M, M', D, F, E**n in argument 130 // 131 const size_t MoonPosition::NL(50), 132 MoonPosition::NB(45), 133 MoonPosition::NP(31); 134 Vector<double> MoonPosition::TL(NL,0.0), 135 MoonPosition::TB(NB,0.0), 136 MoonPosition::TP(NP,0.0); 137 Matrix<int> MoonPosition::ITL(5,NL,0), 138 MoonPosition::ITB(5,NB,0), 139 MoonPosition::ITP(5,NP,0); 140 141 142 143 // Returns the position of Moon ECEF coordinates (meters) at 144 // the indicated time. 145 // 146 // @param[in] t the time to look up 147 // 148 // @return the position of the Moon at time (as a Triple) 149 // 150 // @throw InvalidRequest If the request can not be completed for 151 // any reason, this is thrown. The text may have additional 152 // information as to why the request failed. getPosition(const CommonTime & t) const153 Triple MoonPosition::getPosition(const CommonTime& t) const 154 { 155 // Test if the time interval is correct 156 if ( (t < MoonPosition::initialTime) || 157 (t > MoonPosition::finalTime) ) 158 { 159 InvalidRequest ir("Provided epoch is out of bounds."); 160 GPSTK_THROW(ir); 161 } 162 163 164 // We will store here the results 165 Triple res; 166 167 res = MoonPosition::getPositionCIS(t); 168 res = CIS2CTS(res, t); 169 170 return res; 171 172 } // End MoonPosition::getPosition 173 174 175 176 /* Function to compute Moon position in CIS system (coordinates MJD 177 * in meters) 178 * 179 * @param t Epoch 180 * @throw InvalidRequest 181 */ getPositionCIS(const CommonTime & t) const182 Triple MoonPosition::getPositionCIS(const CommonTime& t) const 183 { 184 // Test if the time interval is correct 185 if ( (t < MoonPosition::initialTime) || 186 (t > MoonPosition::finalTime) ) 187 { 188 InvalidRequest ir("Provided epoch is out of bounds."); 189 GPSTK_THROW(ir); 190 } 191 192 193 194 // Coefficients for Moon position 195 // Tx(N): coefficient of L, B or P term (deg) 196 // ITx(N,0-4): coefficients of M, M', D, F, E**n in argument 197 // 198 199 // Longitude 200 TL( 0) = +6.288750; 201 // M M' D F n 202 203 ITL(0, 0)= +0; ITL(1, 0)= +1; ITL(2, 0)= +0; ITL(3, 0)= +0, ITL(4, 0)= 0; 204 TL( 1) = +1.274018; 205 ITL(0, 1)= +0; ITL(1, 1)= -1; ITL(2, 1)= +2; ITL(3, 1)= +0, ITL(4, 1)= 0; 206 TL( 2) = +0.658309; 207 ITL(0, 2)= +0; ITL(1, 2)= +0; ITL(2, 2)= +2; ITL(3, 2)= +0, ITL(4, 2)= 0; 208 TL( 3) = +0.213616; 209 ITL(0, 3)= +0; ITL(1, 3)= +2; ITL(2, 3)= +0; ITL(3, 3)= +0, ITL(4, 3)= 0; 210 TL( 4) = -0.185596; 211 ITL(0, 4)= +1; ITL(1, 4)= +0; ITL(2, 4)= +0; ITL(3, 4)= +0, ITL(4, 4)= 1; 212 TL( 5) = -0.114336; 213 ITL(0, 5)= +0; ITL(1, 5)= +0; ITL(2, 5)= +0; ITL(3, 5)= +2, ITL(4, 5)= 0; 214 TL( 6) = +0.058793; 215 ITL(0, 6)= +0; ITL(1, 6)= -2; ITL(2, 6)= +2; ITL(3, 6)= +0, ITL(4, 6)= 0; 216 TL( 7) = +0.057212; 217 ITL(0, 7)= -1; ITL(1, 7)= -1; ITL(2, 7)= +2; ITL(3, 7)= +0, ITL(4, 7)= 1; 218 TL( 8) = +0.053320; 219 ITL(0, 8)= +0; ITL(1, 8)= +1; ITL(2, 8)= +2; ITL(3, 8)= +0, ITL(4, 8)= 0; 220 TL( 9) = +0.045874; 221 ITL(0, 9)= -1; ITL(1, 9)= +0; ITL(2, 9)= +2; ITL(3, 9)= +0, ITL(4, 9)= 1; 222 TL(10) = +0.041024; 223 ITL(0,10)= -1; ITL(1,10)= +1; ITL(2,10)= +0; ITL(3,10)= +0, ITL(4,10)= 1; 224 TL(11) = -0.034718; 225 ITL(0,11)= +0; ITL(1,11)= +0; ITL(2,11)= +1; ITL(3,11)= +0, ITL(4,11)= 0; 226 TL(12) = -0.030465; 227 ITL(0,12)= +1; ITL(1,12)= +1; ITL(2,12)= +0; ITL(3,12)= +0, ITL(4,12)= 1; 228 TL(13) = +0.015326; 229 ITL(0,13)= +0; ITL(1,13)= +0; ITL(2,13)= +2; ITL(3,13)= -2, ITL(4,13)= 0; 230 TL(14) = -0.012528; 231 ITL(0,14)= +0; ITL(1,14)= +1; ITL(2,14)= +0; ITL(3,14)= +2, ITL(4,14)= 0; 232 TL(15) = -0.010980; 233 ITL(0,15)= +0; ITL(1,15)= -1; ITL(2,15)= +0; ITL(3,15)= +2, ITL(4,15)= 0; 234 TL(16) = +0.010674; 235 ITL(0,16)= +0; ITL(1,16)= -1; ITL(2,16)= +4; ITL(3,16)= +0, ITL(4,16)= 0; 236 TL(17) = +0.010034; 237 ITL(0,17)= +0; ITL(1,17)= +3; ITL(2,17)= +0; ITL(3,17)= +0, ITL(4,17)= 0; 238 TL(18) = +0.008548; 239 ITL(0,18)= +0; ITL(1,18)= -2; ITL(2,18)= +4; ITL(3,18)= +0, ITL(4,18)= 0; 240 TL(19) = -0.007910; 241 ITL(0,19)= +1; ITL(1,19)= -1; ITL(2,19)= +2; ITL(3,19)= +0, ITL(4,19)= 1; 242 TL(20) = -0.006783; 243 ITL(0,20)= +1; ITL(1,20)= +0; ITL(2,20)= +2; ITL(3,20)= +0, ITL(4,20)= 1; 244 TL(21) = +0.005162; 245 ITL(0,21)= +0; ITL(1,21)= +1; ITL(2,21)= -1; ITL(3,21)= +0, ITL(4,21)= 0; 246 TL(22) = +0.005000; 247 ITL(0,22)= +1; ITL(1,22)= +0; ITL(2,22)= +1; ITL(3,22)= +0, ITL(4,22)= 1; 248 TL(23) = +0.004049; 249 ITL(0,23)= -1; ITL(1,23)= +1; ITL(2,23)= +2; ITL(3,23)= +0, ITL(4,23)= 1; 250 TL(24) = +0.003996; 251 ITL(0,24)= +0; ITL(1,24)= +2; ITL(2,24)= +2; ITL(3,24)= +0, ITL(4,24)= 0; 252 TL(25) = +0.003862; 253 ITL(0,25)= +0; ITL(1,25)= +0; ITL(2,25)= +4; ITL(3,25)= +0, ITL(4,25)= 0; 254 TL(26) = +0.003665; 255 ITL(0,26)= +0; ITL(1,26)= -3; ITL(2,26)= +2; ITL(3,26)= +0, ITL(4,26)= 0; 256 TL(27) = +0.002695; 257 ITL(0,27)= -1; ITL(1,27)= +2; ITL(2,27)= +0; ITL(3,27)= +0, ITL(4,27)= 1; 258 TL(28) = +0.002602; 259 ITL(0,28)= +0; ITL(1,28)= +1; ITL(2,28)= -2; ITL(3,28)= -2, ITL(4,28)= 0; 260 TL(29) = +0.002396; 261 ITL(0,29)= -1; ITL(1,29)= -2; ITL(2,29)= +2; ITL(3,29)= +0, ITL(4,29)= 1; 262 TL(30) = -0.002349; 263 ITL(0,30)= +0; ITL(1,30)= +1; ITL(2,30)= +1; ITL(3,30)= +0, ITL(4,30)= 0; 264 TL(31) = +0.002249; 265 ITL(0,31)= -2; ITL(1,31)= +0; ITL(2,31)= +2; ITL(3,31)= +0, ITL(4,31)= 2; 266 TL(32) = -0.002125; 267 ITL(0,32)= +1; ITL(1,32)= +2; ITL(2,32)= +0; ITL(3,32)= +0, ITL(4,32)= 1; 268 TL(33) = -0.002079; 269 ITL(0,33)= +2; ITL(1,33)= +0; ITL(2,33)= +0; ITL(3,33)= +0, ITL(4,33)= 2; 270 TL(34) = +0.002059; 271 ITL(0,34)= -2; ITL(1,34)= -1; ITL(2,34)= +2; ITL(3,34)= +0, ITL(4,34)= 2; 272 TL(35) = -0.001773; 273 ITL(0,35)= +0; ITL(1,35)= +1; ITL(2,35)= +2; ITL(3,35)= -2, ITL(4,35)= 0; 274 TL(36) = -0.001595; 275 ITL(0,36)= +0; ITL(1,36)= +0; ITL(2,36)= +2; ITL(3,36)= +2, ITL(4,36)= 0; 276 TL(37) = +0.001220; 277 ITL(0,37)= -1; ITL(1,37)= -1; ITL(2,37)= +4; ITL(3,37)= +0, ITL(4,37)= 1; 278 TL(38) = -0.001110; 279 ITL(0,38)= +0; ITL(1,38)= +2; ITL(2,38)= +0; ITL(3,38)= +2, ITL(4,38)= 0; 280 TL(39) = +0.000892; 281 ITL(0,39)= +0; ITL(1,39)= +1; ITL(2,39)= -3; ITL(3,39)= +0, ITL(4,39)= 0; 282 TL(40) = -0.000811; 283 ITL(0,40)= +1; ITL(1,40)= +1; ITL(2,40)= +2; ITL(3,40)= +0, ITL(4,40)= 1; 284 TL(41) = +0.000761; 285 ITL(0,41)= -1; ITL(1,41)= -2; ITL(2,41)= +4; ITL(3,41)= +0, ITL(4,41)= 1; 286 TL(42) = +0.000717; 287 ITL(0,42)= -2; ITL(1,42)= +1; ITL(2,42)= +0; ITL(3,42)= +0, ITL(4,42)= 2; 288 TL(43) = +0.000704; 289 ITL(0,43)= -2; ITL(1,43)= +1; ITL(2,43)= -2; ITL(3,43)= +0, ITL(4,43)= 2; 290 TL(44) = +0.000693; 291 ITL(0,44)= +1; ITL(1,44)= -2; ITL(2,44)= +2; ITL(3,44)= +0, ITL(4,44)= 1; 292 TL(45) = +0.000598; 293 ITL(0,45)= -1; ITL(1,45)= +0; ITL(2,45)= +2; ITL(3,45)= -2, ITL(4,45)= 1; 294 TL(46) = +0.000550; 295 ITL(0,46)= +0; ITL(1,46)= +1; ITL(2,46)= +4; ITL(3,46)= +0, ITL(4,46)= 0; 296 TL(47) = +0.000538; 297 ITL(0,47)= +0; ITL(1,47)= +4; ITL(2,47)= +0; ITL(3,47)= +0, ITL(4,47)= 0; 298 TL(48) = +0.000521; 299 ITL(0,48)= -1; ITL(1,48)= +0; ITL(2,48)= +4; ITL(3,48)= +0, ITL(4,48)= 1; 300 TL(49) = +0.000486; 301 ITL(0,49)= +0; ITL(1,49)= +2; ITL(2,49)= -1; ITL(3,49)= +0, ITL(4,49)= 0; 302 303 304 // Latitude 305 TB( 0) = +5.128189; 306 // M M' D F n 307 ITB(0, 0)= +0; ITB(1, 0)= +0; ITB(2, 0)= +0; ITB(3, 0)= +1, ITB(4, 0)= 0; 308 TB( 1) = +0.280606; 309 ITB(0, 1)= +0; ITB(1, 1)= +1; ITB(2, 1)= +0; ITB(3, 1)= +1, ITB(4, 1)= 0; 310 TB( 2) = +0.277693; 311 ITB(0, 2)= +0; ITB(1, 2)= +1; ITB(2, 2)= +0; ITB(3, 2)= -1, ITB(4, 2)= 0; 312 TB( 3) = +0.173238; 313 ITB(0, 3)= +0; ITB(1, 3)= +0; ITB(2, 3)= +2; ITB(3, 3)= -1, ITB(4, 3)= 0; 314 TB( 4) = +0.055413; 315 ITB(0, 4)= +0; ITB(1, 4)= -1; ITB(2, 4)= +2; ITB(3, 4)= +1, ITB(4, 4)= 0; 316 TB( 5) = +0.046272; 317 ITB(0, 5)= +0; ITB(1, 5)= -1; ITB(2, 5)= +2; ITB(3, 5)= -1, ITB(4, 5)= 0; 318 TB( 6) = +0.032573; 319 ITB(0, 6)= +0; ITB(1, 6)= +0; ITB(2, 6)= +2; ITB(3, 6)= +1, ITB(4, 6)= 0; 320 TB( 7) = +0.017198; 321 ITB(0, 7)= +0; ITB(1, 7)= +2; ITB(2, 7)= +0; ITB(3, 7)= +1, ITB(4, 7)= 0; 322 TB( 8) = +0.009267; 323 ITB(0, 8)= +0; ITB(1, 8)= +1; ITB(2, 8)= +2; ITB(3, 8)= -1, ITB(4, 8)= 0; 324 TB( 9) = +0.008823; 325 ITB(0, 9)= +0; ITB(1, 9)= +2; ITB(2, 9)= +0; ITB(3, 9)= -1, ITB(4, 9)= 0; 326 TB(10) = +0.008247; 327 ITB(0,10)= -1; ITB(1,10)= +0; ITB(2,10)= +2; ITB(3,10)= -1, ITB(4,10)= 1; 328 TB(11) = +0.004323; 329 ITB(0,11)= +0; ITB(1,11)= -2; ITB(2,11)= +2; ITB(3,11)= -1, ITB(4,11)= 0; 330 TB(12) = +0.004200; 331 ITB(0,12)= +0; ITB(1,12)= +1; ITB(2,12)= +2; ITB(3,12)= +1, ITB(4,12)= 0; 332 TB(13) = +0.003372; 333 ITB(0,13)= -1; ITB(1,13)= +0; ITB(2,13)= -2; ITB(3,13)= +1, ITB(4,13)= 1; 334 TB(14) = +0.002472; 335 ITB(0,14)= -1; ITB(1,14)= -1; ITB(2,14)= +2; ITB(3,14)= +1, ITB(4,14)= 1; 336 TB(15) = +0.002222; 337 ITB(0,15)= -1; ITB(1,15)= +0; ITB(2,15)= +2; ITB(3,15)= +1, ITB(4,15)= 1; 338 TB(16) = +0.002072; 339 ITB(0,16)= -1; ITB(1,16)= -1; ITB(2,16)= +2; ITB(3,16)= -1, ITB(4,16)= 1; 340 TB(17) = +0.001877; 341 ITB(0,17)= -1; ITB(1,17)= +1; ITB(2,17)= +0; ITB(3,17)= +1, ITB(4,17)= 1; 342 TB(18) = +0.001828; 343 ITB(0,18)= +0; ITB(1,18)= -1; ITB(2,18)= +4; ITB(3,18)= -1, ITB(4,18)= 0; 344 TB(19) = -0.001803; 345 ITB(0,19)= +1; ITB(1,19)= +0; ITB(2,19)= +0; ITB(3,19)= +1, ITB(4,19)= 1; 346 TB(20) = -0.001750; 347 ITB(0,20)= +0; ITB(1,20)= +0; ITB(2,20)= +0; ITB(3,20)= +3, ITB(4,20)= 0; 348 TB(21) = +0.001570; 349 ITB(0,21)= -1; ITB(1,21)= +1; ITB(2,21)= +0; ITB(3,21)= -1, ITB(4,21)= 1; 350 TB(22) = -0.001487; 351 ITB(0,22)= +0; ITB(1,22)= +0; ITB(2,22)= +1; ITB(3,22)= +1, ITB(4,22)= 0; 352 TB(23) = -0.001481; 353 ITB(0,23)= +1; ITB(1,23)= +1; ITB(2,23)= +0; ITB(3,23)= +1, ITB(4,23)= 1; 354 TB(24) = +0.001417; 355 ITB(0,24)= -1; ITB(1,24)= -1; ITB(2,24)= +0; ITB(3,24)= +1, ITB(4,24)= 1; 356 TB(25) = +0.001350; 357 ITB(0,25)= -1; ITB(1,25)= +0; ITB(2,25)= +0; ITB(3,25)= +1, ITB(4,25)= 1; 358 TB(26) = +0.001330; 359 ITB(0,26)= +0; ITB(1,26)= +0; ITB(2,26)= -1; ITB(3,26)= +1, ITB(4,26)= 0; 360 TB(27) = +0.001106; 361 ITB(0,27)= +0; ITB(1,27)= +3; ITB(2,27)= +0; ITB(3,27)= +1, ITB(4,27)= 0; 362 TB(28) = +0.001020; 363 ITB(0,28)= +0; ITB(1,28)= +0; ITB(2,28)= +4; ITB(3,28)= -1, ITB(4,28)= 0; 364 TB(29) = +0.000833; 365 ITB(0,29)= +0; ITB(1,29)= -1; ITB(2,29)= +4; ITB(3,29)= +1, ITB(4,29)= 0; 366 TB(30) = +0.000781; 367 ITB(0,30)= +0; ITB(1,30)= +1; ITB(2,30)= +0; ITB(3,30)= -3, ITB(4,30)= 0; 368 TB(31) = +0.000670; 369 ITB(0,31)= +0; ITB(1,31)= -2; ITB(2,31)= +4; ITB(3,31)= +1, ITB(4,31)= 0; 370 TB(32) = +0.000606; 371 ITB(0,32)= +0; ITB(1,32)= +0; ITB(2,32)= +2; ITB(3,32)= -3, ITB(4,32)= 0; 372 TB(33) = +0.000597; 373 ITB(0,33)= +0; ITB(1,33)= +2; ITB(2,33)= +2; ITB(3,33)= -1, ITB(4,33)= 0; 374 TB(34) = +0.000492; 375 ITB(0,34)= -1; ITB(1,34)= +1; ITB(2,34)= +2; ITB(3,34)= -1, ITB(4,34)= 1; 376 TB(35) = +0.000450; 377 ITB(0,35)= +0; ITB(1,35)= +2; ITB(2,35)= -2; ITB(3,35)= -1, ITB(4,35)= 0; 378 TB(36) = +0.000439; 379 ITB(0,36)= +0; ITB(1,36)= +3; ITB(2,36)= +0; ITB(3,36)= -1, ITB(4,36)= 0; 380 TB(37) = +0.000423; 381 ITB(0,37)= +0; ITB(1,37)= +2; ITB(2,37)= +2; ITB(3,37)= +1, ITB(4,37)= 0; 382 TB(38) = +0.000422; 383 ITB(0,38)= +0; ITB(1,38)= -3; ITB(2,38)= +2; ITB(3,38)= -1, ITB(4,38)= 0; 384 TB(39) = -0.000367; 385 ITB(0,39)= +1; ITB(1,39)= -1; ITB(2,39)= +2; ITB(3,39)= +1, ITB(4,39)= 1; 386 TB(40) = -0.000353; 387 ITB(0,40)= +1; ITB(1,40)= +0; ITB(2,40)= +2; ITB(3,40)= +1, ITB(4,40)= 1; 388 TB(41) = +0.000331; 389 ITB(0,41)= +0; ITB(1,41)= +0; ITB(2,41)= +4; ITB(3,41)= +1, ITB(4,41)= 0; 390 TB(42) = +0.000317; 391 ITB(0,42)= -1; ITB(1,42)= +1; ITB(2,42)= +2; ITB(3,42)= +1, ITB(4,42)= 1; 392 TB(43) = +0.000306; 393 ITB(0,43)= -2; ITB(1,43)= +0; ITB(2,43)= +2; ITB(3,43)= -1, ITB(4,43)= 2; 394 TB(44) = -0.000283; 395 ITB(0,44)= +0; ITB(1,44)= +1; ITB(2,44)= +0; ITB(3,44)= +3, ITB(4,44)= 0; 396 397 398 // Parallax 399 TP( 0) = +0.950724; 400 // M M' D F n 401 ITP(0, 0)= +0; ITP(1, 0)= +0; ITP(2, 0)= +0; ITP(3, 0)= +0, ITP(4, 0)= 0; 402 TP( 1) = +0.051818; 403 ITP(0, 1)= +0; ITP(1, 1)= +1; ITP(2, 1)= +0; ITP(3, 1)= +0, ITP(4, 1)= 0; 404 TP( 2) = +0.009531; 405 ITP(0, 2)= +0; ITP(1, 2)= -1; ITP(2, 2)= +2; ITP(3, 2)= +0, ITP(4, 2)= 0; 406 TP( 3) = +0.007843; 407 ITP(0, 3)= +0; ITP(1, 3)= +0; ITP(2, 3)= +2; ITP(3, 3)= +0, ITP(4, 3)= 0; 408 TP( 4) = +0.002824; 409 ITP(0, 4)= +0; ITP(1, 4)= +2; ITP(2, 4)= +0; ITP(3, 4)= +0, ITP(4, 4)= 0; 410 TP( 5) = +0.000857; 411 ITP(0, 5)= +0; ITP(1, 5)= +1; ITP(2, 5)= +2; ITP(3, 5)= +0, ITP(4, 5)= 0; 412 TP( 6) = +0.000533; 413 ITP(0, 6)= -1; ITP(1, 6)= +0; ITP(2, 6)= +2; ITP(3, 6)= +0, ITP(4, 6)= 1; 414 TP( 7) = +0.000401; 415 ITP(0, 7)= -1; ITP(1, 7)= -1; ITP(2, 7)= +2; ITP(3, 7)= +0, ITP(4, 7)= 1; 416 TP( 8) = +0.000320; 417 ITP(0, 8)= -1; ITP(1, 8)= +1; ITP(2, 8)= +0; ITP(3, 8)= +0, ITP(4, 8)= 1; 418 TP( 9) = -0.000271; 419 ITP(0, 9)= +0; ITP(1, 9)= +0; ITP(2, 9)= +1; ITP(3, 9)= +0, ITP(4, 9)= 0; 420 TP(10) = -0.000264; 421 ITP(0,10)= +1; ITP(1,10)= +1; ITP(2,10)= +0; ITP(3,10)= +0, ITP(4,10)= 1; 422 TP(11) = -0.000198; 423 ITP(0,11)= +0; ITP(1,11)= -1; ITP(2,11)= +0; ITP(3,11)= +2, ITP(4,11)= 0; 424 TP(12) = +0.000173; 425 ITP(0,12)= +0; ITP(1,12)= +3; ITP(2,12)= +0; ITP(3,12)= +0, ITP(4,12)= 0; 426 TP(13) = +0.000167; 427 ITP(0,13)= +0; ITP(1,13)= -1; ITP(2,13)= +4; ITP(3,13)= +0, ITP(4,13)= 0; 428 TP(14) = -0.000111; 429 ITP(0,14)= +1; ITP(1,14)= +0; ITP(2,14)= +0; ITP(3,14)= +0, ITP(4,14)= 1; 430 TP(15) = +0.000103; 431 ITP(0,15)= +0; ITP(1,15)= -2; ITP(2,15)= +4; ITP(3,15)= +0, ITP(4,15)= 0; 432 TP(16) = -0.000084; 433 ITP(0,16)= +0; ITP(1,16)= +2; ITP(2,16)= -2; ITP(3,16)= +0, ITP(4,16)= 0; 434 TP(17) = -0.000083; 435 ITP(0,17)= +1; ITP(1,17)= +0; ITP(2,17)= +2; ITP(3,17)= +0, ITP(4,17)= 1; 436 TP(18) = +0.000079; 437 ITP(0,18)= +0; ITP(1,18)= +2; ITP(2,18)= +2; ITP(3,18)= +0, ITP(4,18)= 0; 438 TP(19) = +0.000072; 439 ITP(0,19)= +0; ITP(1,19)= +0; ITP(2,19)= +4; ITP(3,19)= +0, ITP(4,19)= 0; 440 TP(20) = +0.000064; 441 ITP(0,20)= -1; ITP(1,20)= +1; ITP(2,20)= +2; ITP(3,20)= +0, ITP(4,20)= 1; 442 TP(21) = -0.000063; 443 ITP(0,21)= +1; ITP(1,21)= -1; ITP(2,21)= +2; ITP(3,21)= +0, ITP(4,21)= 1; 444 TP(22) = +0.000041; 445 ITP(0,22)= +1; ITP(1,22)= +0; ITP(2,22)= +1; ITP(3,22)= +0, ITP(4,22)= 1; 446 TP(23) = +0.000035; 447 ITP(0,23)= -1; ITP(1,23)= +2; ITP(2,23)= +0; ITP(3,23)= +0, ITP(4,23)= 1; 448 TP(24) = -0.000033; 449 ITP(0,24)= +0; ITP(1,24)= +3; ITP(2,24)= -2; ITP(3,24)= +0, ITP(4,24)= 0; 450 TP(25) = -0.000030; 451 ITP(0,25)= +0; ITP(1,25)= +1; ITP(2,25)= +1; ITP(3,25)= +0, ITP(4,25)= 0; 452 TP(26) = -0.000029; 453 ITP(0,26)= +0; ITP(1,26)= +0; ITP(2,26)= -2; ITP(3,26)= +2, ITP(4,26)= 0; 454 TP(27) = -0.000029; 455 ITP(0,27)= +1; ITP(1,27)= +2; ITP(2,27)= +0; ITP(3,27)= +0, ITP(4,27)= 1; 456 TP(28) = +0.000026; 457 ITP(0,28)= -2; ITP(1,28)= +0; ITP(2,28)= +2; ITP(3,28)= +0, ITP(4,28)= 2; 458 TP(29) = -0.000023; 459 ITP(0,29)= +0; ITP(1,29)= +1; ITP(2,29)= -2; ITP(3,29)= +2, ITP(4,29)= 0; 460 TP(30) = +0.000019; 461 ITP(0,30)= -1; ITP(1,30)= -1; ITP(2,30)= +4; ITP(3,30)= +0, ITP(4,30)= 1; 462 463 464 // Centuries since J1900 465 double tt(static_cast<double>((MJD(t).mjd-15019.5)/36525.0)); 466 467 468 // Fundamental arguments (radians) and derivatives (radians per 469 // Julian century) for the current epoch 470 471 // Moon's mean longitude 472 double ELP(D2R*fmod( (ELP0+(ELP1+(ELP2+ELP3*tt)*tt)*tt), 360.0)); 473 double DELP(D2R*(ELP1+(2.0*ELP2+3.0*ELP3*tt)*tt)); 474 475 // Sun's mean anomaly 476 double EM(D2R*fmod( (EM0+(EM1+(EM2+EM3*tt)*tt)*tt), 360.0)); 477 double DEM(D2R*(EM1+(2.0*EM2+3.0*EM3*tt)*tt)); 478 479 // Moon's mean anomaly 480 double EMP(D2R*fmod( (EMP0+(EMP1+(EMP2+EMP3*tt)*tt)*tt), 360.0)); 481 double DEMP(D2R*(EMP1+(2.0*EMP2+3.0*EMP3*tt)*tt)); 482 483 // Moon's mean elongation 484 double D(D2R*fmod( (D0+(D1+(D2+D3*tt)*tt)*tt), 360.0)); 485 double DD(D2R*(D1+(2.0*D2+3.0*D3*tt)*tt)); 486 487 // Mean distance of the Moon from its ascending node 488 double F(D2R*fmod( (F0+(F1+(F2+F3*tt)*tt)*tt), 360.0)); 489 double DF(D2R*(F1+(2.0*F2+3.0*F3*tt)*tt)); 490 491 // Longitude of the Moon's ascending node 492 double OM(D2R*fmod( (OM0+(OM1+(OM2+OM3*tt)*tt)*tt), 360.0)); 493 double DOM(D2R*(OM1+(2.0*OM2+3.0*OM3*tt)*tt)); 494 double SINOM(std::sin(OM)); 495 double COSOM(std::cos(OM)); 496 double DOMCOM(DOM*COSOM); 497 498 499 // Let's add the periodic variations 500 double THETA(D2R*(PA0+PA1*tt)); 501 double WA(std::sin(THETA)); 502 double DWA(D2R*PA1*std::cos(THETA)); 503 THETA = D2R*(PE0+(PE1+PE2*tt)*tt); 504 double WB(PEC*std::sin(THETA)); 505 double DWB(D2R*PEC*(PE1+2.0*PE2*tt)*std::cos(THETA)); 506 ELP = ELP+D2R*(PAC*WA+WB+PFC*SINOM); 507 DELP = DELP+D2R*(PAC*DWA+DWB+PFC*DOMCOM); 508 EM = EM+D2R*PBC*WA; 509 DEM = DEM+D2R*PBC*DWA; 510 EMP = EMP+D2R*(PCC*WA+WB+PGC*SINOM); 511 DEMP = DEMP+D2R*(PCC*DWA+DWB+PGC*DOMCOM); 512 D = D+D2R*(PDC*WA+WB+PHC*SINOM); 513 DD = DD+D2R*(PDC*DWA+DWB+PHC*DOMCOM); 514 double WOM(OM+D2R*(PJ0+PJ1*tt)); 515 double DWOM(DOM+D2R*PJ1); 516 double SINWOM(std::sin(WOM)); 517 double COSWOM(std::cos(WOM)); 518 F = F+D2R*(WB+cPIC*SINOM+PJC*SINWOM); 519 DF = DF+D2R*(DWB+cPIC*DOMCOM+PJC*DWOM*COSWOM); 520 521 522 // E-factor, and square 523 double E(1.0+(E1+E2*tt)*tt); 524 double DE(E1+2.0*E2*tt); 525 double ESQ(E*E); 526 double DESQ(2.0*E*DE); 527 528 529 // Series expansions 530 531 // Longitude 532 double V(0.0); 533 double DV(0.0); 534 for (int n=(NL-1); n>=0; n--) 535 { 536 double EN, DEN; 537 double COEFF(TL(n)); 538 double EMN(static_cast<double>(ITL(0,n))); 539 double EMPN(static_cast<double>(ITL(1,n))); 540 double DN(static_cast<double>(ITL(2,n))); 541 double FN(static_cast<double>(ITL(3,n))); 542 int I=ITL(4,n); 543 if (I == 0) 544 { 545 EN = 1.0; 546 DEN = 0.0; 547 } 548 else 549 { 550 if (I == 1) 551 { 552 EN = E; 553 DEN = DE; 554 } 555 else 556 { 557 EN = ESQ; 558 DEN = DESQ; 559 } 560 } 561 THETA = EMN*EM+EMPN*EMP+DN*D+FN*F; 562 double DTHETA(EMN*DEM+EMPN*DEMP+DN*DD+FN*DF); 563 double FTHETA(std::sin(THETA)); 564 V = V+COEFF*FTHETA*EN; 565 DV = DV+COEFF*(std::cos(THETA)*DTHETA*EN+FTHETA*DEN); 566 } 567 double EL(ELP+D2R*V); 568 569 570 // Latitude 571 V = 0.0; 572 DV = 0.0; 573 for (int n=(NB-1); n>=0; n--) 574 { 575 double EN, DEN; 576 double COEFF(TB(n)); 577 double EMN(static_cast<double>(ITB(0,n))); 578 double EMPN(static_cast<double>(ITB(1,n))); 579 double DN(static_cast<double>(ITB(2,n))); 580 double FN(static_cast<double>(ITB(3,n))); 581 int I=ITB(4,n); 582 if (I == 0) 583 { 584 EN = 1.0; 585 DEN = 0.0; 586 } 587 else 588 { 589 if (I == 1) 590 { 591 EN = E; 592 DEN = DE; 593 } 594 else 595 { 596 EN = ESQ; 597 DEN = DESQ; 598 } 599 } 600 THETA = EMN*EM+EMPN*EMP+DN*D+FN*F; 601 double DTHETA(EMN*DEM+EMPN*DEMP+DN*DD+FN*DF); 602 double FTHETA(std::sin(THETA)); 603 V = V+COEFF*FTHETA*EN; 604 DV = DV+COEFF*(std::cos(THETA)*DTHETA*EN+FTHETA*DEN); 605 } 606 double BF(1.0-CW1*COSOM-CW2*COSWOM); 607 double B(D2R*V*BF); 608 609 610 // Parallax 611 V = 0.0; 612 DV = 0.0; 613 for (int n=(NP-1); n>=0; n--) 614 { 615 double EN, DEN; 616 double COEFF(TP(n)); 617 double EMN(static_cast<double>(ITP(0,n))); 618 double EMPN(static_cast<double>(ITP(1,n))); 619 double DN(static_cast<double>(ITP(2,n))); 620 double FN(static_cast<double>(ITP(3,n))); 621 int I=ITP(4,n); 622 if (I == 0) 623 { 624 EN = 1.0; 625 DEN = 0.0; 626 } 627 else 628 { 629 if (I == 1) 630 { 631 EN = E; 632 DEN = DE; 633 } 634 else 635 { 636 EN = ESQ; 637 DEN = DESQ; 638 } 639 } 640 THETA = EMN*EM+EMPN*EMP+DN*D+FN*F; 641 double DTHETA(EMN*DEM+EMPN*DEMP+DN*DD+FN*DF); 642 double FTHETA(std::cos(THETA)); 643 V = V+COEFF*FTHETA*EN; 644 DV = DV+COEFF*(-std::sin(THETA)*DTHETA*EN+FTHETA*DEN); 645 } 646 double P(D2R*V); 647 648 649 // Transformation into final form 650 651 // Parallax to distance (AU, AU/sec) 652 double SP(std::sin(P)); 653 double R(ERADAU/SP); 654 655 // Longitude, latitude to x,y,z (AU) 656 double SEL(std::sin(EL)); 657 double CEL(std::cos(EL)); 658 double SB(std::sin(B)); 659 double CB(std::cos(B)); 660 double RCB(R*CB); 661 double X(RCB*CEL); 662 double Y(RCB*SEL); 663 double Z(R*SB); 664 665 666 // Julian centuries since J2000 667 tt=(MJD(t).mjd-51544.5)/36525.0; 668 669 // lower-case because 1) upper case is ugly and 2) name 670 // collisions with macros. 671 672 // Fricke equinox correction 673 double epj(2000.0+tt*100.0); 674 double eqcor(DS2R*(0.035+0.00085*(epj-B1950))); 675 676 // Mean obliquity (IAU 1976) 677 double eps(DAS2R*(84381.448 + (-46.8150 + 678 (-0.00059+0.001813*tt)*tt)*tt)); 679 680 // Change to equatorial system, mean of date, FK5 system 681 double sineps(std::sin(eps)); 682 double coseps(std::cos(eps)); 683 double es(eqcor*sineps); 684 double ec(eqcor*coseps); 685 686 Triple res; 687 688 res.theArray[0] = (X-ec*Y+es*Z)*AU_CONST; 689 res.theArray[1] = (eqcor*X+Y*coseps-Z*sineps)*AU_CONST; 690 res.theArray[2] = (Y*sineps+Z*coseps)*AU_CONST; 691 692 693 return res; 694 695 } // End MoonPosition::getPositionCIS() 696 697 698 } // end namespace gpstk 699