Gyoto
GyotoWorldline.h
Go to the documentation of this file.
1 
6 /*
7  Copyright 2011 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 HAVE_BOOST
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 
82 : protected Gyoto::Hook::Listener
83 {
84 
85  // Data :
86  // -----
87  public:
88  int stopcond;
89 
90  protected:
92  double* x0_;
93  double* x1_;
94  double* x2_;
95  double* x3_;
96  double* x0dot_;
97  double* x1dot_;
98  double* x2dot_;
99  double* x3dot_;
100  size_t x_size_;
101  size_t imin_;
102  size_t i0_;
103  size_t imax_;
104  bool adaptive_;
105 
112 
118  double delta_;
119 
120 
130  double tmin_;
131 
132  double * cst_;
133  size_t cst_n_;
134  int wait_pos_;
135  double * init_vel_;
136  size_t maxiter_ ;
137 
145  double delta_min_;
146 
154  double delta_max_;
155 
168 
176  double abstol_;
177 
185  double reltol_;
186 
187  // Constructors - Destructor
188  // -------------------------
189  public:
190  Worldline() ;
191 
192  Worldline(const Worldline& ) ;
193 
195 
208  Worldline(Worldline* orig, size_t i0, int dir, double step_max) ;
209 
210  virtual ~Worldline() ;
211 
212  size_t getImin() const;
213  size_t getImax() const;
214  size_t getI0() const;
215 
216  virtual double getMass() const = 0;
219  virtual void setInitCoord(const double coord[8], int dir = 0);
220 
228  virtual void setInitCoord(double pos[4], double vel[3], int dir=1);
229 
230  virtual void setPosition(double pos[4]);
231  virtual void setVelocity(double vel[3]);
232 
233  void reset() ;
234  void reInit() ;
235 
236  virtual std::string className() const ;
237  virtual std::string className_l() const ;
238 
249  void integrator(std::string type);
250 
254  std::string integrator() const ;
255 
259  double deltaMin() const;
260 
264  void deltaMin(double h1);
265 
269  double deltaMax() const;
270 
271 
272  void absTol(double);
273  double absTol()const;
274  void relTol(double);
275  double relTol()const;
276 
287  virtual double deltaMax(double const pos[8], double delta_max_external) const;
288 
292  void deltaMax(double h1);
293 
294  double deltaMaxOverR() const;
295  void deltaMaxOverR(double t);
296 
297  // Memory management
298  // -----------------
299  protected:
303  virtual void xAllocate();
304 
308  virtual void xAllocate(size_t size);
309 
321  virtual size_t xExpand(int dir);
322 
323 
331  virtual void xExpand(double * &x, int dir);
332 
333  // Mutators / assignment
334  // ---------------------
335  public:
337  void operator=(const Worldline&) ;
338  void delta(const double delta);
339  void delta(double, const std::string &unit);
340  double delta() const ;
341  double delta(const std::string &unit) const ;
342  double tMin() const ;
343  double tMin(const std::string &unit) const ;
344  void tMin(double tlim);
345  void tMin(double, const std::string &unit);
346  void adaptive (bool mode) ;
347  bool adaptive () const ;
348  void secondary (bool sec) ;
349  bool secondary () const ;
350  void maxiter (size_t miter) ;
351  size_t maxiter () const ;
352 
357  double const * getCst() const ;
358 
360 
364  void setCst(double const * cst, size_t const ncsts) ;
365 
367 
375  const double coord[8],
376  const int dir) ;
377 
378  void getInitialCoord(double dest[8]) const;
379  void getCoord(size_t index, double dest[8]) const;
380  void getCartesianPos(size_t index, double dest[4]) const;
381 
382 
383  virtual void xStore(size_t ind, double coord[8]) ;
384  virtual void xFill(double tlim) ;
385 
408  virtual int setParameter(std::string name,
409  std::string content,
410  std::string unit) ;
411 
412 #ifdef GYOTO_USE_XERCES
413 
418  virtual void setParameters(FactoryMessenger *fmp) ;
424  virtual void fillElement(FactoryMessenger *fmp) const ;
426 #endif
427 
428  // Accessors
429  // ---------
430  public:
434  size_t get_nelements() const;
435 
439  void get_t(double *dest) const;
440 
441 
443 
459  void getCartesian(double const * const dates, size_t const n_dates,
460  double * const x, double * const y,
461  double * const z, double * const xprime=NULL,
462  double * const yprime=NULL, double * const zprime=NULL) ;
463 
467  void get_xyz(double* x, double *y, double *z) const;
468 
488  void getCoord(double const * const dates, size_t const n_dates,
489  double * const x1dest,
490  double * const x2dest, double * const x3dest,
491  double * const x0dot=NULL, double * const x1dot=NULL,
492  double * const x2dot=NULL, double * const x3dot=NULL) ;
493 
500  void getCoord(double *x0, double *x1, double *x2, double *x3) const ;
501 
510  void checkPhiTheta(double coord[8]) const;
511 
515  void getSkyPos(SmartPointer<Screen> screen, double *dalpha, double *ddellta, double *dD) const;
516 
520  void get_dot(double *x0dot, double *x1dot, double *x2dot, double *x3dot) const ;
521 
525  void get_prime(double *x1prime, double *x2prime, double *x3prime) const ;
526 
527  // Outputs
528  // -------
529  public:
530  //virtual void sauve(FILE *) const ; ///< Save in a file
531  void save_txyz(char * fichierxyz) const ;
532  void save_txyz(char* const filename, double const t1, double const mass_sun,
533  double const distance_kpc, std::string const unit, SmartPointer<Screen> sc = NULL);
534 
536  friend std::ostream& operator<<(std::ostream& , const Worldline& ) ;
537 
538  protected:
539  virtual void tell(Gyoto::Hook::Teller*);
540 
541  class IntegState {
542  public:
543  class Generic;
544  class Legacy;
545 #ifdef HAVE_BOOST
546  class Boost;
547 #endif
548  };
549 
550 
555 };
556 
557 
564  protected:
566 
571  double delta_;
572  bool adaptive_;
573  double norm_;
574  double normref_;
575 
579  Gyoto::SmartPointer<Gyoto::Metric::Generic> gg_;
580 
581  public:
587  Generic(Worldline *parent);
588 
592  virtual ~Generic();
593 
599  virtual Generic * clone(Worldline*newparent) const =0 ;
600 
607  virtual void init(Worldline * line, const double *coord, const double delta);
608 
617  virtual void init();
618 
625  virtual void checkNorm(double coord[8]);
626 
630  virtual std::string kind()=0;
631 
633 
637  virtual int nextStep(double *coord, double h1max=GYOTO_DEFAULT_DELTA_MAX)=0;
638 
639 };
640 
654 
655  private:
656  double coord_[8];
657 
658  public:
660 
661  Legacy(Worldline *parent);
662  Legacy * clone(Worldline*newparent) const ;
663  using Generic::init;
664  void init(Worldline * line, const double *coord, const double delta);
665  virtual std::string kind();
666 
668 
672  virtual int nextStep(double *coord, double h1max=1e6);
673 
674  virtual ~Legacy();
675 };
676 
677 #ifdef HAVE_BOOST
678 
690  public:
694  enum Kind {runge_kutta_cash_karp54,
695  runge_kutta_fehlberg78,
696  runge_kutta_dopri5,
697  runge_kutta_cash_karp54_classic };
698  private:
701 
703  std::function<boost::numeric::odeint::controlled_step_result
704  (std::array<double,8>&, double&, double&)> try_step_;
705 
707  std::function<void(std::array<double,8>&, double)> do_step_;
708  public:
710 
715  Boost(Worldline* parent, std::string type);
716 
718 
723  Boost(Worldline* parent, Kind type);
724  Boost * clone(Worldline* newparent) const ;
725  virtual ~Boost();
726  virtual void init();
727  virtual void init(Worldline * line, const double *coord, const double delta);
728  virtual int nextStep(double *coord, double h1max=1e6);
729  virtual std::string kind();
730 
731 };
732 #endif
733 
734 #endif
double delta_max_over_r_
Numerical tuning parameter.
Definition: GyotoWorldline.h:167
size_t get_nelements() const
Get number of computed dates.
double * init_vel_
Hack in setParameters()
Definition: GyotoWorldline.h:135
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:704
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:81
SmartPointer< Gyoto::Metric::Generic > metric_
The Gyoto::Metric in this part of the universe.
Definition: GyotoWorldline.h:91
double delta_
Initial integrating step.
Definition: GyotoWorldline.h:118
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:688
void get_t(double *dest) const
Get computed dates.
virtual int setParameter(std::string name, std::string content, std::string unit)
Set parameter by name.
void reInit()
Reset and recompute particle properties.
double tmin_
Time limit for the integration (geometrical units)
Definition: GyotoWorldline.h:130
double tMin() const
Get tmin_.
Kind
Enum to represent the integrator flavour.
Definition: GyotoWorldline.h:694
double delta_min_
Minimum integration step for the adaptive integrator.
Definition: GyotoWorldline.h:145
virtual void checkNorm(double coord[8])
Check norm.
double * cst_
Worldline's csts of motion (if any)
Definition: GyotoWorldline.h:132
friend std::ostream & operator<<(std::ostream &, const Worldline &)
Display.
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:572
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:579
virtual std::string className() const
"Worldline"
size_t i0_
Index of initial condition in array.
Definition: GyotoWorldline.h:102
Home-brewed integrator.
Definition: GyotoWorldline.h:652
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:554
size_t x_size_
Size of x0_, x1_... arrays.
Definition: GyotoWorldline.h:100
double * x0_
t or T
Definition: GyotoWorldline.h:92
double absTol() const
Get abstol_.
bool adaptive_
Whether integration should use adaptive delta.
Definition: GyotoWorldline.h:104
double relTol() const
Get reltol_.
double * x2_
θ or y
Definition: GyotoWorldline.h:94
double * x3dot_
φdot or zdot
Definition: GyotoWorldline.h:99
double * x0dot_
tdot or Tdot
Definition: GyotoWorldline.h:96
Definition: GyotoWorldline.h:541
double norm_
Current norm of the 4-velocity.
Definition: GyotoWorldline.h:573
double reltol_
Absolute tolerance of the integrator.
Definition: GyotoWorldline.h:185
double normref_
Definition: GyotoWorldline.h:574
Namespace for the Gyoto library.
Definition: GyotoAstrobj.h:42
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_.
virtual void setParameters(FactoryMessenger *fmp)
Process XML entity Uses wait_pos_ and init_vel_ to make sure setVelocity() is called after setPositio...
Can be pointed to by a SmartPointer.
Definition: GyotoSmartPointer.h:79
bool secondary_
Experimental: choose 0 to compute only primary image.
Definition: GyotoWorldline.h:111
I might listen to a Teller.
Definition: GyotoHooks.h:64
virtual std::string kind()=0
Return the integrator kind.
virtual ~Worldline()
Destructor.
double delta_
Integration step (current in case of adaptive_).
Definition: GyotoWorldline.h:571
double * x3_
φ or z
Definition: GyotoWorldline.h:95
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:133
size_t imin_
Minimum index for which x0_, x1_... have been computed.
Definition: GyotoWorldline.h:101
double * x1dot_
rdot or xdot
Definition: GyotoWorldline.h:97
double abstol_
Absolute tolerance of the integrator.
Definition: GyotoWorldline.h:176
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:562
size_t getImax() const
Get imax_.
Kind kind_
Integrator flavour.
Definition: GyotoWorldline.h:700
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:570
void operator=(const Worldline &)
Assignment to another Worldline.
double delta_max_
Maximum integration step for the adaptive integrator.
Definition: GyotoWorldline.h:154
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:134
bool adaptive() const
Get adaptive_.
virtual void xAllocate()
Allocate x0, x1 etc. with default size.
double * x2dot_
θdot or ydot
Definition: GyotoWorldline.h:98
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 fillElement(FactoryMessenger *fmp) const
XML output.
virtual void setPosition(double pos[4])
Set initial 4-position.
double * x1_
r or x
Definition: GyotoWorldline.h:93
int stopcond
Whether and why integration is finished.
Definition: GyotoWorldline.h:88
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:136
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:103