25 #ifndef __GyotoWorldline_H_ 26 #define __GyotoWorldline_H_ 35 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS 36 # include <functional> 38 # include <boost/numeric/odeint/stepper/controlled_step_result.hpp> 43 class FactoryMessenger;
58 #define GYOTO_WORLDLINE_PROPERTIES(c) \ 59 Gyoto::Property("Gyoto::Worldline", "Time-like or null geodesic."), \ 60 GYOTO_PROPERTY_BOOL(c, HighOrderImages, PrimaryOnly, _secondary, \ 61 "Whether to stop Photon integration at 180° deflection.") \ 62 GYOTO_PROPERTY_BOOL(c, ParallelTransport, NoParallelTransport, _parallelTransport, \ 63 "Whether to perform parallel transport of a local triad (used for polarization).") \ 64 GYOTO_PROPERTY_DOUBLE(c, RelTol, _relTol, \ 65 "Relative tolerance for the adaptive step integrators.") \ 66 GYOTO_PROPERTY_DOUBLE(c, AbsTol, _absTol, \ 67 "Absolute tolerance for the adaptive step integrators.") \ 68 GYOTO_PROPERTY_DOUBLE(c, DeltaMaxOverR, _deltaMaxOverR, \ 69 "Maximum value of step/distance from center of mass.") \ 70 GYOTO_PROPERTY_DOUBLE(c, DeltaMax, _deltaMax, "Maximum step (geometrical units).") \ 71 GYOTO_PROPERTY_DOUBLE(c, DeltaMin, _deltaMin, "Minimum step (geometrical units).") \ 72 GYOTO_PROPERTY_STRING(c, Integrator, _integrator, \ 73 "Name of integrator (\"runge_kutta_fehlberg78\").") \ 74 GYOTO_PROPERTY_SIZE_T(c, MaxIter, _maxiter, \ 75 "Maximum number of integration steps.") \ 76 GYOTO_PROPERTY_BOOL(c, Adaptive, NonAdaptive, _adaptive, \ 77 "Whether to use an adaptive step.") \ 78 GYOTO_PROPERTY_DOUBLE_UNIT(c, MinimumTime, _tMin, \ 79 "Do not integrate earlier than this date (geometrical_time).") \ 80 GYOTO_PROPERTY_DOUBLE_UNIT(c, Delta, _delta, \ 81 "Initial integration step (geometrical units).") \ 82 GYOTO_PROPERTY_VECTOR_DOUBLE(c, InitCoord, _initCoord, \ 83 "Initial 8-coordinate.") \ 84 GYOTO_PROPERTY_METRIC(c, Metric, _metric, \ 85 "The geometry of space-time at this end of the Universe.") 102 #define GYOTO_WORLDLINE_ACCESSORS(c) \ 103 void c::_secondary(bool s) {secondary(s);} \ 104 bool c::_secondary() const {return secondary();} \ 105 void c::_parallelTransport(bool s) {parallelTransport(s);} \ 106 bool c::_parallelTransport() const {return parallelTransport();} \ 107 void c::_adaptive(bool s) {adaptive(s);} \ 108 bool c::_adaptive() const {return adaptive();} \ 109 void c::_relTol(double f){relTol(f);} \ 110 double c::_relTol()const{return relTol();} \ 111 void c::_absTol(double f){absTol(f);} \ 112 double c::_absTol()const{return absTol();} \ 113 void c::_deltaMin(double f){deltaMin(f);} \ 114 double c::_deltaMin()const{return deltaMin();} \ 115 void c::_deltaMax(double f){deltaMax(f);} \ 116 double c::_deltaMax()const{return deltaMax();} \ 117 void c::_deltaMaxOverR(double f){deltaMaxOverR(f);} \ 118 double c::_deltaMaxOverR()const{return deltaMaxOverR();} \ 119 void c::_delta(double f){delta(f);} \ 120 double c::_delta()const{return delta();} \ 121 void c::_delta(double f, std::string const &u){delta(f, u);} \ 122 double c::_delta(std::string const &u)const{return delta(u);} \ 123 void c::_tMin(double f){tMin(f);} \ 124 double c::_tMin()const{return tMin();} \ 125 void c::_tMin(double f, std::string const &u){tMin(f, u);} \ 126 double c::_tMin(std::string const &u)const{return tMin(u);} \ 127 void c::_maxiter(size_t f){maxiter(f);} \ 128 size_t c::_maxiter()const{return maxiter();} \ 129 void c::_integrator(std::string const &f){integrator(f);} \ 130 std::string c::_integrator() const {return integrator();} \ 131 std::vector<double> c::_initCoord()const{return initCoord();} \ 132 void c::_initCoord(std::vector<double> const &f){initCoord(f);} \ 133 void c::_metric(SmartPointer<Metric::Generic>f){metric(f);} \ 134 SmartPointer<Metric::Generic> c::_metric() const{return metric();} 143 #define GYOTO_WORLDLINE_PROPERTY_END(c, a) \ 144 GYOTO_WORLDLINE_PROPERTIES(c) \ 145 GYOTO_PROPERTY_END(c, a) \ 146 GYOTO_WORLDLINE_ACCESSORS(c) 156 #define GYOTO_WORLDLINE \ 157 void _delta(const double delta); \ 158 void _delta(double, const std::string &unit); \ 159 double _delta() const ; \ 160 double _delta(const std::string &unit) const ; \ 161 void _tMin(const double tmin); \ 162 void _tMin(double, const std::string &unit); \ 163 double _tMin() const ; \ 164 double _tMin(const std::string &unit) const ; \ 165 void _adaptive (bool mode) ; \ 166 bool _adaptive () const ; \ 167 void _secondary (bool sec) ; \ 168 bool _secondary () const ; \ 169 void _parallelTransport (bool sec) ; \ 170 bool _parallelTransport () const ; \ 171 void _maxiter (size_t miter) ; \ 172 size_t _maxiter () const ; \ 173 void _integrator(std::string const & type); \ 174 std::string _integrator() const ; \ 175 double _deltaMin() const; \ 176 void _deltaMin(double h1); \ 177 void _absTol(double); \ 178 double _absTol()const; \ 179 void _relTol(double); \ 180 double _relTol()const; \ 181 void _deltaMax(double h1); \ 182 double _deltaMax()const; \ 183 double _deltaMaxOverR() const; \ 184 void _deltaMaxOverR(double t); \ 185 std::vector<double> _initCoord()const; \ 186 void _initCoord(std::vector<double> const&f); \ 187 void _metric(SmartPointer<Metric::Generic>); \ 188 SmartPointer<Metric::Generic> _metric() const; 376 size_t getI0()
const;
378 virtual double getMass()
const = 0;
382 void initCoord(std::vector<double>
const&);
383 std::vector<double> initCoord()
const;
407 virtual void setInitCoord(
const double coord[8],
int dir,
408 double const Ephi[4],
409 double const Etheta[4]);
430 virtual void setInitCoord(
const double coord[8],
int dir = 0);
439 virtual void setInitCoord(
double const pos[4],
double const vel[3],
int dir=0);
498 virtual double deltaMax(
double const pos[8],
double delta_max_external)
const;
532 virtual size_t xExpand(
int dir);
541 virtual void xExpand(
double * &x,
int dir);
563 void delta(
double,
const std::string &unit);
564 double delta()
const ;
565 double delta(
const std::string &unit)
const ;
566 double tMin()
const ;
567 double tMin(
const std::string &unit)
const ;
568 void tMin(
double tlim);
569 void tMin(
double,
const std::string &unit);
583 double const *
getCst()
const ;
590 void setCst(
double const * cst,
size_t const ncsts) ;
605 const double coord[8],
607 double const Ephi[4],
double const Etheta[4]) ;
618 const double coord[8],
633 void getCoord(
size_t index, Gyoto::state_t &dest)
const;
639 void getCoord(
size_t index, Gyoto::state_t &dest) ;
646 void getCoord(
double date, Gyoto::state_t &dest,
bool proper=
false);
651 virtual void xStore(
size_t ind, state_t
const &coord,
double tau) ;
652 virtual void xStore(
size_t ind, state_t
const &coord) =
delete;
653 virtual void xStore(
size_t ind,
double const coord[8]) =
delete;
654 virtual void xFill(
double tlim,
bool proper=
false) ;
669 void get_t(
double *dest)
const;
674 void get_tau(
double *dest)
const;
693 void getCartesian(
double const *
const dates,
size_t const n_dates,
694 double *
const x,
double *
const y,
695 double *
const z,
double *
const xprime=NULL,
696 double *
const yprime=NULL,
double *
const zprime=NULL) ;
701 void get_xyz(
double* x,
double *y,
double *z)
const;
734 void getCoord(
double const *
const dates,
size_t const n_dates,
735 double *
const x1dest,
736 double *
const x2dest,
double *
const x3dest,
737 double *
const x0dot=NULL,
double *
const x1dot=NULL,
738 double *
const x2dot=NULL,
double *
const x3dot=NULL,
739 double * ep0=NULL,
double * ep1=NULL,
double * ep2=NULL,
double * ep3=NULL,
740 double * et0=NULL,
double * et1=NULL,
double * et2=NULL,
double * et3=NULL,
741 double * otime=NULL,
bool proper=
false) ;
749 void getCoord(
double *x0,
double *x1,
double *x2,
double *x3)
const ;
769 void get_dot(
double *x0dot,
double *x1dot,
double *x2dot,
double *x3dot)
const ;
774 void get_prime(
double *x1prime,
double *x2prime,
double *x3prime)
const ;
780 void save_txyz(
char * fichierxyz)
const ;
781 void save_txyz(
char*
const filename,
double const t1,
double const mass_sun,
791 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS 803 #ifndef GYOTO_SWIGIMPORTED 855 virtual void
init(Worldline * line, const state_t &coord, const double delta);
857 virtual void
init(Worldline * line, const double *coord, const double delta) = delete;
880 virtual std::string
kind()=0;
888 virtual int
nextStep(state_t &coord, double &tau, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0;
890 virtual int
nextStep(double *coord, double h1max=GYOTO_DEFAULT_DELTA_MAX) = delete;
902 virtual void
doStep(state_t const &coordin,
904 state_t &coordout)=0;
906 virtual void
doStep(double const coordin[8],
908 double coordout[8]) = delete;
934 Legacy(Worldline *parent);
937 void
init(Worldline * line, const state_t &coord, const double delta);
938 virtual std::string
kind();
940 virtual int
nextStep(state_t &coord, double &tau, double h1max=1e6);
942 virtual void
doStep(state_t const &coordin,
949 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS 966 enum
Kind {runge_kutta_cash_karp54,
967 runge_kutta_fehlberg78,
969 runge_kutta_cash_karp54_classic };
974 typedef std::function<boost::numeric::odeint::controlled_step_result
975 (state_t&,
double&,
double&)> try_step_t;
976 typedef std::function<void(state_t&, double)> do_step_t;
977 typedef std::function<void(
const state_t &,
979 const double )> system_t;
1005 virtual void init();
1007 virtual int nextStep(state_t &coord,
double &tau,
double h1max=1e6);
1008 virtual void doStep(state_t
const &coordin,
1011 virtual std::string
kind();
double delta_max_over_r_
Numerical tuning parameter.
Definition: GyotoWorldline.h:329
double * init_vel_
Hack in setParameters()
Definition: GyotoWorldline.h:297
void getSkyPos(SmartPointer< Screen > screen, double *dalpha, double *ddellta, double *dD) const
Get computed positions in sky coordinates.
size_t get_nelements() const
Get number of computed dates.
void reset()
Forget integration, keeping initial contition.
virtual void setPosition(double const pos[4])
Set initial 4-position.
double deltaMax() const
Get delta_max_.
Timelike or null geodesics.
Definition: GyotoWorldline.h:226
SmartPointer< Gyoto::Metric::Generic > metric_
The Gyoto::Metric in this part of the universe.
Definition: GyotoWorldline.h:236
double delta_
Initial integrating step.
Definition: GyotoWorldline.h:280
void get_tau(double *dest) const
Get computed proper times or values of the affine parameter.
virtual void xFill(double tlim, bool proper=false)
Fill x0, x1... by integrating the Worldline from previously set inittial condition to time tlim...
void save_txyz(char *fichierxyz) const
Save in a file.
Tellers tell Listeners when they mutate.
virtual void xStore(size_t ind, state_t const &coord, double tau)
Store coord at index ind.
bool parallel_transport_
Whether to parallel-transport base vectors.
Definition: GyotoWorldline.h:820
Reference-counting pointers.
virtual std::string className_l() const
"worldline"
Boost integrator.
Definition: GyotoWorldline.h:960
do_step_t do_step_
Stepper used by the non-adaptive-step integrator.
Definition: GyotoWorldline.h:985
void setInitialCondition(SmartPointer< Metric::Generic > gg, const double coord[8], const int dir, double const Ephi[4], double const Etheta[4])
Set or re-set the initial condition prior to integration.
void reInit()
Reset and recompute particle properties.
double tmin_
Time limit for the integration (geometrical units)
Definition: GyotoWorldline.h:292
Kind
Enum to represent the integrator flavour.
Definition: GyotoWorldline.h:966
double delta_min_
Minimum integration step for the adaptive integrator.
Definition: GyotoWorldline.h:307
double * tau_
proper time or affine parameter
Definition: GyotoWorldline.h:237
size_t getImax() const
Get imax_.
virtual void eExpand(int dir)
virtual void checkNorm(double coord[8])
Check norm.
double * cst_
Worldline's csts of motion (if any)
Definition: GyotoWorldline.h:294
double * et1_
Coordinate of Second base vector to parallel transport.
Definition: GyotoWorldline.h:251
double relTol() const
Get reltol_.
double delta() const
Get delta_.
void getInitialCoord(std::vector< double > &dest) const
Get initial coordinates + base vectors.
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)
Get the 6 Cartesian coordinates for specific dates.
virtual int nextStep(state_t &coord, double &tau, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0
Make one step.
bool adaptive_
Whether to use an adaptive step.
Definition: GyotoWorldline.h:819
virtual void tell(Gyoto::Hook::Teller *)
This is how a Teller tells.
bool secondary() const
Get secondary_.
Gyoto::SmartPointer< Gyoto::Metric::Generic > gg_
The Metric in this end of the Universe.
Definition: GyotoWorldline.h:827
double * et2_
Coordinate of Second base vector to parallel transport.
Definition: GyotoWorldline.h:252
void get_xyz(double *x, double *y, double *z) const
Get 3-position in cartesian coordinates for computed dates.
size_t i0_
Index of initial condition in array.
Definition: GyotoWorldline.h:256
void getCartesianPos(size_t index, double dest[4]) const
Get Cartesian expression of 4-position at index.
Obsolete: Home-brewed integrator.
Definition: GyotoWorldline.h:925
double tMin() const
Get tmin_.
double const * getCst() const
Returns the worldline's cst of motion (if any)
Gyoto ubiquitous macros and typedefs.
Base class for metric description.
SmartPointer< Worldline::IntegState::Generic > state_
An object to hold the integration state.
Definition: GyotoWorldline.h:800
size_t x_size_
Coordinate of Second base vector to parallel transport.
Definition: GyotoWorldline.h:254
double * x0_
t or T
Definition: GyotoWorldline.h:238
bool adaptive_
Whether integration should use adaptive delta.
Definition: GyotoWorldline.h:258
double * ep3_
Coordinate of first base vector to parallel transport.
Definition: GyotoWorldline.h:249
size_t getImin() const
Get imin_.
size_t maxiter() const
Get maxiter_.
double * x2_
θ or y
Definition: GyotoWorldline.h:240
double * x3dot_
φdot or zdot
Definition: GyotoWorldline.h:245
double * x0dot_
tdot or Tdot
Definition: GyotoWorldline.h:242
bool parallel_transport_
Whether to parallel transport a local triad.
Definition: GyotoWorldline.h:273
Definition: GyotoWorldline.h:787
double norm_
Current norm of the 4-velocity.
Definition: GyotoWorldline.h:821
double reltol_
Absolute tolerance of the integrator.
Definition: GyotoWorldline.h:347
double normref_
Definition: GyotoWorldline.h:822
double deltaMaxOverR() const
Get delta_max_over_r_.
double * et3_
Coordinate of Second base vector to parallel transport.
Definition: GyotoWorldline.h:253
Namespace for the Gyoto library.
Definition: GyotoAstrobj.h:43
virtual double getMass() const =0
Get mass of particule.
virtual Generic * clone(Worldline *newparent) const =0
Deep copy.
virtual void init()
Cache whatever needs to be cached.
virtual void doStep(state_t const &coordin, double step, state_t &coordout)=0
Make one step of exactly this size.
virtual void eAllocate()
Allocate ep0_ ... et3_.
void get_t(double *dest) const
Get computed dates.
double * ep1_
Coordinate of first base vector to parallel transport.
Definition: GyotoWorldline.h:247
bool adaptive() const
Get adaptive_.
void checkPhiTheta(double coord[8]) const
Bring θ in [0,Π] and φ in [0,2Π].
double deltaMin() const
Get delta_min_.
Can be pointed to by a SmartPointer.
Definition: GyotoSmartPointer.h:80
bool secondary_
Experimental: choose 0 to compute only primary image.
Definition: GyotoWorldline.h:265
try_step_t try_step_
Stepper used by the adaptive-step integrator.
Definition: GyotoWorldline.h:982
I might listen to a Teller.
Definition: GyotoHooks.h:64
void get_prime(double *x1prime, double *x2prime, double *x3prime) const
Get computed 3-velocities.
bool parallelTransport() const
Get parallel_transport_.
virtual std::string kind()=0
Return the integrator kind.
virtual ~Worldline()
Destructor.
double delta_
Integration step (current in case of adaptive_).
Definition: GyotoWorldline.h:818
double * x3_
φ or z
Definition: GyotoWorldline.h:241
void getCoord(size_t index, Gyoto::state_t &dest) const
Get coordinates+base vectors corresponding to index.
double * ep2_
Coordinate of first base vector to parallel transport.
Definition: GyotoWorldline.h:248
virtual std::string className() const
"Worldline"
size_t cst_n_
Number of constants of motion.
Definition: GyotoWorldline.h:295
size_t imin_
Minimum index for which x0_, x1_... have been computed.
Definition: GyotoWorldline.h:255
double * x1dot_
rdot or xdot
Definition: GyotoWorldline.h:243
double abstol_
Absolute tolerance of the integrator.
Definition: GyotoWorldline.h:338
size_t getI0() const
Get i0_.
virtual size_t xExpand(int dir)
Expand x0, x1 etc... to hold more elements.
Current state of a geodesic integration.
Definition: GyotoWorldline.h:808
virtual void setVelocity(double const vel[3])
Set initial 3-velocity.
Kind kind_
Integrator flavour.
Definition: GyotoWorldline.h:972
Description of the observer screen.
virtual void setInitCoord(const double coord[8], int dir, double const Ephi[4], double const Etheta[4])
Set Initial coordinate.
SmartPointer< Metric::Generic > metric() const
Get metric.
virtual void eDeallocate()
Deallocate ep0_ ... et3_.
Worldline * line_
Worldline that we are integrating.
Definition: GyotoWorldline.h:817
double delta_max_
Maximum integration step for the adaptive integrator.
Definition: GyotoWorldline.h:316
double * et0_
Coordinate of first base vector to parallel transport.
Definition: GyotoWorldline.h:250
int wait_pos_
Hack in setParameters()
Definition: GyotoWorldline.h:296
std::string integrator() const
Describe the integrator used by state_.
virtual void xAllocate()
Allocate x0, x1 etc. with default size.
double * x2dot_
θdot or ydot
Definition: GyotoWorldline.h:244
Listen to me and I'll warn you when I change.
Definition: GyotoHooks.h:82
void get_dot(double *x0dot, double *x1dot, double *x2dot, double *x3dot) const
Get computed 4-velocities.
double * x1_
r or x
Definition: GyotoWorldline.h:239
int stopcond
Whether and why integration is finished.
Definition: GyotoWorldline.h:233
void setCst(double const *cst, size_t const ncsts)
Set Metric-specific constants of motion.
size_t maxiter_
Maximum number of iterations when integrating.
Definition: GyotoWorldline.h:298
Worldline()
Default constructor.
double absTol() const
Get abstol_.
size_t imax_
Maximum index for which x0_, x1_... have been computed.
Definition: GyotoWorldline.h:257