GyotoAstrobj.h

Go to the documentation of this file.

 
/*
    Copyright 2011-2020 Thibaut Paumard, Frederic Vincent
 
    This file is part of Gyoto.
 
    Gyoto is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    Gyoto is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with Gyoto.  If not, see <http://www.gnu.org/licenses/>.
 */
 
 
#ifndef __GyotoAstrobj_H_ 
#define __GyotoAstrobj_H_ 
 
#include "GyotoConfig.h"
 
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
 
#include <GyotoDefs.h>
#include <GyotoSmartPointer.h>
#include <GyotoConverters.h>
#include <GyotoObject.h>
#include <eigen3/Eigen/Dense>
 
namespace Gyoto{
  class Photon;
  namespace Register { class Entry; }
  namespace Metric { class Generic; }
  class FactoryMessenger;
  namespace Astrobj {
    class Generic;
    class Properties;
 
    typedef SmartPointer<Gyoto::Astrobj::Generic>
      Subcontractor_t(Gyoto::FactoryMessenger*, std::vector<std::string> const &);
 
    template<typename T> SmartPointer<Astrobj::Generic> Subcontractor
      (FactoryMessenger* fmp, std::vector<std::string> const &plugin) {
      SmartPointer<T> ao = new T();
      ao -> plugins(plugin) ;
#ifdef GYOTO_USE_XERCES
      if (fmp) ao -> setParameters(fmp);
#endif
      return ao;
    }
 
 
    Gyoto::Astrobj::Subcontractor_t* getSubcontractor(std::string name,
                                                      std::vector<std::string> &plugin,
                                                      int errmode = 0);
 
    extern Gyoto::Register::Entry * Register_;
 
    void initRegister(); 
 
    void Register(std::string name, Gyoto::Astrobj::Subcontractor_t* scp);
  }
}
 
