C++ MJD类代码示例

void dsp::TimeDivide::set_start_time (MJD _start_time)
{
  if (start_time == _start_time)
    return;

  if (Operation::verbose)
    cerr << "dsp::TimeDivide::set_start_time start_time=" 
         << _start_time.printall() << endl;

  start_time = _start_time;
  start_phase = Pulsar::Phase::zero;
  is_valid = false;

  if( reference_epoch != MJD::zero )
  {
    if (Operation::verbose)
      cerr << "dsp::TimeDivide::set_start_time set to reference_epoch=" 
           << reference_epoch.printall() << endl;

    start_time = reference_epoch;
  }
  else if( division_seconds && division_seconds == unsigned(division_seconds) )
  {
    unsigned integer_seconds = unsigned(division_seconds);
    unsigned seconds = start_time.get_secs();
    unsigned divisions = seconds / integer_seconds;
    start_time = MJD (start_time.intday(), divisions * integer_seconds, 0.0);

    if (Operation::verbose)
      cerr << "dsp::TimeDivide::set_start_time rounded start_time=" 
           << _start_time.printall() << endl;
  }
}

//! Get the mid-time of the integration
MJD dsp::PhaseSeries::get_mid_time (bool phased) const
{
  if (verbose)
  {
    cerr << "PhaseSeries::get_mid_time start=" 
      << start_time << " end=" << end_time << endl;
  }

  MJD midtime = 0.5 * (start_time + end_time);

  if (!phased)
    return midtime;

  if (folding_predictor)
  {
    // truncate midtime to the nearest pulse phase = reference_phase
    Pulsar::Phase phase = folding_predictor->phase(midtime).Floor() + reference_phase;
    midtime = folding_predictor->iphase (phase, &midtime);
  }

  if (folding_period)
  {
    double phase = reference_phase + 
      fmod (midtime.in_seconds(), folding_period)/folding_period;
    midtime -= phase * folding_period;
  }

  return midtime;
}

//! convert an MJD to long double
long double from_MJD (const MJD& t)
{
  const long double secs_in_day = 86400.0L;

  return 
    (long double) (t.intday()) +
    (long double) (t.get_secs()) / secs_in_day +
    (long double) (t.get_fracsec()) / secs_in_day;
}

void qt_MJD::value_Entered_CB ()
{
  MJD newval;
  if (newval.Construct (value.text().ascii()) < 0) {
    if (MJD::verbose)
      std::cerr << "qt_MJD:: invalid mjd:" << value.text() << std::endl;
    newval = valset;
  }
  setMJD (newval);
}

//! Return the phase, given the epoch
Phase Tempo2::Predictor::phase (const MJD& t) const
{
  if (verbose)
    cerr << "Tempo2::Predictor::phase epoch=" << t << " frequency=" 
	 << observing_frequency << endl;

  if (observing_frequency <= 0)
    Error error (InvalidState, "Tempo2::Predictor::phase",
		 "observing_frequency=%lf", (double) observing_frequency);

  long double p = T2Predictor_GetPhase ( &predictor, from_MJD (t),
					 observing_frequency );

  if (ChebyModelSet_OutOfRange)
    throw Error (InvalidParam, "Tempo2::Predictor::phase",
		 "epoch %s not spanned by ChebyModelSet",
		 t.printdays(20).c_str());

  if (!finite(p)) {
    Error error (InvalidState, "Tempo2::Predictor::phase",
		 "T2Predictor_GetPhase result = ");
    error << p;
    throw error;
  }

  return to_Phase( p );
}

//! Return the spin frequency, given the epoch
long double Pulsar::SimplePredictor::frequency (const MJD& t) const
{
  if (reference_epoch == MJD::zero)
     const_cast<SimplePredictor*>(this)->reference_epoch = t;

  MJD diff = t - reference_epoch;
  long double seconds = diff.get_secs();
  seconds += diff.get_fracsec();

  long double power_of_t = 1.0;
  long double result = 0;

  for (unsigned i=0; i<coefs.size(); i++)
  {
    result += coefs[i] * (i+1) * power_of_t;
    power_of_t *= seconds;
  }

  return result;
}

