ComplexRVMFit.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2009 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// psrchive/More/Polarimetry/Pulsar/ComplexRVMFit.h
 
#ifndef __Pulsar_ComplexRVMFit_h
#define __Pulsar_ComplexRVMFit_h
 
#include "MEAL/Axis.h"
#include "Estimate.h"
 
#include <complex>
 
// forward declarations
namespace MEAL
{
  class ComplexRVM;
  class OrthoRVM;
  class RotatingVectorModel;  
}
 
typedef std::pair<double,double> range;
 
namespace Pulsar
{
  class PolnProfile;

  class ComplexRVMFit : public Reference::Able
  {
    
  public:

    static bool verbose;

    ComplexRVMFit ();

    double get_gate_duty_cycle() const { return gate; }
    
    virtual void set_gate_duty_cycle (double turns) { gate = turns; }


    void add_opm (const range& radians);

    bool is_opm (double phase_in_radians) const;


    void add_include (const range& radians);

    bool is_included (double phase_in_radians) const;


    void add_exclude (const range& radians);

    bool is_excluded (double phase_in_radians) const;

    void set_observation (const PolnProfile*);
    const PolnProfile* get_observation () const;

    void set_threshold (float sigma);
    float get_threshold () const;
    
    void set_guess_smooth (unsigned phase_bins);
    unsigned get_guess_smooth () const;
    
    MEAL::ComplexRVM* get_model ();

    void solve ();
 
    void set_chisq_map (bool flag) { chisq_map = flag; }
 
    void set_range_alpha (const range& r) { range_alpha = r; }
    void set_range_zeta (const range& r) { range_zeta = r; }
    void set_range_beta (const range& r) { range_beta = r; }
 
    void search_2D (unsigned nalpha, unsigned nzeta);
    void search_1D (unsigned nzeta);

    double evaluate (double phi_radians);
 
    float get_chisq () const { return chisq; }
    unsigned get_nfree () const { return nfree; }
 
    // get the post-fit residuals of Stokes Q and U
    void get_residuals (std::vector<double>& phases,
                        std::vector< std::complex< Estimate<double> > >& res);
 
    // get the post-fit residuals of position angle
    void get_psi_residuals (std::vector<double>& phases,
                            std::vector< Estimate<double> >& res);
 
  protected:
 
    // fraction of pulse period stored in pulse profile
    double gate;
    
    void check_parameters ();
 
    void init (MEAL::RotatingVectorModel*);
    void init (MEAL::OrthoRVM*);
 
    void check_parameters (MEAL::RotatingVectorModel*);
    
    Reference::To<const PolnProfile> data;
    Reference::To<MEAL::ComplexRVM> model;
 
    float threshold;
 
    // ranges of pulse longitude to include in model
    std::vector<range> range_include;
    // ranges of pulse longitude to exclude from model
    std::vector<range> range_exclude;
 
    MEAL::Axis<unsigned> state;
    std::vector< MEAL::Axis<unsigned>::Value > data_x;
    std::vector< std::complex< Estimate<double> > > data_y;
 
    // orthogonally polarized mode phase ranges
    std::vector<range> opm;
 
    range range_alpha;
    range range_beta;
    range range_zeta;
 
    // length of the first guess ellipse along the sin(alpha) axis
    float guess_alpha;
    // length of the first guess ellipse along the sin(beta) axis
    float guess_beta;
 
    // number of phase bins in smoothing window used to stabilize first guess
    unsigned guess_smooth;
 
  private:
 
    // used by set_observation method to find the maximum in delpsi/delphi
    void find_delpsi_delphi_max ();
 
    // set L of each bin to least-squares best fit
    void renormalize();
 
    std::vector< std::complex< Estimate<double> > > linear;
    double delpsi_delphi;
 
    float peak_phase;
    float peak_pa;
 
    float chisq;
    unsigned nfree;
    bool chisq_map;
 
    double max_L;
 
    bool auto_detect_opm;
  };
}
 
#endif