C++ SolarSystem::JPLNumber方法代码示例

// -----------------------------------------------------------------------------------
// Version with JPL SolarSystem ephemeris. Throw if the SolarSystem is not valid
Matrix<double> SatelliteAttitude(const CommonTime& tt, const Position& SV,
                                 const SolarSystem& SSEph, const EarthOrientation& EO,
                                 double& sf)
   throw(Exception)
{
   if(SSEph.JPLNumber() == -1 || SSEph.startTime()-tt > 1.e-8
                              || tt - SSEph.endTime() > 1.e-8) {
      Exception e("Solar system ephemeris invalid");
      GPSTK_THROW(e);
   }

   return doSatAtt(tt,SV,SSEph,EO,sf);
}

//.........这里部分代码省略.........
   }
   if(dataFilenames.size() == 0) {
      LOG(ERROR) << "Must specify data file name(s)";
      ok = false;
   }
   if(outputFilename.empty()) {
      LOG(ERROR) << "Must specify an output file name";
      ok = false;
   }
   if(!ok) return -1;

   // set up the log file
   if(!logFilename.empty()) {
      // choose the log file
      logstrm.open(logFilename.c_str(),ios_base::out);
      ConfigureLOG::Stream() = &logstrm;
      // if not the above, output is to stderr

      cout << Title << endl;
      cout << "Output logged in file " << logFilename << endl;
   }
      // set the maximum level to be logged
   if(debug)
      ConfigureLOG::ReportingLevel() = ConfigureLOG::FromString("DEBUG");
   else if(verbose)
      ConfigureLOG::ReportingLevel() = ConfigureLOG::FromString("VERBOSE");
   else
      ConfigureLOG::ReportingLevel() = ConfigureLOG::FromString("INFO");
      // = any of ERROR,WARNING,INFO,DEBUG,DEBUGn (n=1,2,3,4,5,6,7)
   //cout << "Reporting in main is "
   //   << ConfigureLOG::ToString(ConfigureLOG::ReportingLevel()) << endl;
   ConfigureLOG::ReportLevels() = ConfigureLOG::ReportTimeTags() = true;

   // display title in log file
   LOG(INFO) << Title;

   // read header file
   eph.readASCIIheader(headerFilename);
   LOG(VERBOSE) << "Finished reading ASCII header " << headerFilename;
   LOG(INFO) << "Ephemeris number from header is " << eph.JPLNumber();

   // read the data files
   eph.readASCIIdata(dataFilenames);
   for(j=0; j<dataFilenames.size(); j++)
      LOG(VERBOSE) << "Finished reading ASCII data " << dataFilenames[j];
   LOG(INFO) << "Ephemeris number from data is " << eph.JPLNumber();

   // dump to a file
   LOG(INFO) << "Dump ASCII header to csse.header.asc";
   ofstream ofs;
   ofs.open("csse.header.asc",ios_base::out);
   eph.writeASCIIheader(ofs);
   ofs.close();

   LOG(INFO) << "Dump ASCII data to csse.data.asc";
   ofs.open("csse.data.asc",ios_base::out);
   eph.writeASCIIdata(ofs);
   ofs.close();

   // write the whole thing out to a binary file
   LOG(INFO) << "Write to binary file " << outputFilename;
   eph.writeBinaryFile(outputFilename);
   LOG(INFO) << "Finished writing binary file.";

   // read it back in
   LOG(INFO) << "Read from binary file " << outputFilename;
   eph.readBinaryFile(outputFilename);
   LOG(INFO) << "Finished reading binary file " << outputFilename;

   // dump to a file
   LOG(INFO) << "Dump ASCII header to csse.header.bin.asc";
   ofs.open("csse.header.bin.asc",ios_base::out);
   eph.writeASCIIheader(ofs);
   ofs.close();

   LOG(INFO) << "Dump ASCII data to csse.data.bin.asc";
   ofs.open("csse.data.bin.asc",ios_base::out);
   eph.writeASCIIdata(ofs);
   ofs.close();
   LOG(INFO) << "Now compare the outputs by differencing";
   LOG(INFO) << " Try 'diff csse.data.asc csse.data.bin.asc'";
   LOG(INFO) << " and 'diff csse.data.asc csse.data.bin.asc'";

      // compute run time
   totaltime = clock()-totaltime;
   LOG(INFO) << PrgmName << " timing: " << fixed << setprecision(9)
      << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds.";
   //if(LOGstrm != cout) cout << PrgmName << " timing: " << fixed << setprecision(9)
   //   << double(totaltime)/double(CLOCKS_PER_SEC) << " seconds." << endl;

   return iret;
}
catch(Exception& e) {
   LOG(ERROR) << "GPSTk Exception : " << e.what();
}
catch (...) {
   LOG(ERROR) << "Unknown error in " << PrgmName << ".  Abort." << endl;
}
   return -1;
}   // end main()