//! Return the spin frequency, given the epoch
long double Tempo2::Predictor::frequency (const MJD& t) const
{
  long double f = T2Predictor_GetFrequency (&predictor, from_MJD (t),
					    observing_frequency);

  if (ChebyModelSet_OutOfRange)
    throw Error (InvalidParam, "Tempo2::Predictor::frequency",
		 "epoch %s not spanned by ChebyModelSet",
		 t.printdays(20).c_str());

  return f;
}

//! Return the phase, given the epoch
Pulsar::Phase Pulsar::SimplePredictor::phase (const MJD& t) const
{
  if (reference_epoch == MJD::zero)
    const_cast<SimplePredictor*>(this)->reference_epoch = t;

  MJD diff = t - reference_epoch;
  long double seconds = diff.get_secs();
  seconds += diff.get_fracsec();

  long double power_of_t = seconds;
  long double result = 0;

  for (unsigned i=0; i<coefs.size(); i++)
  {
    result += coefs[i] * power_of_t;
    power_of_t *= seconds;
  }

  int64_t iphase = int64_t(result);

  return Pulsar::Phase(iphase,result-iphase);
}

MJD dsp::LoadToFold::parse_epoch (const std::string& epoch_string)
{
  MJD epoch;

  if (epoch_string == "start")
  {
    epoch = manager->get_info()->get_start_time();
    epoch += manager->get_input()->tell_seconds();

    if (Operation::verbose)
      cerr << "dsp::LoadToFold::parse reference epoch=start_time=" 
	   << epoch.printdays(13) << endl;
  }
  else if (!epoch_string.empty())
  {
    epoch = MJD( epoch_string );
    if (Operation::verbose)
      cerr << "dsp::LoadToFold::parse reference epoch="
	   << epoch.printdays(13) << endl;
  }

  return epoch;
}

//! Seek to a sample close to the specified MJD
void dsp::Input::seek(MJD mjd)
{
  int misplacement = 0;

  if (mjd+1.0/get_info()->get_rate() < get_info()->get_start_time())
    misplacement = -1;

  if (mjd-1.0/get_info()->get_rate() > get_info()->get_end_time())
    misplacement = 1;

  double seek_seconds = (mjd-get_info()->get_start_time()).in_seconds();

  if (misplacement)
  {
    string msg = "The given MJD (" + mjd.printall() + ") is ";
    if (misplacement < 0)
      msg += "before the start time";
    else
      msg += "after the end time";
    msg += "of the input data "
           "(" + get_info()->get_start_time().printall() + "); "
           "difference is %lf seconds";

    throw Error (InvalidParam, "dsp::Input::seek", msg.c_str(), seek_seconds);
  }
 
  double seek_samples = seek_seconds*get_info()->get_rate();
  uint64_t actual_seek = 0;

  if( seek_samples<0.0 )
    actual_seek = 0;
  else if( uint64_t(seek_samples) > get_info()->get_ndat() )
    actual_seek = get_info()->get_ndat();
  else
    actual_seek = uint64_t(seek_samples);

  if( verbose )
    fprintf(stderr,"dsp::Input::seek(MJD) will seek %f = "UI64" samples\n",
	    seek_samples, actual_seek);

  seek( actual_seek, SEEK_SET );
}

std::string dsp::FilenameEpoch::get_filename (const PhaseSeries* data)
{
  MJD epoch = data->get_start_time();

  if (Observation::verbose)
    cerr << "dsp::FilenameEpoch::get_filename epoch=" 
         << epoch.printall() << endl;

  if (integer_seconds)
  {
    // ensure that the epoch is rounded up into the current division
    epoch = data->get_mid_time (false);

    DEBUG("dsp::FilenameEpoch::get_filename mid_time=" << epoch.printall());

    unsigned seconds = epoch.get_secs();
    unsigned divisions = seconds / integer_seconds;
    epoch = MJD (epoch.intday(), divisions * integer_seconds, 0.0);

    if (Observation::verbose)
      cerr << "dsp::FilenameEpoch::get_filename division_start=" 
           << epoch.printall() << endl;
  }

  vector<char> fname (FILENAME_MAX);
  char* filename = &fname[0];

  if (!epoch.datestr( filename, FILENAME_MAX, datestr_pattern.c_str() ))
    throw Error (FailedSys, "dsp::PhaseSeriesUnloader::get_filename",
       "error MJD::datestr(" + datestr_pattern + ")");

  if (report_unload)
    cerr << "unloading " << tostring(data->get_integration_length(),2)
	 << " seconds: " << filename << endl;

  return filename;
}

