00001 #pragma ident "$Id: GeodeticFrames.cpp 2741 2011-06-22 16:37:02Z nwu $" 00002 00003 //============================================================================ 00004 // 00005 // This file is part of GPSTk, the GPS Toolkit. 00006 // 00007 // The GPSTk is free software; you can redistribute it and/or modify 00008 // it under the terms of the GNU Lesser General Public License as published 00009 // by the Free Software Foundation; either version 2.1 of the License, or 00010 // any later version. 00011 // 00012 // The GPSTk is distributed in the hope that it will be useful, 00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 // GNU Lesser General Public License for more details. 00016 // 00017 // You should have received a copy of the GNU Lesser General Public 00018 // License along with GPSTk; if not, write to the Free Software Foundation, 00019 // Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00020 // 00021 // Copyright 2004, The University of Texas at Austin 00022 // 00023 //============================================================================ 00024 00025 //============================================================================ 00026 // 00027 //This software developed by Applied Research Laboratories at the University of 00028 //Texas at Austin, under contract to an agency or agencies within the U.S. 00029 //Department of Defense. The U.S. Government retains all rights to use, 00030 //duplicate, distribute, disclose, or release this software. 00031 // 00032 //Pursuant to DoD Directive 523024 00033 // 00034 // DISTRIBUTION STATEMENT A: This software has been approved for public 00035 // release, distribution is unlimited. 00036 // 00037 //============================================================================= 00038 00051 //------------------------------------------------------------------------------------ 00052 // GPSTk includes 00053 #include "StringUtils.hpp" 00054 #include "geometry.hpp" // for DEG_TO_RAD 00055 #include "icd_200_constants.hpp" // for PI and TWO_PI 00056 #include "GeodeticFrames.hpp" 00057 00058 using namespace std; 00059 00060 namespace gpstk 00061 { 00062 //--------------------------------------------------------------------------------- 00063 //--------------------------------------------------------------------------------- 00064 // constants 00065 // epoch for CoordTransTime 00066 const long GeodeticFrames::JulianEpoch=2451545; 00067 00068 //--------------------------------------------------------------------------------- 00069 //--------------------------------------------------------------------------------- 00070 // functions used internally 00071 00072 //--------------------------------------------------------------------------------- 00073 // Compute the 'coordinate transformation time' which is the time since 00074 // epoch J2000 = January 1 2000 12h UT = 2451545.0JD, divided by 36525 days. 00075 // This quantity is used throughout the terrestrial / inertial coordinate 00076 // transformations. 00077 double GeodeticFrames::CoordTransTime(DayTime t) 00078 throw() 00079 { 00080 return (t.JD()-double(JulianEpoch))/36525.0; 00081 } 00082 00083 //--------------------------------------------------------------------------------- 00084 // Return mean longitude of lunar ascending node, in degrees, 00085 // given T, the CoordTransTime at the Epoch of interest. (Ref: F5 pg 23) 00086 double GeodeticFrames::Omega(double T) 00087 throw() 00088 { 00089 double om; 00090 om = -0.00005939*T; // seconds of arc 00091 om = (om + 0.007702)*T; 00092 om = (om + 7.4722)*T; 00093 om /= 3600.0; // convert to degrees 00094 om = (om - 1934.136185139)*T; // 1934.136185139 = 6962890.2665/3600.0 00095 om += 125.04455501; 00096 //om = fmod(om,360.0); 00097 return om; 00098 } 00099 00100 //--------------------------------------------------------------------------------- 00101 // Return mean longitude of the moon - Omega, in degrees, 00102 // given T, the CoordTransTime at the Epoch of interest. (Ref: F3 pg 23) 00103 double GeodeticFrames::F(double T) 00104 throw() 00105 { 00106 double f; 00107 f = 0.00000417*T; // seconds of arc 00108 f = (f - 0.001037)*T; 00109 f = (f - 12.7512)*T; 00110 f /= 3600.0; // convert to degrees 00111 f = (f + 483202.01745772)*T; // 483202.01745772 = 1739527262.8478/3600.0 00112 f += 93.27209062; 00113 //f = fmod(f,360.0); 00114 return 0.0; 00115 } 00116 00117 //--------------------------------------------------------------------------------- 00118 // Return mean elongation of the moon from the sun, in degrees, 00119 // given T, the CoordTransTime at the Epoch of interest. (Ref: F4 pg 23) 00120 double GeodeticFrames::D(double T) 00121 throw() 00122 { 00123 double d; 00124 d = -0.00003169*T; // seconds of arc 00125 d = (d + 0.006593)*T; 00126 d = (d - 6.3706)*T; 00127 d /= 3600.0; // convert to degrees 00128 d = (d + 445267.111446944)*T; // 445267.111446944 = 1602961601.2090 / 3600.0 00129 d += 297.85019547; 00130 //d = fmod(d,360.0); 00131 return d; 00132 } 00133 00134 //--------------------------------------------------------------------------------- 00135 // Return mean anomaly of the moon, in degrees, 00136 // given T, the CoordTransTime at the Epoch of interest. (Ref: F1 pg 23) 00137 double GeodeticFrames::L(double T) 00138 throw() 00139 { 00140 double l; 00141 l = -0.00024470*T; // seconds of arc 00142 l = (l + 0.051635)*T; 00143 l = (l + 31.8792)*T; 00144 l /= 3600.0; // convert to degrees 00145 l = (l + 477198.8675605)*T; // 477198.8675605 = 1717915923.2178 / 3600.0 00146 l += 134.96340251; 00147 //l = fmod(l,360.0); 00148 return l; 00149 } 00150 00151 //--------------------------------------------------------------------------------- 00152 // Return mean anomaly of the sun, in degrees, 00153 // given T, the CoordTransTime at the Epoch of interest. (Ref: F2 pg 23) 00154 double GeodeticFrames::Lp(double T) 00155 throw() 00156 { 00157 double lp; 00158 lp = -0.00001149*T; // seconds of arc 00159 lp = (lp + 0.000136)*T; 00160 lp = (lp - 0.5532)*T; 00161 lp /= 3600.0; // convert to degrees 00162 lp = (lp + 35999.050291139)*T;// 35999.050291139 = 129596581.0481 / 3600.0 00163 lp += 357.52910918; 00164 //lp = fmod(lp,360.0); 00165 return lp; 00166 } 00167 00168 //--------------------------------------------------------------------------------- 00169 // Compute eps, the obliquity of the ecliptic, in degrees, 00170 // given T, the CoordTransTime at the time of interest. IAU76 00171 double GeodeticFrames::Obliquity(double T) 00172 throw() 00173 { 00174 double ep; 00175 // seconds of arc 00176 ep = T*(-46.8150 + T*(-0.00059 + T*0.001813)); 00177 ep /= 3600.0; // convert to degrees 00178 // degrees 00179 ep += 23.43929111; // = 84381.448/3600.0 00180 00181 return ep; 00182 } 00183 00184 //--------------------------------------------------------------------------------- 00185 // Nutation of the obliquity (deps) and of the longitude (dpsi), IERS 1996 00186 // model (ref pg 26), given 00187 // @param T, the coordinate transformation time at the time of interest 00188 // @param deps, nutation of the obliquity, in arc seconds (output) 00189 // @param dpsi, nutation of the longitude, in arc seconds (output) 00190 void GeodeticFrames::NutationAngles(double T, 00191 double& deps, 00192 double& dpsi) 00193 throw() 00194 { 00195 //----------------------------------------------------------------------- 00196 // Code to implement Table 5.2 of IERS Conventions 1996 series for 00197 // nutation in longitude dpsi and obliquity deps in arc seconds. 00198 // (Generated by perl script Table52.pl) 00199 double arg; 00200 dpsi = deps = 0.0; 00201 00202 // double T must be defined as time in Julian centuries from epoch J2000.0 00203 // Also define (all doubles): 00204 double o = Omega(T); // = mean longitude of lunar ascending node, in degrees, 00205 double f = F(T) ; // = mean longitude of the moon - Omega, in degrees, 00206 double d = D(T) ; // = mean elongation of the moon from the sun, in degrees, 00207 double l = L(T) ; // = mean anomaly of the moon, in degrees, and 00208 double lp = Lp(T) ; // = mean anomaly of the sun, in degrees. 00209 o *= DEG_TO_RAD; 00210 f *= DEG_TO_RAD; 00211 d *= DEG_TO_RAD; 00212 l *= DEG_TO_RAD; 00213 lp *= DEG_TO_RAD; 00214 00215 // line 1 of Table 5.2, period = -6798.38 days 00216 arg = o; 00217 dpsi += (-17.206277 - 0.017419*T) * ::sin(arg) + 0.003645 * ::cos(arg); 00218 deps += (+9.205356 + 0.000886*T) * ::cos(arg) + 0.001553 * ::sin(arg); 00219 // line 2 of Table 5.2, period = 182.62 days 00220 arg = 2*f-2*d+2*o; 00221 dpsi += (-1.317014 - 0.000156*T) * ::sin(arg) - 0.001400 * ::cos(arg); 00222 deps += (+0.573058 - 0.000306*T) * ::cos(arg) - 0.000464 * ::sin(arg); 00223 // line 3 of Table 5.2, period = 13.66 days 00224 arg = 2*f+2*o; 00225 dpsi += (-0.227720 - 0.000023*T) * ::sin(arg) + 0.000269 * ::cos(arg); 00226 deps += (+0.097864 - 0.000048*T) * ::cos(arg) + 0.000136 * ::sin(arg); 00227 // line 4 of Table 5.2, period = -3399.18 days 00228 arg = 2*o; 00229 dpsi += (+0.207429 + 0.000021*T) * ::sin(arg) - 0.000071 * ::cos(arg); 00230 deps += (-0.089747 + 0.000047*T) * ::cos(arg) - 0.000029 * ::sin(arg); 00231 // line 5 of Table 5.2, period = -365.26 days 00232 arg = -lp; 00233 dpsi += (-0.147538 + 0.000364*T) * ::sin(arg) + 0.001121 * ::cos(arg); 00234 deps += (+0.007388 - 0.000019*T) * ::cos(arg) + 0.000198 * ::sin(arg); 00235 // line 6 of Table 5.2, period = 121.75 days 00236 arg = lp+2*f-2*d+2*o; 00237 dpsi += (-0.051687 + 0.000123*T) * ::sin(arg) - 0.000054 * ::cos(arg); 00238 deps += (+0.022440 - 0.000068*T) * ::cos(arg) - 0.000018 * ::sin(arg); 00239 // line 7 of Table 5.2, period = 27.55 days 00240 arg = l; 00241 dpsi += (+0.071118 + 0.000007*T) * ::sin(arg) - 0.000094 * ::cos(arg); 00242 deps -= 0.000687 * ::cos(arg) + 0.000039 * ::sin(arg); 00243 // line 8 of Table 5.2, period = 13.63 days 00244 arg = 2*f+o; 00245 dpsi += (-0.038752 - 0.000037*T) * ::sin(arg) + 0.000034 * ::cos(arg); 00246 deps += (+0.020076 + 0.000002*T) * ::cos(arg) + 0.000032 * ::sin(arg); 00247 // line 9 of Table 5.2, period = 9.13 days 00248 arg = l+2*f+2*o; 00249 dpsi += (-0.030137 - 0.000004*T) * ::sin(arg) + 0.000077 * ::cos(arg); 00250 deps += (+0.012896 - 0.000006*T) * ::cos(arg) + 0.000035 * ::sin(arg); 00251 // line 10 of Table 5.2, period = 365.22 days 00252 arg = -lp+2*f-2*d+2*o; 00253 dpsi += (+0.021583 - 0.000049*T) * ::sin(arg) + 0.000006 * ::cos(arg); 00254 deps += (-0.009591 + 0.000030*T) * ::cos(arg) + 0.000012 * ::sin(arg); 00255 // line 11 of Table 5.2, period = 177.84 days 00256 arg = 2*f-2*d+o; 00257 dpsi += (+0.012820 + 0.000014*T) * ::sin(arg) + 0.000018 * ::cos(arg); 00258 deps += (-0.006897 - 0.000001*T) * ::cos(arg) + 0.000004 * ::sin(arg); 00259 // line 12 of Table 5.2, period = 27.09 days 00260 arg = -l+2*f+2*o; 00261 dpsi += (+0.012353 + 0.000001*T) * ::sin(arg) + 0.000002 * ::cos(arg); 00262 deps += (-0.005334 + 0.000003*T) * ::cos(arg); 00263 // line 13 of Table 5.2, period = 31.81 days 00264 arg = -l+2*d; 00265 dpsi += (+0.015699 + 0.000001*T) * ::sin(arg) - 0.000018 * ::cos(arg); 00266 deps -= 0.000127 * ::cos(arg) + 0.000009 * ::sin(arg); 00267 // line 14 of Table 5.2, period = 27.67 days 00268 arg = l+o; 00269 dpsi += (+0.006314 + 0.000006*T) * ::sin(arg) + 0.000003 * ::cos(arg); 00270 deps -= 0.003323 * ::cos(arg) - 0.000001 * ::sin(arg); 00271 // line 15 of Table 5.2, period = -27.44 days 00272 arg = -l+o; 00273 dpsi += (-0.005797 - 0.000006*T) * ::sin(arg) - 0.000019 * ::cos(arg); 00274 deps += 0.003141 * ::cos(arg) - 0.000008 * ::sin(arg); 00275 // line 16 of Table 5.2, period = 9.56 days 00276 arg = -l+2*f+2*d+2*o; 00277 dpsi += (-0.005965 - 0.000001*T) * ::sin(arg) + 0.000014 * ::cos(arg); 00278 deps += (+0.002554 - 0.000001*T) * ::cos(arg) + 0.000007 * ::sin(arg); 00279 // line 17 of Table 5.2, period = 9.12 days 00280 arg = l+2*f+o; 00281 dpsi += (-0.005163 - 0.000004*T) * ::sin(arg) + 0.000012 * ::cos(arg); 00282 deps += 0.002635 * ::cos(arg) + 0.000008 * ::sin(arg); 00283 // line 18 of Table 5.2, period = 1305.48 days 00284 arg = -2*l+2*f+o; 00285 dpsi += (+0.004590 + 0.000005*T) * ::sin(arg) + 0.000001 * ::cos(arg); 00286 deps += (-0.002424 - 0.000001*T) * ::cos(arg) + 0.000001 * ::sin(arg); 00287 // line 19 of Table 5.2, period = 14.77 days 00288 arg = 2*d; 00289 dpsi += (+0.006336 + 0.000001*T) * ::sin(arg) - 0.000015 * ::cos(arg); 00290 deps -= 0.000125 * ::cos(arg) + 0.000003 * ::sin(arg); 00291 // line 20 of Table 5.2, period = 7.10 days 00292 arg = 2*f+2*d+2*o; 00293 dpsi -= 0.003854 * ::sin(arg) + 0.000015 * ::cos(arg); 00294 deps += 0.001643 * ::cos(arg) + 0.000006 * ::sin(arg); 00295 // line 21 of Table 5.2, period = -205.89 days 00296 arg = -2*l+2*d; 00297 dpsi -= 0.004774 * ::sin(arg) - 0.000002 * ::cos(arg); 00298 deps += 0.000048 * ::cos(arg) - 0.000003 * ::sin(arg); 00299 // line 22 of Table 5.2, period = 6.86 days 00300 arg = 2*l+2*f+2*o; 00301 dpsi -= 0.003102 * ::sin(arg) + 0.000012 * ::cos(arg); 00302 deps += (+0.001323 - 0.000001*T) * ::cos(arg) + 0.000005 * ::sin(arg); 00303 // line 23 of Table 5.2, period = 23.94 days 00304 arg = l+2*f-2*d+2*o; 00305 dpsi += 0.002863 * ::sin(arg); 00306 deps += (-0.001235 + 0.000001*T) * ::cos(arg); 00307 // line 24 of Table 5.2, period = 26.98 days 00308 arg = -l+2*f+o; 00309 dpsi += (+0.002044 + 0.000002*T) * ::sin(arg) + 0.000001 * ::cos(arg); 00310 deps -= 0.001076 * ::cos(arg); 00311 // line 25 of Table 5.2, period = 13.78 days 00312 arg = 2*l; 00313 dpsi += 0.002923 * ::sin(arg) - 0.000008 * ::cos(arg); 00314 deps -= 0.000062 * ::cos(arg) + 0.000001 * ::sin(arg); 00315 // line 26 of Table 5.2, period = 13.61 days 00316 arg = 2*f; 00317 dpsi += 0.002585 * ::sin(arg) - 0.000007 * ::cos(arg); 00318 deps -= 0.000056 * ::cos(arg) + 0.000001 * ::sin(arg); 00319 // line 27 of Table 5.2, period = 386.00 days 00320 arg = lp+o; 00321 dpsi += (-0.001406 - 0.000003*T) * ::sin(arg) + 0.000008 * ::cos(arg); 00322 deps += 0.000857 * ::cos(arg) - 0.000004 * ::sin(arg); 00323 // line 28 of Table 5.2, period = 31.96 days 00324 arg = -l+2*d+o; 00325 dpsi += (+0.001517 + 0.000001*T) * ::sin(arg) + 0.000001 * ::cos(arg); 00326 deps -= 0.000801 * ::cos(arg); 00327 // line 29 of Table 5.2, period = 91.31 days 00328 arg = 2*lp+2*f-2*d+2*o; 00329 dpsi += (-0.001578 + 0.000007*T) * ::sin(arg) - 0.000002 * ::cos(arg); 00330 deps += (+0.000685 - 0.000004*T) * ::cos(arg) - 0.000001 * ::sin(arg); 00331 // line 30 of Table 5.2, period = -173.31 days 00332 arg = -2*f+2*d; 00333 dpsi += 0.002178 * ::sin(arg) + 0.000001 * ::cos(arg); 00334 deps -= 0.000015 * ::cos(arg) + 0.000001 * ::sin(arg); 00335 // line 31 of Table 5.2, period = -31.66 days 00336 arg = l-2*d+o; 00337 dpsi += (-0.001286 - 0.000001*T) * ::sin(arg) - 0.000004 * ::cos(arg); 00338 deps += 0.000694 * ::cos(arg) - 0.000002 * ::sin(arg); 00339 // line 32 of Table 5.2, period = -346.64 days 00340 arg = -lp+o; 00341 dpsi += (-0.001269 + 0.000001*T) * ::sin(arg) + 0.000006 * ::cos(arg); 00342 deps += (+0.000642 + 0.000001*T) * ::cos(arg) + 0.000002 * ::sin(arg); 00343 // line 33 of Table 5.2, period = 9.54 days 00344 arg = -l+2*f+2*d+o; 00345 dpsi += (-0.001022 - 0.000001*T) * ::sin(arg) + 0.000002 * ::cos(arg); 00346 deps += 0.000522 * ::cos(arg) + 0.000001 * ::sin(arg); 00347 // line 34 of Table 5.2, period = -182.63 days 00348 arg = -2*lp; 00349 dpsi += (-0.001671 + 0.000008*T) * ::sin(arg) - 0.000001 * ::cos(arg); 00350 deps += 0.000014 * ::cos(arg) - 0.000001 * ::sin(arg); 00351 // line 35 of Table 5.2, period = 5.64 days 00352 arg = l+2*f+2*d+2*o; 00353 dpsi -= 0.000768 * ::sin(arg) + 0.000004 * ::cos(arg); 00354 deps += 0.000325 * ::cos(arg) + 0.000002 * ::sin(arg); 00355 // line 36 of Table 5.2, period = 1095.18 days 00356 arg = -2*l+2*f; 00357 dpsi -= 0.001102 * ::sin(arg) - 0.000001 * ::cos(arg); 00358 deps += 0.000010 * ::cos(arg); 00359 /* 00360 // line 37 of Table 5.2, period = 13.17 days 00361 arg = lp+2*f+2*o; 00362 dpsi += (+0.000757 - 0.000002*T) * ::sin(arg) - 0.000001 * ::cos(arg); 00363 deps += (-0.000326 - 0.000002*T) * ::cos(arg); 00364 // line 38 of Table 5.2, period = 7.09 days 00365 arg = 2*f+2*d+o; 00366 dpsi += (-0.000664 - 0.000001*T) * ::sin(arg) + 0.000002 * ::cos(arg); 00367 deps += (+0.000335 - 0.000001*T) * ::cos(arg) + 0.000001 * ::sin(arg); 00368 // line 39 of Table 5.2, period = 14.19 days 00369 arg = -lp+2*f+2*o; 00370 dpsi += (-0.000714 + 0.000002*T) * ::sin(arg) + 0.000001 * ::cos(arg); 00371 deps += (+0.000307 + 0.000002*T) * ::cos(arg); 00372 // line 40 of Table 5.2, period = 14.80 days 00373 arg = 2*d+o; 00374 dpsi += (-0.000631 - 0.000001*T) * ::sin(arg); 00375 deps += 0.000327 * ::cos(arg); 00376 // line 41 of Table 5.2, period = 23.86 days 00377 arg = l+2*f-2*d+o; 00378 dpsi += (+0.000580 + 0.000001*T) * ::sin(arg); 00379 deps -= 0.000307 * ::cos(arg); 00380 // line 42 of Table 5.2, period = 12.81 days 00381 arg = 2*l+2*f-2*d+2*o; 00382 dpsi += 0.000643 * ::sin(arg) - 0.000001 * ::cos(arg); 00383 deps -= 0.000277 * ::cos(arg); 00384 // line 43 of Table 5.2, period = -199.84 days 00385 arg = -2*l+2*d+o; 00386 dpsi += (-0.000579 - 0.000001*T) * ::sin(arg) - 0.000001 * ::cos(arg); 00387 deps += 0.000304 * ::cos(arg); 00388 // line 44 of Table 5.2, period = 6.85 days 00389 arg = 2*l+2*f+o; 00390 dpsi -= 0.000533 * ::sin(arg) + 0.000002 * ::cos(arg); 00391 deps += 0.000269 * ::cos(arg) + 0.000001 * ::sin(arg); 00392 // line 45 of Table 5.2, period = 346.60 days 00393 arg = -lp+2*f-2*d+o; 00394 dpsi += (-0.000477 - 0.000001*T) * ::sin(arg); 00395 deps += (+0.000271 - 0.000001*T) * ::cos(arg); 00396 // line 46 of Table 5.2, period = -14.73 days 00397 arg = -2*d+o; 00398 dpsi += (-0.000493 - 0.000001*T) * ::sin(arg) - 0.000002 * ::cos(arg); 00399 deps += 0.000272 * ::cos(arg) - 0.000001 * ::sin(arg); 00400 // line 47 of Table 5.2, period = 34.85 days 00401 arg = -l-lp+2*d; 00402 dpsi += 0.000735 * ::sin(arg) - 0.000001 * ::cos(arg); 00403 deps -= 0.000005 * ::cos(arg); 00404 // line 48 of Table 5.2, period = 212.32 days 00405 arg = 2*l-2*d+o; 00406 dpsi += 0.000405 * ::sin(arg) + 0.000001 * ::cos(arg); 00407 deps -= 0.000220 * ::cos(arg); 00408 // line 49 of Table 5.2, period = 9.61 days 00409 arg = l+2*d; 00410 dpsi += 0.000657 * ::sin(arg) - 0.000002 * ::cos(arg); 00411 deps -= 0.000020 * ::cos(arg); 00412 // line 50 of Table 5.2, period = 119.61 days 00413 arg = lp+2*f-2*d+o; 00414 dpsi += 0.000361 * ::sin(arg) + 0.000001 * ::cos(arg); 00415 deps -= 0.000194 * ::cos(arg); 00416 // line 51 of Table 5.2, period = 29.80 days 00417 arg = l-lp; 00418 dpsi += 0.000471 * ::sin(arg) - 0.000001 * ::cos(arg); 00419 deps -= 0.000004 * ::cos(arg); 00420 // line 52 of Table 5.2, period = 1615.76 days 00421 arg = -2*l+2*f+2*o; 00422 dpsi -= 0.000311 * ::sin(arg); 00423 deps += 0.000131 * ::cos(arg); 00424 // line 53 of Table 5.2, period = 5.49 days 00425 arg = 3*l+2*f+2*o; 00426 dpsi -= 0.000289 * ::sin(arg) + 0.000002 * ::cos(arg); 00427 deps += 0.000124 * ::cos(arg) + 0.000001 * ::sin(arg); 00428 // line 54 of Table 5.2, period = 15.39 days 00429 arg = -lp+2*d; 00430 dpsi += 0.000435 * ::sin(arg) - 0.000001 * ::cos(arg); 00431 deps -= 0.000009 * ::cos(arg); 00432 // line 55 of Table 5.2, period = 9.37 days 00433 arg = l-lp+2*f+2*o; 00434 dpsi -= 0.000287 * ::sin(arg) + 0.000001 * ::cos(arg); 00435 deps += 0.000123 * ::cos(arg); 00436 // line 56 of Table 5.2, period = 9.81 days 00437 arg = -l-lp+2*f+2*d+2*o; 00438 dpsi -= 0.000282 * ::sin(arg) + 0.000001 * ::cos(arg); 00439 deps += 0.000122 * ::cos(arg); 00440 // line 57 of Table 5.2, period = 29.53 days 00441 arg = d; 00442 dpsi -= 0.000422 * ::sin(arg) + 0.000001 * ::cos(arg); 00443 deps += 0.000003 * ::cos(arg); 00444 // line 58 of Table 5.2, period = 26.88 days 00445 arg = -l+2*f; 00446 dpsi -= 0.000404 * ::sin(arg) + 0.000001 * ::cos(arg); 00447 deps += 0.000004 * ::cos(arg); 00448 // line 59 of Table 5.2, period = 7.24 days 00449 arg = -lp+2*f+2*d+2*o; 00450 dpsi -= 0.000264 * ::sin(arg) + 0.000001 * ::cos(arg); 00451 deps += 0.000114 * ::cos(arg); 00452 // line 60 of Table 5.2, period = -13.75 days 00453 arg = -2*l+o; 00454 dpsi -= 0.000228 * ::sin(arg) - 0.000001 * ::cos(arg); 00455 deps += 0.000126 * ::cos(arg); 00456 // line 61 of Table 5.2, period = 8.91 days 00457 arg = l+lp+2*f+2*o; 00458 dpsi += 0.000246 * ::sin(arg) - 0.000001 * ::cos(arg); 00459 deps -= 0.000106 * ::cos(arg); 00460 // line 62 of Table 5.2, period = 13.81 days 00461 arg = 2*l+o; 00462 dpsi += 0.000218 * ::sin(arg); 00463 deps -= 0.000114 * ::cos(arg); 00464 // line 63 of Table 5.2, period = 3232.87 days 00465 arg = -l+lp+d; 00466 dpsi += 0.000327 * ::sin(arg); 00467 deps -= 0.000001 * ::cos(arg); 00468 // line 64 of Table 5.2, period = 25.62 days 00469 arg = l+lp; 00470 dpsi -= 0.000338 * ::sin(arg); 00471 deps += 0.000004 * ::cos(arg); 00472 // line 65 of Table 5.2, period = 9.11 days 00473 arg = l+2*f; 00474 dpsi += 0.000334 * ::sin(arg) - 0.000001 * ::cos(arg); 00475 deps -= 0.000011 * ::cos(arg); 00476 // line 66 of Table 5.2, period = -32.61 days 00477 arg = -l+2*f-2*d+o; 00478 dpsi -= 0.000199 * ::sin(arg) - 0.000001 * ::cos(arg); 00479 deps += 0.000107 * ::cos(arg); 00480 // line 67 of Table 5.2, period = 27.78 days 00481 arg = l+2*o; 00482 dpsi -= 0.000197 * ::sin(arg); 00483 deps += 0.000085 * ::cos(arg); 00484 // line 68 of Table 5.2, period = -411.78 days 00485 arg = -l+d; 00486 dpsi += 0.000405 * ::sin(arg) - 0.000035 * ::cos(arg); 00487 deps -= 0.000055 * ::cos(arg) - 0.000014 * ::sin(arg); 00488 // line 69 of Table 5.2, period = 9.34 days 00489 arg = 2*f+d+2*o; 00490 dpsi += 0.000165 * ::sin(arg); 00491 deps -= 0.000072 * ::cos(arg); 00492 // line 70 of Table 5.2, period = 5.80 days 00493 arg = -l+2*f+4*d+2*o; 00494 dpsi -= 0.000151 * ::sin(arg) + 0.000001 * ::cos(arg); 00495 deps += 0.000066 * ::cos(arg); 00496 // line 71 of Table 5.2, period = 6786.31 days 00497 arg = -2*lp+2*f-2*d+o; 00498 dpsi -= 0.000130 * ::sin(arg); 00499 deps += 0.000069 * ::cos(arg); 00500 // line 72 of Table 5.2, period = 6164.17 days 00501 arg = -l+lp+d+o; 00502 dpsi += 0.000132 * ::sin(arg); 00503 deps -= 0.000068 * ::cos(arg); 00504 // line 73 of Table 5.2, period = 5.64 days 00505 arg = l+2*f+2*d+o; 00506 dpsi -= 0.000133 * ::sin(arg) + 0.000001 * ::cos(arg); 00507 deps += 0.000066 * ::cos(arg); 00508 // line 74 of Table 5.2, period = 14.63 days 00509 arg = -2*l+2*f+2*d+2*o; 00510 dpsi += 0.000139 * ::sin(arg); 00511 deps -= 0.000060 * ::cos(arg); 00512 // line 75 of Table 5.2, period = -27.33 days 00513 arg = -l+2*o; 00514 dpsi += 0.000139 * ::sin(arg); 00515 deps -= 0.000060 * ::cos(arg); 00516 // line 76 of Table 5.2, period = 22.47 days 00517 arg = l+lp+2*f-2*d+2*o; 00518 dpsi += 0.000128 * ::sin(arg); 00519 deps -= 0.000055 * ::cos(arg); 00520 // line 77 of Table 5.2, period = 7.35 days 00521 arg = -2*l+2*f+4*d+2*o; 00522 dpsi -= 0.000121 * ::sin(arg); 00523 deps += 0.000052 * ::cos(arg); 00524 // line 78 of Table 5.2, period = 9.06 days 00525 arg = -l+4*f+2*o; 00526 dpsi += 0.000115 * ::sin(arg); 00527 deps -= 0.000049 * ::cos(arg); 00528 // line 79 of Table 5.2, period = 12.79 days 00529 arg = 2*l+2*f-2*d+o; 00530 dpsi += 0.000101 * ::sin(arg); 00531 deps -= 0.000054 * ::cos(arg); 00532 // line 80 of Table 5.2, period = 4.68 days 00533 arg = 2*l+2*f+2*d+2*o; 00534 dpsi -= 0.000108 * ::sin(arg) + 0.000001 * ::cos(arg); 00535 deps += 0.000047 * ::cos(arg); 00536 // line 81 of Table 5.2, period = 9.63 days 00537 arg = l+2*d+o; 00538 dpsi -= 0.000095 * ::sin(arg); 00539 deps += 0.000049 * ::cos(arg); 00540 // line 82 of Table 5.2, period = 9.18 days 00541 arg = 3*l; 00542 dpsi += 0.000157 * ::sin(arg) - 0.000001 * ::cos(arg); 00543 deps -= 0.000005 * ::cos(arg); 00544 // line 83 of Table 5.2, period = 8.75 days 00545 arg = 3*l+2*f-2*d+2*o; 00546 dpsi += 0.000094 * ::sin(arg); 00547 deps -= 0.000040 * ::cos(arg); 00548 // line 84 of Table 5.2, period = 12.66 days 00549 arg = 4*f-2*d+2*o; 00550 dpsi += 0.000091 * ::sin(arg); 00551 deps -= 0.000039 * ::cos(arg); 00552 // line 85 of Table 5.2, period = -169.00 days 00553 arg = -2*f+2*d+o; 00554 dpsi += 0.000087 * ::sin(arg); 00555 deps -= 0.000044 * ::cos(arg); 00556 // line 86 of Table 5.2, period = 13.14 days 00557 arg = lp+2*f+o; 00558 dpsi += 0.000081 * ::sin(arg); 00559 deps -= 0.000042 * ::cos(arg); 00560 // line 87 of Table 5.2, period = 187.66 days 00561 arg = 2*f-2*d+3*o; 00562 dpsi += 0.000123 * ::sin(arg); 00563 deps -= 0.000020 * ::cos(arg); 00564 // line 88 of Table 5.2, period = 10.08 days 00565 arg = -l+4*d; 00566 dpsi += 0.000133 * ::sin(arg); 00567 deps -= 0.000004 * ::cos(arg); 00568 // line 89 of Table 5.2, period = -943.23 days 00569 arg = 2*l-2*f+o; 00570 dpsi += 0.000071 * ::sin(arg); 00571 deps -= 0.000038 * ::cos(arg); 00572 // line 90 of Table 5.2, period = -15.91 days 00573 arg = 2*l-4*d; 00574 dpsi -= 0.000128 * ::sin(arg); 00575 deps += 0.000001 * ::cos(arg); 00576 // line 91 of Table 5.2, period = 35.03 days 00577 arg = -l-lp+2*d+o; 00578 dpsi += 0.000075 * ::sin(arg); 00579 deps -= 0.000039 * ::cos(arg); 00580 // line 92 of Table 5.2, period = -131.67 days 00581 arg = -2*l-lp+2*d; 00582 dpsi -= 0.000115 * ::sin(arg); 00583 deps += 0.000001 * ::cos(arg); 00584 // line 93 of Table 5.2, period = 14.16 days 00585 arg = -lp+2*f+o; 00586 dpsi -= 0.000066 * ::sin(arg); 00587 deps += 0.000035 * ::cos(arg); 00588 // line 94 of Table 5.2, period = -388.27 days 00589 arg = -l+d+o; 00590 dpsi += 0.000101 * ::sin(arg) - 0.000003 * ::cos(arg); 00591 deps -= 0.000049 * ::cos(arg) - 0.000001 * ::sin(arg); 00592 // line 95 of Table 5.2, period = -13.58 days 00593 arg = -2*f+o; 00594 dpsi -= 0.000068 * ::sin(arg); 00595 deps += 0.000036 * ::cos(arg); 00596 // line 96 of Table 5.2, period = 409.23 days 00597 arg = lp+2*o; 00598 dpsi += 0.000069 * ::sin(arg) - 0.000001 * ::cos(arg); 00599 deps -= 0.000033 * ::cos(arg); 00600 // line 97 of Table 5.2, period = 25.42 days 00601 arg = 2*f-d+2*o; 00602 dpsi -= 0.000074 * ::sin(arg); 00603 deps += 0.000031 * ::cos(arg); 00604 // line 98 of Table 5.2, period = 4.79 days 00605 arg = 2*f+4*d+2*o; 00606 dpsi -= 0.000069 * ::sin(arg); 00607 deps += 0.000029 * ::cos(arg); 00608 // line 99 of Table 5.2, period = -34.67 days 00609 arg = l+lp-2*d+o; 00610 dpsi -= 0.000061 * ::sin(arg); 00611 deps += 0.000032 * ::cos(arg); 00612 // line 100 of Table 5.2, period = 29.26 days 00613 arg = -l+lp+2*d; 00614 dpsi -= 0.000094 * ::sin(arg); 00615 // line 101 of Table 5.2, period = 5.73 days 00616 arg = l-lp+2*f+2*d+2*o; 00617 dpsi -= 0.000059 * ::sin(arg); 00618 deps += 0.000025 * ::cos(arg); 00619 // line 102 of Table 5.2, period = 29.93 days 00620 arg = l-lp+o; 00621 dpsi += 0.000051 * ::sin(arg); 00622 deps -= 0.000027 * ::cos(arg); 00623 // line 103 of Table 5.2, period = -329.79 days 00624 arg = lp-2*f+2*d; 00625 dpsi -= 0.000090 * ::sin(arg); 00626 deps += 0.000003 * ::cos(arg); 00627 // line 104 of Table 5.2, period = 5.49 days 00628 arg = 3*l+2*f+o; 00629 dpsi -= 0.000050 * ::sin(arg); 00630 deps += 0.000025 * ::cos(arg); 00631 // line 105 of Table 5.2, period = 9.31 days 00632 arg = -l+lp+2*f+2*d+2*o; 00633 dpsi += 0.000056 * ::sin(arg); 00634 deps -= 0.000024 * ::cos(arg); 00635 // line 106 of Table 5.2, period = 6.96 days 00636 arg = lp+2*f+2*d+2*o; 00637 dpsi += 0.000054 * ::sin(arg); 00638 deps -= 0.000022 * ::cos(arg); 00639 // line 107 of Table 5.2, period = -9.60 days 00640 arg = -l-2*d+o; 00641 dpsi -= 0.000050 * ::sin(arg); 00642 deps += 0.000027 * ::cos(arg); 00643 // line 108 of Table 5.2, period = 66079.30 days 00644 arg = -l+lp+d+2*o; 00645 dpsi -= 0.000052 * ::sin(arg); 00646 deps += 0.000023 * ::cos(arg); 00647 // line 109 of Table 5.2, period = 7.23 days 00648 arg = -lp+2*f+2*d+o; 00649 dpsi -= 0.000044 * ::sin(arg); 00650 deps += 0.000024 * ::cos(arg); 00651 // line 110 of Table 5.2, period = -38.74 days 00652 arg = l+2*f-4*d+o; 00653 dpsi -= 0.000047 * ::sin(arg); 00654 deps += 0.000024 * ::cos(arg); 00655 // line 111 of Table 5.2, period = -23.77 days 00656 arg = -l-2*f+2*d; 00657 dpsi += 0.000077 * ::sin(arg); 00658 // line 112 of Table 5.2, period = 9.80 days 00659 arg = -l-lp+2*f+2*d+o; 00660 dpsi -= 0.000046 * ::sin(arg); 00661 deps += 0.000024 * ::cos(arg); 00662 // line 113 of Table 5.2, period = -329.82 days 00663 arg = -lp+2*o; 00664 dpsi += 0.000059 * ::sin(arg); 00665 deps -= 0.000025 * ::cos(arg); 00666 // line 114 of Table 5.2, period = 6.99 days 00667 arg = 2*l-lp+2*f+2*o; 00668 dpsi -= 0.000048 * ::sin(arg); 00669 deps += 0.000021 * ::cos(arg); 00670 // line 115 of Table 5.2, period = 9.35 days 00671 arg = l-lp+2*f+o; 00672 dpsi -= 0.000042 * ::sin(arg); 00673 deps += 0.000022 * ::cos(arg); 00674 // line 116 of Table 5.2, period = 14.83 days 00675 arg = 2*d+2*o; 00676 dpsi -= 0.000046 * ::sin(arg); 00677 deps += 0.000020 * ::cos(arg); 00678 // line 117 of Table 5.2, period = 14.19 days 00679 arg = lp+2*d; 00680 dpsi -= 0.000067 * ::sin(arg); 00681 // line 118 of Table 5.2, period = 25.22 days 00682 arg = -l+lp+2*f+2*o; 00683 dpsi += 0.000047 * ::sin(arg); 00684 deps -= 0.000020 * ::cos(arg); 00685 // line 119 of Table 5.2, period = 73.05 days 00686 arg = 3*lp+2*f-2*d+2*o; 00687 dpsi -= 0.000044 * ::sin(arg); 00688 deps += 0.000019 * ::cos(arg); 00689 // line 120 of Table 5.2, period = -117.54 days 00690 arg = -lp-2*f+2*d; 00691 dpsi += 0.000066 * ::sin(arg); 00692 // line 121 of Table 5.2, period = 29.66 days 00693 arg = d+o; 00694 dpsi -= 0.000037 * ::sin(arg); 00695 deps += 0.000020 * ::cos(arg); 00696 // line 122 of Table 5.2, period = -9.53 days 00697 arg = l-2*f-2*d; 00698 dpsi -= 0.000064 * ::sin(arg); 00699 deps += 0.000001 * ::cos(arg); 00700 // line 123 of Table 5.2, period = 8.90 days 00701 arg = l+lp+2*f+o; 00702 dpsi += 0.000036 * ::sin(arg); 00703 deps -= 0.000018 * ::cos(arg); 00704 // line 124 of Table 5.2, period = 6.73 days 00705 arg = 2*l+lp+2*f+2*o; 00706 dpsi += 0.000040 * ::sin(arg); 00707 deps -= 0.000017 * ::cos(arg); 00708 // line 125 of Table 5.2, period = 27.32 days 00709 arg = lp+d; 00710 dpsi += 0.000057 * ::sin(arg); 00711 // line 126 of Table 5.2, period = 32.76 days 00712 arg = l-2*f+2*d; 00713 dpsi -= 0.000058 * ::sin(arg); 00714 // line 127 of Table 5.2, period = 25.72 days 00715 arg = l+lp+o; 00716 dpsi -= 0.000034 * ::sin(arg); 00717 deps += 0.000019 * ::cos(arg); 00718 // line 128 of Table 5.2, period = -7.13 days 00719 arg = -2*l-2*d; 00720 dpsi -= 0.000059 * ::sin(arg); 00721 deps += 0.000001 * ::cos(arg); 00722 // line 129 of Table 5.2, period = 32.11 days 00723 arg = -l+2*d+2*o; 00724 dpsi -= 0.000038 * ::sin(arg); 00725 deps += 0.000017 * ::cos(arg); 00726 // line 130 of Table 5.2, period = -29.40 days 00727 arg = -d+o; 00728 dpsi += 0.000033 * ::sin(arg); 00729 deps -= 0.000018 * ::cos(arg); 00730 // line 131 of Table 5.2, period = -15.35 days 00731 arg = lp-2*d+o; 00732 dpsi -= 0.000033 * ::sin(arg); 00733 deps += 0.000018 * ::cos(arg); 00734 // line 132 of Table 5.2, period = -32.45 days 00735 arg = -l+2*f-2*d+2*o; 00736 dpsi += 0.000036 * ::sin(arg); 00737 deps -= 0.000016 * ::cos(arg); 00738 // line 133 of Table 5.2, period = -29.67 days 00739 arg = -l+lp+o; 00740 dpsi -= 0.000031 * ::sin(arg); 00741 deps += 0.000017 * ::cos(arg); 00742 // line 134 of Table 5.2, period = 6.98 days 00743 arg = l+2*f+d+2*o; 00744 dpsi += 0.000033 * ::sin(arg); 00745 deps -= 0.000014 * ::cos(arg); 00746 // line 135 of Table 5.2, period = -7.38 days 00747 arg = -4*d; 00748 dpsi -= 0.000048 * ::sin(arg); 00749 deps += 0.000001 * ::cos(arg); 00750 // line 136 of Table 5.2, period = 9.33 days 00751 arg = 2*f+d+o; 00752 dpsi += 0.000027 * ::sin(arg); 00753 deps -= 0.000014 * ::cos(arg); 00754 // line 137 of Table 5.2, period = -31.52 days 00755 arg = l-2*d+2*o; 00756 dpsi += 0.000032 * ::sin(arg); 00757 deps -= 0.000014 * ::cos(arg); 00758 // line 138 of Table 5.2, period = 13.22 days 00759 arg = l+2*f-d+2*o; 00760 dpsi -= 0.000033 * ::sin(arg); 00761 deps += 0.000013 * ::cos(arg); 00762 // line 139 of Table 5.2, period = 9.87 days 00763 arg = l-lp+2*d; 00764 dpsi += 0.000048 * ::sin(arg); 00765 // line 140 of Table 5.2, period = 5.80 days 00766 arg = -l+2*f+4*d+o; 00767 dpsi -= 0.000026 * ::sin(arg); 00768 deps += 0.000013 * ::cos(arg); 00769 // line 141 of Table 5.2, period = -7.08 days 00770 arg = -2*f-2*d; 00771 dpsi -= 0.000041 * ::sin(arg); 00772 deps += 0.000001 * ::cos(arg); 00773 // line 142 of Table 5.2, period = -26.77 days 00774 arg = l-2*f+o; 00775 dpsi += 0.000027 * ::sin(arg); 00776 deps -= 0.000014 * ::cos(arg); 00777 // line 143 of Table 5.2, period = 313.04 days 00778 arg = -l+2*f-d+o; 00779 dpsi -= 0.000023 * ::sin(arg); 00780 deps += 0.000014 * ::cos(arg); 00781 // line 144 of Table 5.2, period = 22.40 days 00782 arg = l+lp+2*f-2*d+o; 00783 dpsi += 0.000023 * ::sin(arg); 00784 deps -= 0.000012 * ::cos(arg); 00785 // line 145 of Table 5.2, period = 4.58 days 00786 arg = 4*l+2*f+2*o; 00787 dpsi -= 0.000026 * ::sin(arg); 00788 deps += 0.000011 * ::cos(arg); 00789 // line 146 of Table 5.2, period = 9.11 days 00790 arg = lp+2*f+d+2*o; 00791 dpsi -= 0.000024 * ::sin(arg); 00792 deps += 0.000010 * ::cos(arg); 00793 // line 147 of Table 5.2, period = -6.85 days 00794 arg = -2*l-2*f; 00795 dpsi -= 0.000036 * ::sin(arg); 00796 deps += 0.000001 * ::cos(arg); 00797 // line 148 of Table 5.2, period = 12.38 days 00798 arg = 2*l+lp+2*f-2*d+2*o; 00799 dpsi += 0.000025 * ::sin(arg); 00800 deps -= 0.000010 * ::cos(arg); 00801 // line 149 of Table 5.2, period = 14.32 days 00802 arg = 2*l-lp; 00803 dpsi += 0.000038 * ::sin(arg); 00804 // line 150 of Table 5.2, period = -25.53 days 00805 arg = -l-lp+o; 00806 dpsi += 0.000021 * ::sin(arg); 00807 deps -= 0.000012 * ::cos(arg); 00808 // line 151 of Table 5.2, period = 14.60 days 00809 arg = -2*l+2*f+2*d+o; 00810 dpsi += 0.000022 * ::sin(arg); 00811 deps -= 0.000011 * ::cos(arg); 00812 // line 152 of Table 5.2, period = -2266.12 days 00813 arg = 3*o; 00814 dpsi -= 0.000022 * ::sin(arg); 00815 deps += 0.000010 * ::cos(arg); 00816 // line 153 of Table 5.2, period = 8.68 days 00817 arg = l+4*f-2*d+2*o; 00818 dpsi += 0.000023 * ::sin(arg); 00819 deps -= 0.000009 * ::cos(arg); 00820 // line 154 of Table 5.2, period = 4.68 days 00821 arg = 2*l+2*f+2*d+o; 00822 dpsi -= 0.000019 * ::sin(arg); 00823 deps += 0.000010 * ::cos(arg); 00824 // line 155 of Table 5.2, period = 7.34 days 00825 arg = -2*l+2*f+4*d+o; 00826 dpsi -= 0.000020 * ::sin(arg); 00827 deps += 0.000010 * ::cos(arg); 00828 // line 156 of Table 5.2, period = 14.22 days 00829 arg = lp+2*d+o; 00830 dpsi += 0.000018 * ::sin(arg); 00831 deps -= 0.000009 * ::cos(arg); 00832 // line 157 of Table 5.2, period = 14.25 days 00833 arg = l+d; 00834 dpsi -= 0.000033 * ::sin(arg); 00835 // line 158 of Table 5.2, period = 10.10 days 00836 arg = -l+4*d+o; 00837 dpsi -= 0.000018 * ::sin(arg); 00838 deps += 0.000009 * ::cos(arg); 00839 // line 159 of Table 5.2, period = 9.05 days 00840 arg = -l+4*f+o; 00841 dpsi += 0.000019 * ::sin(arg); 00842 deps -= 0.000009 * ::cos(arg); 00843 // line 160 of Table 5.2, period = -35.23 days 00844 arg = 2*f-3*d+2*o; 00845 dpsi -= 0.000020 * ::sin(arg); 00846 deps += 0.000008 * ::cos(arg); 00847 // line 161 of Table 5.2, period = 6.82 days 00848 arg = 4*f+2*o; 00849 dpsi += 0.000019 * ::sin(arg); 00850 deps -= 0.000008 * ::cos(arg); 00851 // line 162 of Table 5.2, period = 13.28 days 00852 arg = 2*l+lp; 00853 dpsi -= 0.000028 * ::sin(arg); 00854 // line 163 of Table 5.2, period = -16.10 days 00855 arg = 2*f-4*d+o; 00856 dpsi -= 0.000016 * ::sin(arg); 00857 deps += 0.000009 * ::cos(arg); 00858 // line 164 of Table 5.2, period = 5.90 days 00859 arg = -l-lp+2*f+4*d+2*o; 00860 dpsi -= 0.000017 * ::sin(arg); 00861 deps += 0.000007 * ::cos(arg); 00862 // line 165 of Table 5.2, period = 38.52 days 00863 arg = -l-2*lp+2*d; 00864 dpsi += 0.000027 * ::sin(arg); 00865 // line 166 of Table 5.2, period = 7.39 days 00866 arg = 4*d+o; 00867 dpsi -= 0.000016 * ::sin(arg); 00868 deps += 0.000007 * ::cos(arg); 00869 // line 167 of Table 5.2, period = 15.42 days 00870 arg = -lp+2*d+o; 00871 dpsi -= 0.000014 * ::sin(arg); 00872 deps += 0.000007 * ::cos(arg); 00873 // line 168 of Table 5.2, period = 4.08 days 00874 arg = l+2*f+4*d+2*o; 00875 dpsi -= 0.000016 * ::sin(arg); 00876 deps += 0.000007 * ::cos(arg); 00877 // line 169 of Table 5.2, period = -194.13 days 00878 arg = -2*l+2*d+2*o; 00879 dpsi += 0.000018 * ::sin(arg); 00880 deps -= 0.000008 * ::cos(arg); 00881 // line 170 of Table 5.2, period = 1616.44 days 00882 arg = -2*l+2*lp+2*d; 00883 dpsi -= 0.000022 * ::sin(arg); 00884 // line 171 of Table 5.2, period = -507.16 days 00885 arg = -2*l-lp+2*f+o; 00886 dpsi += 0.000009 * ::sin(arg); 00887 deps -= 0.000005 * ::cos(arg); 00888 // line 172 of Table 5.2, period = -9.17 days 00889 arg = -3*l+o; 00890 dpsi -= 0.000014 * ::sin(arg); 00891 deps += 0.000007 * ::cos(arg); 00892 // line 173 of Table 5.2, period = 13.69 days 00893 arg = 2*f+3*o; 00894 dpsi += 0.000020 * ::sin(arg); 00895 // line 174 of Table 5.2, period = 4.79 days 00896 arg = 2*f+4*d+o; 00897 dpsi -= 0.000012 * ::sin(arg); 00898 deps += 0.000006 * ::cos(arg); 00899 // line 175 of Table 5.2, period = 12.64 days 00900 arg = 4*f-2*d+o; 00901 dpsi += 0.000012 * ::sin(arg); 00902 deps -= 0.000007 * ::cos(arg); 00903 // line 176 of Table 5.2, period = 16.06 days 00904 arg = -2*lp+2*d; 00905 dpsi += 0.000021 * ::sin(arg); 00906 // line 177 of Table 5.2, period = 438.33 days 00907 arg = l-d+o; 00908 dpsi += 0.000017 * ::sin(arg) - 0.000003 * ::cos(arg); 00909 deps -= 0.000005 * ::cos(arg) + 0.000001 * ::sin(arg); 00910 // line 178 of Table 5.2, period = 5.56 days 00911 arg = l+lp+2*f+2*d+2*o; 00912 dpsi += 0.000015 * ::sin(arg); 00913 deps -= 0.000006 * ::cos(arg); 00914 // line 179 of Table 5.2, period = 8.73 days 00915 arg = 3*l+2*f-2*d+o; 00916 dpsi += 0.000012 * ::sin(arg); 00917 deps -= 0.000007 * ::cos(arg); 00918 // line 180 of Table 5.2, period = 29.26 days 00919 arg = -l-lp+2*f+2*o; 00920 dpsi -= 0.000016 * ::sin(arg); 00921 deps += 0.000006 * ::cos(arg); 00922 // line 181 of Table 5.2, period = -129.17 days 00923 arg = -2*l-lp+2*d+o; 00924 dpsi -= 0.000013 * ::sin(arg); 00925 deps += 0.000007 * ::cos(arg); 00926 // line 182 of Table 5.2, period = -14.70 days 00927 arg = -2*d+2*o; 00928 dpsi += 0.000013 * ::sin(arg); 00929 deps -= 0.000005 * ::cos(arg); 00930 // line 183 of Table 5.2, period = 7.38 days 00931 arg = -2*lp+2*f+2*d+2*o; 00932 dpsi -= 0.000013 * ::sin(arg); 00933 deps += 0.000005 * ::cos(arg); 00934 // line 184 of Table 5.2, period = -10.07 days 00935 arg = l-4*d+o; 00936 dpsi -= 0.000012 * ::sin(arg); 00937 deps += 0.000006 * ::cos(arg); 00938 // line 185 of Table 5.2, period = 29.39 days 00939 arg = -l+lp+2*d+o; 00940 dpsi -= 0.000010 * ::sin(arg); 00941 deps += 0.000006 * ::cos(arg); 00942 // line 186 of Table 5.2, period = 15.94 days 00943 arg = -2*l+4*d+o; 00944 dpsi += 0.000011 * ::sin(arg); 00945 deps -= 0.000006 * ::cos(arg); 00946 // line 187 of Table 5.2, period = 25.33 days 00947 arg = 2*f-d+o; 00948 dpsi -= 0.000010 * ::sin(arg); 00949 deps += 0.000005 * ::cos(arg); 00950 // line 188 of Table 5.2, period = 187.67 days 00951 arg = 2*lp+o; 00952 dpsi -= 0.000009 * ::sin(arg); 00953 deps += 0.000005 * ::cos(arg); 00954 // line 189 of Table 5.2, period = 90.10 days 00955 arg = 2*lp+2*f-2*d+o; 00956 dpsi += 0.000008 * ::sin(arg); 00957 deps -= 0.000005 * ::cos(arg); 00958 // line 190 of Table 5.2, period = 7.13 days 00959 arg = 2*l+2*d+o; 00960 dpsi -= 0.000009 * ::sin(arg); 00961 deps += 0.000005 * ::cos(arg); 00962 // line 191 of Table 5.2, period = -15.87 days 00963 arg = 2*l-4*d+o; 00964 dpsi -= 0.000011 * ::sin(arg); 00965 deps += 0.000005 * ::cos(arg); 00966 // line 192 of Table 5.2, period = 95.42 days 00967 arg = 2*l+2*f-4*d+o; 00968 dpsi += 0.000010 * ::sin(arg); 00969 deps -= 0.000005 * ::cos(arg); 00970 // line 193 of Table 5.2, period = -9.10 days 00971 arg = -l-2*f+o; 00972 dpsi -= 0.000010 * ::sin(arg); 00973 deps += 0.000005 * ::cos(arg); 00974 // line 194 of Table 5.2, period = 25.13 days 00975 arg = -l+lp+2*f+o; 00976 dpsi += 0.000009 * ::sin(arg); 00977 deps -= 0.000005 * ::cos(arg); 00978 // line 195 of Table 5.2, period = -35.80 days 00979 arg = -l+lp+2*f-2*d+o; 00980 dpsi -= 0.000011 * ::sin(arg); 00981 deps += 0.000005 * ::cos(arg); 00982 // line 196 of Table 5.2, period = 10.37 days 00983 arg = -l-lp+4*d; 00984 dpsi += 0.000015 * ::sin(arg); 00985 // line 197 of Table 5.2, period = 37.63 days 00986 arg = -3*l+4*d; 00987 dpsi += 0.000016 * ::sin(arg); 00988 // line 198 of Table 5.2, period = 4.00 days 00989 arg = 3*l+2*f+2*d+2*o; 00990 dpsi -= 0.000014 * ::sin(arg); 00991 // line 199 of Table 5.2, period = 471.89 days 00992 arg = 2*l-lp-2*d; 00993 dpsi -= 0.000009 * ::sin(arg) - 0.000001 * ::cos(arg); 00994 deps += 0.000001 * ::cos(arg); 00995 // line 200 of Table 5.2, period = -3396.16 days 00996 arg = 2*lp-2*f+2*d; 00997 dpsi -= 0.000009 * ::sin(arg); 00998 // line 201 of Table 5.2, period = 4.86 days 00999 arg = -lp+2*f+4*d+2*o; 01000 dpsi -= 0.000009 * ::sin(arg); 01001 // line 202 of Table 5.2, period = 27.32 days 01002 arg = -lp+2*f-d+2*o; 01003 dpsi += 0.000009 * ::sin(arg); 01004 // line 203 of Table 5.2, period = 9.37 days 01005 arg = l+lp+2*d; 01006 dpsi -= 0.000010 * ::sin(arg); 01007 // line 204 of Table 5.2, period = 219.17 days 01008 arg = 2*l-2*d+2*o; 01009 dpsi -= 0.000011 * ::sin(arg); 01010 // line 205 of Table 5.2, period = 4.74 days 01011 arg = 2*l-lp+2*f+2*d+2*o; 01012 dpsi -= 0.000009 * ::sin(arg); 01013 // line 206 of Table 5.2, period = 6.89 days 01014 arg = 4*l; 01015 dpsi += 0.000009 * ::sin(arg); 01016 // line 207 of Table 5.2, period = 6.64 days 01017 arg = 4*l+2*f-2*d+2*o; 01018 dpsi += 0.000012 * ::sin(arg); 01019 // line 208 of Table 5.2, period = 15.31 days 01020 arg = -l+3*d; 01021 dpsi -= 0.000010 * ::sin(arg); 01022 // line 209 of Table 5.2, period = 23.43 days 01023 arg = -l+4*f-2*d+2*o; 01024 dpsi -= 0.000009 * ::sin(arg); 01025 // line 210 of Table 5.2, period = 10.08 days 01026 arg = -l-2*lp+2*f+2*d+2*o; 01027 dpsi -= 0.000009 * ::sin(arg); 01028 // line 211 of Table 5.2, period = 16.63 days 01029 arg = -2*l-lp+4*d; 01030 dpsi += 0.000012 * ::sin(arg); 01031 // line 212 of Table 5.2, period = 7.50 days 01032 arg = -2*l-lp+2*f+4*d+2*o; 01033 dpsi -= 0.000012 * ::sin(arg); 01034 // line 213 of Table 5.2, period = 6.95 days 01035 arg = lp+2*f+2*d+o; 01036 dpsi += 0.000007 * ::sin(arg); 01037 // line 214 of Table 5.2, period = 12.71 days 01038 arg = 2*lp+2*f+2*o; 01039 dpsi += 0.000007 * ::sin(arg); 01040 // line 215 of Table 5.2, period = 14.77 days 01041 arg = -2*lp+2*f+2*o; 01042 dpsi -= 0.000008 * ::sin(arg); 01043 // line 216 of Table 5.2, period = 5.82 days 01044 arg = l+4*d; 01045 dpsi += 0.000008 * ::sin(arg); 01046 // line 217 of Table 5.2, period = 5.63 days 01047 arg = l+2*f+2*d; 01048 dpsi += 0.000008 * ::sin(arg); 01049 // line 218 of Table 5.2, period = -38.52 days 01050 arg = l+2*f-4*d+2*o; 01051 dpsi += 0.000007 * ::sin(arg); 01052 // line 219 of Table 5.2, period = 5.73 days 01053 arg = l-lp+2*f+2*d+o; 01054 dpsi -= 0.000008 * ::sin(arg); 01055 // line 220 of Table 5.2, period = 25.62 days 01056 arg = l-lp+2*f-2*d+2*o; 01057 dpsi -= 0.000007 * ::sin(arg); 01058 // line 221 of Table 5.2, period = 32.45 days 01059 arg = l-2*lp; 01060 dpsi += 0.000008 * ::sin(arg); 01061 // line 222 of Table 5.2, period = 13.83 days 01062 arg = 2*l+2*o; 01063 dpsi -= 0.000008 * ::sin(arg); 01064 // line 223 of Table 5.2, period = 134.27 days 01065 arg = 2*l+lp-2*d+o; 01066 dpsi += 0.000008 * ::sin(arg); 01067 // line 224 of Table 5.2, period = 9.20 days 01068 arg = 3*l+o; 01069 dpsi += 0.000007 * ::sin(arg); 01070 // line 225 of Table 5.2, period = 14.13 days 01071 arg = -l+2*f+d+2*o; 01072 dpsi += 0.000008 * ::sin(arg); 01073 // line 226 of Table 5.2, period = 7.22 days 01074 arg = -l+2*f+3*d+2*o; 01075 dpsi += 0.000008 * ::sin(arg); 01076 // line 227 of Table 5.2, period = 38.96 days 01077 arg = -l-2*f+4*d; 01078 dpsi -= 0.000007 * ::sin(arg); 01079 // line 228 of Table 5.2, period = 9.30 days 01080 arg = -l+lp+2*f+2*d+o; 01081 dpsi += 0.000007 * ::sin(arg); 01082 // line 229 of Table 5.2, period = 27.09 days 01083 arg = -l+2*lp+2*d; 01084 dpsi -= 0.000008 * ::sin(arg); 01085 // line 230 of Table 5.2, period = 2189.73 days 01086 arg = -l-lp+2*f-d+o; 01087 dpsi += 0.000007 * ::sin(arg); 01088 // line 231 of Table 5.2, period = -14.93 days 01089 arg = -2*l+2*f-2*d+o; 01090 dpsi -= 0.000008 * ::sin(arg); 01091 // line 232 of Table 5.2, period = 13.49 days 01092 arg = -2*l+4*f+2*o; 01093 dpsi -= 0.000007 * ::sin(arg); 01094 // line 233 of Table 5.2, period = -12.76 days 01095 arg = -2*l-2*f+2*d; 01096 dpsi += 0.000008 * ::sin(arg); 01097 // line 234 of Table 5.2, period = 285.41 days 01098 arg = -2*l+lp+2*f+o; 01099 dpsi += 0.000009 * ::sin(arg); 01100 // line 235 of Table 5.2, period = -28.15 days 01101 arg = -3*l+2*f+o; 01102 dpsi -= 0.000008 * ::sin(arg); 01103 // line 236 of Table 5.2, period = 27.43 days 01104 arg = lp+d+o; 01105 dpsi += 0.000005 * ::sin(arg); 01106 // line 237 of Table 5.2, period = 7.53 days 01107 arg = -lp+4*d; 01108 dpsi += 0.000006 * ::sin(arg); 01109 // line 238 of Table 5.2, period = -14.16 days 01110 arg = -lp-2*d+o; 01111 dpsi += 0.000005 * ::sin(arg); 01112 // line 239 of Table 5.2, period = -177.85 days 01113 arg = -2*lp+o; 01114 dpsi -= 0.000006 * ::sin(arg); 01115 // line 240 of Table 5.2, period = 6.97 days 01116 arg = l+2*f+d+o; 01117 dpsi += 0.000005 * ::sin(arg); 01118 // line 241 of Table 5.2, period = 126.51 days 01119 arg = l+2*f-3*d+2*o; 01120 dpsi -= 0.000006 * ::sin(arg); 01121 // line 242 of Table 5.2, period = -299.26 days 01122 arg = l-2*f+d; 01123 dpsi -= 0.000007 * ::sin(arg); 01124 // line 243 of Table 5.2, period = 13.72 days 01125 arg = l+lp+d; 01126 dpsi += 0.000005 * ::sin(arg); 01127 // line 244 of Table 5.2, period = -29.14 days 01128 arg = l-lp-2*d+o; 01129 dpsi += 0.000006 * ::sin(arg); 01130 // line 245 of Table 5.2, period = 8.93 days 01131 arg = 2*l+2*f-d+2*o; 01132 dpsi -= 0.000006 * ::sin(arg); 01133 // line 246 of Table 5.2, period = 6.73 days 01134 arg = 2*l+lp+2*f+o; 01135 dpsi += 0.000005 * ::sin(arg); 01136 // line 247 of Table 5.2, period = 6.98 days 01137 arg = 2*l-lp+2*f+o; 01138 dpsi -= 0.000006 * ::sin(arg); 01139 // line 248 of Table 5.2, period = 13.28 days 01140 arg = 2*l-lp+2*f-2*d+2*o; 01141 dpsi += 0.000005 * ::sin(arg); 01142 // line 249 of Table 5.2, period = 5.66 days 01143 arg = 3*l+2*d; 01144 dpsi += 0.000005 * ::sin(arg); 01145 // line 250 of Table 5.2, period = 5.58 days 01146 arg = 3*l-lp+2*f+2*o; 01147 dpsi -= 0.000005 * ::sin(arg); 01148 // line 251 of Table 5.2, period = 29.14 days 01149 arg = -l-lp+2*f+o; 01150 dpsi -= 0.000006 * ::sin(arg); 01151 // line 252 of Table 5.2, period = -13.72 days 01152 arg = -2*l+2*o; 01153 dpsi += 0.000006 * ::sin(arg); 01154 // line 253 of Table 5.2, period = 34.48 days 01155 arg = -2*l+3*d; 01156 dpsi -= 0.000005 * ::sin(arg); 01157 // line 254 of Table 5.2, period = -7.12 days 01158 arg = -2*l-2*d+o; 01159 dpsi -= 0.000005 * ::sin(arg); 01160 // line 255 of Table 5.2, period = 14.57 days 01161 arg = -2*l+2*f+2*d; 01162 dpsi -= 0.000006 * ::sin(arg); 01163 // line 256 of Table 5.2, period = -548.04 days 01164 arg = -2*l-lp+2*f; 01165 dpsi -= 0.000005 * ::sin(arg); 01166 // line 257 of Table 5.2, period = 15.24 days 01167 arg = -2*l-lp+2*f+2*d+2*o; 01168 dpsi += 0.000006 * ::sin(arg); 01169 // line 258 of Table 5.2, period = 27.21 days 01170 arg = f; 01171 dpsi += 0.000008 * ::cos(arg); 01172 // line 259 of Table 5.2, period = 27.32 days 01173 arg = f+o; 01174 dpsi -= 0.000016 * ::cos(arg); 01175 deps -= 0.000014 * ::sin(arg); 01176 // line 260 of Table 5.2, period = 2190.35 days 01177 arg = -l+f; 01178 dpsi += 0.000033 * ::cos(arg); 01179 // line 261 of Table 5.2, period = 3231.51 days 01180 arg = -l+f+o; 01181 dpsi -= 0.000105 * ::cos(arg); 01182 deps -= 0.000089 * ::sin(arg); 01183 // line 262 of Table 5.2, period = 6159.22 days 01184 arg = -l+f+2*o; 01185 dpsi += 0.000036 * ::cos(arg); 01186 deps += 0.000018 * ::sin(arg); 01187 // line 263 of Table 5.2, period = 65514.10 days 01188 arg = -l+f+3*o; 01189 dpsi -= 0.000006 * ::cos(arg); 01190 */ 01191 //End Code implementing Table 5.2 IERS Conventions 1996 nutation series. 01192 //----------------------------------------------------------------------- 01193 } 01194 01195 //------------------------------------------------------------------------------ 01196 // Zonal tide terms for corrections of UT1mUTC when that quantity does not 01197 // include tides (e.g. NGA EOP), ref. IERS 1996 Ch. 8, table 8.1 pg 74. 01198 // @param T, the coordinate transformation time at the time of interest 01199 // @param UT1mUT1R, the correction to UT1mUTC (seconds) 01200 // @param dlodR, the correction to the length of day (seconds) 01201 // @param domegaR, the correction to the Earth rotation rate, rad/second. 01202 void GeodeticFrames::UT1mUTCTidalCorrections(double T, 01203 double& UT1mUT1R, 01204 double& dlodR, 01205 double& domegaR) 01206 throw() 01207 { 01208 //----------------------------------------------------------------------- 01209 // Code to implement Table 8.1 of IERS Conventions 1996 series for 01210 // Zonal Tide terms for UT1R, length of the day and omega(Earth). 01211 // Units for UT1R, length of the day and omega(Earth) are: 01212 // seconds, seconds and radians/sec. 01213 // (Generated by perl script Table81.pl) 01214 double arg; 01215 UT1mUT1R = dlodR = domegaR = 0.0; 01216 01217 // Define (all doubles) and all in RADIANS: 01218 double o = Omega(T); // = mean longitude of lunar ascending node, in degrees, 01219 double f = F(T) ; // = mean longitude of the moon - Omega, in degrees, 01220 double d = D(T) ; // = mean elongation of the moon from the sun, in degrees, 01221 double l = L(T) ; // = mean anomaly of the moon, in degrees, and 01222 double lp = Lp(T) ; // = mean anomaly of the sun, in degrees. 01223 o *= DEG_TO_RAD; 01224 f *= DEG_TO_RAD; 01225 d *= DEG_TO_RAD; 01226 l *= DEG_TO_RAD; 01227 lp *= DEG_TO_RAD; 01228 01229 // line 1 of Table 8.1, period = 5.64 days 01230 arg = l+2*f+2*d+2*o; 01231 UT1mUT1R -= 0.02e-4 * ::sin(arg); 01232 dlodR += 0.3e-5 * ::cos(arg); 01233 domegaR -= 0.2e-14 * ::cos(arg); 01234 // line 2 of Table 8.1, period = 6.85 days 01235 arg = 2*l+2*f+o; 01236 UT1mUT1R -= 0.04e-4 * ::sin(arg); 01237 dlodR += 0.4e-5 * ::cos(arg); 01238 domegaR -= 0.3e-14 * ::cos(arg); 01239 // line 3 of Table 8.1, period = 6.86 days 01240 arg = 2*l+2*f+2*o; 01241 UT1mUT1R -= 0.10e-4 * ::sin(arg); 01242 dlodR += 0.9e-5 * ::cos(arg); 01243 domegaR -= 0.8e-14 * ::cos(arg); 01244 // line 4 of Table 8.1, period = 7.09 days 01245 arg = 2*f+2*d+o; 01246 UT1mUT1R -= 0.05e-4 * ::sin(arg); 01247 dlodR += 0.4e-5 * ::cos(arg); 01248 domegaR -= 0.4e-14 * ::cos(arg); 01249 // line 5 of Table 8.1, period = 7.10 days 01250 arg = 2*f+2*d+2*o; 01251 UT1mUT1R -= 0.12e-4 * ::sin(arg); 01252 dlodR += 1.1e-5 * ::cos(arg); 01253 domegaR -= 0.9e-14 * ::cos(arg); 01254 // line 6 of Table 8.1, period = 9.11 days 01255 arg = l+2*f; 01256 UT1mUT1R -= 0.04e-4 * ::sin(arg); 01257 dlodR += 0.3e-5 * ::cos(arg); 01258 domegaR -= 0.2e-14 * ::cos(arg); 01259 // line 7 of Table 8.1, period = 9.12 days 01260 arg = l+2*f+o; 01261 UT1mUT1R -= 0.41e-4 * ::sin(arg); 01262 dlodR += 2.8e-5 * ::cos(arg); 01263 domegaR -= 2.4e-14 * ::cos(arg); 01264 // line 8 of Table 8.1, period = 9.13 days 01265 arg = l+2*f+2*o; 01266 UT1mUT1R -= 0.99e-4 * ::sin(arg); 01267 dlodR += 6.8e-5 * ::cos(arg); 01268 domegaR -= 5.8e-14 * ::cos(arg); 01269 // line 9 of Table 8.1, period = 9.18 days 01270 arg = 3*l; 01271 UT1mUT1R -= 0.02e-4 * ::sin(arg); 01272 dlodR += 0.1e-5 * ::cos(arg); 01273 domegaR -= 0.1e-14 * ::cos(arg); 01274 // line 10 of Table 8.1, period = 9.54 days 01275 arg = -l+2*f+2*d+o; 01276 UT1mUT1R -= 0.08e-4 * ::sin(arg); 01277 dlodR += 0.5e-5 * ::cos(arg); 01278 domegaR -= 0.5e-14 * ::cos(arg); 01279 // line 11 of Table 8.1, period = 9.56 days 01280 arg = -l+2*f+2*d+2*o; 01281 UT1mUT1R -= 0.20e-4 * ::sin(arg); 01282 dlodR += 1.3e-5 * ::cos(arg); 01283 domegaR -= 1.1e-14 * ::cos(arg); 01284 // line 12 of Table 8.1, period = 9.61 days 01285 arg = l+2*d; 01286 UT1mUT1R -= 0.08e-4 * ::sin(arg); 01287 dlodR += 0.5e-5 * ::cos(arg); 01288 domegaR -= 0.4e-14 * ::cos(arg); 01289 // line 13 of Table 8.1, period = 12.81 days 01290 arg = 2*l+2*f-2*d+2*o; 01291 UT1mUT1R += 0.02e-4 * ::sin(arg); 01292 dlodR -= 0.1e-5 * ::cos(arg); 01293 domegaR += 0.1e-14 * ::cos(arg); 01294 // line 14 of Table 8.1, period = 13.17 days 01295 arg = lp+2*f+2*o; 01296 UT1mUT1R += 0.03e-4 * ::sin(arg); 01297 dlodR -= 0.1e-5 * ::cos(arg); 01298 domegaR += 0.1e-14 * ::cos(arg); 01299 // line 15 of Table 8.1, period = 13.61 days 01300 arg = 2*f; 01301 UT1mUT1R -= 0.30e-4 * ::sin(arg); 01302 dlodR += 1.4e-5 * ::cos(arg); 01303 domegaR -= 1.2e-14 * ::cos(arg); 01304 // line 16 of Table 8.1, period = 13.63 days 01305 arg = 2*f+o; 01306 UT1mUT1R -= 3.21e-4 * ::sin(arg); 01307 dlodR += 14.8e-5 * ::cos(arg); 01308 domegaR -= 12.5e-14 * ::cos(arg); 01309 // line 17 of Table 8.1, period = 13.66 days 01310 arg = 2*f+2*o; 01311 UT1mUT1R -= 7.76e-4 * ::sin(arg); 01312 dlodR += 35.7e-5 * ::cos(arg); 01313 domegaR -= 30.1e-14 * ::cos(arg); 01314 // line 18 of Table 8.1, period = 13.75 days 01315 arg = 2*l-o; 01316 UT1mUT1R += 0.02e-4 * ::sin(arg); 01317 dlodR -= 0.1e-5 * ::cos(arg); 01318 domegaR += 0.1e-14 * ::cos(arg); 01319 // line 19 of Table 8.1, period = 13.78 days 01320 arg = 2*l; 01321 UT1mUT1R -= 0.34e-4 * ::sin(arg); 01322 dlodR += 1.5e-5 * ::cos(arg); 01323 domegaR -= 1.3e-14 * ::cos(arg); 01324 // line 20 of Table 8.1, period = 13.81 days 01325 arg = 2*l+o; 01326 UT1mUT1R += 0.02e-4 * ::sin(arg); 01327 dlodR -= 0.1e-5 * ::cos(arg); 01328 domegaR += 0.1e-14 * ::cos(arg); 01329 // line 21 of Table 8.1, period = 14.19 days 01330 arg = -lp+2*f+2*o; 01331 UT1mUT1R -= 0.02e-4 * ::sin(arg); 01332 dlodR += 0.1e-5 * ::cos(arg); 01333 domegaR -= 0.1e-14 * ::cos(arg); 01334 // line 22 of Table 8.1, period = 14.73 days 01335 arg = 2*d-o; 01336 UT1mUT1R += 0.05e-4 * ::sin(arg); 01337 dlodR -= 0.2e-5 * ::cos(arg); 01338 domegaR += 0.2e-14 * ::cos(arg); 01339 // line 23 of Table 8.1, period = 14.77 days 01340 arg = 2*d; 01341 UT1mUT1R -= 0.73e-4 * ::sin(arg); 01342 dlodR += 3.1e-5 * ::cos(arg); 01343 domegaR -= 2.6e-14 * ::cos(arg); 01344 // line 24 of Table 8.1, period = 14.80 days 01345 arg = 2*d+o; 01346 UT1mUT1R -= 0.05e-4 * ::sin(arg); 01347 dlodR += 0.2e-5 * ::cos(arg); 01348 domegaR -= 0.2e-14 * ::cos(arg); 01349 // line 25 of Table 8.1, period = 15.39 days 01350 arg = -lp+2*d; 01351 UT1mUT1R -= 0.05e-4 * ::sin(arg); 01352 dlodR += 0.2e-5 * ::cos(arg); 01353 domegaR -= 0.2e-14 * ::cos(arg); 01354 // line 26 of Table 8.1, period = 23.86 days 01355 arg = l+2*f-2*d+o; 01356 UT1mUT1R += 0.05e-4 * ::sin(arg); 01357 dlodR -= 0.1e-5 * ::cos(arg); 01358 domegaR += 0.1e-14 * ::cos(arg); 01359 // line 27 of Table 8.1, period = 23.94 days 01360 arg = l+2*f-2*d+2*o; 01361 UT1mUT1R += 0.10e-4 * ::sin(arg); 01362 dlodR -= 0.3e-5 * ::cos(arg); 01363 domegaR += 0.2e-14 * ::cos(arg); 01364 // line 28 of Table 8.1, period = 25.62 days 01365 arg = l+lp; 01366 UT1mUT1R += 0.04e-4 * ::sin(arg); 01367 dlodR -= 0.1e-5 * ::cos(arg); 01368 domegaR += 0.1e-14 * ::cos(arg); 01369 // line 29 of Table 8.1, period = 26.88 days 01370 arg = -l+2*f; 01371 UT1mUT1R += 0.05e-4 * ::sin(arg); 01372 dlodR -= 0.1e-5 * ::cos(arg); 01373 domegaR += 0.1e-14 * ::cos(arg); 01374 // line 30 of Table 8.1, period = 26.98 days 01375 arg = -l+2*f+o; 01376 UT1mUT1R += 0.18e-4 * ::sin(arg); 01377 dlodR -= 0.4e-5 * ::cos(arg); 01378 domegaR += 0.3e-14 * ::cos(arg); 01379 // line 31 of Table 8.1, period = 27.09 days 01380 arg = -l+2*f+2*o; 01381 UT1mUT1R += 0.44e-4 * ::sin(arg); 01382 dlodR -= 1.0e-5 * ::cos(arg); 01383 domegaR += 0.9e-14 * ::cos(arg); 01384 // line 32 of Table 8.1, period = 27.44 days 01385 arg = l-o; 01386 UT1mUT1R += 0.53e-4 * ::sin(arg); 01387 dlodR -= 1.2e-5 * ::cos(arg); 01388 domegaR += 1.0e-14 * ::cos(arg); 01389 // line 33 of Table 8.1, period = 27.56 days 01390 arg = l; 01391 UT1mUT1R -= 8.26e-4 * ::sin(arg); 01392 dlodR += 18.8e-5 * ::cos(arg); 01393 domegaR -= 15.9e-14 * ::cos(arg); 01394 // line 34 of Table 8.1, period = 27.67 days 01395 arg = l+o; 01396 UT1mUT1R += 0.54e-4 * ::sin(arg); 01397 dlodR -= 1.2e-5 * ::cos(arg); 01398 domegaR += 1.0e-14 * ::cos(arg); 01399 // line 35 of Table 8.1, period = 29.53 days 01400 arg = d; 01401 UT1mUT1R += 0.05e-4 * ::sin(arg); 01402 dlodR -= 0.1e-5 * ::cos(arg); 01403 domegaR += 0.1e-14 * ::cos(arg); 01404 // line 36 of Table 8.1, period = 29.80 days 01405 arg = l-lp; 01406 UT1mUT1R -= 0.06e-4 * ::sin(arg); 01407 dlodR += 0.1e-5 * ::cos(arg); 01408 domegaR -= 0.1e-14 * ::cos(arg); 01409 // line 37 of Table 8.1, period = 31.66 days 01410 arg = -l+2*d-o; 01411 UT1mUT1R += 0.12e-4 * ::sin(arg); 01412 dlodR -= 0.2e-5 * ::cos(arg); 01413 domegaR += 0.2e-14 * ::cos(arg); 01414 // line 38 of Table 8.1, period = 31.81 days 01415 arg = -l+2*d; 01416 UT1mUT1R -= 1.82e-4 * ::sin(arg); 01417 dlodR += 3.6e-5 * ::cos(arg); 01418 domegaR -= 3.0e-14 * ::cos(arg); 01419 // line 39 of Table 8.1, period = 31.96 days 01420 arg = -l+2*d+o; 01421 UT1mUT1R += 0.13e-4 * ::sin(arg); 01422 dlodR -= 0.3e-5 * ::cos(arg); 01423 domegaR += 0.2e-14 * ::cos(arg); 01424 // line 40 of Table 8.1, period = 32.61 days 01425 arg = l-2*f+2*d-o; 01426 UT1mUT1R += 0.02e-4 * ::sin(arg); 01427 // line 41 of Table 8.1, period = 34.85 days 01428 arg = -l-lp+2*d; 01429 UT1mUT1R -= 0.09e-4 * ::sin(arg); 01430 dlodR += 0.2e-5 * ::cos(arg); 01431 domegaR -= 0.1e-14 * ::cos(arg); 01432 01433 //End Code implementing Table 8.1 IERS Conventions 1996 UT1R tide series. 01434 //----------------------------------------------------------------------- 01435 } 01436 01437 //--------------------------------------------------------------------------------- 01438 // Compute the Greenwich hour angle of the true vernal equinox, or 01439 // Greenwich Apparent Sidereal Time (GAST) in radians, 01440 // given the (UT) time of interest t, and, where T = CoordTransTime(t), 01441 // o = Omega(T) = mean longitude of lunar ascending node, in degrees, 01442 // eps = the obliquity of the ecliptic, in degrees, 01443 // dpsi = nutation in longitude (counted in the ecliptic), 01444 // in seconds of arc. 01445 // 01446 // GAST = Greenwich hour angle of the true vernal equinox 01447 // GAST = GMST + dpsi*cos(eps) + 0.00264" * sin(Omega) +0.000063" * sin(2*Omega) 01448 // (these terms account for the accumulated precession and nutation in 01449 // right ascension and minimize any discontinuity in UT1) 01450 // 01451 // GMST = Greenwich hour angle of the mean vernal equinox 01452 // = Greenwich Mean Sideral Time 01453 // = GMST0 + r*[UTC + (UT1-UTC)] 01454 // r = is the ratio of universal to sidereal time 01455 // = 1.002737909350795 + 5.9006E-11*T' - 5.9e-15*T'^2 01456 // T' = days'/36525 01457 // days'= number of days elapsed since the Julian Epoch t0 (J2000) 01458 // = +/-(integer+0.5) 01459 // and 01460 // (UT1-UTC) (seconds) is taken from the IERS bulletin 01461 // 01462 // GMST0 = GMST at 0h UT1 01463 // = 6h 41min (50.54841+8640184.812866*T'+0.093104*T'^2-6.2E-6*T'^3)sec 01464 // 01465 // see pg 21 of the Reference (IERS 1996). 01466 double GeodeticFrames::gast(DayTime t, 01467 double om, 01468 double eps, 01469 double dpsi, 01470 double UT1mUTC) 01471 throw() 01472 { 01473 double G = GMST(t,UT1mUTC); 01474 // add dpsi, eps and Omega terms 01475 om *= DEG_TO_RAD; 01476 eps *= DEG_TO_RAD; 01477 G += ( dpsi * ::cos(eps) 01478 + 0.00264 * ::sin(om) 01479 + 0.000063 * ::sin(2.0*om) ) * DEG_TO_RAD / 3600.0; 01480 01481 return G; 01482 } 01483 01484 //--------------------------------------------------------------------------------- 01485 // Compute the precession matrix, a 3x3 rotation matrix, given T, 01486 // the coordinate transformation time at the time of interest 01487 Matrix<double> GeodeticFrames::PrecessionMatrix(double T) 01488 throw(InvalidRequest) 01489 { 01490 try { 01491 // IAU76 - ref McCarthy - seconds of arc 01492 double zeta = T*(2306.2181 + T*(0.30188 + T*0.017998)); 01493 double theta = T*(2004.3109 - T*(0.42665 + T*0.041833)); 01494 double z = T*(2306.2181 + T*(1.09468 + T*0.018203)); 01495 01496 // convert to degrees 01497 zeta /= 3600.0; 01498 theta /= 3600.0; 01499 z /= 3600.0; 01500 01501 Matrix<double> R1 = rotation(-zeta*DEG_TO_RAD, 3); 01502 Matrix<double> R2 = rotation(theta*DEG_TO_RAD, 2); 01503 Matrix<double> R3 = rotation(-z*DEG_TO_RAD, 3); 01504 Matrix<double> P = R3*R2*R1; 01505 01506 return P; 01507 } 01508 catch(InvalidRequest& ire) { 01509 GPSTK_RETHROW(ire); 01510 } 01511 } 01512 01513 //--------------------------------------------------------------------------------- 01514 // Compute the nutation matrix, given 01515 // eps, the obliquity of the ecliptic, in degrees, 01516 // dpsi, the nutation in longitude (counted in the ecliptic), 01517 // in seconds of arc, and 01518 // deps, the nutation in obliquity, in seconds of arc. 01519 Matrix<double> GeodeticFrames::NutationMatrix(double eps, 01520 double dpsi, 01521 double deps) 01522 throw(InvalidRequest) 01523 { 01524 Matrix<double> N; 01525 try { 01526 Matrix<double> R1 = rotation(-eps*DEG_TO_RAD, 1); 01527 Matrix<double> R2 = rotation(dpsi*DEG_TO_RAD/3600.0, 3); 01528 Matrix<double> R3 = rotation((eps+deps/3600.0)*DEG_TO_RAD, 1); 01529 N = R1*R2*R3; 01530 return N; 01531 } 01532 catch(InvalidRequest& ire) { 01533 GPSTK_RETHROW(ire); 01534 } 01535 } 01536 01537 //--------------------------------------------------------------------------------- 01538 // public functions 01539 01540 //--------------------------------------------------------------------------------- 01541 // Compute Greenwich Mean Sidereal Time, or the Greenwich hour angle of 01542 // the mean vernal equinox (radians), given the coordinate time of interest, 01543 // and UT1-UTC (sec), which comes from the IERS bulletin. 01544 // @param t DayTime epoch of the rotation. 01545 // @param UT1mUTC, UT1-UTC in seconds, as found in the IERS bulletin. 01546 // @param reduced, bool true when UT1mUTC is 'reduced', meaning assumes 01547 // 'no tides', as is the case with the NGA EOPs (default=F). 01548 double GeodeticFrames::GMST(DayTime t, 01549 double UT1mUTC, 01550 bool reduced) 01551 throw() 01552 { 01553 // days since epoch 01554 double days = t.JD() - JulianEpoch; // days 01555 if(days <= 0.0) days -= 1.0; 01556 double Tp = days/36525.0; // dim-less 01557 01558 // Compute GMST in radians 01559 double G; 01560 // seconds (24060s = 6h 41min) 01561 //G = 24110.54841 + (8640184.812866 + (0.093104 - 6.2e-6*Tp)*Tp)*Tp; // sec 01562 //G /= 86400.0; // instead, divide the numbers above manually 01563 G = 0.279057273264 + 100.0021390378009*Tp // seconds/86400 = days 01564 + (0.093104 - 6.2e-6*Tp)*Tp*Tp/86400.0; 01565 01566 // if reduced, compute tidal terms 01567 if(reduced) { 01568 double dlodR,domegaR,UT1mUT1R; 01569 UT1mUTCTidalCorrections(CoordTransTime(t), UT1mUT1R, dlodR, domegaR); 01570 UT1mUTC = UT1mUT1R-UT1mUTC; 01571 } 01572 01573 G += (UT1mUTC + t.secOfDay())/86400.0; // days 01574 01575 // put answer between 0 and 2pi 01576 G = fmod(G,1.0); 01577 while(G < 0.0) G += 1.0; 01578 G *= TWO_PI; // radians 01579 01580 return G; 01581 } 01582 01583 //--------------------------------------------------------------------------------- 01584 // Compute Greenwich Apparent Sidereal Time, or the Greenwich hour angle of 01585 // the true vernal equinox (radians), given the coordinate time of interest, 01586 // and UT1-UTC (sec), which comes from the IERS bulletin. 01587 // @param t DayTime epoch of the rotation. 01588 // @param UT1mUTC, UT1-UTC in seconds, as found in the IERS bulletin. 01589 double GeodeticFrames::GAST(DayTime t, 01590 double UT1mUTC, 01591 bool reduced) 01592 throw() 01593 { 01594 double T = CoordTransTime(t); 01595 double o = Omega(T); 01596 double eps = Obliquity(T); 01597 double deps,dpsi; 01598 01599 NutationAngles(T,deps,dpsi); // deps is not used... 01600 01601 // if reduced (NGA), correct for tides 01602 double UT1mUT1R,dlodR,domegaR; 01603 if(reduced) 01604 UT1mUTCTidalCorrections(T, UT1mUT1R, dlodR, domegaR); 01605 01606 return gast(t, o, eps, dpsi, reduced ? UT1mUT1R-UT1mUTC : UT1mUTC); 01607 } 01608 01609 //--------------------------------------------------------------------------------- 01610 // Generate transformation matrix (3X3 rotation) due to polar motion (xp,yp) 01611 // xp and yp are in arc seconds, as found in the IERS bulletin 01612 Matrix<double> GeodeticFrames::PolarMotion(double xp, 01613 double yp) 01614 throw(InvalidRequest) 01615 { 01616 try { 01617 xp *= DEG_TO_RAD/3600.0; 01618 yp *= DEG_TO_RAD/3600.0; 01619 Matrix<double> R1,R2; 01620 R1 = rotation(yp,1); 01621 R2 = rotation(xp,2); 01622 return (R1*R2); 01623 } 01624 catch(InvalidRequest& ire) { 01625 GPSTK_RETHROW(ire); 01626 } 01627 } 01628 01629 //--------------------------------------------------------------------------------- 01630 // Generate precise transformation matrix (3X3 rotation) due to Earth rotation 01631 // at Greenwich hour angle of the true vernal equinox and which accounts for 01632 // precession and nutation in right ascension, given the UT time of interest 01633 // and the UT1-UTC correction (in sec), obtained from the IERS bulletin. 01634 Matrix<double> GeodeticFrames::PreciseEarthRotation(DayTime t, 01635 double UT1mUTC, 01636 bool reduced) 01637 throw(InvalidRequest) 01638 { 01639 try { 01640 return (rotation(-GAST(t,UT1mUTC,reduced),3)); 01641 } 01642 catch(InvalidRequest& ire) { 01643 GPSTK_RETHROW(ire); 01644 } 01645 } 01646 01647 //--------------------------------------------------------------------------------- 01648 // Generate an Earth Nutation Matrix (3X3 rotation) at the given DayTime. 01649 // @param t DayTime epoch of the rotation. 01650 Matrix<double> GeodeticFrames::Nutation(DayTime t) 01651 throw(InvalidRequest) 01652 { 01653 try { 01654 double T=CoordTransTime(t); 01655 double eps=Obliquity(T); // degrees 01656 double deps,dpsi; 01657 01658 NutationAngles(T,deps,dpsi); 01659 01660 return NutationMatrix(eps,dpsi,deps); 01661 } 01662 catch(InvalidRequest& ire) { 01663 GPSTK_RETHROW(ire); 01664 } 01665 } 01666 01667 //--------------------------------------------------------------------------------- 01668 // Generate the full transformation matrix (3x3 rotation) relating the ECEF 01669 // frame to the conventional inertial frame. 01670 // throw(); Input is the time of interest, 01671 // the polar motion angles xp and yp (arcseconds), and UT1-UTC (seconds) 01672 // (xp,yp and UT1-UTC are just as found in the IERS bulletin). 01673 Matrix<double> GeodeticFrames::ECEFtoInertial(DayTime t, 01674 double xp, 01675 double yp, 01676 double UT1mUTC, 01677 bool reduced) 01678 throw(InvalidRequest) 01679 { 01680 try { 01681 Matrix<double> P,N,W,S; 01682 01683 double T=CoordTransTime(t); 01684 P = PrecessionMatrix(T); 01685 01686 double eps=Obliquity(T); // degrees 01687 double deps,dpsi; 01688 NutationAngles(T,deps,dpsi); 01689 N = NutationMatrix(eps,dpsi,deps); 01690 01691 // PolarMotion converts xp, yp to radians 01692 W = PolarMotion(xp, yp); 01693 01694 // if reduced (NGA), correct UT1mUTC for tides 01695 double UT1mUT1R,dlodR,domegaR; 01696 if(reduced) 01697 UT1mUTCTidalCorrections(T, UT1mUT1R, dlodR, domegaR); 01698 01699 double omega = Omega(T); 01700 double g = gast(t, omega, eps, dpsi, reduced ? UT1mUT1R-UT1mUTC : UT1mUTC); 01701 S = rotation(-g,3); 01702 01703 return (P*N*W*S); 01704 } 01705 catch(InvalidRequest& ire) { 01706 GPSTK_RETHROW(ire); 01707 } 01708 } 01709 01710 //--------------------------------------------------------------------------------- 01711 // Given a rotation matrix R (3x3), inverse(R)=transpose(R), 01712 // find the Euler angles (theta,phi,psi) which produce this rotation, 01713 // and also determine the magnitude (alpha) and direction (nhat = unit 3-vector) 01714 // of the net rotation. 01715 // Throw InvalidRequest if the matrix is not a rotation matrix. 01716 // 01717 // Euler angles (this is one convention - there are others): 01718 // Let R1 = rotation about z through angle phi 01719 // R2 = rotation about x through angle theta ( 0 <= theta <= pi) 01720 // R3 = rotation about z through angle psi 01721 // Any rotation matrix can be expressed as the product of these rotations: 01722 // R = R3*R2*R1. In particular, by definition 01723 // 01724 // [ cos(phi) sin(phi) 0 ] 01725 // R1 = [ -sin(phi) cos(phi) 0 ] 01726 // [ 0 0 1 ] 01727 // 01728 // [ cos(theta) 0 -sin(theta) ] 01729 // R2 = [ 0 1 0 ] 01730 // [ sin(theta) 0 cos(theta) ] 01731 // 01732 // [ cos(psi) sin(psi) 0 ] 01733 // R3 = [ -sin(psi) cos(psi) 0 ] 01734 // [ 0 0 1 ] 01735 // 01736 // and if we define 01737 // [ r11 r12 r13 ] 01738 // R = [ r21 r22 r23 ] 01739 // [ r31 r32 r33 ] 01740 // 01741 // then 01742 // r11 = cos(phi)cos(psi) - cos(theta)sin(phi)sin(psi) 01743 // r12 = sin(phi)cos(psi) + cos(theta)cos(phi)sin(psi) 01744 // r13 = sin(psi)sin(theta) 01745 // r21 = -cos(phi)sin(psi) - cos(theta)sin(phi)cos(psi) 01746 // r22 = -sin(phi)sin(psi) + cos(theta)cos(phi)cos(psi) 01747 // r23 = cos(psi)sin(theta) 01748 // r31 = sin(phi)sin(theta) 01749 // r32 = -cos(phi)sin(theta) 01750 // r33 = cos(theta) 01751 // 01752 // This can be inverted to get theta,phi,psi by 01753 // cos(theta) = r33 01754 // sin(theta) = +sqrt(1-cos^2(theta)) [recall 0 <= theta <= pi] 01755 // cos(phi) = -r32/sin(theta) [if sin(theta) != 0] 01756 // sin(phi) = r31/sin(theta) 01757 // cos(psi) = r23/sin(theta) 01758 // sin(psi) = r13/sin(theta) 01759 // or better 01760 // theta = acos(r33) 01761 // phi = atan2(r31,-r32) 01762 // psi = atan2(r13,r23) 01763 // 01764 // If sin(theta) == 0 then cos(theta) = {+1 OR -1}, theta = {0 OR pi}, 01765 // and (refer to R3*R2*R1 above), R becomes 01766 // [ cos(phi+psi) sin(phi+psi) 0 ] OR [ cos(phi-psi) sin(phi-psi) 0 ] 01767 // R = [ -sin(phi+psi) cos(phi+psi) 0 ] OR [ sin(phi-psi) -cos(phi-psi) 0 ] 01768 // [ 0 0 1 ] OR [ 0 0 -1 ] 01769 // and thus 01770 // cos(phi+psi) = r11 = r22 OR cos(phi-psi) = r11 = -r22 01771 // sin(phi+psi) = r12 = -r21 OR sin(phi-psi) = r12 = r21 01772 // 01773 // Now let E = e0,e1,e2,e3 = (e0,e) = quaternion [scalar (e0) and 3-vector (e)]. 01774 // E describes a rotation (through angle alpha) about axis (unit vector) nhat as 01775 // cos(alpha) = 2e0^2-1 = e0^2-dot(e,e) = e0^2 - sum_squares(ei,i=1,2,3). 01776 // nhat * sin(alpha) = 2e*e0 or 01777 // nhat = (e1,e2,e3)/dot(e,e) 01778 // [e0^2 + dot(e,e) = sum_squares(ei,i=1,2,3,4) = 1 by identity]. 01779 // It can be shown that the Euler rotation is equal to the quaternion 01780 // e0 = cos[(phi+psi)/2]*cos(theta/2) 01781 // e1 = sin[(phi-psi)/2]*sin(theta/2) 01782 // e2 = cos[(phi-psi)/2]*sin(theta/2) 01783 // e3 = sin[(phi+psi)/2]*cos(theta/2) 01784 // and 01785 // rij = dij(e0^2-ekek)+2eiej+2(epsijk)e0ek 01786 // where dij is the Kroncker delta and epsijk is the permutation symbol; 01787 // r11 = e0^2 + e1^2 - e2^2 - e3^2 01788 // r12 = 2(e1e2 + e0e3) 01789 // r13 = 2(e1e3 - e0e2) 01790 // r21 = 2(e1e2 - e0e3) 01791 // r22 = e0^2 - e1^2 + e2^2 - e3^2 01792 // r23 = 2(e2e3 + e0e1) 01793 // r31 = 2(e1e3 + e0e2) 01794 // r32 = 2(e2e3 - e0e1) 01795 // r33 = e0^2 - e1^2 - e2^2 + e3^2 01796 // 01797 // If theta=0 this reduces trivially to a simple rotation about z thru phi+psi 01798 // e0 = cos[(phi+psi)/2] 01799 // e1 = e2 = 0 01800 // e3 = sin[(phi+psi)/2] 01801 // alpha = phi+psi (phi and psi cannot be separated) 01802 // nhat = (0,0,1) 01803 // while if theta=pi it reduces to 01804 // e0 = e3 = 0 01805 // e1 = sin[(phi-psi)/2] 01806 // e2 = cos[(phi-psi)/2] 01807 // alpha = pi 01808 // nhat = (e1,e2,0) (dot(e,e)=1) 01809 void GeodeticFrames::ResolveRotation(const Matrix<double>& R, 01810 double& theta, 01811 double& phi, 01812 double& psi, 01813 double& alpha, 01814 Vector<double>& nhat) 01815 throw(InvalidRequest) 01816 { 01817 if(R.rows() != 3 || R.cols() != 3) { 01818 using namespace StringUtils; 01819 InvalidRequest ir(string("Input matrix has dimension ") 01820 + asString<int>(R.rows()) + string(",") 01821 + asString<int>(R.cols())); 01822 GPSTK_THROW(ir); 01823 } 01824 01825 const double tol=1.e-12; // tolerance TD use limits 01826 Matrix<double> T=transpose(R)*R-ident<double>(3); 01827 if(normF(T) > tol) { // RSS of elements 01828 InvalidRequest ir(string("Input matrix is not a rotation")); 01829 GPSTK_THROW(ir); 01830 } 01831 01832 // first find the Euler angles 01833 double st = SQRT(1.0-R(2,2)*R(2,2)); 01834 if(st < tol) { // theta is 0 or pi 01835 if(R(2,2) > 0.0) theta = 0.0; 01836 else theta = PI; 01837 psi = 0.0; // arbitrary, since only phi +/- psi can be known 01838 // tan(phi) = r12/r11 = sin(phi+/-psi) / cos(phi+/-psi) 01839 phi = atan2(R(0,1),R(0,0)); 01840 } 01841 else { 01842 theta = ::acos(R(2,2)); 01843 // tan(psi) = r13/r23 = sin(psi)sin(theta)/cos(psi)sin(theta) 01844 psi = ::atan2(R(0,2),R(1,2)); 01845 // tan(phi) = r31/-r32 = sin(phi)sin(theta)/cos(phi)sin(theta) 01846 phi = ::atan2(R(2,0),-R(2,1)); 01847 } 01848 01849 // now find the rotation angle alpha and the axis of rotation 01850 nhat.resize(3); 01851 if(theta == 0.0) { 01852 alpha = phi; 01853 nhat(0) = nhat(1) = 0.0; 01854 nhat(2) = 1.0; 01855 } 01856 else if(theta == PI) { 01857 alpha = PI; 01858 nhat(0) = ::sin(phi/2.0); 01859 nhat(1) = ::cos(phi/2.0); 01860 nhat(2) = 0.0; 01861 } 01862 else { 01863 double e0 = ::cos((phi+psi)/2.0) * ::cos(theta/2.0); 01864 alpha = ::acos(2.0*e0*e0-1.0); 01865 // construct e, then normalize 01866 nhat(0) = ::sin((phi-psi)/2.0) * ::sin(theta/2.0); 01867 nhat(1) = ::cos((phi-psi)/2.0) * ::sin(theta/2.0); 01868 nhat(2) = ::sin((phi+psi)/2.0) * ::cos(theta/2.0); 01869 e0 = norm(nhat); 01870 nhat /= e0; 01871 } 01872 } 01873 01874 } // end namespace gpstk 01875 /* 01876 # Table 5.2. IERS 1996 series for nutation in longitude Dpsi and obliquity Deps, 01877 # referred to the mean equator and equinox of date, with t measured in Julian 01878 # centuries from epoch J2000.0. The signs of the fundamental arguments, periods, 01879 # and coefficients may differ from the original publication. These 01880 # have been changed to be consistent with other portions of this chapter. 01881 # 263 01882 # Dpsi = sum { (A_i+A'_i*t)*sin(ARGUMENT) + A''_i*cos(ARGUMENT) } 01883 # i=1 01884 # 01885 # 263 01886 # Deps = sum { (B_i+B'_i*t)*cos(ARGUMENT) + B''_i*sin(ARGUMENT) } 01887 # i=1 01888 # LONGITUDE OBLIQUITY 01889 #MULTIPLIERS OF: PERIOD (0.001 mas) (0.001 mas) 01890 # l l' F D Om (days) A A' B B' A'' B'' 01891 0 0 0 0 1 -6798.38 -17206277 -17419 9205356 886 3645 1553 01892 0 0 2 -2 2 182.62 -1317014 -156 573058 -306 -1400 -464 01893 0 0 2 0 2 13.66 -227720 -23 97864 -48 269 136 01894 0 0 0 0 2 -3399.18 207429 21 -89747 47 -71 -29 01895 0 -1 0 0 0 -365.26 -147538 364 7388 -19 1121 198 01896 0 1 2 -2 2 121.75 -51687 123 22440 -68 -54 -18 01897 1 0 0 0 0 27.55 71118 7 -687 0 -94 39 01898 0 0 2 0 1 13.63 -38752 -37 20076 2 34 32 01899 1 0 2 0 2 9.13 -30137 -4 12896 -6 77 35 01900 0 -1 2 -2 2 365.22 21583 -49 -9591 30 6 12 01901 0 0 2 -2 1 177.84 12820 14 -6897 -1 18 4 01902 -1 0 2 0 2 27.09 12353 1 -5334 3 2 0 01903 -1 0 0 2 0 31.81 15699 1 -127 0 -18 9 01904 1 0 0 0 1 27.67 6314 6 -3323 0 3 -1 01905 -1 0 0 0 1 -27.44 -5797 -6 3141 0 -19 -8 01906 -1 0 2 2 2 9.56 -5965 -1 2554 -1 14 7 01907 1 0 2 0 1 9.12 -5163 -4 2635 0 12 8 01908 -2 0 2 0 1 1305.48 4590 5 -2424 -1 1 1 01909 0 0 0 2 0 14.77 6336 1 -125 0 -15 3 01910 0 0 2 2 2 7.10 -3854 0 1643 0 15 6 01911 -2 0 0 2 0 -205.89 -4774 0 48 0 -2 -3 01912 2 0 2 0 2 6.86 -3102 0 1323 -1 12 5 01913 1 0 2 -2 2 23.94 2863 0 -1235 1 0 0 01914 -1 0 2 0 1 26.98 2044 2 -1076 0 1 0 01915 2 0 0 0 0 13.78 2923 0 -62 0 -8 1 01916 0 0 2 0 0 13.61 2585 0 -56 0 -7 1 01917 0 1 0 0 1 386.00 -1406 -3 857 0 8 -4 01918 -1 0 0 2 1 31.96 1517 1 -801 0 1 0 01919 0 2 2 -2 2 91.31 -1578 7 685 -4 -2 -1 01920 0 0 -2 2 0 -173.31 2178 0 -15 0 1 1 01921 1 0 0 -2 1 -31.66 -1286 -1 694 0 -4 -2 01922 0 -1 0 0 1 -346.64 -1269 1 642 1 6 2 01923 -1 0 2 2 1 9.54 -1022 -1 522 0 2 1 01924 0 -2 0 0 0 -182.63 -1671 8 14 0 -1 -1 01925 1 0 2 2 2 5.64 -768 0 325 0 4 2 01926 -2 0 2 0 0 1095.18 -1102 0 10 0 -1 0 01927 0 1 2 0 2 13.17 757 -2 -326 -2 -1 0 01928 0 0 2 2 1 7.09 -664 -1 335 -1 2 1 01929 0 -1 2 0 2 14.19 -714 2 307 2 1 0 01930 0 0 0 2 1 14.80 -631 -1 327 0 0 0 01931 1 0 2 -2 1 23.86 580 1 -307 0 0 0 01932 2 0 2 -2 2 12.81 643 0 -277 0 -1 0 01933 -2 0 0 2 1 -199.84 -579 -1 304 0 -1 0 01934 2 0 2 0 1 6.85 -533 0 269 0 2 1 01935 0 -1 2 -2 1 346.60 -477 -1 271 -1 0 0 01936 0 0 0 -2 1 -14.73 -493 -1 272 0 -2 -1 01937 -1 -1 0 2 0 34.85 735 0 -5 0 -1 0 01938 2 0 0 -2 1 212.32 405 0 -220 0 1 0 01939 1 0 0 2 0 9.61 657 0 -20 0 -2 0 01940 0 1 2 -2 1 119.61 361 0 -194 0 1 0 01941 1 -1 0 0 0 29.80 471 0 -4 0 -1 0 01942 -2 0 2 0 2 1615.76 -311 0 131 0 0 0 01943 3 0 2 0 2 5.49 -289 0 124 0 2 1 01944 0 -1 0 2 0 15.39 435 0 -9 0 -1 0 01945 1 -1 2 0 2 9.37 -287 0 123 0 1 0 01946 -1 -1 2 2 2 9.81 -282 0 122 0 1 0 01947 0 0 0 1 0 29.53 -422 0 3 0 1 0 01948 -1 0 2 0 0 26.88 -404 0 4 0 1 0 01949 0 -1 2 2 2 7.24 -264 0 114 0 1 0 01950 -2 0 0 0 1 -13.75 -228 0 126 0 -1 0 01951 1 1 2 0 2 8.91 246 0 -106 0 -1 0 01952 2 0 0 0 1 13.81 218 0 -114 0 0 0 01953 -1 1 0 1 0 3232.87 327 0 -1 0 0 0 01954 1 1 0 0 0 25.62 -338 0 4 0 0 0 01955 1 0 2 0 0 9.11 334 0 -11 0 -1 0 01956 -1 0 2 -2 1 -32.61 -199 0 107 0 -1 0 01957 1 0 0 0 2 27.78 -197 0 85 0 0 0 01958 -1 0 0 1 0 -411.78 405 0 -55 0 -35 -14 01959 0 0 2 1 2 9.34 165 0 -72 0 0 0 01960 -1 0 2 4 2 5.80 -151 0 66 0 1 0 01961 0 -2 2 -2 1 6786.31 -130 0 69 0 0 0 01962 -1 1 0 1 1 6164.17 132 0 -68 0 0 0 01963 1 0 2 2 1 5.64 -133 0 66 0 1 0 01964 -2 0 2 2 2 14.63 139 0 -60 0 0 0 01965 -1 0 0 0 2 -27.33 139 0 -60 0 0 0 01966 1 1 2 -2 2 22.47 128 0 -55 0 0 0 01967 -2 0 2 4 2 7.35 -121 0 52 0 0 0 01968 -1 0 4 0 2 9.06 115 0 -49 0 0 0 01969 2 0 2 -2 1 12.79 101 0 -54 0 0 0 01970 2 0 2 2 2 4.68 -108 0 47 0 1 0 01971 1 0 0 2 1 9.63 -95 0 49 0 0 0 01972 3 0 0 0 0 9.18 157 0 -5 0 -1 0 01973 3 0 2 -2 2 8.75 94 0 -40 0 0 0 01974 0 0 4 -2 2 12.66 91 0 -39 0 0 0 01975 0 0 -2 2 1 -169.00 87 0 -44 0 0 0 01976 0 1 2 0 1 13.14 81 0 -42 0 0 0 01977 0 0 2 -2 3 187.66 123 0 -20 0 0 0 01978 -1 0 0 4 0 10.08 133 0 -4 0 0 0 01979 2 0 -2 0 1 -943.23 71 0 -38 0 0 0 01980 2 0 0 -4 0 -15.91 -128 0 1 0 0 0 01981 -1 -1 0 2 1 35.03 75 0 -39 0 0 0 01982 -2 -1 0 2 0 -131.67 -115 0 1 0 0 0 01983 0 -1 2 0 1 14.16 -66 0 35 0 0 0 01984 -1 0 0 1 1 -388.27 101 0 -49 0 -3 -1 01985 0 0 -2 0 1 -13.58 -68 0 36 0 0 0 01986 0 1 0 0 2 409.23 69 0 -33 0 -1 0 01987 0 0 2 -1 2 25.42 -74 0 31 0 0 0 01988 0 0 2 4 2 4.79 -69 0 29 0 0 0 01989 1 1 0 -2 1 -34.67 -61 0 32 0 0 0 01990 -1 1 0 2 0 29.26 -94 0 0 0 0 0 01991 1 -1 2 2 2 5.73 -59 0 25 0 0 0 01992 1 -1 0 0 1 29.93 51 0 -27 0 0 0 01993 0 1 -2 2 0 -329.79 -90 0 3 0 0 0 01994 3 0 2 0 1 5.49 -50 0 25 0 0 0 01995 -1 1 2 2 2 9.31 56 0 -24 0 0 0 01996 0 1 2 2 2 6.96 54 0 -22 0 0 0 01997 -1 0 0 -2 1 -9.60 -50 0 27 0 0 0 01998 -1 1 0 1 2 66079.30 -52 0 23 0 0 0 01999 0 -1 2 2 1 7.23 -44 0 24 0 0 0 02000 1 0 2 -4 1 -38.74 -47 0 24 0 0 0 02001 -1 0 -2 2 0 -23.77 77 0 0 0 0 0 02002 -1 -1 2 2 1 9.80 -46 0 24 0 0 0 02003 0 -1 0 0 2 -329.82 59 0 -25 0 0 0 02004 2 -1 2 0 2 6.99 -48 0 21 0 0 0 02005 1 -1 2 0 1 9.35 -42 0 22 0 0 0 02006 0 0 0 2 2 14.83 -46 0 20 0 0 0 02007 0 1 0 2 0 14.19 -67 0 0 0 0 0 02008 -1 1 2 0 2 25.22 47 0 -20 0 0 0 02009 0 3 2 -2 2 73.05 -44 0 19 0 0 0 02010 0 -1 -2 2 0 -117.54 66 0 0 0 0 0 02011 0 0 0 1 1 29.66 -37 0 20 0 0 0 02012 1 0 -2 -2 0 -9.53 -64 0 1 0 0 0 02013 1 1 2 0 1 8.90 36 0 -18 0 0 0 02014 2 1 2 0 2 6.73 40 0 -17 0 0 0 02015 0 1 0 1 0 27.32 57 0 0 0 0 0 02016 1 0 -2 2 0 32.76 -58 0 0 0 0 0 02017 1 1 0 0 1 25.72 -34 0 19 0 0 0 02018 -2 0 0 -2 0 -7.13 -59 0 1 0 0 0 02019 -1 0 0 2 2 32.11 -38 0 17 0 0 0 02020 0 0 0 -1 1 -29.40 33 0 -18 0 0 0 02021 0 1 0 -2 1 -15.35 -33 0 18 0 0 0 02022 -1 0 2 -2 2 -32.45 36 0 -16 0 0 0 02023 -1 1 0 0 1 -29.67 -31 0 17 0 0 0 02024 1 0 2 1 2 6.98 33 0 -14 0 0 0 02025 0 0 0 -4 0 -7.38 -48 0 1 0 0 0 02026 0 0 2 1 1 9.33 27 0 -14 0 0 0 02027 1 0 0 -2 2 -31.52 32 0 -14 0 0 0 02028 1 0 2 -1 2 13.22 -33 0 13 0 0 0 02029 1 -1 0 2 0 9.87 48 0 0 0 0 0 02030 -1 0 2 4 1 5.80 -26 0 13 0 0 0 02031 0 0 -2 -2 0 -7.08 -41 0 1 0 0 0 02032 1 0 -2 0 1 -26.77 27 0 -14 0 0 0 02033 -1 0 2 -1 1 313.04 -23 0 14 0 0 0 02034 1 1 2 -2 1 22.40 23 0 -12 0 0 0 02035 4 0 2 0 2 4.58 -26 0 11 0 0 0 02036 0 1 2 1 2 9.11 -24 0 10 0 0 0 02037 -2 0 -2 0 0 -6.85 -36 0 1 0 0 0 02038 2 1 2 -2 2 12.38 25 0 -10 0 0 0 02039 2 -1 0 0 0 14.32 38 0 0 0 0 0 02040 -1 -1 0 0 1 -25.53 21 0 -12 0 0 0 02041 -2 0 2 2 1 14.60 22 0 -11 0 0 0 02042 0 0 0 0 3 -2266.12 -22 0 10 0 0 0 02043 1 0 4 -2 2 8.68 23 0 -9 0 0 0 02044 2 0 2 2 1 4.68 -19 0 10 0 0 0 02045 -2 0 2 4 1 7.34 -20 0 10 0 0 0 02046 0 1 0 2 1 14.22 18 0 -9 0 0 0 02047 1 0 0 1 0 14.25 -33 0 0 0 0 0 02048 -1 0 0 4 1 10.10 -18 0 9 0 0 0 02049 -1 0 4 0 1 9.05 19 0 -9 0 0 0 02050 0 0 2 -3 2 -35.23 -20 0 8 0 0 0 02051 0 0 4 0 2 6.82 19 0 -8 0 0 0 02052 2 1 0 0 0 13.28 -28 0 0 0 0 0 02053 0 0 2 -4 1 -16.10 -16 0 9 0 0 0 02054 -1 -1 2 4 2 5.90 -17 0 7 0 0 0 02055 -1 -2 0 2 0 38.52 27 0 0 0 0 0 02056 0 0 0 4 1 7.39 -16 0 7 0 0 0 02057 0 -1 0 2 1 15.42 -14 0 7 0 0 0 02058 1 0 2 4 2 4.08 -16 0 7 0 0 0 02059 -2 0 0 2 2 -194.13 18 0 -8 0 0 0 02060 -2 2 0 2 0 1616.44 -22 0 0 0 0 0 02061 -2 -1 2 0 1 -507.16 9 0 -5 0 0 0 02062 -3 0 0 0 1 -9.17 -14 0 7 0 0 0 02063 0 0 2 0 3 13.69 20 0 0 0 0 0 02064 0 0 2 4 1 4.79 -12 0 6 0 0 0 02065 0 0 4 -2 1 12.64 12 0 -7 0 0 0 02066 0 -2 0 2 0 16.06 21 0 0 0 0 0 02067 1 0 0 -1 1 438.33 17 0 -5 0 -3 1 02068 1 1 2 2 2 5.56 15 0 -6 0 0 0 02069 3 0 2 -2 1 8.73 12 0 -7 0 0 0 02070 -1 -1 2 0 2 29.26 -16 0 6 0 0 0 02071 -2 -1 0 2 1 -129.17 -13 0 7 0 0 0 02072 0 0 0 -2 2 -14.70 13 0 -5 0 0 0 02073 0 -2 2 2 2 7.38 -13 0 5 0 0 0 02074 1 0 0 -4 1 -10.07 -12 0 6 0 0 0 02075 -1 1 0 2 1 29.39 -10 0 6 0 0 0 02076 -2 0 0 4 1 15.94 11 0 -6 0 0 0 02077 0 0 2 -1 1 25.33 -10 0 5 0 0 0 02078 0 2 0 0 1 187.67 -9 0 5 0 0 0 02079 0 2 2 -2 1 90.10 8 0 -5 0 0 0 02080 2 0 0 2 1 7.13 -9 0 5 0 0 0 02081 2 0 0 -4 1 -15.87 -11 0 5 0 0 0 02082 2 0 2 -4 1 95.42 10 0 -5 0 0 0 02083 -1 0 -2 0 1 -9.10 -10 0 5 0 0 0 02084 -1 1 2 0 1 25.13 9 0 -5 0 0 0 02085 -1 1 2 -2 1 -35.80 -11 0 5 0 0 0 02086 -1 -1 0 4 0 10.37 15 0 0 0 0 0 02087 -3 0 0 4 0 37.63 16 0 0 0 0 0 02088 3 0 2 2 2 4.00 -14 0 0 0 0 0 02089 2 -1 0 -2 0 471.89 -9 0 1 0 -1 0 02090 0 2 -2 2 0 -3396.16 -9 0 0 0 0 0 02091 0 -1 2 4 2 4.86 -9 0 0 0 0 0 02092 0 -1 2 -1 2 27.32 9 0 0 0 0 0 02093 1 1 0 2 0 9.37 -10 0 0 0 0 0 02094 2 0 0 -2 2 219.17 -11 0 0 0 0 0 02095 2 -1 2 2 2 4.74 -9 0 0 0 0 0 02096 4 0 0 0 0 6.89 9 0 0 0 0 0 02097 4 0 2 -2 2 6.64 12 0 0 0 0 0 02098 -1 0 0 3 0 15.31 -10 0 0 0 0 0 02099 -1 0 4 -2 2 23.43 -9 0 0 0 0 0 02100 -1 -2 2 2 2 10.08 -9 0 0 0 0 0 02101 -2 -1 0 4 0 16.63 12 0 0 0 0 0 02102 -2 -1 2 4 2 7.50 -12 0 0 0 0 0 02103 0 1 2 2 1 6.95 7 0 0 0 0 0 02104 0 2 2 0 2 12.71 7 0 0 0 0 0 02105 0 -2 2 0 2 14.77 -8 0 0 0 0 0 02106 1 0 0 4 0 5.82 8 0 0 0 0 0 02107 1 0 2 2 0 5.63 8 0 0 0 0 0 02108 1 0 2 -4 2 -38.52 7 0 0 0 0 0 02109 1 -1 2 2 1 5.73 -8 0 0 0 0 0 02110 1 -1 2 -2 2 25.62 -7 0 0 0 0 0 02111 1 -2 0 0 0 32.45 8 0 0 0 0 0 02112 2 0 0 0 2 13.83 -8 0 0 0 0 0 02113 2 1 0 -2 1 134.27 8 0 0 0 0 0 02114 3 0 0 0 1 9.20 7 0 0 0 0 0 02115 -1 0 2 1 2 14.13 8 0 0 0 0 0 02116 -1 0 2 3 2 7.22 8 0 0 0 0 0 02117 -1 0 -2 4 0 38.96 -7 0 0 0 0 0 02118 -1 1 2 2 1 9.30 7 0 0 0 0 0 02119 -1 2 0 2 0 27.09 -8 0 0 0 0 0 02120 -1 -1 2 -1 1 2189.73 7 0 0 0 0 0 02121 -2 0 2 -2 1 -14.93 -8 0 0 0 0 0 02122 -2 0 4 0 2 13.49 -7 0 0 0 0 0 02123 -2 0 -2 2 0 -12.76 8 0 0 0 0 0 02124 -2 1 2 0 1 285.41 9 0 0 0 0 0 02125 -3 0 2 0 1 -28.15 -8 0 0 0 0 0 02126 0 1 0 1 1 27.43 5 0 0 0 0 0 02127 0 -1 0 4 0 7.53 6 0 0 0 0 0 02128 0 -1 0 -2 1 -14.16 5 0 0 0 0 0 02129 0 -2 0 0 1 -177.85 -6 0 0 0 0 0 02130 1 0 2 1 1 6.97 5 0 0 0 0 0 02131 1 0 2 -3 2 126.51 -6 0 0 0 0 0 02132 1 0 -2 1 0 -299.26 -7 0 0 0 0 0 02133 1 1 0 1 0 13.72 5 0 0 0 0 0 02134 1 -1 0 -2 1 -29.14 6 0 0 0 0 0 02135 2 0 2 -1 2 8.93 -6 0 0 0 0 0 02136 2 1 2 0 1 6.73 5 0 0 0 0 0 02137 2 -1 2 0 1 6.98 -6 0 0 0 0 0 02138 2 -1 2 -2 2 13.28 5 0 0 0 0 0 02139 3 0 0 2 0 5.66 5 0 0 0 0 0 02140 3 -1 2 0 2 5.58 -5 0 0 0 0 0 02141 -1 -1 2 0 1 29.14 -6 0 0 0 0 0 02142 -2 0 0 0 2 -13.72 6 0 0 0 0 0 02143 -2 0 0 3 0 34.48 -5 0 0 0 0 0 02144 -2 0 0 -2 1 -7.12 -5 0 0 0 0 0 02145 -2 0 2 2 0 14.57 -6 0 0 0 0 0 02146 -2 -1 2 0 0 -548.04 -5 0 0 0 0 0 02147 -2 -1 2 2 2 15.24 6 0 0 0 0 0 02148 0 0 1 0 0 27.21 0 0 0 0 8 0 02149 0 0 1 0 1 27.32 0 0 0 0 -16 -14 02150 -1 0 1 0 0 2190.35 0 0 0 0 33 0 02151 -1 0 1 0 1 3231.51 0 0 0 0 -105 -89 02152 -1 0 1 0 2 6159.22 0 0 0 0 36 18 02153 -1 0 1 0 3 65514.10 0 0 0 0 -6 0 02154 # end of table 5.2 02155 */
1.3.9.1