1 //============================================================================== 2 // 3 // This file is part of GPSTk, the GPS Toolkit. 4 // 5 // The GPSTk is free software; you can redistribute it and/or modify 6 // it under the terms of the GNU Lesser General Public License as published 7 // by the Free Software Foundation; either version 3.0 of the License, or 8 // any later version. 9 // 10 // The GPSTk is distributed in the hope that it will be useful, 11 // but WITHOUT ANY WARRANTY; without even the implied warranty of 12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 // GNU Lesser General Public License for more details. 14 // 15 // You should have received a copy of the GNU Lesser General Public 16 // License along with GPSTk; if not, write to the Free Software Foundation, 17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA 18 // 19 // This software was developed by Applied Research Laboratories at the 20 // University of Texas at Austin. 21 // Copyright 2004-2020, The Board of Regents of The University of Texas System 22 // 23 //============================================================================== 24 25 //============================================================================== 26 // 27 // This software was developed by Applied Research Laboratories at the 28 // University of Texas at Austin, under contract to an agency or agencies 29 // within the U.S. Department of Defense. The U.S. Government retains all 30 // rights to use, duplicate, distribute, disclose, or release this software. 31 // 32 // Pursuant to DoD Directive 523024 33 // 34 // DISTRIBUTION STATEMENT A: This software has been approved for public 35 // release, distribution is unlimited. 36 // 37 //============================================================================== 38 39 40 /** 41 * @file PlanetEphemeris.hpp 42 * 43 */ 44 45 #ifndef GPSTK_PLANETEPHEMERIS_HPP 46 #define GPSTK_PLANETEPHEMERIS_HPP 47 48 #include <iostream> 49 #include <fstream> 50 #include <string> 51 #include <vector> 52 #include <map> 53 54 #include "Exception.hpp" 55 #include "CommonTime.hpp" 56 #include "MJD.hpp" 57 //#include "Position.hpp" 58 59 namespace gpstk 60 { 61 /** This class Handle planet ephemeris from JPL. 62 * 63 * @WARNING It's just a copy class of SolarSystem. 64 */ 65 class PlanetEphemeris 66 { 67 public: 68 /// These are indexes used by the caller of computeState(). 69 enum Planet { 70 // the following are relative to the solar system barycenter, except MOON 71 None=0, ///< 0 place holder 72 Mercury, ///< 1 Mercury 73 Venus, ///< 2 Venus 74 Earth, ///< 3 Earth 75 Mars, ///< 4 Mars 76 Jupiter, ///< 5 Jupiter 77 Saturn, ///< 6 Saturn 78 Uranus, ///< 7 Uranus 79 Neptune, ///< 8 Neptune 80 Pluto, ///< 9 Pluto 81 Moon, ///< 10 Moon (Geocentric coordinates) 82 Sun, ///< 11 Sun 83 SolarSystemBarycenter, ///< 12 Solar system barycenter 84 EarthMoonBarycenter, ///< 13 Earth-moon barycenter 85 Nutations, ///< 14 Nutations (psi, epsilon and their rates) 86 Librations ///< 15 Lunar Librations (3 euler angles) 87 }; 88 89 /// Constructor. Set EphemerisNumber to -1 to indicate that nothing has been 90 /// read yet. PlanetEphemeris(void)91 PlanetEphemeris(void) throw() : EphemerisNumber(-1) {}; 92 93 /// Read the header from a JPL ASCII planetary ephemeris file. Note that this 94 /// routine clears the 'store' map and defines the 'constants' hash. It also 95 /// sets EphemerisNumber to the constant "DENUM" if successful. 96 /// @param filename the name of the ASCII header file. 97 /// @throw if the file cannot be opened. 98 /// @throw if the header ends prematurely or if it is not properly formatted. 99 /// @throw Exception if any stream error occurs. 100 void readASCIIheader(std::string filename); 101 102 /// Read one or more ASCII data files. Call only after having read the header, 103 /// and call only with data files for the same ephemeris as the header (the JPL 104 /// files are named 'header.NNN' and 'ascSYYYY.NNN' where NNN is the ephemeris 105 /// number (appears inside the header file, but not inside the data files), 106 /// S is either 'p' or 'n' as the year is positive or negative (AD or BC), and 107 /// YYYY is the year of the first record in the file. 108 /// @param filenames vector containting the names of the ASCII data files 109 /// (downloaded from JPL), in any order. 110 /// @return 0 ok, -1 if a stream error occurred. 111 /// @throw Exception if the header has not yet been read. 112 /// @throw Exception if any file could not be opened. 113 /// @throw Exception if any record in any file has a 'number 114 /// of coefficients' that differs from the header 115 /// value. 116 int readASCIIdata(std::vector<std::string>& filenames); 117 118 /// Read only one ASCII data file. Also see the documentation for the 119 /// other version of this routine. 120 /// @param filename name of an ASCII data files (downloaded from JPL) 121 /// @return 0 ok, -1 if a stream error occurred, -4 header has not been read. 122 /// @throw Exception if the header has not yet been read. 123 /// @throw Exception if the file could not be opened. 124 /// @throw Exception if any record has a 'number of coefficients' that differs 125 /// from the header value. 126 int readASCIIdata(std::string filename); 127 128 /// Write the header (ASCII) to an output stream 129 /// @throw Exception if any stream error occurs 130 /// @return 0 success, 131 /// -4 header has not yet been read. 132 int writeASCIIheader(std::ostream& os); 133 134 /// Write the stored data (ASCII) to an output stream 135 /// NB. This routine does NOT clear the store - use clearStorage() 136 /// @throw Exception if any stream error occurs 137 /// @return 0 success, 138 /// -4 header has not yet been read. 139 int writeASCIIdata(std::ostream& os); 140 141 /// Write the header and the stored data to a binary output file. 142 /// NB. This routine does NOT clear the store - use clearStorage() 143 /// @param filename name of binary file to be read. 144 /// @return 0 ok 145 /// -4 if data has not been read into the object 146 /// @throw Exception if any stream error occurs. 147 int writeBinaryFile(std::string filename); 148 149 /// clear the store map containing all the data read by 150 /// readASCIIdata() or readBinaryData(true). clearStorage(void)151 void clearStorage(void) throw() { store.clear(); } 152 153 /// Read header and data from a binary file, storing ALL the data in store. 154 /// For use with copying, merging or editing data files. Closes the stream before 155 /// returning. 156 /// @param filename name of binary file to be read. 157 /// @return 0 success, 158 /// -3 input stream is not open or not valid 159 /// -4 header has not yet been read. 160 /// @throw Exception if a gap in time is found between consecutive records. 161 int readBinaryFile(std::string filename); 162 163 /// Open the given binary file, read the header and prepare for reading data 164 /// records at random using seekToJD() and computing positions and velocities 165 /// with computeState(). Does not store the data. 166 /// @param filename name of binary file to be read. 167 /// @return 0 success, 168 /// -3 input stream is not open or not valid 169 /// -4 header has not yet been read. 170 /// @throw Exception if a gap in time is found between consecutive records. 171 int initializeWithBinaryFile(std::string filename); 172 173 /// Compute position and velocity of given 'target' body, relative to the 'center' 174 /// body, at the given time. On successful return, PV contains position 175 /// (in components 0-2) and velocity (components 3-5) (units: see param km) for 176 /// regular bodies; for nutations and librations the units are radians and 177 /// radians/day; nutations (components 0-3 only) are longitude or psi (component 0) 178 /// and obliquity or epsilon (component 1) and their rates (components 2,3); 179 /// librations (components 0-5) are the 3 euler angles in radians and their 180 /// rates in radians/day. 181 /// @param tt Time (Julian Date) of interest. 182 /// @param target Body for which position and velocity are to be computed. 183 /// @param center Body relative to which the results apply. However, center 184 /// may == SolarSystem::None, in which case the results are 185 /// relative to the solar system barycenter. 186 /// If target == Nutations or Librations, center is ignored. 187 /// @param PV Double array of length 6 containing output position and velocity 188 /// components of target relative to center, in the order 189 /// X,Y,Z,Vx,Vy,Vz. Units are determined by parameter km. 190 /// If target == Nutations, PV contains 4 results, psi, eps, 191 /// psi dot, eps dot in units radians and radians/day. 192 /// If target == Librations, PV contains 3 euler angles in radians 193 /// and their rates in radians/day. 194 /// @param km boolean: if true (default), units are km, km/day; else AU, AU/day 195 /// (but not Nutations or Librations - see above). 196 /// @return 0 success, or (same as seekToJD()) 197 /// -1 given time is before the first record in the file, 198 /// -2 given time is after the last record, or in a gap between records, 199 /// -3 input stream is not open or not valid, or EOF was found prematurely, 200 /// -4 ephemeris is not initialized 201 /// -3 or -4 => initializeWithBinaryFile() has not been called, or reading failed. 202 int computeState(double tt, 203 Planet target, 204 Planet center, 205 double PV[6], 206 bool kilometers = true); 207 208 /// Return the value of 1 AU (Astronomical Unit) in km. If the file header has not 209 /// been read, return -1.0. 210 /// @return the value of 1 AU in km; 211 /// return -1 if ephemeris has not been initialized. AU(void)212 double AU(void) throw() 213 { if(EphemerisNumber == -1) return -1.0; return constants["AU"]; } 214 215 /// Return the ephemeris number. 216 /// @return JPL ephemeris number, e.g. 403, 217 /// -1 if ephemeris has not been initialized. JPLNumber(void) const218 int JPLNumber(void) const throw() 219 { return EphemerisNumber; } 220 221 /// @return the value of the contant with the given name. If the header 222 /// has not been read, return -1. Return zero if the constant is not found. getConstant(std::string name)223 double getConstant(std::string name) throw() { 224 if(EphemerisNumber == -1) return -1.0; 225 if(constants.find(name) != constants.end()) return constants[name]; 226 return 0.0; 227 } 228 229 /// Return the start time of the data startTime(void) const230 gpstk::CommonTime startTime(void) const 231 { return gpstk::MJD(startJD - gpstk::MJD_TO_JD); } 232 233 /// Return the end time of the data endTime(void) const234 gpstk::CommonTime endTime(void) const 235 { return gpstk::MJD(endJD - gpstk::MJD_TO_JD); } 236 237 238 private: 239 /// Helper routine for binary writing. 240 /// @throw Exception if there is any stream error. 241 void writeBinary(std::ofstream& strm, const char *ptr, size_t size); 242 243 /// Helper routine for binary reading. 244 /// @throw Exception if there is any error or if EOF is found. 245 void readBinary(char *ptr, size_t size); 246 247 /// Read header from a binary file. 248 /// @param filename name of binary file, probably written by writeBinaryFile(). 249 /// @throw Exception if read error or premature EOF if found. 250 void readBinaryHeader(std::string filename); 251 252 /// Read data from a binary file, already opened by readBinaryHeader. 253 /// Build the file position map, and store the first set of coefficients. 254 /// If calling argument is true, save all the coefficient data in a map. 255 /// @param save if true, save all the data in store, else clear the store. 256 /// @return 0 success, 257 /// -3 input stream is not open or not valid 258 /// -4 header has not yet been read. 259 /// @throw Exception if a gap in time is found between consecutive records. 260 int readBinaryData(bool save); 261 262 /// Read a single binary record (not a header record) at the current file 263 /// position, into the given vector. For use by readBinaryData() and seekToJD(). 264 /// @param data_vector vector<double> to hold coefficients. 265 /// @return 0 success, 266 /// -2 EOF was reached 267 /// -3 input stream is not open or not valid 268 int readBinaryRecord(std::vector<double>& data_vector); 269 270 /// These are indexes used in the actual computation, and correspond to indexes 271 /// in the ephemeris file; for example computation for the SUN is done using 272 /// c_offset[SUN], c_ncoeff[SUN] and c_nsets[SUN]. 273 enum computeID { 274 // the following are relative to the solar system barycenter, except MOON 275 NONE=-1, ///< -1 Place holder 276 MERCURY, ///< 0 Mercury 277 VENUS, ///< 1 Venus 278 EMBARY, ///< 2 Earth-Moon barycenter 279 MARS, ///< 3 Mars 280 JUPITER, ///< 4 Jupiter 281 SATURN, ///< 5 Saturn 282 URANUS, ///< 6 Uranus 283 NEPTUNE, ///< 7 Neptune 284 PLUTO, ///< 8 Pluto 285 MOON, ///< 9 Moon (Geocentric coordinates) 286 SUN, ///< 10 Sun 287 NUTATIONS, ///< 11 Nutations (psi, epsilon and their rates) 288 LIBRATIONS ///< 12 Lunar Librations (3 euler angles) 289 }; 290 291 /// Search the data records of the file opened by initializeWithBinaryFile() and 292 /// read the one whose time limits include the given time. May be called only 293 /// after initializeWithBinaryFile(). 294 /// @param JD the time (Julian Date) of interest 295 /// @return 0 success, or 296 /// -1 given time is before the first record in the file, 297 /// -2 given time is after the last record, or in a gap between records, 298 /// -3 input stream is not open or not valid, or EOF was found prematurely, 299 /// -4 ephemeris (binary file) is not initialized 300 /// -3 or -4 => initializeWithBinaryFile() has not been called, or reading failed. 301 int seekToJD(double JD); 302 303 /// Compute position and velocity of given body at given time, using the current 304 /// coefficient array. NB caller MUST call seekToJD(time) BEFORE calling this. 305 /// On successful return, PV[0-2] contains the three position components, in km, 306 /// and PV[3-5] the velocity components in km/day (for regular bodies), relative 307 /// to the solar system barycenter, except for the moon, which is relative to 308 /// Earth. For nutations and librations the units are radians and radians/day; 309 /// nutations (components 0-3 only) are longitude and obliquity, and librations 310 /// are the three euler angles. 311 /// @param tt Time (Julian Date) of interest. 312 /// @param which computeID of the body of interest. 313 /// @param PV double(6) array containing the output position and velocity. 314 void computeState(double tt, computeID which, double PV[6]); 315 316 // member data --------------------------------------------------------- 317 318 /// input stream, for use by readBinary...() 319 std::ifstream istrm; ///< input stream for binary files 320 321 /// header information 322 /// -1 if the header has not been filled; also, for binary file input, 0 if 323 /// the file position map has not yet been filled; otherwise it equals the 324 /// number JPL assigns the ephemeris, e.g. 403, 405, which is identical to 325 /// constants["DENUM"]. 326 int EphemerisNumber; 327 328 /// The number of coefficients in a single record. This will determine the 329 /// binary record size 330 int Ncoeff; 331 int Nconst; ///< number of constants in the header (see map constants) 332 std::string label[3]; ///< lines under group 1010 333 double startJD; ///< JD of the start of the first record in the file 334 double endJD; ///< JD of the end of the last record in the file 335 double interval; ///< number of days covered by each block of coeff.s 336 int c_offset[13]; ///< starting index in the coefficients array for each planet 337 int c_ncoeff[13]; ///< number of coefficients per component for each planet 338 int c_nsets[13]; ///< number of sets of coefficients for each planet 339 340 /// Hash of labels and values of constants read from the header. 341 /// This is taken directly from the JPL documentation: 342 /// <pre> 343 /// The following is a partial list of constants found on the ephemeris file: 344 /// DENUM Planetary ephemeris number. 345 /// LENUM Lunar ephemeris number. 346 /// TDATEF, TDATEB Dates of the Forward and Backward Integrations 347 /// CLIGHT Speed of light (km/s). 348 /// AU Number of kilometers per astronomical unit. 349 /// EMRAT Earth-Moon mass ratio. 350 /// GMi GM for ith planet [au**3/day**2]. 351 /// GMB GM for the Earth-Moon Barycenter [au**3/day**2]. 352 /// GMS Sun (= k**2) [au**3/day**2]. 353 /// X1, ..., ZD9 Initial conditions for the numerical integration, 354 /// given at "JDEPOC", with respect to "CENTER". 355 /// JDEPOC Epoch (JED) of initial conditions, normally JED 2440400.5. 356 /// CENTER Reference center for the initial conditions. 357 /// (Sun: 11, Solar System Barycenter: 12) 358 /// RADi Radius of ith planet [km]. 359 /// MA0001...MA0324 GM's of asteroid number 0001 ... 0234 [au**3/day**2]. 360 /// PHASE The phase angle of the moon's rotation. 361 /// LOVENO The Love Number, k2, for the moon. 362 /// PHI, THT, PSI Euler angles of the orientation of the lunar mantle. 363 /// OMEGAX, ... Rotational velocities of the lunar mantle. 364 /// PHIC,THTC,PSIC Euler angles of the orientation of the lunar core. 365 /// OMGCX, ... Rotational velocities of the lunar core. 366 /// </pre> 367 std::map<std::string,double> constants; 368 369 /// Hash of start times (JD) and coefficient arrays. 370 /// This object is not intended to store the entire dataset, except temporarily 371 /// for the purpose of reading/writing files, NOT for ephemeris computation. 372 std::map<double, std::vector<double> > store; 373 374 /// Hash of start times (JD) and file positions. This object is filled by 375 /// readBinaryData(), which is called by initializeWithBinaryFile(), and is 376 /// used by seekToJD() to read records in random order. 377 std::map<double, long> fileposMap; 378 379 /// One complete data record (Ncoeff doubles) consisting of times and coefficients. 380 /// seekToJD() stores the current record here, and computeState() makes use of it. 381 std::vector<double> coefficients; 382 383 }; // end class PlanetEphemeris 384 385 } // End of namespace gpstk 386 387 388 #endif //GPSTK_PLANETEPHEMERIS_HPP 389