double double_cast (const MJD& mjd)
{
  return mjd.in_seconds();
}

//.........这里部分代码省略.........
  case 3: /* timing mode data - not relevant, but needs to work! */
    break;
  case 4:
    get_info()->set_nbit (8);
    break;
  default:
    throw Error (InvalidState, "dsp::WAPPFile::open_file",
		 "lagformat variable in header should be 0, 1 or 4");
    break;
  }

  cerr << "WAPP nbit=" << get_info()->get_nbit() << endl;

  // ////////////////////////////////////////////////////////////////////
  //
  // start_time
  //

  /* built by WAPP from yyyymmdd */  
  string utc = head->obs_date;  utc += "-";

  /* UT seconds after midnight (start on 1-sec tick) [hh:mm:ss] */
  utc += head->start_time;

#if WVS_FIXES_STR2TM

  struct tm time;

  /* the str2tm function has been deprecated in favour of the standard C strptime */
  if (str2tm (&time, utc.c_str()) < 0)
    throw Error (InvalidState, "dsp::WAPPFile::open_file",
		 "Could not parse UTC from " + utc);

  MJD mjd (time);

#else

  // copied from WAPPArchive

  struct tm obs_date_greg;
  struct WAPP_HEADER* hdr = head;

  int rv = sscanf(hdr->obs_date, "%4d%2d%2d", 
      &obs_date_greg.tm_year, &obs_date_greg.tm_mon,
      &obs_date_greg.tm_mday);
  obs_date_greg.tm_year -= 1900;
  obs_date_greg.tm_mon -= 1;
  if (rv!=3) 
    throw Error (InvalidState, "dsp::WAPPFile::open_file",
        "Error converting obs_date string (rv=%d, obs_date=%s)", 
        rv, hdr->obs_date);
  rv = sscanf(hdr->start_time, "%2d:%2d:%2d", 
      &obs_date_greg.tm_hour, &obs_date_greg.tm_min, 
      &obs_date_greg.tm_sec);
  if (rv!=3) 
    throw Error (InvalidState, "dsp::WAPPFile::open_file",
        "Error converting start_time string (rv=%d, start_time=%s)", 
        rv, hdr->start_time);

  MJD mjd (obs_date_greg); 

#endif

  char buff[64];
  cerr << "UTC=" << utc << " MJD=" << mjd << " -> "
       << mjd.datestr (buff, 64, "%Y-%m-%d %H:%M:%S") << endl;

/*! This function assumes that the Integration will have the global
  attributes of the file. */
