GyotoNumericalMetricLorene.h

Go to the documentation of this file.

/*
 *   Copyright (c) 2014-2018, 2020 Frederic Vincent, Thibaut Paumard
 *
 *
 */

#ifndef __GyotoNumericalMetricLoreneMetric_H_
#define __GyotoNumericalMetricLoreneMetric_H_ 

#include <iostream>
#include <fstream>

namespace Gyoto {
  namespace Metric { class NumericalMetricLorene; }
  class FactoryMessenger;
}

// Forward declarations of Lorene classes
namespace Lorene {
  class Scalar;
  class Vector;
  class Sym_tensor;
  class Valeur;
}

#include <GyotoMetric.h>
#include <GyotoWorldline.h>
#include <GyotoSmartPointer.h>

#ifdef GYOTO_USE_XERCES
#include <GyotoRegister.h>
#endif

class Gyoto::Metric::NumericalMetricLorene
: public Gyoto::Metric::Generic
{
  friend class Gyoto::SmartPointer<Gyoto::Metric::NumericalMetricLorene>;

 private:
  char* filename_; 
  bool mapet_; 
  bool axisymCirc_; 
  bool bosonstarcircular_; 
  int has_surface_; 
  int has_acceleration_vector_; 
  int specify_marginalorbits_; 
  double horizon_; 
  double initial_time_; 
  Lorene::Scalar** lapse_tab_;
  Lorene::Vector** shift_tab_;
  Lorene::Sym_tensor** gamcov_tab_;
  Lorene::Sym_tensor** gamcon_tab_;
  Lorene::Sym_tensor** kij_tab_;
  double* times_; 
  int nb_times_; 
  Lorene::Valeur** nssurf_tab_; 
  Lorene::Vector** vsurf_tab_; 
  Lorene::Vector** accel_tab_; 
  Lorene::Scalar** lorentz_tab_; 
  Lorene::Valeur** hor_tab_; 
  double risco_; 
  double rico_; 
  double rmb_; 

  void free(); 

 public:
  GYOTO_OBJECT;
  GYOTO_OBJECT_THREAD_SAFETY;
  NumericalMetricLorene(); 
  NumericalMetricLorene(const NumericalMetricLorene&); 
  virtual NumericalMetricLorene* clone() const ;
  virtual ~NumericalMetricLorene() ;        

  virtual void setMetricSource();

  void        directory(std::string const &dir) ;
  std::string directory() const ;
  double initialTime() const ;
  void   initialTime(double t0);
  double horizon() const ;
  void   horizon(double t0);
  double rico() const ;
  void   rico(double r0);
  bool hasSurface() const;
  void hasSurface(bool s);
  bool hasAccelerationVector() const;
  void hasAccelerationVector(bool aa);
  bool bosonstarcircular() const;
  void bosonstarcircular(bool);
  bool specifyMarginalOrbits() const;
  void specifyMarginalOrbits(bool s);
  bool mapEt() const;
  void mapEt(bool s);
  bool axisymCirc() const;
  void axisymCirc(bool s);

  Lorene::Vector** getShift_tab() const;
  Lorene::Scalar** getLapse_tab() const;
  Lorene::Sym_tensor** getGamcon_tab() const;
  Lorene::Sym_tensor** getGamcov_tab() const;
  double* getTimes() const;
  int getNbtimes() const;
  Lorene::Valeur** getNssurf_tab() const;
  Lorene::Vector** getVsurf_tab() const;
  Lorene::Vector** getAccel_tab() const;
  Lorene::Scalar** getLorentz_tab() const;
  Lorene::Valeur** getHor_tab() const;
  double getRms() const;
  double getRmb() const;
  void setLapse_tab(Lorene::Scalar* lapse, int ii);
  void setShift_tab(Lorene::Vector* shift, int ii);
  void setGamcov_tab(Lorene::Sym_tensor* gamcov, int ii);
  void setGamcon_tab(Lorene::Sym_tensor* gamcon, int ii);
  void setKij_tab(Lorene::Sym_tensor* kij, int ii);
  void setTimes(double time,int ii);

  virtual double getSpecificAngularMomentum(double rr) const;
  virtual double getPotential(double const pos[4], double l_cst) const;
  
  virtual void computeNBeta(const double coord[4],double &NN,double beta[3]) const;//Compute lapse and shift at coord

  using Generic::gmunu;
  double gmunu(const double x[4], int mu, int nu) const ;

  double gmunu(const double x[3], int indice_time, int mu, int nu) const ;

  virtual void gmunu_up(double ARGOUT_ARRAY2[4][4], const double IN_ARRAY1[4]) const ;

  void gmunu_up(double gup[4][4], const double x[4], int indice_time) const ;

  void gmunu_di(const double pos[4],
                double gmunudr[4][4],
                double gmunudth[4][4]) const ;
  
  void gmunu_di(const double pos[4],
                int indice_time,
                double gmunudr[4][4],
                double gmunudth[4][4]) const ;

  virtual void jacobian(double ARGOUT_ARRAY3[4][4][4], const double IN_ARRAY1[4]) const ;

  double gmunu_up_dr(const double x[4], int mu, int nu) const ;

  double gmunu_up_dr(const double x[3], int indice_time, int mu, int nu) const ;
  
  double christoffel(const double coord[4], const int alpha, const int mu,
                     const int nu) const ;
  double christoffel(const double coord[4],
                     const int alpha, 
                     const int mu, const int nu,
                     const int indice_time) const;
  virtual int christoffel(double dst[4][4][4],
                          const double coord[4]) const;
  int christoffel(double dst[4][4][4],
                  const double coord[4],
                  const int indice_time) const;


  double Interpol3rdOrder(double tt, int indice_time, double values[4]) const;

  double computeHorizon(const double* pos) const;
  double computeHorizon(const double* pos, int indice) const;

  //using Generic::diff;
  virtual int diff(state_t const &coord, state_t &res, double mass) const;
  int diff(double tt, const double y[7], double res[7]) const ;
  virtual int diff(const double y[7], double res[7], int indice_time) const ;
  virtual int diff31(state_t const &x, state_t &dxdt, double mass) const ;

  void circularVelocity(double const coor[4], double vel[3],
                        double dir) const ;
  void circularVelocity(double const coor[4], double vel[3],
                        double dir, int indice_time) const ;

};

#endif