00001 //# MeasTable.h: MeasTable provides Measure computing database data 00002 //# Copyright (C) 1995-1999,2000-2004 00003 //# Associated Universities, Inc. Washington DC, USA. 00004 //# 00005 //# This library is free software; you can redistribute it and/or modify it 00006 //# under the terms of the GNU Library General Public License as published by 00007 //# the Free Software Foundation; either version 2 of the License, or (at your 00008 //# option) any later version. 00009 //# 00010 //# This library is distributed in the hope that it will be useful, but WITHOUT 00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00013 //# License for more details. 00014 //# 00015 //# You should have received a copy of the GNU Library General Public License 00016 //# along with this library; if not, write to the Free Software Foundation, 00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00018 //# 00019 //# Correspondence concerning AIPS++ should be addressed as follows: 00020 //# Internet email: aips2-request@nrao.edu. 00021 //# Postal address: AIPS++ Project Office 00022 //# National Radio Astronomy Observatory 00023 //# 520 Edgemont Road 00024 //# Charlottesville, VA 22903-2475 USA 00025 //# 00026 //# 00027 //# $Id$ 00028 00029 #ifndef MEASURES_MEASTABLE_H 00030 #define MEASURES_MEASTABLE_H 00031 00032 //# Includes 00033 #include <casacore/casa/aips.h> 00034 #include <casacore/measures/Measures/MeasTableMul.h> 00035 #include <casacore/measures/Measures/MeasData.h> 00036 #include <casacore/measures/Measures/MPosition.h> 00037 #include <casacore/measures/Measures/MDirection.h> 00038 #include <casacore/measures/Measures/MFrequency.h> 00039 #include <casacore/scimath/Functionals/Polynomial.h> 00040 #include <casacore/casa/OS/Mutex.h> 00041 #include <vector> 00042 00043 namespace casacore { //# NAMESPACE CASACORE - BEGIN 00044 00045 //# Forward Declarations 00046 class RotMatrix; 00047 class Euler; 00048 00049 // <summary> 00050 // MeasTable provides Measure computing database data 00051 // </summary> 00052 00053 // <use visibility=local> 00054 00055 // <reviewed reviewer="UNKNOWN" date="before2004/08/25" tests="tMeasMath" demos=""> 00056 // </reviewed> 00057 00058 // <prerequisite> 00059 // <li> <linkto class=Measure>Measure</linkto> class 00060 // <li> <linkto class=MeasData>MeasData</linkto> class for constant data 00061 // <li> <linkto class=Aipsrc>Aipsrc</linkto> class for data placement 00062 // </prerequisite> 00063 // 00064 // <etymology> 00065 // MeasTable from Measure and Table 00066 // </etymology> 00067 // 00068 // <synopsis> 00069 // MeasTable contains the database interface for all 00070 // data necessary for precession, nutation and other 00071 // <linkto class=Measure>Measure</linkto> related calculations.<br> 00072 // All data are obtained by calls to a method. E.g. 00073 // <src> fundArg(1) </src> will provide the first fundamental argument for 00074 // nutation calculations, i.e. 'l'. <br> 00075 // This class contains no constructors or destructors, only static 00076 // methods and (static) constants. 00077 // <br> References:<br> Explanatory supplements to the Astronomical Almanac 00078 // <br> C. Ron and J. Vondrak, Bull. Astron. Inst. Czechosl. 37, p96, 1986 00079 // <br> M. Soma, Th. Hirayama and H. Kinoshita, Celest. Mech. 41, p389, 1988 00080 // <br> V.S. Gubanov, Astron. Zh. 49, p1112, 1972 00081 // 00082 // Where strings are passed in as arguments (observatory names, sources), they 00083 // will be case insensitive, and minimum match. 00084 // </synopsis> 00085 // 00086 // <example> 00087 // Usage examples can be found in <linkto class=Precession>Precession</linkto> 00088 // </example> 00089 // 00090 // <motivation> 00091 // To create a clean interface between the actual calculations and the 00092 // methods to obtain the parameters for these calculations. Note that the 00093 // tables are in general in the format and units found in the literature. This 00094 // is to be able to easy check and change them. However, in the future 00095 // re-arrangement could produce faster and more compact code. 00096 // </motivation> 00097 // 00098 // <todo asof="1997/09/02"> 00099 // <li> more database interfaces, rather than constants 00100 // </todo> 00101 00102 class MeasTable { 00103 00104 public: 00105 00106 //# Enumerations 00107 // Types to be used in different calls 00108 enum Types { 00109 // Planetary information 00110 MERCURY = 1, 00111 VENUS = 2, 00112 EARTH = 3, 00113 MARS = 4, 00114 JUPITER = 5, 00115 SATURN = 6, 00116 URANUS = 7, 00117 NEPTUNE = 8, 00118 PLUTO = 9, 00119 MOON = 10, 00120 SUN = 11, 00121 // Solar system barycentre 00122 BARYSOLAR = 12, 00123 // Earth-Moon system barycentre 00124 BARYEARTH = 13, 00125 // Nutations 00126 NUTATION = 14, 00127 // Librations 00128 LIBRATION = 15, 00129 // Number of types 00130 N_Types }; 00131 00132 // Codes for JPL constants: order should be same as in MeasJPL, length less 00133 // than or equal 00134 enum JPLconst { 00135 // Light velocity used in AU/d 00136 CAU, 00137 // Solar mass (GM0)/c<sup>2</sup> in AU 00138 GMS, 00139 // AU in km 00140 AU, 00141 // Solar radius in AU 00142 RADS, 00143 // # of codes 00144 N_JPLconst }; 00145 00146 //# General Member Functions 00147 // Selection related data 00148 // <group> 00149 // Are the IAU2000 precession/nutation to be used or not (IAU1984) 00150 static Bool useIAU2000(); 00151 // If IAU2000 model, do we use the high precision 2000A model? 00152 static Bool useIAU2000A(); 00153 // </group> 00154 00155 // Precession related data 00156 // <group> 00157 // Get the precession-rate part of the IAU2000 precession-nutation models 00158 // (which 0=dpsi (long) and 1=deps (obliquity) and 2 =0) 00159 static Double precRate00(const uInt which); 00160 00161 // Get the frame bias matrix for IAU2000 model. 00162 static RotMatrix frameBias00(); 00163 00164 // Generate the precession calculation polynomials for a fixed Epoch T 00165 // in the result area specified. 00166 // T is given in Julian centuries since J2000.0. 00167 static void 00168 precessionCoef(Double T, Polynomial<Double> result[3]); 00169 00170 // Generate the precession polynomials for IAU2000 system. 00171 static void 00172 precessionCoef2000(Polynomial<Double> result[3]); 00173 00174 // Generate the precession polynomials for 1950 system for a fixed Epoch T 00175 // in the area specified. T is given in Tropical centuries since B1850.0 00176 static void 00177 precessionCoef1950(Double T, Polynomial<Double> result[3]); 00178 // </group> 00179 00180 // Nutation related data 00181 // <group> 00182 // Generate the polynomial for the fundamental arguments (eps, l, l', 00183 // F, D, omega) as a function of Julian centuries 00184 // <group> 00185 static const Polynomial<Double> &fundArg(uInt which); 00186 static const Polynomial<Double> &fundArg1950(uInt which); 00187 static const Polynomial<Double> &fundArg2000(uInt which); 00188 // </group> 00189 00190 // Get the planetary arguments (L, L', F, D, Om, Me, Ve, E, Ma, Ju Sa, 00191 // Ur, Ne, pre) 00192 static const Polynomial<Double> &planetaryArg2000(uInt which); 00193 00194 // Generate the which' vector of the nutation series arguments 00195 // <group> 00196 static const Double* mulArg(uInt which); 00197 static const Double* mulArg1950(uInt which); 00198 static const Double* mulArg2000A(uInt which); 00199 static const Double* mulArg2000B(uInt which); 00200 static const Double* mulPlanArg2000A(uInt which); 00201 // </group> 00202 00203 // Generate the which' vector of the equation of equinoxes (IAU2000) 00204 // complementary terms series arguments 00205 static const Double* mulArgEqEqCT2000(uInt which); 00206 00207 // Generate the which' vector of the nutation series multipliers 00208 // at T, measured in Julian centuries since J2000.0, respectively B1900.0 00209 // <group> 00210 static CountedPtr<Matrix<Double> > mulSC(Double time, Double epsilon); 00211 static CountedPtr<Matrix<Double> > mulSC1950(Double time, Double epsilon); 00212 static CountedPtr<Matrix<Double> > mulSC2000A(Double time, Double epsilon); 00213 static CountedPtr<Matrix<Double> > mulSC2000B(Double time, Double epsilon); 00214 static const Double* mulPlanSC2000A(uInt which); 00215 // </group> 00216 00217 // Generate the which' vector of the equation of equinoxes (IAU2000) 00218 // complementary terms series multipliers 00219 // at T, measured in Julian centuries since J2000.0, respectively B1900.0 00220 static const Double* mulSCEqEqCT2000(uInt which); 00221 00222 // Get nutation angles corrections for UTC T in rad. 00223 // which = 0 : dPsi as given by IERS for IAU nutation theory; 00224 // = 1: dEps as same. 00225 static Double dPsiEps(uInt which, Double T); 00226 // </group> 00227 00228 // Planetary (JPL DE) related data 00229 // <group> 00230 // Get the position (AU or rad) and velocity (AU/d or rad/d) for specified 00231 // code at TDB T. The ephemeris to use (now DE200 or DE405) can be selected 00232 // with the 'measures.jpl.ephemeris' aipsrc resource (default DE200). 00233 static Vector<Double> Planetary(MeasTable::Types which, 00234 Double T); 00235 // Get the JPL DE constant indicated 00236 static Double Planetary(MeasTable::JPLconst what); 00237 // </group> 00238 00239 // Observatory positions 00240 // <group> 00241 // Initialise list of all observatories from Observatories table 00242 static void initObservatories(); 00243 // Get list of all observatories 00244 static const Vector<String> &Observatories(); 00245 // Get position of observatory nam (False if not present) 00246 static Bool Observatory(MPosition &obs, const String &nam); 00247 00248 // Get _absolute_ path to AntennaResponses table of observatory 00249 // <src>nam</src>. It returns False if no _valid_ path can be found or the 00250 // observatory is unknown. If the observatory is known, antRespPath will 00251 // be set to the entry in the AntennaResponses column of the 00252 // Observatories table even if it doesn't describe a valid path; if the 00253 // entry is not an absolute path, the data directory name will be 00254 // prepended and validity verified. 00255 static Bool AntennaResponsesPath(String &antRespPath, const String &nam); 00256 // </group> 00257 00258 // Source list positions 00259 // <group> 00260 // Initialise list of all source from Sources table 00261 static void initSources(); 00262 // Get list of all sources 00263 static const Vector<String> &Sources(); 00264 // Get position of source <src>nam</src> (False if not present) 00265 static Bool Source(MDirection &obs, const String &nam); 00266 // </group> 00267 00268 // Rest frequencies 00269 // <group> 00270 // Initialise list from internal Table for now 00271 static void initLines(); 00272 // Get list of all frequencies 00273 static const Vector<String> &Lines(); 00274 // Get frequency of line name (False if not present) 00275 static Bool Line(MFrequency &obs, const String &nam); 00276 // </group> 00277 00278 // Initialise list of IGRF data 00279 static void initIGRF(); 00280 // Earth magnetic field (IGRF) data 00281 // Get the harmonic terms for specified time (mjd) 00282 static Vector<Double> IGRF(Double t); 00283 00284 // Aberration related data 00285 // <group> 00286 // Generate the polynomial for the fundamental arguments (l1-l8, w, D, l, 00287 // l', F) for the Ron/Vondrak aberration calculations as a function of 00288 // Julian centuries(J2000), or the comparable ones for the Gubanov expansion 00289 // (B1950). 00290 // <group> 00291 static const Polynomial<Double> &aberArg(uInt which); 00292 static const Polynomial<Double> &aberArgDeriv(uInt which); 00293 static const Polynomial<Double> &aber1950Arg(uInt which); 00294 static const Polynomial<Double> &aber1950ArgDeriv(uInt which); 00295 // </group> 00296 00297 // Generate the 'which' vector of the aberration series arguments 00298 // <group> 00299 static const Double* mulAberArg(uInt which); 00300 static const Double* mulAber1950Arg(uInt which); 00301 static const Double* mulAberSunArg(uInt which); 00302 static const Double* mulAberEarthArg(uInt which); 00303 // </group> 00304 00305 // Generate the 'which' vector of the aberration series multipliers 00306 // at T, measured in Julian centuries since J2000.0 (or J1900.0, yes, 00307 // J1900.0, for B1950). 00308 // <group> 00309 static CountedPtr<Matrix<Double> > mulAber(Double time, Double epsilon); 00310 static CountedPtr<Matrix<Double> > mulAber1950(Double time, Double epsilon); 00311 static const Vector<Double> &mulSunAber(uInt which); 00312 static const Vector<Double> &mulEarthAber(uInt which); 00313 // </group> 00314 00315 // Get the E-terms of Aberration correction (0 for position, 1 for velocity) 00316 // <group> 00317 static const Vector<Double> &AberETerm(uInt which); 00318 // </group> 00319 00320 // </group> 00321 00322 // Diurnal aberration factor 00323 static Double diurnalAber(Double radius, Double T); 00324 00325 // LSR (kinematical) velocity conversion: 0 gives J2000; 1 gives B1950. 00326 // In both cases a velocity of 20.0 km/s is assumed, and a B1900 RA/Dec 00327 // direction of (270,30) degrees. This value has been defined between 00328 // the groups doing HI radio work in the mid 1950s. 00329 static const Vector<Double> &velocityLSRK(uInt which); 00330 // LSR (dynamical, IAU definition). Velocity (9,12,7) km/s in galactic 00331 // coordinates. Or 16.552945 towards l,b = 53.13, +25.02 deg. 00332 // 0 gives J2000, 1 gives B1950 velocities. 00333 static const Vector<Double> &velocityLSR(uInt which); 00334 // Velocity of LSR with respect to galactic centre. 220 km/s in direction 00335 // l,b = 270, +0 deg. 0 returns J2000, 1 B1950 00336 static const Vector<Double> &velocityLSRGal(uInt which); 00337 // Velocity of Local Group wrt bary center (F.Ghigo): 308km/s towards 00338 // l,b = 105,-7. 0 for J2000, 1 for B1950 00339 static const Vector<Double> &velocityCMB(uInt which); 00340 // Velocity of CMB wrt bary center (F.Ghigo): 369.5km/s towards 00341 // l,b = 264.4,48.4. 0 for J2000, 1 for B1950 00342 00343 static const Vector<Double> &velocityLGROUP(uInt which); 00344 // Earth and Sun position related data 00345 // <group> 00346 // Fundamental arguments for Soma et al. methods 00347 // <group> 00348 static const Polynomial<Double> &posArg(uInt which); 00349 // Precomputed derivative of PosArg 00350 static const Polynomial<Double> &posArgDeriv(uInt which); 00351 // </group> 00352 // Generate the which' vector of the position series arguments 00353 // <group> 00354 static const Double* mulPosEarthXYArg(uInt which); 00355 static const Double* mulPosEarthZArg(uInt which); 00356 static const Double* mulPosSunXYArg(uInt which); 00357 static const Double* mulPosSunZArg(uInt which); 00358 // </group> 00359 00360 // Generate the which' vector of the position series multipliers 00361 // at T, measured in Julian centuries since J2000.0 00362 // <group> 00363 static CountedPtr<Matrix<Double> > mulPosEarthXY(Double time, Double epsilon); 00364 static CountedPtr<Matrix<Double> > mulPosEarthZ (Double time, Double epsilon); 00365 static CountedPtr<Matrix<Double> > mulPosSunXY (Double time, Double epsilon); 00366 static CountedPtr<Matrix<Double> > mulPosSunZ (Double time, Double epsilon); 00367 // </group> 00368 // Get the rotation matrix to change position from ecliptic to rectangular 00369 // for Soma et al. analytical expression 00370 static const RotMatrix &posToRect(); 00371 // Get the rotation matrix to change position from rectangular to ecliptic 00372 // for Soma et al. analytical expression 00373 static const RotMatrix &rectToPos(); 00374 // Get the rotation matrix from galactic to supergalactic. 00375 // Based on De Vaucouleurs 1976: Pole at 47.37/6.32 deg; 137.37 l0 00376 // Euler angles: 90, 83.68, 47.37 degrees 00377 static const RotMatrix &galToSupergal(); 00378 // Get the rotation matrix from ICRS to J2000/FK5. 00379 // Based on the IAU 2000 resolutions (the bias matrix) 00380 static const RotMatrix &ICRSToJ2000(); 00381 // </group> 00382 00383 // Position related routines 00384 // <group> 00385 // Equatorial radius (0) and flattening(1) of geodetic reference spheroids 00386 static Double WGS84(uInt which); 00387 // </group> 00388 00389 // Polar motion related routines 00390 // <group> 00391 // Get the polar motion (-x,-y,0)(2,1,3) angles at the given epoch 00392 static Euler polarMotion(Double ut); 00393 // </group> 00394 00395 // Time related routines 00396 // <note> 00397 // WARNING given if correction not obtainable 00398 // </note> 00399 // <thrown> 00400 // <li> AipsError if table seems to be corrupted 00401 // </thrown> 00402 // <group> 00403 // Give TAI-UTC (in s) for MJD utc UTC 00404 static Double dUTC(Double utc); 00405 // UT1-UTC (in s) for MJD tai TAI 00406 static Double dUT1(Double utc); 00407 // TDT-TAI (in s) for MJD tai TAI. Note this is equal to TT2000-TAI 00408 static Double dTAI(Double tai=0.0); 00409 // TDB-TDT (in s) for MJD ut1 UT1 00410 static Double dTDT(Double ut1); 00411 // TCB-TDB (in s) for MJD tai TAI 00412 static Double dTDB(Double tai); 00413 // TCG-TT (in s) for MJD tai TAI 00414 static Double dTCG(Double tai); 00415 // GMST1 at MJD ut1 UT1 00416 static Double GMST0(Double ut1); 00417 // GMST (IAU2000) including the ERA (IAU2000 Earth Rotation Angle) in rad 00418 static Double GMST00(Double ut1, Double tt); 00419 // Earth Rotation Angle (IAU2000) in rad 00420 static Double ERA00(Double ut1); 00421 // s' (IAU2000) in rad (approximate value) 00422 static Double sprime00(Double tt); 00423 // UT1 at GMSD gmst1 GMST1 00424 static Double GMUT0(Double gmst1); 00425 // Ratio UT1/MST at MJD ut1 UT1 00426 static Double UTtoST(Double ut1); 00427 // </group> 00428 00429 private: 00430 00431 // Copy assign, NOT defined 00432 MeasTable &operator=(const MeasTable &other); 00433 00434 //# General member functions 00435 00436 static void doInitObservatories (void*); 00437 static void doInitLines (void*); 00438 static void doInitSources (void*); 00439 static void doInitIGRF (void*); 00440 00441 // Calculate precessionCoef 00442 // <group> 00443 static void calcPrecesCoef(Double T, Polynomial<Double> result[3], 00444 const Double coeff[3][6]); 00445 static void calcPrecesCoef2000(Polynomial<Double> result[3], 00446 const Double coeff[3][6]); 00447 // </group> 00448 00449 // Calculate fundArg 00450 // <group> 00451 static void calcFundArg(volatile Bool &need, Polynomial<Double> result[6], 00452 const Double coeff[6][4]); 00453 static void calcFundArg00(volatile Bool &need, Polynomial<Double> result[6], 00454 const Double coeff[6][5]); 00455 static void calcPlanArg00(volatile Bool &need, 00456 Polynomial<Double> result[14], 00457 const Double coeff[8][2]); 00458 // </group> 00459 00460 00461 //# Data 00462 // Multipliers for nutation, etc. 00463 // <group> 00464 static MeasTableMulSC theirMulSC; 00465 static MeasTableMulSC1950 theirMulSC1950; 00466 static MeasTableMulSC2000A theirMulSC2000A; 00467 static MeasTableMulSC2000B theirMulSC2000B; 00468 static MeasTableMulAber theirMulAber; 00469 static MeasTableMulAber1950 theirMulAber1950; 00470 static MeasTableMulPosSunXY theirMulPosSunXY; 00471 static MeasTableMulPosSunZ theirMulPosSunZ; 00472 static MeasTableMulPosEarthXY theirMulPosEarthXY; 00473 static MeasTableMulPosEarthZ theirMulPosEarthZ; 00474 // </group> 00475 // Observatories table data 00476 // <group> 00477 static MutexedInit obsMutexedInit; 00478 static Vector<String> obsNams; 00479 static Vector<MPosition> obsPos; 00480 static Vector<String> antResponsesPath; 00481 // </group> 00482 // Spectral line table data 00483 // <group> 00484 static MutexedInit lineMutexedInit; 00485 static Vector<String> lineNams; 00486 static Vector<MFrequency> linePos; 00487 // </group> 00488 // Sources table data 00489 // <group> 00490 static MutexedInit srcMutexedInit; 00491 static Vector<String> srcNams; 00492 static Vector<MDirection> srcPos; 00493 // </group> 00494 // IGRF data 00495 // <group> 00496 static MutexedInit igrfMutexedInit; 00497 static Double dtimeIGRF; 00498 static Double firstIGRF; 00499 static Double lastIGRF; 00500 static Double time0IGRF; 00501 static Double timeIGRF; 00502 static std::vector<Vector<Double> > coefIGRF; 00503 static std::vector<Vector<Double> > dIGRF; 00504 // </group> 00505 // Aipsrc registration (for speed) of use of iau2000 and if so 00506 // the 2000a version 00507 // <group> 00508 static uInt iau2000_reg; 00509 static uInt iau2000a_reg; 00510 // </group> 00511 // Mutex for thread-safety. 00512 static Mutex theirMutex; 00513 }; 00514 00515 00516 } //# NAMESPACE CASACORE - END 00517 00518 #endif