GyotoMetric.h

Go to the documentation of this file.

/*
    Copyright 2011-2016, 2018-2020 Frederic Vincent, Thibaut Paumard

    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 __GyotoMetric_H_
#define __GyotoMetric_H_ 

#include <iostream>
#include <fstream>
#include <string>
#include <vector>

#include <GyotoSmartPointer.h>
#include <GyotoObject.h>
#include <GyotoAstrobj.h>
#include <GyotoRegister.h>
#include <GyotoHooks.h>
#include <GyotoDefs.h>

namespace Gyoto {
  namespace Metric {

    class Generic;


    typedef SmartPointer<Metric::Generic> Subcontractor_t(FactoryMessenger*, std::vector<std::string> const &);


    template<typename T> SmartPointer<Metric::Generic> Subcontractor
      (FactoryMessenger* fmp, std::vector<std::string> const &plugins) {
      SmartPointer<T> gg = new T();
      gg -> plugins(plugins);
#ifdef GYOTO_USE_XERCES
      if (fmp) gg -> setParameters(fmp);
#endif
      return gg;
    }


    Gyoto::Metric::Subcontractor_t* getSubcontractor(std::string name,
                                                     std::vector<std::string> &plugin,
                                                     int errmode=0);


    extern Register::Entry * Register_;


     void Register(std::string kind, Gyoto::Metric::Subcontractor_t* scp);


     void initRegister();

  }

  /* Documented elswhere */
  class Worldline;
}

class Gyoto::Metric::Generic
: public Gyoto::SmartPointee,
  public Gyoto::Object,
  public Gyoto::Hook::Teller
{
  friend class Gyoto::SmartPointer<Gyoto::Metric::Generic>;

 private:
  double mass_;     
  int coordkind_; 

  int __defaultfeatures;

 protected:
  double delta_min_; 
  double delta_max_; 

  double delta_max_over_r_;

  bool keplerian_; 

 protected:
  void kind(const std::string); 

  void coordKind(int coordkind); 


 public:
  GYOTO_OBJECT;

  const std::string kind() const; 
  int getRefCount();
  
  // Constructors - Destructor
  // -------------------------
  Generic(const int coordkind, const std::string &name); 
  Generic(Generic const &o); 
  virtual ~Generic() ;                        
  
  // Mutators / assignment
  // ---------------------
  virtual Generic * clone() const ; 

  virtual void mass(const double);        
  virtual void mass(const double, const std::string &unit);        

  // Accessors

  int coordKind() const; 

  double mass() const;        
  double mass(const std::string &unit) const; 

  double unitLength() const ; 
  double unitLength(const std::string &unit) const ; 

  virtual double getRmb() const;

  virtual double getRms() const;

  virtual double getSpecificAngularMomentum(double rr) const;

  virtual double getPotential(double const pos[4], double l_cst) const;

  double deltaMin() const;

  void deltaMin(double h1);

  double deltaMax() const;

  virtual double deltaMax(double const pos[8], double delta_max_external) const;

  void deltaMax(double h1);

  double deltaMaxOverR() const; 
  void deltaMaxOverR(double t); 

  bool keplerian() const; 
  void keplerian(bool); 

  virtual void cartesianVelocity(double const coord[8], double vel[3]);

  virtual double SysPrimeToTdot(const double coord[4], const double v[3]) const;

  virtual void circularVelocity(double const pos[4], double vel[4],
                                double dir=1.) const ;

  virtual void zamoVelocity(double const pos[4], double vel[4]) const ;

  virtual void nullifyCoord(double coord[8]) const;

  virtual void nullifyCoord(double coord[8], double& tdot2) const;

  virtual void normalizeFourVel(double coord[8]) const;

  virtual void normalizeFourVel(double const pos[4], double fourvel[4]) const;


  virtual double ScalarProd(const double pos[4],
                    const double u1[4], const double u2[4]) const; 

  double norm(const double pos[4],
                      const double u1[4]) const; 

  void multiplyFourVect(double vect[4], double a) const;

  void addFourVect(double u1[4], double const u2[4]) const;

  void projectFourVect(double const pos[4], double u1[4], double const u2[4]) const;

  void dualOneForm(double const IN_ARRAY1_1[4], double const IN_ARRAY1_2[4], double ARGOUT_ARRAY1[4]) const ;

  virtual void observerTetrad(obskind_t obskind,
                              double const pos[4], double fourvel[4],
                              double screen1[4], double screen2[4],
                              double screen3[4]) const ;

  virtual void observerTetrad(double const pos[4], double fourvel[4],
                              double screen1[4], double screen2[4],
                              double screen3[4]) const ;

  void GramSchmidt(double const pos[4], double u0[4],
                   double u1[4], double u2[4], double u3[4]) const;
  // Outputs

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

  // Keep argument names for swig!
  virtual void gmunu(double ARGOUT_ARRAY2[4][4], double const IN_ARRAY1[4]) const;

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

  // Keep argument names for swig!
  virtual void gmunu_up(double ARGOUT_ARRAY2[4][4], const double IN_ARRAY1[4]) const;

  // Keep argument names for swig!
  virtual void jacobian(double ARGOUT_ARRAY3[4][4][4], const double IN_ARRAY1[4]) const;

  // Keep argument names for swig!
  virtual void gmunu_up_and_jacobian(double ARGOUT_ARRAY2[4][4], double ARGOUT_ARRAY3[4][4][4], const double IN_ARRAY1[4]) const;

  virtual void computeNBeta(const double coord[4],double &NN,double beta[3]) const;

  virtual double christoffel(const double coord[4],
                             const int alpha, const int mu, const int nu) const;

  virtual int christoffel(double dst[4][4][4], const double coord[4]) const ;



  virtual int myrk4(Worldline * line, state_t const &coord, double h, state_t &res) const;
  virtual int myrk4(Worldline * line, const double coord[8], double h, double res[8]) const = delete;
  
  virtual int myrk4_adaptive(Gyoto::Worldline* line, state_t const &coord,
                             double lastnorm, double normref,
                             state_t &coordnew, double h0, double& h1,
                             double deltamax=GYOTO_DEFAULT_DELTA_MAX) const;
  virtual int myrk4_adaptive(Gyoto::Worldline* line, const double coord[8],
                             double lastnorm, double normref,
                             double coordnew[8], double h0, double& h1,
                             double deltamax=GYOTO_DEFAULT_DELTA_MAX) const = delete;


  virtual int isStopCondition(double const coord[8]) const;

  virtual int diff(state_t const &x, state_t &dxdt, double mass) const ;
  virtual int diff(state_t const &x, state_t &dxdt)  const = delete;
  virtual int diff(const double y[8], double res[8]) const = delete ;

  virtual int diff31(state_t const &x, state_t &dxdt, double mass) const ;

  virtual void setParticleProperties(Gyoto::Worldline* line,
                                     double const coord[8]) const;
  

};

#endif