25 #ifndef __GyotoWorldline_H_
26 #define __GyotoWorldline_H_
34 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS
35 # include <functional>
37 # include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
42 class FactoryMessenger;
57 #define GYOTO_WORLDLINE_PROPERTIES(c) \
58 GYOTO_PROPERTY_BOOL(c, HighOrderImages, PrimaryOnly, _secondary) \
59 GYOTO_PROPERTY_DOUBLE(c, RelTol, _relTol) \
60 GYOTO_PROPERTY_DOUBLE(c, AbsTol, _absTol) \
61 GYOTO_PROPERTY_DOUBLE(c, DeltaMaxOverR, _deltaMaxOverR) \
62 GYOTO_PROPERTY_DOUBLE(c, DeltaMax, _deltaMax) \
63 GYOTO_PROPERTY_DOUBLE(c, DeltaMin, _deltaMin) \
64 GYOTO_PROPERTY_STRING(c, Integrator, _integrator) \
65 GYOTO_PROPERTY_SIZE_T(c, MaxIter, _maxiter) \
66 GYOTO_PROPERTY_BOOL(c, Adaptive, NonAdaptive, _adaptive) \
67 GYOTO_PROPERTY_DOUBLE_UNIT(c, MinimumTime, _tMin) \
68 GYOTO_PROPERTY_DOUBLE_UNIT(c, Delta, _delta) \
69 GYOTO_PROPERTY_VECTOR_DOUBLE(c, InitCoord, _initCoord) \
70 GYOTO_PROPERTY_METRIC(c, Metric, _metric)
87 #define GYOTO_WORLDLINE_ACCESSORS(c) \
88 void c::_secondary(bool s) {secondary(s);} \
89 bool c::_secondary() const {return secondary();} \
90 void c::_adaptive(bool s) {adaptive(s);} \
91 bool c::_adaptive() const {return adaptive();} \
92 void c::_relTol(double f){relTol(f);} \
93 double c::_relTol()const{return relTol();} \
94 void c::_absTol(double f){absTol(f);} \
95 double c::_absTol()const{return absTol();} \
96 void c::_deltaMin(double f){deltaMin(f);} \
97 double c::_deltaMin()const{return deltaMin();} \
98 void c::_deltaMax(double f){deltaMax(f);} \
99 double c::_deltaMax()const{return deltaMax();} \
100 void c::_deltaMaxOverR(double f){deltaMaxOverR(f);} \
101 double c::_deltaMaxOverR()const{return deltaMaxOverR();} \
102 void c::_delta(double f){delta(f);} \
103 double c::_delta()const{return delta();} \
104 void c::_delta(double f, std::string const &u){delta(f, u);} \
105 double c::_delta(std::string const &u)const{return delta(u);} \
106 void c::_tMin(double f){tMin(f);} \
107 double c::_tMin()const{return tMin();} \
108 void c::_tMin(double f, std::string const &u){tMin(f, u);} \
109 double c::_tMin(std::string const &u)const{return tMin(u);} \
110 void c::_maxiter(size_t f){maxiter(f);} \
111 size_t c::_maxiter()const{return maxiter();} \
112 void c::_integrator(std::string const &f){integrator(f);} \
113 std::string c::_integrator() const {return integrator();} \
114 std::vector<double> c::_initCoord()const{return initCoord();} \
115 void c::_initCoord(std::vector<double> const &f){initCoord(f);} \
116 void c::_metric(SmartPointer<Metric::Generic>f){metric(f);} \
117 SmartPointer<Metric::Generic> c::_metric() const{return metric();}
126 #define GYOTO_WORLDLINE_PROPERTY_END(c, a) \
127 GYOTO_WORLDLINE_PROPERTIES(c) \
128 GYOTO_PROPERTY_END(c, a) \
129 GYOTO_WORLDLINE_ACCESSORS(c)
139 #define GYOTO_WORLDLINE \
140 void _delta(const double delta); \
141 void _delta(double, const std::string &unit); \
142 double _delta() const ; \
143 double _delta(const std::string &unit) const ; \
144 void _tMin(const double tmin); \
145 void _tMin(double, const std::string &unit); \
146 double _tMin() const ; \
147 double _tMin(const std::string &unit) const ; \
148 void _adaptive (bool mode) ; \
149 bool _adaptive () const ; \
150 void _secondary (bool sec) ; \
151 bool _secondary () const ; \
152 void _maxiter (size_t miter) ; \
153 size_t _maxiter () const ; \
154 void _integrator(std::string const & type); \
155 std::string _integrator() const ; \
156 double _deltaMin() const; \
157 void _deltaMin(double h1); \
158 void _absTol(double); \
159 double _absTol()const; \
160 void _relTol(double); \
161 double _relTol()const; \
162 void _deltaMax(double h1); \
163 double _deltaMax()const; \
164 double _deltaMaxOverR() const; \
165 void _deltaMaxOverR(double t); \
166 std::vector<double> _initCoord()const; \
167 void _initCoord(std::vector<double> const&f); \
168 void _metric(SmartPointer<Metric::Generic>); \
169 SmartPointer<Metric::Generic> _metric() const;
340 size_t getI0()
const;
342 virtual double getMass()
const = 0;
346 void initCoord(std::vector<double>
const&);
347 std::vector<double> initCoord()
const;
349 virtual void setInitCoord(
const double coord[8],
int dir = 0);
358 virtual void setInitCoord(
double pos[4],
double vel[3],
int dir=1);
417 virtual double deltaMax(
double const pos[8],
double delta_max_external)
const;
451 virtual size_t xExpand(
int dir);
461 virtual void xExpand(
double * &x,
int dir);
468 void delta(
double,
const std::string &unit);
469 double delta()
const ;
470 double delta(
const std::string &unit)
const ;
471 double tMin()
const ;
472 double tMin(
const std::string &unit)
const ;
473 void tMin(
double tlim);
474 void tMin(
double,
const std::string &unit);
486 double const *
getCst()
const ;
493 void setCst(
double const * cst,
size_t const ncsts) ;
504 const double coord[8],
508 void getCoord(
size_t index,
double dest[8])
const;
512 virtual void xStore(
size_t ind,
double coord[8]) ;
513 virtual void xFill(
double tlim) ;
528 void get_t(
double *dest)
const;
548 void getCartesian(
double const *
const dates,
size_t const n_dates,
549 double *
const x,
double *
const y,
550 double *
const z,
double *
const xprime=NULL,
551 double *
const yprime=NULL,
double *
const zprime=NULL) ;
556 void get_xyz(
double* x,
double *y,
double *z)
const;
577 void getCoord(
double const *
const dates,
size_t const n_dates,
578 double *
const x1dest,
579 double *
const x2dest,
double *
const x3dest,
580 double *
const x0dot=NULL,
double *
const x1dot=NULL,
581 double *
const x2dot=NULL,
double *
const x3dot=NULL) ;
589 void getCoord(
double *x0,
double *x1,
double *x2,
double *x3)
const ;
609 void get_dot(
double *x0dot,
double *x1dot,
double *x2dot,
double *x3dot)
const ;
614 void get_prime(
double *x1prime,
double *x2prime,
double *x3prime)
const ;
620 void save_txyz(
char * fichierxyz)
const ;
621 void save_txyz(
char*
const filename,
double const t1,
double const mass_sun,
631 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS
643 #ifndef GYOTO_SWIGIMPORTED
693 virtual void
init(Worldline * line, const double *coord, const double delta);
716 virtual std::string
kind()=0;
723 virtual int
nextStep(double *coord, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0;
735 virtual void
doStep(double const coordin[8],
737 double coordout[8])=0;
760 Legacy(Worldline *parent);
763 void
init(Worldline * line, const double *coord, const double delta);
764 virtual std::string
kind();
766 virtual int
nextStep(double *coord, double h1max=1e6);
768 virtual void
doStep(double const coordin[8],
775 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS
792 enum
Kind {runge_kutta_cash_karp54,
793 runge_kutta_fehlberg78,
795 runge_kutta_cash_karp54_classic };
801 std::function<boost::numeric::odeint::controlled_step_result
805 std::function<void(std::array<double,8>&, double)> do_step_;
826 virtual int nextStep(
double *coord,
double h1max=1e6);
827 virtual void doStep(
double const coordin[8],
830 virtual std::string
kind();
double delta_max_over_r_
Numerical tuning parameter.
Definition: GyotoWorldline.h:293
size_t get_nelements() const
Get number of computed dates.
double * init_vel_
Hack in setParameters()
Definition: GyotoWorldline.h:261
std::function< boost::numeric::odeint::controlled_step_result(std::array< double, 8 > &, double &, double &)> try_step_
Stepper used by the adaptive-step integrator.
Definition: GyotoWorldline.h:802
size_t maxiter() const
Get maxiter_.
void getInitialCoord(double dest[8]) const
Get initial coordinate.
virtual void xStore(size_t ind, double coord[8])
Store coord at index ind.
void reset()
Forget integration, keeping initial contition.
void get_dot(double *x0dot, double *x1dot, double *x2dot, double *x3dot) const
Get computed 4-velocities.
Timelike or null geodesics.
Definition: GyotoWorldline.h:207
SmartPointer< Gyoto::Metric::Generic > metric_
The Gyoto::Metric in this part of the universe.
Definition: GyotoWorldline.h:217
double delta_
Initial integrating step.
Definition: GyotoWorldline.h:244
bool secondary() const
Get secondary_.
std::string integrator() const
Describe the integrator used by state_.
double const * getCst() const
Returns the worldline's cst of motion (if any)
size_t getImin() const
Get imin_.
Tellers tell Listeners when they mutate.
Reference-counting pointers.
Boost integrator.
Definition: GyotoWorldline.h:786
void get_t(double *dest) const
Get computed dates.
void reInit()
Reset and recompute particle properties.
double tmin_
Time limit for the integration (geometrical units)
Definition: GyotoWorldline.h:256
double tMin() const
Get tmin_.
Kind
Enum to represent the integrator flavour.
Definition: GyotoWorldline.h:792
double delta_min_
Minimum integration step for the adaptive integrator.
Definition: GyotoWorldline.h:271
virtual void checkNorm(double coord[8])
Check norm.
double * cst_
Worldline's csts of motion (if any)
Definition: GyotoWorldline.h:258
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.
double deltaMin() const
Get delta_min_.
bool adaptive_
Whether to use an adaptive step.
Definition: GyotoWorldline.h:658
SmartPointer< Metric::Generic > metric() const
Get metric.
virtual void tell(Gyoto::Hook::Teller *)
This is how a Teller tells.
void getSkyPos(SmartPointer< Screen > screen, double *dalpha, double *ddellta, double *dD) const
Get computed positions in sky coordinates.
double deltaMax() const
Get delta_max_.
void checkPhiTheta(double coord[8]) const
Bring θ in [0,Π] and φ in [0,2Π].
Gyoto::SmartPointer< Gyoto::Metric::Generic > gg_
The Metric in this end of the Universe.
Definition: GyotoWorldline.h:665
virtual std::string className() const
"Worldline"
size_t i0_
Index of initial condition in array.
Definition: GyotoWorldline.h:228
Home-brewed integrator.
Definition: GyotoWorldline.h:751
Gyoto ubiquitous macros and typedefs.
void setInitialCondition(SmartPointer< Metric::Generic > gg, const double coord[8], const int dir)
Set or re-set the initial condition prior to integration.
Base class for metric description.
void save_txyz(char *fichierxyz) const
Save in a file.
SmartPointer< Worldline::IntegState::Generic > state_
An object to hold the integration state.
Definition: GyotoWorldline.h:640
size_t x_size_
Size of x0_, x1_... arrays.
Definition: GyotoWorldline.h:226
double * x0_
t or T
Definition: GyotoWorldline.h:218
double absTol() const
Get abstol_.
bool adaptive_
Whether integration should use adaptive delta.
Definition: GyotoWorldline.h:230
double relTol() const
Get reltol_.
double * x2_
θ or y
Definition: GyotoWorldline.h:220
double * x3dot_
φdot or zdot
Definition: GyotoWorldline.h:225
double * x0dot_
tdot or Tdot
Definition: GyotoWorldline.h:222
Definition: GyotoWorldline.h:627
double norm_
Current norm of the 4-velocity.
Definition: GyotoWorldline.h:659
double reltol_
Absolute tolerance of the integrator.
Definition: GyotoWorldline.h:311
double normref_
Definition: GyotoWorldline.h:660
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.
size_t getI0() const
Get i0_.
Can be pointed to by a SmartPointer.
Definition: GyotoSmartPointer.h:78
bool secondary_
Experimental: choose 0 to compute only primary image.
Definition: GyotoWorldline.h:237
I might listen to a Teller.
Definition: GyotoHooks.h:64
virtual void doStep(double const coordin[8], double step, double coordout[8])=0
Make one step of exactly this size.
virtual std::string kind()=0
Return the integrator kind.
virtual ~Worldline()
Destructor.
double delta_
Integration step (current in case of adaptive_).
Definition: GyotoWorldline.h:657
double * x3_
φ or z
Definition: GyotoWorldline.h:221
void get_xyz(double *x, double *y, double *z) const
Get 3-position in cartesian coordinates for computed dates.
void get_prime(double *x1prime, double *x2prime, double *x3prime) const
Get computed 3-velocities.
size_t cst_n_
Number of constants of motion.
Definition: GyotoWorldline.h:259
size_t imin_
Minimum index for which x0_, x1_... have been computed.
Definition: GyotoWorldline.h:227
double * x1dot_
rdot or xdot
Definition: GyotoWorldline.h:223
double abstol_
Absolute tolerance of the integrator.
Definition: GyotoWorldline.h:302
double delta() const
Get delta_.
void getCartesianPos(size_t index, double dest[4]) const
Get Cartesian expression of 4-position at index.
virtual size_t xExpand(int dir)
Expand x0, x1 etc... to hold more elements.
Current state of a geodesic integration.
Definition: GyotoWorldline.h:648
size_t getImax() const
Get imax_.
Kind kind_
Integrator flavour.
Definition: GyotoWorldline.h:798
Description of the observer screen.
virtual void setInitCoord(const double coord[8], int dir=0)
Set Initial coordinate.
Worldline * line_
Worldline that we are integrating.
Definition: GyotoWorldline.h:656
double delta_max_
Maximum integration step for the adaptive integrator.
Definition: GyotoWorldline.h:280
virtual void xFill(double tlim)
Fill x0, x1... by integrating the Worldline from previously set inittial condition to time tlim...
virtual std::string className_l() const
"worldline"
virtual int nextStep(double *coord, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0
Make one step.
double deltaMaxOverR() const
Get delta_max_over_r_.
int wait_pos_
Hack in setParameters()
Definition: GyotoWorldline.h:260
bool adaptive() const
Get adaptive_.
virtual void xAllocate()
Allocate x0, x1 etc. with default size.
double * x2dot_
θdot or ydot
Definition: GyotoWorldline.h:224
void getCoord(size_t index, double dest[8]) const
Get coordinates corresponding to index.
Listen to me and I'll warn you when I change.
Definition: GyotoHooks.h:82
virtual void setPosition(double pos[4])
Set initial 4-position.
double * x1_
r or x
Definition: GyotoWorldline.h:219
int stopcond
Whether and why integration is finished.
Definition: GyotoWorldline.h:214
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:262
Worldline()
Default constructor.
virtual void setVelocity(double vel[3])
Set initial 3-velocity.
size_t imax_
Maximum index for which x0_, x1_... have been computed.
Definition: GyotoWorldline.h:229