// -----------------------------------------------------------------------------------
// Compute the satellite attitude, given the time and the satellite position SV.
// If the SolarSystem is valid, use it; otherwise use SolarPosition.
// See two versions of SatelliteAttitude() for the user interface.
// Return a 3x3 Matrix which contains, as rows, the unit (ECEF) vectors X,Y,Z in the
// body frame of the satellite, namely
//    Z = along the boresight (i.e. towards Earth center),
//    Y = perpendicular to both Z and the satellite-sun direction, and
//    X = completing the orthonormal triad. X will generally point toward the sun.
// Thus this rotation matrix R * (ECEF XYZ vector) = components in body frame, and
// R.transpose() * (sat. body. frame vector) = ECEF XYZ components.
// Also return the shadow factor = fraction of sun's area not visible to satellite.
Matrix<double> doSatAtt(const CommonTime& tt, const Position& SV,
                        const SolarSystem& SSEph, const EarthOrientation& EO,
                        double& sf)
   throw(Exception)
{
   try {
      int i;
      double d,svrange,DistSun,AngRadSun,AngRadEarth,AngSeparation;
      Position X,Y,Z,T,S,Sun;
      Matrix<double> R(3,3);

      // Z points from satellite to Earth center - along the antenna boresight
      Z = SV;
      Z.transformTo(Position::Cartesian);
      svrange = Z.mag();
      d = -1.0/Z.mag();
      Z = d * Z;                                // reverse and normalize Z

      // get the Sun's position
      if(SSEph.JPLNumber() > -1) {
         //SolarSystem& mySSEph=const_cast<SolarSystem&>(SSEph);
         Sun = const_cast<SolarSystem&>(SSEph).WGS84Position(SolarSystem::Sun,tt,EO);
      }
      else {
         double AR;
         Sun = SolarPosition(tt, AR);
      }
      DistSun = Sun.radius();

      // apparent angular radius of sun = 0.2666/distance in AU (deg)
      AngRadSun = 0.2666/(DistSun/149598.0e6);
      AngRadSun *= DEG_TO_RAD;

      // angular radius of earth at sat
      AngRadEarth = ::asin(6378137.0/svrange);

      // T points from satellite to sun
      T = Sun;                                  // vector earth to sun
      T.transformTo(Position::Cartesian);
      S = SV;
      S.transformTo(Position::Cartesian);
      T = T - S;                                // sat to sun=(E to sun)-(E to sat)
      d = 1.0/T.mag();
      T = d * T;                                // normalize T

      AngSeparation = ::acos(Z.dot(T));         // apparent angular distance, earth
                                                // to sun, as seen at satellite
      // is satellite in eclipse?
      try { sf = ShadowFactor(AngRadEarth, AngRadSun, AngSeparation); }
      catch(Exception& e) { GPSTK_RETHROW(e); }

      // Y is perpendicular to Z and T, such that ...
      Y = Z.cross(T);
      d = 1.0/Y.mag();
      Y = d * Y;                                // normalize Y

      // ... X points generally in the direction of the sun
      X = Y.cross(Z);                           // X will be unit vector
      if(X.dot(T) < 0) {                        // need to reverse X, hence Y also
         X = -1.0 * X;
         Y = -1.0 * Y;
      }

      // fill the matrix and return it
      for(i=0; i<3; i++) {
         R(0,i) = X[i];
         R(1,i) = Y[i];
         R(2,i) = Z[i];
      }

      return R;
   }
   catch(Exception& e) { GPSTK_RETHROW(e); }
   catch(exception& e) {Exception E("std except: "+string(e.what())); GPSTK_THROW(E);}
   catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); }
}