class Gyoto::Astrobj::Generic
: public Gyoto::SmartPointee,
  public Gyoto::Object
{
  friend class Gyoto::SmartPointer<Gyoto::Astrobj::Generic>;
 
 
  // Data : 
  // -----
private:
 
  int __defaultfeatures;
 
 protected:
 
  SmartPointer<Gyoto::Metric::Generic> gg_;
 
 
  double rmax_; 
 
  double deltamaxinsidermax_; 
 
  bool flag_radtransf_; 
 
  int shadow_; 
  int noredshift_; 
  // Constructors - Destructor
  // -------------------------
 public:
  GYOTO_OBJECT;
 
  Generic(); 
 
  Generic(double radmax); 
 
  Generic(std::string kind); 
 
  Generic(const Generic& ) ; 
 
  virtual Generic* clone() const = 0 ; 
  
  virtual ~Generic() ; 
 
  // Accessors
  // ---------
 public:
  virtual SmartPointer<Metric::Generic> metric() const;
 
  virtual void metric(SmartPointer<Metric::Generic>) ;
 
  virtual double rMax(); 
  virtual double rMax() const; 
 
  virtual double rMax(std::string const &unit); 
  virtual double rMax(std::string const &unit) const; 
 
 
  virtual double deltaMax(double coord[8]);
 
  virtual void rMax(double val); 
 
  virtual void rMax(double val, std::string const &unit); 
 
  GYOTO_OBJECT_ACCESSORS_UNIT(deltaMaxInsideRMax);
 
  void opticallyThin(bool flag);
 
  bool opticallyThin() const ;
 
  void showshadow(bool flag);
  bool showshadow() const ;
 
  void redshift(bool flag);
  bool redshift() const ;
 
  virtual Gyoto::Quantity_t getDefaultQuantities();
 
  //XML I/O
 public:
 
#ifdef GYOTO_USE_XERCES
  virtual void setParameters(FactoryMessenger *fmp);
 
 
#endif
  
  // Outputs
  // -------
 public:
  virtual int Impact(Gyoto::Photon* ph, size_t index,
                     Astrobj::Properties *data=NULL) = 0 ;
  
  virtual void processHitQuantities(Photon * ph, state_t const &coord_ph_hit,
                                   double const * coord_obj_hit, double dt,
                                   Astrobj::Properties* data) const;
  virtual void processHitQuantities(Photon* ph, double* coord_ph_hit,
                                   double* coord_obj_hit, double dt,
                                   Astrobj::Properties* data) const = delete ;
 
  virtual double emission(double nu_em, double dsem, state_t const &coord_ph,
                          double const coord_obj[8]=NULL)
    const ;
  virtual double emission(double nu_em, double dsem, double coord_ph[8],
                          double coord_obj[8]=NULL)
    const = delete; 
 
  virtual void emission(double Inu[], double const nu_em[], size_t nbnu,
                        double dsem, state_t const &coord_ph,
                        double const coord_obj[8]=NULL) const ;
  virtual void emission(double Inu[], double nu_em[], size_t nbnu,
                        double dsem, double coord_ph[8],
                        double coord_obj[8]=NULL) const = delete ; 
 
  virtual void radiativeQ(double Inu[], double Taunu[], 
                          double const nu_em[], size_t nbnu,
                          double dsem, state_t const &coord_ph,
                          double const coord_obj[8]=NULL) const ;
  virtual void radiativeQ(double Inu[], double Taunu[], 
                          double nu_em[], size_t nbnu,
                          double dsem, double coord_ph[8],
                          double coord_obj[8]=NULL) const = delete ;
  virtual void radiativeQ(double *Inu, double *Qnu, double *Unu, double *Vnu,
                          Eigen::Matrix4d *Onu, double const *nuem , size_t nbnu, double dsem,
                          state_t const &cph, double const *co) const ;
 
  virtual double integrateEmission(double nu1, double nu2, double dsem,
                                   state_t const &c_ph, double const c_obj[8]=NULL) const;
  virtual double integrateEmission(double nu1, double nu2, double dsem,
                                   double c_ph[8], double c_obj[8]=NULL) const=delete; 
 
  virtual void integrateEmission(double * I, double const * boundaries,
                                 size_t const * chaninds, size_t nbnu,
                                 double dsem, state_t const &cph, double const *co) const;
  virtual void integrateEmission(double * I, double const * boundaries,
                                 size_t const * chaninds, size_t nbnu,
                                 double dsem, double *cph, double *co) const = delete;
 
  virtual double transmission(double nuem, double dsem, state_t const &coord_ph, double const coord_obj[8]) const ;
  virtual double transmission(double nuem, double dsem, state_t const &coord) const = delete;  
  virtual double transmission(double nuem, double dsem, double coord[8]) const = delete;
 
  Eigen::Matrix4d Omatrix(double alphanu[4], double rnu[3], double Chi, double dsem) const;
  Eigen::Matrix4d Omatrix(double alphaInu, double alphaQnu, double alphaUnu, double alphaVnu,
        double rQnu, double rUnu, double rVnu, double Chi, double dsem) const;
  Eigen::Matrix4d Omatrix(double alphanu[4], double rnu[3], double sin2Chi, double cos2Chi, double dsem) const;
  Eigen::Matrix4d Omatrix(double alphaInu, double alphaQnu, double alphaUnu, double alphaVnu,
        double rQnu, double rUnu, double rVnu, double sin2Chi, double cos2Chi, double dsem) const;
 
  Eigen::Matrix4d Pmatrix(double alphaInu, double alphaQnu, double alphaUnu, double alphaVnu,
        double rQnu, double rUnu, double rVnu, double sin2Chi, double cos2Chi, double dsem) const;
  
  Eigen::Vector4d rotateJs(double jInu, double jQnu, double jUnu, double jVnu, double sin2Chi, double cos2Chi) const;
  Eigen::Vector4d rotateJs(double jInu, double jQnu, double jUnu, double jVnu, double Chi) const;
 
  double getChi(double const fourvect[4], state_t const &cph, double const vel[4], bool elec=false) const;
  void getSinCos2Chi(double const fourvect[4], state_t const &cph, double const vel[4], double* sin2Chi, double* cos2Chi, bool elec=false) const;
  void computeB4vect(double B4vect[4], std::string const magneticConfig, double const co[8], state_t const &cph) const;
 
  void computeB4vect_ipole(double B4vect[4], std::string const magneticConfig, double const co[8], state_t const &cph, double spin) const;
  double interpolate(int const N, double* const array, double* const Xq, double** const X, int* const X_params, std::string const *cond_limits) const;
  double interpolate(int const N, double* const array, double* const Xq, double** const X_params, std::string const *cond_limits) const;
 
  private :
  double interp1d(double const x, double const y0, double const y1) const;
  double interpNd(int const N, double* const Xq, double** const X, double* const Y, std::string const *cond_limit) const;
  int getIndice(double &xq, std::string const cond_limit, double const X_params[3], double* const X=NULL) const;
 
};
 
class Gyoto::Astrobj::Properties : public Gyoto::SmartPointee {
  friend class Gyoto::SmartPointer<Gyoto::Astrobj::Properties>;
 public:
  double *intensity; 
  double *time; 
 
  double *distance; 
 
  double * first_dmin; 
 
  int first_dmin_found; 
 
  double *redshift; 
 
  double *nbcrosseqplane; 
 
  double *spectrum; 
 
  double *stokesQ;
  double *stokesU;
  double *stokesV;
 
  double *binspectrum; 
 
  ptrdiff_t offset; 
 
  double * impactcoords; 
 
  double *user1;
 
  double *user2;
 
  double *user3;
 
  double *user4;
 
  double *user5;
# ifdef HAVE_UDUNITS
  Gyoto::SmartPointer<Gyoto::Units::Converter> intensity_converter_ ;
  Gyoto::SmartPointer<Gyoto::Units::Converter> spectrum_converter_ ;
  Gyoto::SmartPointer<Gyoto::Units::Converter> binspectrum_converter_ ;
# endif
 
  bool alloc;
 
 public:
  Properties(); 
  Properties (double*, double*); 
 
  void init(size_t nbnuobs=0);
 
  Properties& operator++();
 
  Properties& operator+=(ptrdiff_t offset);
 
  operator Gyoto::Quantity_t () const;
 
# ifdef HAVE_UDUNITS
  void intensityConverter(Gyoto::SmartPointer<Gyoto::Units::Converter>);
  void intensityConverter(std::string);
  void spectrumConverter(Gyoto::SmartPointer<Gyoto::Units::Converter>);
  void spectrumConverter(std::string);
  void binSpectrumConverter(Gyoto::SmartPointer<Gyoto::Units::Converter>);
  void binSpectrumConverter(std::string);
# endif
};
 
#endif