Gyoto
GyotoWorldline.h
Go to the documentation of this file.
1 
6 /*
7  Copyright 2011-2015 Frederic Vincent, Thibaut Paumard
8 
9  This file is part of Gyoto.
10 
11  Gyoto is free software: you can redistribute it and/or modify
12  it under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  Gyoto is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  GNU General Public License for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with Gyoto. If not, see <http://www.gnu.org/licenses/>.
23  */
24 
25 #ifndef __GyotoWorldline_H_
26 #define __GyotoWorldline_H_
27 
28 #include <iostream>
29 #include <fstream>
30 #include <string>
31 
32 #include <GyotoDefs.h>
33 
34 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS
35 # include <functional>
36 # include <array>
37 # include <boost/numeric/odeint/stepper/controlled_step_result.hpp>
38 #endif
39 
40 namespace Gyoto {
41  class Worldline;
42  class FactoryMessenger;
43 }
44 
45 #include <GyotoSmartPointer.h>
46 #include <GyotoMetric.h>
47 #include <GyotoScreen.h>
48 #include <GyotoHooks.h>
49 
51 
57 #define GYOTO_WORLDLINE_PROPERTIES(c) \
58  Gyoto::Property("Gyoto::Worldline", "Time-like or null geodesic."), \
59  GYOTO_PROPERTY_BOOL(c, HighOrderImages, PrimaryOnly, _secondary, \
60  "Whether to stop Photon integration at 180° deflection.") \
61  GYOTO_PROPERTY_DOUBLE(c, RelTol, _relTol, \
62  "Relative tolerance for the adaptive step integrators.") \
63  GYOTO_PROPERTY_DOUBLE(c, AbsTol, _absTol, \
64  "Absolute tolerance for the adaptive step integrators.") \
65  GYOTO_PROPERTY_DOUBLE(c, DeltaMaxOverR, _deltaMaxOverR, \
66  "Maximum value of step/distance from center of mass.") \
67  GYOTO_PROPERTY_DOUBLE(c, DeltaMax, _deltaMax, "Maximum step (geometrical units).") \
68  GYOTO_PROPERTY_DOUBLE(c, DeltaMin, _deltaMin, "Minimum step (geometrical units).") \
69  GYOTO_PROPERTY_STRING(c, Integrator, _integrator, \
70  "Name of integrator (\"runge_kutta_fehlberg78\").") \
71  GYOTO_PROPERTY_SIZE_T(c, MaxIter, _maxiter, \
72  "Maximum number of integration steps.") \
73  GYOTO_PROPERTY_BOOL(c, Adaptive, NonAdaptive, _adaptive, \
74  "Whether to use an adaptive step.") \
75  GYOTO_PROPERTY_DOUBLE_UNIT(c, MinimumTime, _tMin, \
76  "Do not integrate earlier than this date (geometrical_time).") \
77  GYOTO_PROPERTY_DOUBLE_UNIT(c, Delta, _delta, \
78  "Initial integration step (geometrical units).") \
79  GYOTO_PROPERTY_VECTOR_DOUBLE(c, InitCoord, _initCoord, \
80  "Initial 8-coordinate.") \
81  GYOTO_PROPERTY_METRIC(c, Metric, _metric, \
82  "The geometry of space-time at this end of the Universe.")
83 
85 
99 #define GYOTO_WORLDLINE_ACCESSORS(c) \
100  void c::_secondary(bool s) {secondary(s);} \
101  bool c::_secondary() const {return secondary();} \
102  void c::_adaptive(bool s) {adaptive(s);} \
103  bool c::_adaptive() const {return adaptive();} \
104  void c::_relTol(double f){relTol(f);} \
105  double c::_relTol()const{return relTol();} \
106  void c::_absTol(double f){absTol(f);} \
107  double c::_absTol()const{return absTol();} \
108  void c::_deltaMin(double f){deltaMin(f);} \
109  double c::_deltaMin()const{return deltaMin();} \
110  void c::_deltaMax(double f){deltaMax(f);} \
111  double c::_deltaMax()const{return deltaMax();} \
112  void c::_deltaMaxOverR(double f){deltaMaxOverR(f);} \
113  double c::_deltaMaxOverR()const{return deltaMaxOverR();} \
114  void c::_delta(double f){delta(f);} \
115  double c::_delta()const{return delta();} \
116  void c::_delta(double f, std::string const &u){delta(f, u);} \
117  double c::_delta(std::string const &u)const{return delta(u);} \
118  void c::_tMin(double f){tMin(f);} \
119  double c::_tMin()const{return tMin();} \
120  void c::_tMin(double f, std::string const &u){tMin(f, u);} \
121  double c::_tMin(std::string const &u)const{return tMin(u);} \
122  void c::_maxiter(size_t f){maxiter(f);} \
123  size_t c::_maxiter()const{return maxiter();} \
124  void c::_integrator(std::string const &f){integrator(f);} \
125  std::string c::_integrator() const {return integrator();} \
126  std::vector<double> c::_initCoord()const{return initCoord();} \
127  void c::_initCoord(std::vector<double> const &f){initCoord(f);} \
128  void c::_metric(SmartPointer<Metric::Generic>f){metric(f);} \
129  SmartPointer<Metric::Generic> c::_metric() const{return metric();}
130 
132 
138 #define GYOTO_WORLDLINE_PROPERTY_END(c, a) \
139  GYOTO_WORLDLINE_PROPERTIES(c) \
140  GYOTO_PROPERTY_END(c, a) \
141  GYOTO_WORLDLINE_ACCESSORS(c)
142 
144 
151 #define GYOTO_WORLDLINE \
152  void _delta(const double delta); \
153  void _delta(double, const std::string &unit); \
154  double _delta() const ; \
155  double _delta(const std::string &unit) const ; \
156  void _tMin(const double tmin); \
157  void _tMin(double, const std::string &unit); \
158  double _tMin() const ; \
159  double _tMin(const std::string &unit) const ; \
160  void _adaptive (bool mode) ; \
161  bool _adaptive () const ; \
162  void _secondary (bool sec) ; \
163  bool _secondary () const ; \
164  void _maxiter (size_t miter) ; \
165  size_t _maxiter () const ; \
166  void _integrator(std::string const & type); \
167  std::string _integrator() const ; \
168  double _deltaMin() const; \
169  void _deltaMin(double h1); \
170  void _absTol(double); \
171  double _absTol()const; \
172  void _relTol(double); \
173  double _relTol()const; \
174  void _deltaMax(double h1); \
175  double _deltaMax()const; \
176  double _deltaMaxOverR() const; \
177  void _deltaMaxOverR(double t); \
178  std::vector<double> _initCoord()const; \
179  void _initCoord(std::vector<double> const&f); \
180  void _metric(SmartPointer<Metric::Generic>); \
181  SmartPointer<Metric::Generic> _metric() const;
182 
220 : protected Gyoto::Hook::Listener
221 {
222 
223  // Data :
224  // -----
225  public:
226  int stopcond;
227 
228  protected:
230  double* x0_;
231  double* x1_;
232  double* x2_;
233  double* x3_;
234  double* x0dot_;
235  double* x1dot_;
236  double* x2dot_;
237  double* x3dot_;
238  size_t x_size_;
239  size_t imin_;
240  size_t i0_;
241  size_t imax_;
242  bool adaptive_;
243 
250 
256  double delta_;
257 
258 
268  double tmin_;
269 
270  double * cst_;
271  size_t cst_n_;
272  int wait_pos_;
273  double * init_vel_;
274  size_t maxiter_ ;
275 
283  double delta_min_;
284 
292  double delta_max_;
293 
306 
314  double abstol_;
315 
323  double reltol_;
324 
325  // Constructors - Destructor
326  // -------------------------
327  public:
328  Worldline() ;
329 
330  Worldline(const Worldline& ) ;
331 
333 
346  Worldline(Worldline* orig, size_t i0, int dir, double step_max) ;
347 
348  virtual ~Worldline() ;
349 
350  size_t getImin() const;
351  size_t getImax() const;
352  size_t getI0() const;
353 
354  virtual double getMass() const = 0;
357 
358  void initCoord(std::vector<double> const&);
359  std::vector<double> initCoord() const;
360 
379  virtual void setInitCoord(const double coord[8], int dir = 0);
380 
388  virtual void setInitCoord(double const pos[4], double const vel[3], int dir=0);
389 
390  virtual void setPosition(double const pos[4]);
391  virtual void setVelocity(double const vel[3]);
392 
393  void reset() ;
394  void reInit() ;
395 
396  virtual std::string className() const ;
397  virtual std::string className_l() const ;
398 
409  void integrator(std::string const & type);
410 
414  std::string integrator() const ;
415 
419  double deltaMin() const;
420 
424  void deltaMin(double h1);
425 
429  double deltaMax() const;
430 
431 
432  void absTol(double);
433  double absTol()const;
434  void relTol(double);
435  double relTol()const;
436 
447  virtual double deltaMax(double const pos[8], double delta_max_external) const;
448 
452  void deltaMax(double h1);
453 
454  double deltaMaxOverR() const;
455  void deltaMaxOverR(double t);
456 
457  // Memory management
458  // -----------------
459  protected:
463  virtual void xAllocate();
464 
468  virtual void xAllocate(size_t size);
469 
481  virtual size_t xExpand(int dir);
482 
483 
491  virtual void xExpand(double * &x, int dir);
492 
493  // Mutators / assignment
494  // ---------------------
495  public:
497  void delta(const double delta);
498  void delta(double, const std::string &unit);
499  double delta() const ;
500  double delta(const std::string &unit) const ;
501  double tMin() const ;
502  double tMin(const std::string &unit) const ;
503  void tMin(double tlim);
504  void tMin(double, const std::string &unit);
505  void adaptive (bool mode) ;
506  bool adaptive () const ;
507  void secondary (bool sec) ;
508  bool secondary () const ;
509  void maxiter (size_t miter) ;
510  size_t maxiter () const ;
511 
516  double const * getCst() const ;
517 
519 
523  void setCst(double const * cst, size_t const ncsts) ;
524 
526 
534  const double coord[8],
535  const int dir) ;
536 
537  void getInitialCoord(double dest[8]) const;
538  void getCoord(size_t index, double dest[8]) const;
539  void getCartesianPos(size_t index, double dest[4]) const;
540 
541 
542  virtual void xStore(size_t ind, double const coord[8]) ;
543  virtual void xFill(double tlim) ;
544 
545 
546 
547  // Accessors
548  // ---------
549  public:
553  size_t get_nelements() const;
554 
558  void get_t(double *dest) const;
559 
560 
562 
578  void getCartesian(double const * const dates, size_t const n_dates,
579  double * const x, double * const y,
580  double * const z, double * const xprime=NULL,
581  double * const yprime=NULL, double * const zprime=NULL) ;
582 
586  void get_xyz(double* x, double *y, double *z) const;
587 
607  void getCoord(double const * const dates, size_t const n_dates,
608  double * const x1dest,
609  double * const x2dest, double * const x3dest,
610  double * const x0dot=NULL, double * const x1dot=NULL,
611  double * const x2dot=NULL, double * const x3dot=NULL) ;
612 
619  void getCoord(double *x0, double *x1, double *x2, double *x3) const ;
620 
629  void checkPhiTheta(double coord[8]) const;
630 
634  void getSkyPos(SmartPointer<Screen> screen, double *dalpha, double *ddellta, double *dD) const;
635 
639  void get_dot(double *x0dot, double *x1dot, double *x2dot, double *x3dot) const ;
640 
644  void get_prime(double *x1prime, double *x2prime, double *x3prime) const ;
645 
646  // Outputs
647  // -------
648  public:
649  //virtual void sauve(FILE *) const ; ///< Save in a file
650  void save_txyz(char * fichierxyz) const ;
651  void save_txyz(char* const filename, double const t1, double const mass_sun,
652  double const distance_kpc, std::string const unit, SmartPointer<Screen> sc = NULL);
653 
654  protected:
655  virtual void tell(Gyoto::Hook::Teller*);
656 
657  class IntegState {
658  public:
659  class Generic;
660  class Legacy;
661 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS
662  class Boost;
663 #endif
664  };
665 
666 
671 };
672 
673 #ifndef GYOTO_SWIGIMPORTED
674 
680  protected:
682 
687  double delta_;
688  bool adaptive_;
689  double norm_;
690  double normref_;
691 
695  Gyoto::SmartPointer<Gyoto::Metric::Generic> gg_;
696 
697  public:
703  Generic(Worldline *parent);
704 
708  virtual ~Generic();
709 
715  virtual Generic * clone(Worldline*newparent) const =0 ;
716 
723  virtual void init(Worldline * line, const double *coord, const double delta);
724 
733  virtual void init();
734 
741  virtual void checkNorm(double coord[8]);
742 
746  virtual std::string kind()=0;
747 
749 
753  virtual int nextStep(double *coord, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0;
754 
756 
765  virtual void doStep(double const coordin[8],
766  double step,
767  double coordout[8])=0;
768 };
769 
783 
784  private:
785  double coord_[8];
786 
787  public:
789 
790  Legacy(Worldline *parent);
791  Legacy * clone(Worldline*newparent) const ;
792  using Generic::init;
793  void init(Worldline * line, const double *coord, const double delta);
794  virtual std::string kind();
795 
796  virtual int nextStep(double *coord, double h1max=1e6);
797 
798  virtual void doStep(double const coordin[8],
799  double step,
800  double coordout[8]);
801 
802  virtual ~Legacy();
803 };
804 
805 #ifdef GYOTO_HAVE_BOOST_INTEGRATORS
806 
818  public:
822  enum Kind {runge_kutta_cash_karp54,
823  runge_kutta_fehlberg78,
824  runge_kutta_dopri5,
825  runge_kutta_cash_karp54_classic };
826  private:
829 
831  std::function<boost::numeric::odeint::controlled_step_result
832  (std::array<double,8>&, double&, double&)> try_step_;
833 
835  std::function<void(std::array<double,8>&, double)> do_step_;
836  public:
838 
843  Boost(Worldline* parent, std::string type);
844 
846 
851  Boost(Worldline* parent, Kind type);
852  Boost * clone(Worldline* newparent) const ;
853  virtual ~Boost();
854  virtual void init();
855  virtual void init(Worldline * line, const double *coord, const double delta);
856  virtual int nextStep(double *coord, double h1max=1e6);
857  virtual void doStep(double const coordin[8],
858  double step,
859  double coordout[8]);
860  virtual std::string kind();
861 
862 };
863 #endif
864 #endif
865 #endif
double delta_max_over_r_
Numerical tuning parameter.
Definition: GyotoWorldline.h:305
double * init_vel_
Hack in setParameters()
Definition: GyotoWorldline.h:273
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:832
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:219
SmartPointer< Gyoto::Metric::Generic > metric_
The Gyoto::Metric in this part of the universe.
Definition: GyotoWorldline.h:229
double delta_
Initial integrating step.
Definition: GyotoWorldline.h:256
void getInitialCoord(double dest[8]) const
Get initial coordinate.
void save_txyz(char *fichierxyz) const
Save in a file.
Tellers tell Listeners when they mutate.
Reference-counting pointers.
virtual std::string className_l() const
"worldline"
Boost integrator.
Definition: GyotoWorldline.h:816
virtual void xStore(size_t ind, double const coord[8])
Store coord at index ind.
void reInit()
Reset and recompute particle properties.
double tmin_
Time limit for the integration (geometrical units)
Definition: GyotoWorldline.h:268
Kind
Enum to represent the integrator flavour.
Definition: GyotoWorldline.h:822
double delta_min_
Minimum integration step for the adaptive integrator.
Definition: GyotoWorldline.h:283
size_t getImax() const
Get imax_.
virtual void checkNorm(double coord[8])
Check norm.
double * cst_
Worldline&#39;s csts of motion (if any)
Definition: GyotoWorldline.h:270
double relTol() const
Get reltol_.
double delta() const
Get delta_.
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.
void getCoord(size_t index, double dest[8]) const
Get coordinates corresponding to index.
bool adaptive_
Whether to use an adaptive step.
Definition: GyotoWorldline.h:688
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:695
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:240
void getCartesianPos(size_t index, double dest[4]) const
Get Cartesian expression of 4-position at index.
Home-brewed integrator.
Definition: GyotoWorldline.h:781
double tMin() const
Get tmin_.
double const * getCst() const
Returns the worldline&#39;s cst of motion (if any)
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.
SmartPointer< Worldline::IntegState::Generic > state_
An object to hold the integration state.
Definition: GyotoWorldline.h:670
size_t x_size_
Size of x0_, x1_... arrays.
Definition: GyotoWorldline.h:238
double * x0_
t or T
Definition: GyotoWorldline.h:230
bool adaptive_
Whether integration should use adaptive delta.
Definition: GyotoWorldline.h:242
size_t getImin() const
Get imin_.
size_t maxiter() const
Get maxiter_.
double * x2_
θ or y
Definition: GyotoWorldline.h:232
double * x3dot_
φdot or zdot
Definition: GyotoWorldline.h:237
double * x0dot_
tdot or Tdot
Definition: GyotoWorldline.h:234
Definition: GyotoWorldline.h:657
double norm_
Current norm of the 4-velocity.
Definition: GyotoWorldline.h:689
double reltol_
Absolute tolerance of the integrator.
Definition: GyotoWorldline.h:323
double normref_
Definition: GyotoWorldline.h:690
double deltaMaxOverR() const
Get delta_max_over_r_.
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.
void get_t(double *dest) const
Get computed dates.
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:78
bool secondary_
Experimental: choose 0 to compute only primary image.
Definition: GyotoWorldline.h:249
I might listen to a Teller.
Definition: GyotoHooks.h:64
void get_prime(double *x1prime, double *x2prime, double *x3prime) const
Get computed 3-velocities.
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:687
double * x3_
φ or z
Definition: GyotoWorldline.h:233
virtual std::string className() const
"Worldline"
size_t cst_n_
Number of constants of motion.
Definition: GyotoWorldline.h:271
size_t imin_
Minimum index for which x0_, x1_... have been computed.
Definition: GyotoWorldline.h:239
double * x1dot_
rdot or xdot
Definition: GyotoWorldline.h:235
double abstol_
Absolute tolerance of the integrator.
Definition: GyotoWorldline.h:314
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:678
virtual void setVelocity(double const vel[3])
Set initial 3-velocity.
Kind kind_
Integrator flavour.
Definition: GyotoWorldline.h:828
Description of the observer screen.
SmartPointer< Metric::Generic > metric() const
Get metric.
virtual void setInitCoord(const double coord[8], int dir=0)
Set Initial coordinate.
Worldline * line_
Worldline that we are integrating.
Definition: GyotoWorldline.h:686
double delta_max_
Maximum integration step for the adaptive integrator.
Definition: GyotoWorldline.h:292
virtual void xFill(double tlim)
Fill x0, x1... by integrating the Worldline from previously set inittial condition to time tlim...
virtual int nextStep(double *coord, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0
Make one step.
int wait_pos_
Hack in setParameters()
Definition: GyotoWorldline.h:272
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:236
Listen to me and I&#39;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:231
int stopcond
Whether and why integration is finished.
Definition: GyotoWorldline.h:226
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:274
Worldline()
Default constructor.
double absTol() const
Get abstol_.
size_t imax_
Maximum index for which x0_, x1_... have been computed.
Definition: GyotoWorldline.h:241