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