inverse_phase.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2007 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// psrchive/Util/tempo/inverse_phase.h
 
#ifndef __P_inverse_phase_h
#define __P_inverse_phase_h
 
#include "true_math.h"
#include "MJD.h"
#include "Phase.h"
 
#include <iostream>
 
namespace Pulsar {

  extern unsigned inverse_phase_calls;
  extern unsigned inverse_phase_iterations;


  template<typename P>
  MJD inverse_phase (const P& predictor, const Phase& p, MJD guess)
  {
 
    MJD dt;
    int gi = 0;
 
    double precision = std::max (P::precision, MJD::precision);
    if (Pulsar::Predictor::verbose)
      std::cerr << "inverse_phase: precision=" << precision*1e6 
                << " us" << std::endl;
 
    inverse_phase_calls ++;
 
    for (gi=0; gi<10000; gi++) {
 
      inverse_phase_iterations ++;
 
      dt = (predictor.phase(guess) - p) / predictor.frequency(guess);
 
      if (!::true_math::finite(dt.in_seconds()))
        throw Error (InvalidState, "inverse_phase",
                     "dt not finite; freq=%lf", predictor.frequency(guess));
 
      guess -= dt; // * converge_faster;
 
      if (Pulsar::Predictor::verbose)
        std::cerr << "inverse_phase: guess=" << guess.printdays(20)
                  << " dt=" << dt.in_seconds()*1e6 << " us" << std::endl;
 
      if (fabs (dt.in_seconds()) < precision)
        return guess;
    }
    
    std::cerr << "inverse_phase maximum iterations exceeded - error="
              << dt.in_seconds() * 1e6 << "us" << std::endl;
    
    return guess;
  }


 
  template<typename P>
  MJD inverse_phase (const P& predictor, const Phase& p)
  {
    MJD guess = predictor.get_reftime()
      + (p - predictor.get_refphase()) / predictor.get_reffrequency();
 
    if (Pulsar::Predictor::verbose)
      std::cerr << "inverse_phase: reftime=" << predictor.get_reftime() 
                << " refphase=" << predictor.get_refphase() 
                << " reffreq=" << predictor.get_reffrequency() 
                << " first guess = " << guess << std::endl;
 
    return inverse_phase (predictor, p, guess);
  }
 
 
 
}
 
#endif