//------------------------------------------------------------------------------------
int main(int argc, char **argv)
{
try {
   clock_t totaltime;            // for timing
   totaltime = clock();

   // locals
   int i,iret=0;
   DayTime CurrEpoch;         // defaults to current local time
   SolarSystem SSEphemeris;

   // program name, title and version
   PrgmName = string("testSSEph");
   Title = PrgmName + ", test program for JPL ephemeris, version "
      + Version + ", " + CurrEpoch.printf("Run %04Y/%02m/%02d at %02H:%02M:%02S");

   // default command line input
   bool verbose=false,debug=false;
   string inputFilename,testFilename,logFilename;

   // parse the command line input
   for(i=0; i<argc; i++) {
      string word = string(argv[i]);
      if(argc == 1 || word == "-h" || word == "--help") {
         cout << "Program " << PrgmName
            << " reads a binary JPL planetary ephemeris file, created by "
            << "convertSSEph,\n  and a test file, downloaded from the JPL ftp site, "
            << "containing times and planets\n  with JPL-generated ephemeris "
            << "coordinate values. The coordinates are computed using\n  the binary "
            << "file and the SolarSystem class, and compared with the JPL values;\n  "
            << "any difference larger than 10^-13 is noted with the word 'Failure' "
            << "at EOL.\n  Note that some large coordinate values may differ at the "
            << "level of 10^-13 because the\n  size of double precision is barely "
            << "able to hold that much precision; compare the\n  computed value "
            << "with the JPL value (copied as a string) in the output file.\n"
            << "\n Usage: " << PrgmName << " [options]\n Options are:\n"
            << "   --log <file>   name of optional log file (otherwise stderr)\n"
            << "   --file <file>  name of binary SS ephemeris file\n"
            << "   --test <file>  name of JPL test file (e.g. testpo.403)\n"
            << "   --verbose      print info to the log file.\n"
            << "   --debug        print debugging info to the log file.\n"
            << "   --help         print this and quit.\n"
            ;
         return 0;
      }
      else if(i == 0) continue;
      else if(word == "-d" || word == "--debug") debug = true;
      else if(word == "-v" || word == "--verbose") verbose = true;
      else if(word == "--log") logFilename = string(argv[++i]);
      else if(word == "--file") inputFilename = string(argv[++i]);
      else if(word == "--test") testFilename = string(argv[++i]);
   }

   // test input
   if(inputFilename.empty()) {
      LOG(ERROR) << "Must specify an input binary file name";
      return -1;
   }
   if(testFilename.empty()) {
      LOG(ERROR) << "Must specify an input test file name";
      return -1;
   }

   if(!logFilename.empty()) {
      // choose the log file
      logstrm.open(logFilename.c_str(),ios_base::out);
      ConfigureLOG::Stream() = &logstrm;
      // if not the above, output is to stderr

      cout << Title << endl;
      cout << "Output is logged to file " << logFilename << endl;
   }

      // set the maximum level to be logged
   ConfigureLOG::ReportLevels() = ConfigureLOG::ReportTimeTags() = true;
   if(debug)
      ConfigureLOG::ReportingLevel() = ConfigureLOG::FromString("DEBUG");
   else if(verbose)
      ConfigureLOG::ReportingLevel() = ConfigureLOG::FromString("VERBOSE");
   else
      ConfigureLOG::ReportingLevel() = ConfigureLOG::FromString("INFO");
      // = any of ERROR,WARNING,INFO,VERBOSE,DEBUG,DEBUGn (n=1,2,3,4,5,6,7)
   //cout << "Reporting in main is "
   //   << ConfigureLOG::ToString(ConfigureLOG::ReportingLevel()) << endl;

   // display title in the log file
   LOG(INFO) << Title;

   // now read the binary file, and read selected records
   // use the binary to test using the JPL file testpo.<EPH#>
   LOG(VERBOSE) << "Initialize with file " << inputFilename;
   SSEphemeris.initializeWithBinaryFile(inputFilename);
   LOG(VERBOSE) << "End Initialize";
   LOG(INFO) << "Ephemeris number is " << SSEphemeris.JPLNumber();

   bool foundEOT=false;
   int target,center,coord;
   double JD,PV[6],value,diff;
   SolarSystem::Planet Target,Center;
//.........这里部分代码省略.........