Pulsar::Integration* 
Pulsar::FITSArchive::load_Integration (const char* filename, unsigned isubint)
try {

  if (!filename)
    throw Error (InvalidParam, "FITSArchive::load_Integration",
		 "filename unspecified");

  if (search_mode)
    throw Error (InvalidParam, "FITSArchive::load_Integration",
		 "SEARCH mode data -- no Integrations to load");

  Reference::To<Pulsar::Integration> integ = new_Integration();
  init_Integration (integ);

  int row = isubint + 1;
  
  int status = 0;  

  if (verbose > 2)
    cerr << "FITSArchive::load_Integration number " << isubint << endl;
  
  double nulldouble = 0.0;
  float nullfloat = 0.0;
  
  int initflag = 0;
  int colnum = 0;
  
  // Open the file

  fitsfile* fptr = 0;
  
  if (verbose > 2)
    cerr << "FITSArchive::load_Integration"
      " fits_open_file (" << filename << ")" << endl;
  
  fits_open_file (&fptr, filename, READONLY, &status);
  
  if (status != 0)
    throw FITSError (status, "FITSArchive::load_Integration", 
		     "fits_open_file(%s)", filename);
  

  // Move to the SUBINT Header Data Unit
  
  fits_movnam_hdu (fptr, BINARY_TBL, "SUBINT", 0, &status);
  if (status != 0)
    throw FITSError (status, "FITSArchive::load_Integration", 
		     "fits_movnam_hdu SUBINT");


  // Load the convention for the epoch definition
  string epoch_def;
  string default_def = "STT_MJD";
  psrfits_read_key (fptr, "EPOCHS", &epoch_def,
		    default_def, verbose > 2);

  if (verbose > 2)
    cerr << "FITSArchive::load_Integration epochs are " << epoch_def << endl;

  // By default, correct epochs using the phase model
  bool correct_epoch_phase = true;

  // By default, correct epochs such that phase(epoch)=phase(start)
  bool phase_match_start_time = true;

  if (epoch_def == "VALID")
    correct_epoch_phase = phase_match_start_time = false;

  if (epoch_def == "MIDTIME")
    phase_match_start_time = false;

  if (get<Pulsar::IntegrationOrder>())
  {
    colnum = 0;
    fits_get_colnum (fptr, CASEINSEN, "INDEXVAL", &colnum, &status);
    
    double value = 0.0;

    fits_read_col (fptr, TDOUBLE, colnum, row, 1, 1, &nulldouble,
		   &value, &initflag, &status);
    
    if (status != 0)
      throw FITSError (status, "FITSArchive::load_Integration", 
		       "fits_read_col INDEXVAL");
    
    get<Pulsar::IntegrationOrder>()->set_Index(row-1,value);
  }
  
  // Get the reference epoch from the primary header
  const Pulsar::FITSHdrExtension* hdr_ext = get<Pulsar::FITSHdrExtension>();
  
  if (!hdr_ext)
  {
    throw Error (InvalidParam, "FITSArchive::load_Integration",
		 "No FITSHdrExtension found");
  }
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
               {
                  convertSODtoTime( asDouble( info['s'] ), 
                                    tt.hour, tt.minute, tt.second );
               }
               t = tt.convertToCommonTime();
               return;
            }
            else  // use YDSTime as default
            {
               YDSTime tt;
               tt.setFromInfo( info );
               if( hhour && hmin && hsec )
               {
                  tt.sod = convertTimeToSOD( asInt( info['H'] ), 
                                             asInt( info['M'] ), 
                                             asDouble( info['S'] ) );
               }
               t = tt.convertToCommonTime();
               return;
            }

         } // end of if( hyear )

         if( hzcount32 ||
             (hfullweek && hzcount) ||
             (hepoch && (hzcount29 || 
                         (hweek && hzcount))) )
         {
            GPSWeekZcount tt;
            tt.setFromInfo( info );
            t = tt.convertToCommonTime();
            return;
         }

         if ( (hepoch && hweek) || hfullweek )
         {
            WeekSecond* ptt;
            if(hbdse || hbdsfw || hbdsw) ptt = new BDSWeekSecond();
            else if(hqzse || hqzsfw || hqzsw) ptt = new QZSWeekSecond();
            else if(hgale || hgalfw || hgalw) ptt = new GALWeekSecond();
            else ptt = new GPSWeekSecond();
            ptt->setFromInfo(info);
            if( hdow && !hsow )
            {
               ptt->sow = asInt( info['w'] ) * SEC_PER_DAY;
               if( hsod )
               {
                  ptt->sow += asDouble( info['s'] );
               }
               else if( hhour && hmin && hsec )
               {
                  ptt->sow += convertTimeToSOD( asInt( info['H'] ), 
                                              asInt( info['M'] ), 
                                              asDouble( info['S'] ) );
               }
            }
            t = ptt->convertToCommonTime();
            return;
         }

         if( hmjd )
         {
            MJD tt;
            tt.setFromInfo( info );
            t = tt.convertToCommonTime();
            return;
         }

         if( hjulian )
         {
            JulianDate tt;
            tt.setFromInfo( info );
            t = tt.convertToCommonTime();
            return;
         }

         if( hansi )
         {
            ANSITime tt;
            tt.setFromInfo( info );
            t = tt.convertToCommonTime();
            return;
         } 
         
         if( hunixsec || hunixusec )
         {
            UnixTime tt;
            tt.setFromInfo( info );
            t = tt.convertToCommonTime();
            return;
         }

         InvalidRequest ir("Incomplete time specification for readTime");
         GPSTK_THROW( ir );
      }
      catch( gpstk::StringUtils::StringException& se )
      {
         GPSTK_RETHROW( se );
      }
   }