GyotoWorldline.h

Go to the documentation of this file.

/*
    Copyright 2011-2025 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 __GyotoWorldline_H_
#define __GyotoWorldline_H_ 

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

#include <GyotoDefs.h>

#ifdef GYOTO_HAVE_BOOST_INTEGRATORS
# include <functional>
# include <array>
# include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
# include <boost/numeric/odeint/algebra/vector_space_algebra.hpp>
# include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>


#endif

namespace Gyoto {
  class Worldline;
  class FactoryMessenger;
}

#include <GyotoSmartPointer.h>
#include <GyotoMetric.h>
#include <GyotoScreen.h>
#include <GyotoHooks.h>


#define GYOTO_WORLDLINE_PROPERTIES(c)                                   \
  Gyoto::Property("Gyoto::Worldline", "Time-like or null geodesic."),   \
    GYOTO_PROPERTY_BOOL(c, Integ31, Integ4D, _integ31,                  \
                        "Whether to integrate the 3+1 or 4D equation of geodesics.") \
    GYOTO_PROPERTY_BOOL(c, HighOrderImages, PrimaryOnly, _secondary,    \
                        "Whether to stop Photon integration at 180° deflection.") \
    GYOTO_PROPERTY_BOOL(c, ParallelTransport, NoParallelTransport, _parallelTransport, \
                        "Whether to perform parallel transport of a local triad (used for polarization).") \
    GYOTO_PROPERTY_BOOL(c, CheckBasis, NoCheckBasis, _checkBasis, \
                        "Whether to check basis orthogonality during parallel transport.") \
    GYOTO_PROPERTY_DOUBLE(c, MaxCrossEqplane, _maxCrossEqplane,         \
                          "Maximum number of crossings of the equatorial plane allowed for this worldline") \
    GYOTO_PROPERTY_DOUBLE(c, RelTol, _relTol,                           \
                          "Relative tolerance for the adaptive step integrators.") \
    GYOTO_PROPERTY_DOUBLE(c, AbsTol, _absTol,                           \
                          "Absolute tolerance for the adaptive step integrators.") \
    GYOTO_PROPERTY_DOUBLE(c, NormTol, _normTol,                         \
                          "Tolerance on 4-velocity norm variations.")   \
    GYOTO_PROPERTY_DOUBLE(c, DeltaMaxOverR, _deltaMaxOverR,             \
                          "Maximum value of step/distance from center of mass.") \
    GYOTO_PROPERTY_DOUBLE(c, DeltaMax, _deltaMax, "Maximum step (geometrical units).")  \
    GYOTO_PROPERTY_DOUBLE(c, DeltaMin, _deltaMin, "Minimum step (geometrical units).")  \
    GYOTO_PROPERTY_STRING(c, Integrator, _integrator,                   \
                          "Name of integrator (\"runge_kutta_fehlberg78\").")                   \
    GYOTO_PROPERTY_SIZE_T(c, MaxIter, _maxiter,                         \
                          "Maximum number of integration steps.")       \
    GYOTO_PROPERTY_BOOL(c, Adaptive, NonAdaptive, _adaptive,            \
                        "Whether to use an adaptive step.")             \
    GYOTO_PROPERTY_DOUBLE_UNIT(c, MinimumTime, _tMin,                   \
                               "Do not integrate earlier than this date (geometrical_time).") \
    GYOTO_PROPERTY_DOUBLE_UNIT(c, Delta, _delta,                        \
                               "Initial integration step (geometrical units).")         \
    GYOTO_PROPERTY_VECTOR_DOUBLE(c, InitCoord, _initCoord,              \
                                 "Initial 8-coordinate.")               \
    GYOTO_PROPERTY_METRIC(c, Metric, _metric,                           \
                          "The geometry of space-time at this end of the Universe.")


#define GYOTO_WORLDLINE_ACCESSORS(c)                                    \
  void c::_integ31(bool s) {integ31(s);}                                \
  bool c::_integ31() const {return integ31();}                          \
  void c::_secondary(bool s) {secondary(s);}                            \
  bool c::_secondary() const {return secondary();}                      \
  void c::_parallelTransport(bool s) {parallelTransport(s);}            \
  bool c::_parallelTransport() const {return parallelTransport();}      \
  void c::_checkBasis(bool s) {checkBasis(s);}                          \
  bool c::_checkBasis() const {return checkBasis();}                    \
  void c::_adaptive(bool s) {adaptive(s);}                              \
  bool c::_adaptive() const {return adaptive();}                        \
  void c::_maxCrossEqplane(double max){maxCrossEqplane(max);}           \
  double c::_maxCrossEqplane()const{return maxCrossEqplane();}          \
  void c::_relTol(double f){relTol(f);}                                 \
  double c::_relTol()const{return relTol();}                            \
  void c::_absTol(double f){absTol(f);}                                 \
  double c::_absTol()const{return absTol();}                            \
  void c::_deltaMin(double f){deltaMin(f);}                             \
  double c::_deltaMin()const{return deltaMin();}                        \
  void c::_deltaMax(double f){deltaMax(f);}                             \
  double c::_deltaMax()const{return deltaMax();}                        \
  void c::_normTol(double f){normTol(f);}                               \
  double c::_normTol()const{return normTol();}                          \
  void c::_deltaMaxOverR(double f){deltaMaxOverR(f);}                   \
  double c::_deltaMaxOverR()const{return deltaMaxOverR();}              \
  void c::_delta(double f){delta(f);}                                   \
  double c::_delta()const{return delta();}                              \
  void c::_delta(double f, std::string const &u){delta(f, u);}          \
  double c::_delta(std::string const &u)const{return delta(u);}         \
  void c::_tMin(double f){tMin(f);}                                     \
  double c::_tMin()const{return tMin();}                                \
  void c::_tMin(double f, std::string const &u){tMin(f, u);}            \
  double c::_tMin(std::string const &u)const{return tMin(u);}           \
  void c::_maxiter(size_t f){maxiter(f);}                               \
  size_t c::_maxiter()const{return maxiter();}                          \
  void c::_integrator(std::string const &f){integrator(f);}             \
  std::string c::_integrator() const {return integrator();}             \
  std::vector<double> c::_initCoord()const{return initCoord();}         \
  void c::_initCoord(std::vector<double> const &f){initCoord(f);}       \
  void c::_metric(SmartPointer<Metric::Generic>f){metric(f);}           \
  SmartPointer<Metric::Generic> c::_metric() const{return metric();}


#define GYOTO_WORLDLINE_PROPERTY_END(c, a) \
  GYOTO_WORLDLINE_PROPERTIES(c)            \
  GYOTO_PROPERTY_END(c, a)                 \
  GYOTO_WORLDLINE_ACCESSORS(c)


#define GYOTO_WORLDLINE                                 \
  void _delta(const double delta);                      \
  void _delta(double, const std::string &unit);         \
  double _delta() const ;                               \
  double _delta(const std::string &unit) const ;        \
  void _tMin(const double tmin);                        \
  void _tMin(double, const std::string &unit);          \
  double _tMin() const ;                                \
  double _tMin(const std::string &unit) const ;         \
  void _adaptive (bool mode) ;                          \
  bool _adaptive () const ;                             \
  void _secondary (bool sec) ;                          \
  bool _secondary () const ;                            \
  void _integ31 (bool sec) ;                            \
  bool _integ31 () const ;                              \
  void _parallelTransport (bool sec) ;                  \
  bool _parallelTransport () const ;                    \
  void _checkBasis (bool sec) ;                         \
  bool _checkBasis () const ;                           \
  void _maxiter (size_t miter) ;                        \
  size_t _maxiter () const ;                            \
  void _integrator(std::string const & type);           \
  std::string _integrator() const ;                     \
  double _deltaMin() const;                             \
  void _deltaMin(double h1);                            \
  double _normTol() const;                              \
  void _normTol(double normtol);                        \
  void _absTol(double);                                 \
  double _absTol()const;                                \
  void _maxCrossEqplane(double);                        \
  double _maxCrossEqplane()const;                       \
  void _relTol(double);                                 \
  double _relTol()const;                                \
  void _deltaMax(double h1);                            \
  double _deltaMax()const;                              \
  double _deltaMaxOverR() const;                        \
  void _deltaMaxOverR(double t);                        \
  std::vector<double> _initCoord()const;                \
  void _initCoord(std::vector<double> const&f);         \
  void _metric(SmartPointer<Metric::Generic>);          \
  SmartPointer<Metric::Generic> _metric() const;

class Gyoto::Worldline
:  protected Gyoto::Hook::Listener
{

  // Data : 
  // -----
 public:
  int stopcond; 

 protected:
  SmartPointer<Gyoto::Metric::Generic> metric_ ; 
  double* tau_; 
  double* x0_;
  double* x1_;
  double* x2_;
  double* x3_;
  double* x0dot_;
  double* x1dot_;
  double* x2dot_;
  double* x3dot_;
  double* ep0_;
  double* ep1_;
  double* ep2_;
  double* ep3_;
  double* et0_;
  double* et1_;
  double* et2_;
  double* et3_;
  size_t x_size_;
  size_t imin_;
  size_t i0_;  
  size_t imax_;
  bool   adaptive_; 

  bool secondary_;

  bool parallel_transport_;

  bool check_basis_;

  double delta_;


  double tmin_;

  double * cst_; 
  size_t cst_n_; 
  int wait_pos_; 
  double * init_vel_; 
  size_t maxiter_ ; 

  double delta_min_;

  double delta_max_;

  double delta_max_over_r_;

  double abstol_;

  double reltol_;

  double normtol_;

  double maxCrossEqplane_;

  // Constructors - Destructor
  // -------------------------
 public: 
  Worldline() ; 
  
  Worldline(const Worldline& ) ;                
  

  Worldline(Worldline* orig, size_t i0, int dir, double step_max) ;

  virtual ~Worldline() ;                        

  size_t getImin() const; 
  size_t getImax() const; 
  size_t getI0() const; 

  virtual double getMass() const = 0; 
  void   metric(SmartPointer<Metric::Generic>); 
  SmartPointer<Metric::Generic> metric() const; 

  void initCoord(std::vector<double> const&); 
  std::vector<double> initCoord() const; 

  virtual void   setInitCoord(const double coord[8], int dir,
                              double const Ephi[4],
                              double const Etheta[4]);


  virtual void   setInitCoord(const double coord[8], int dir = 0);

  virtual void setInitCoord(double const pos[4], double const vel[3], int dir=0);

  virtual void setPosition(double const pos[4]); 
  virtual void setVelocity(double const vel[3]); 

  void reset() ; 
  void reInit() ; 

  virtual std::string className() const ; 
  virtual std::string className_l() const ; 

  void integrator(std::string const & type);

  std::string integrator() const ;

  void integ31(bool integ);

  bool integ31() const ;

  double deltaMin() const;

  void deltaMin(double h1);

  double deltaMax() const;

  double normTol() const; 
  void normTol(double normtol); 
  void absTol(double); 
  double absTol()const; 
  void relTol(double); 
  double relTol()const; 

  void maxCrossEqplane(double); 
  double maxCrossEqplane()const; 

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

  void deltaMax(double h1);

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

  // Memory management
  // ----------------- 
 protected:
  virtual void xAllocate(); 

  virtual void xAllocate(size_t size); 

  virtual size_t xExpand(int dir); 

  virtual void xExpand(double * &x, int dir); 

  virtual void eAllocate (); 

  virtual void eDeallocate (); 

  virtual void eExpand(int dir); 

  // Mutators / assignment
  // ---------------------
 public:
  void delta(const double delta); 
  void delta(double, const std::string &unit);   
  double delta() const ; 
  double delta(const std::string &unit) const ;  
  double tMin() const ; 
  double tMin(const std::string &unit) const ;  
  void tMin(double tlim); 
  void tMin(double, const std::string &unit);   
  void adaptive (bool mode) ; 
  bool adaptive () const ; 
  void secondary (bool sec) ; 
  bool secondary () const ; 
  void parallelTransport (bool pt) ; 
  bool parallelTransport () const ; 
  void checkBasis (bool cb) ; 
  bool checkBasis () const ; 
  void maxiter (size_t miter) ; 
  size_t maxiter () const ; 

  double const * getCst() const ; 


  void setCst(double const * cst, size_t const ncsts) ;


  void constantsOfMotion(std::vector<double> const cstv) ;


  std::vector<double> constantsOfMotion() const ;


  void setInitialCondition(SmartPointer<Metric::Generic> gg, 
                           const double coord[8],
                           const int dir,
                           double const Ephi[4], double const Etheta[4]) ;


  void setInitialCondition(SmartPointer<Metric::Generic> gg,
                           const double coord[8],
                           const int dir) ;

  void getInitialCoord(std::vector<double> &dest) const; 

  void getCoord(size_t index, Gyoto::state_t &dest) const; 

  void getCoord(size_t index, Gyoto::state_t &dest) ; 

  void getCoord(double date, Gyoto::state_t &dest, bool proper=false); 

  void getCartesianPos(size_t index, double dest[4]) const; 


  virtual void xStore(size_t ind, state_t const &coord, double tau) ; 
  virtual void xStore(size_t ind, state_t const &coord) =delete; 
  virtual void xStore(size_t ind, double const coord[8]) = delete; 
  virtual void xFill(double tlim, bool proper=false) ; 



  // Accessors
  // ---------
 public:
  size_t get_nelements() const;

  void get_t(double *dest) const;

  void get_tau(double *dest) const;
  

  void getCartesian(double const * const dates, size_t const n_dates,
                double * const x, double * const y,
                double * const z, double * const xprime=NULL,
                double * const yprime=NULL,  double * const zprime=NULL) ;

  void get_xyz(double* x, double *y, double *z) const;

  void getCoord(double const * const dates, size_t const n_dates,
                double * const x1dest,
                double * const x2dest, double * const x3dest,
                double * const x0dot=NULL,  double * const x1dot=NULL,
                double * const x2dot=NULL,  double * const x3dot=NULL,
                double * ep0=NULL, double * ep1=NULL, double * ep2=NULL, double * ep3=NULL,
                double * et0=NULL, double * et1=NULL, double * et2=NULL, double * et3=NULL,
                double * otime=NULL, bool proper=false) ;

  void getCoord(double *x0, double *x1, double *x2, double *x3) const ;

  void checkPhiTheta(double coord[8]) const;

  void getSkyPos(SmartPointer<Screen> screen, double *dalpha, double *ddellta, double *dD) const;

  void get_dot(double *x0dot, double *x1dot, double *x2dot, double *x3dot) const ;

  void get_prime(double *x1prime, double *x2prime, double *x3prime) const ;
  
  // Outputs
  // -------
 public:
  //virtual void sauve(FILE *) const ;            ///< Save in a file
  void save_txyz(char * fichierxyz) const ;            
  void save_txyz(char* const filename, double const t1, double const  mass_sun,
                 double const distance_kpc, std::string const unit, SmartPointer<Screen> sc = NULL);

 protected:
  virtual void tell(Gyoto::Hook::Teller*);
  void checkBasis(state_t &coord) const;

  class IntegState {
  public:
    class Generic;
    class Legacy;
#ifdef GYOTO_HAVE_BOOST_INTEGRATORS
    class Boost;
#endif
  };


  SmartPointer<Worldline::IntegState::Generic> state_;
};

#ifndef GYOTO_SWIGIMPORTED

class Gyoto::Worldline::IntegState::Generic : public SmartPointee {
  friend class Gyoto::SmartPointer<Gyoto::Worldline::IntegState::Generic>;

 protected:

  Worldline * line_;
  // add 31 flag
  bool integ_31_; 
  double delta_; 
  bool adaptive_; 
  bool parallel_transport_; 
  bool check_basis_; 
  double norm_; 
  double normref_; 

  Gyoto::SmartPointer<Gyoto::Metric::Generic> gg_;

 public:
  Generic(Worldline *parent);

  virtual ~Generic();

  virtual Generic * clone(Worldline*newparent) const =0 ;

  virtual void init(Worldline * line, const state_t &coord, const double delta);
  virtual void init(Worldline * line, const double *coord, const double delta) = delete;

  virtual void init();

  virtual void checkNorm(double coord[8]);

  virtual std::string kind()=0;

  void integ31(bool integ);
  bool integ31() const ;


  virtual int nextStep(state_t &coord, double &tau, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0;
  virtual int nextStep(double *coord, double h1max=GYOTO_DEFAULT_DELTA_MAX) = delete;


  virtual void doStep(state_t const &coordin, 
                      double step,
                      state_t &coordout)=0;
  virtual void doStep(double const coordin[8], 
                      double step,
                      double coordout[8]) = delete;

};

class Gyoto::Worldline::IntegState::Legacy : public Generic {
  friend class Gyoto::SmartPointer<Gyoto::Worldline::IntegState::Legacy>;

 private:
  state_t coord_; 

 public:

  Legacy(Worldline *parent);
  Legacy * clone(Worldline*newparent) const ;
  using Generic::init;
  void init(Worldline * line, const state_t &coord, const double delta);
  virtual std::string kind();

  virtual int nextStep(state_t &coord, double &tau, double h1max=1e6);

  virtual void doStep(state_t const &coordin, 
                      double step,
                      state_t &coordout);

  virtual ~Legacy();
};

#ifdef GYOTO_HAVE_BOOST_INTEGRATORS

class Gyoto::Worldline::IntegState::Boost : public Generic {
  friend class Gyoto::SmartPointer<Gyoto::Worldline::IntegState::Boost>;
  friend class Gyoto::SmartPointer<Gyoto::Worldline::IntegState>;
 public:
  enum Kind {runge_kutta_cash_karp54,
             runge_kutta_fehlberg78,
             runge_kutta_dopri5,
             runge_kutta_cash_karp54_classic };
 private:
  Kind kind_;

  typedef std::function<boost::numeric::odeint::controlled_step_result
    (state_t&, double&, double&)> try_step_t;
  typedef std::function<void(state_t&, double)> do_step_t;
  typedef std::function<void(const state_t &/*x*/,
                             state_t & /*dxdt*/,
                             const double /* t*/ )> system_t;

  try_step_t try_step_;

  do_step_t do_step_;

  void setup_stepper();
 public:
  system_t system_;
  template<class Value, class Algebra = boost::numeric::odeint::range_algebra, class Operations = boost::numeric::odeint::default_operations>
  struct my_error_checker {
    using value_type = Value;
    using default_checker = boost::numeric::odeint::default_error_checker<Value, Algebra, Operations>;

    default_checker checker;
    const Metric::Generic* met; // Metric gg_
    Value normTol;              // tolerance on the norm
    const double& ref_norm;     // pointer for getting normref_

    my_error_checker(Value abs_tol,
                     Value rel_tol,
                     const Metric::Generic* met,
                     const double& ref_norm,
                     Value normTol_,
                     const do_step_t& do_step)
      : checker(abs_tol, rel_tol), met(met), normTol(normTol_), ref_norm(ref_norm), do_step_(do_step) {}

    const do_step_t& do_step_;
    // Signature for Boost.Odeint
    template<class State, class Deriv, class Err, class Time>
    Value error(Algebra &algebra,
                const State &x_old,
                const Deriv &dxdt_old,
                Err &x_err,
                const Time &dt) const
    {
      // Standard Boost numerical error
      Value err_default = checker.error(algebra, x_old, dxdt_old, x_err, dt);
      
      // Compute new vector (4-position, 4-velocity, polar basis if parallel_transport)
      state_t x_new=x_old;
      do_step_(x_new, dt);

      // Compute norm of 4-velocity
      GYOTO_DEBUG << "Computing norm" << std::endl;
      Value u_norm = met->ScalarProd(&x_new[0], &x_new[4], &x_new[4]);
      GYOTO_DEBUG_EXPR(u_norm);

      // Compute drifting of norm
      Value norm_drift = fabs(u_norm - ref_norm) / normTol;
      GYOTO_DEBUG_EXPR(ref_norm);

      GYOTO_DEBUG << "err_default = " << err_default
          << " norm_drift = " << norm_drift
          << " max = " << std::max(err_default, norm_drift)
          << std::endl;
      
      // Return max between numerical error and drift of the norm if drift > 1
      // This avoid to the norm_drift to affects too much the length's computation of the next step
      return std::max(err_default, (norm_drift > 1.0 ? norm_drift : 0.0));


    }
  };


  Boost(Worldline* parent, std::string type);


  Boost(Worldline* parent, Kind type);
  Boost * clone(Worldline* newparent) const ;
  virtual ~Boost();
  virtual void init();
  virtual void init(Worldline * line, const state_t &coord, const double delta);
  virtual int nextStep(state_t &coord, double &tau, double h1max=1e6);
  virtual void doStep(state_t const &coordin, 
                      double step,
                      state_t &coordout);
  virtual std::string kind();

  private:
    double checkBasis(state_t const &coord) const;
    void reorthonormalizeBasis(state_t &coord) const;
  
};
#endif 
#endif 
#endif