polyco.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 1999 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// psrchive/Util/tempo/polyco.h
 
#ifndef __POLY_H
#define __POLY_H
 
#include "Pulsar/Predictor.h"
#include "Pulsar/Generator.h"
 
#include <string>
#include <vector>

class polynomial {
  
  friend class polyco;
  
protected:
  
  std::string psrname;
  
  std::string date;

  std::string utc;

  Pulsar::Phase ref_phase;

  double ref_freq;

  MJD ref_time;

  char telescope;

  double freq;

  bool binary;

  double binph;

  double binfreq;

  float dm;

  double nspan_mins;

  bool tempov11;

  double doppler_shift;

  double log_rms_resid;

  std::vector<double> coefs;

  void init();

  mutable bool keep;
 
public:

  class Expert;
  friend class Expert;

  static double precision;

  polynomial ();

  polynomial (const polynomial& poly);

  ~polynomial ();
 
  const polynomial& operator = (const polynomial& poly);
 
  int load (std::string* instr);
  int unload (std::string *outstr) const;
  int unload (FILE* fptr) const;
 
  void prettyprint  () const;

  double period (const MJD &t) const;

  Pulsar::Phase phase (const MJD &t) const;

  MJD iphase (const Pulsar::Phase& p, const MJD* = 0) const;

  long double frequency (const MJD &t) const;

  Pulsar::Phase dispersion (const MJD &t, long double MHz) const;

  double chirp (const MJD &t) const;

  double accel (const MJD &t) const { return chirp(t)/frequency(t)*2.9979e8; }
 
  bool   is_tempov11       () const { return tempov11; }
  char   get_telescope     () const { return telescope; }
  double get_freq          () const { return freq; }
  MJD    get_reftime       () const { return ref_time; }
  Pulsar::Phase  get_refphase      () const { return ref_phase; }
  double get_reffrequency  () const { return ref_freq; }
  double get_nspan         () const { return nspan_mins; }
  float  get_dm            () const { return dm; }
  int    get_ncoeff        () const { return (int) coefs.size(); }
  double get_doppler_shift () const { return doppler_shift / 1e4; }
  std::string get_psrname  () const { return psrname; }
  bool   get_binary        () const { return binary; } 
  double get_binph         () const { return binph; }
  double get_log_rms_resid () const { return log_rms_resid; }
 
  static double flexibility;
 
  MJD start_time( double f = flexibility ) const
  { return ref_time - nspan_mins * (1.0+f) * 60.0/2.0; };
  MJD end_time( double f = flexibility ) const 
  { return ref_time + nspan_mins * (1.0+f) * 60.0/2.0; };
 
  Pulsar::Phase start_phase( double f = flexibility ) const
  { return phase (start_time(f)); };
  Pulsar::Phase end_phase( double f = flexibility ) const
  { return phase (end_time(f)); };
  
  friend bool operator == (const polynomial &, const polynomial &);
  friend bool operator != (const polynomial &, const polynomial &);
};

class polynomial::Expert {
 
public:
  Expert (polynomial* p) { instance = p; }
 
  void set_tempov11      (bool v) { instance->tempov11 = v; }
  void set_telescope     (char v) { instance->telescope = v; }
  void set_freq          (double v) { instance->freq = v; }
  void set_reftime       (const MJD& v) { instance->ref_time = v; }
  void set_refphase      (const Pulsar::Phase& v) { instance->ref_phase = v; }
  void set_reffrequency  (double v) { instance->ref_freq = v; }
  void set_nspan         (double v) { instance->nspan_mins = v; }
  void set_dm            (float v) { instance->dm = v; }
  void set_ncoef         (unsigned v) { instance->coefs.resize(v); }
  void set_doppler_shift (double v) { instance->doppler_shift = v * 1e4; }
  void set_psrname       (const std::string& v) { instance->psrname = v; }
  void set_binary        (bool v) { instance->binary = v; }
  void set_binph         (double v) { instance->binph = v; }
 
  double* get_coefs () { return &(instance->coefs[0]); }
 
  polynomial* instance;
};

class polyco : public Pulsar::Predictor {
 
public:
  static bool debug;

  std::vector<polynomial> pollys;

  polyco ();

  polyco (const polyco&);

  const polyco& operator = (const polyco&);

  ~polyco ();
 
  //
  // Pulsar::Predictor implementation
  //

  Predictor* clone () const;

  Pulsar::Generator* generator () const;

  void insert (const Predictor*);

  void keep (const std::vector<MJD>& epochs);

  bool matches (const Predictor*) const;

  void set_observing_frequency (long double MHz);

  long double get_observing_frequency () const;

  Pulsar::Phase phase (const MJD& t) const
  { return best(t).phase(t); }

  Pulsar::Phase phase (const MJD& t, long double MHz) const
  { const polynomial& b = best(t); return b.phase(t) + b.dispersion(t,MHz); }

  MJD iphase (const Pulsar::Phase& phase, const MJD* guess = 0) const
  { return best(phase).iphase(phase, guess); }

  long double frequency (const MJD& t) const
  { return best(t).frequency(t); }

  Pulsar::Phase dispersion (const MJD &t, long double MHz) const
  { return best(t).dispersion(t,MHz); }

  void load (FILE*);

  void unload (FILE*) const;
 
  //
  // the rest
  //

  polyco (const std::string& id);

  int load (const std::string& filename);
  int load (std::string* instr);
 
  // these functions return -1 upon error
  int unload (const std::string& filename) const;
 
  // these functions return the number of bytes unloaded (-1 on error)
  int unload (std::string *outstr) const;
 
  void append (const polyco& poly);
 
  void  prettyprint () const;
 
  const polynomial* nearest (const MJD &t) const;
 
  const polynomial& best (const MJD &t) const;
  const polynomial& best (const Pulsar::Phase &p) const;
 
  int i_nearest (const MJD &t, bool throw_exception = false) const;
  int i_nearest (const Pulsar::Phase &p, bool throw_exception = false) const;
 
  double doppler_shift (const MJD& t) const
  { return best(t).get_doppler_shift(); };
 
  double period(const MJD& t) const
  { return best(t).period(t); };
 
  double chirp(const MJD& t) const
  { return best(t).chirp(t); };
 
  double accel(const MJD& t) const
  { return best(t).accel(t); };
 
  char   get_telescope () const;
  double get_freq      () const;
  MJD    get_reftime   () const;
  double get_refperiod () const;
  double get_nspan     () const;
  float  get_dm        () const;
  int    get_ncoeff    () const;
  std::string get_psrname   () const;
 
  bool is_tempov11() const;
 
  MJD  start_time () const { return pollys.front().start_time(); };
  MJD  end_time ()   const { return pollys.back().end_time(); };
 
  Pulsar::Phase start_phase () const { return pollys.front().start_phase(); };
  Pulsar::Phase end_phase ()   const { return pollys.back().end_phase(); };
 
  friend bool operator == (const polyco &, const polyco &);
  friend bool operator != (const polyco &, const polyco &);
 
private:
 
  // these attributes enable optimized i_nearest
  mutable int last_index;
 
  mutable MJD last_epoch;
  mutable double last_span_epoch;
 
  mutable Pulsar::Phase last_phase;
  mutable double last_span_phase;
 
  void init();
  void set_last (int i) const;
};
 
inline std::ostream& operator<< (std::ostream& ostr, const polyco& p)
{
  std::string out; p.unload(&out); return ostr << out;
}
 
#endif