Karim Sabri / gammalib

Commit 571bc6fee609aff595e2ed23106dc2e22b803075

Authored by Jürgen Knödlseder
1 parent 659bf459

Do not longer use GModelTemporalPhaseCurve for phase evaluation and implement ha… 

…ndling of different time systems
Inline Side-by-side
Showing 3 changed files with 142 additions and 49 deletions
include/GPulsarEphemeris.hpp
  View file @571bc6f
... ... @@ -32,7 +32,6 @@
32 32 #include "GBase.hpp"
33 33 #include "GTime.hpp"
34 34 #include "GSkyDir.hpp"
35   -#include "GModelTemporalPhaseCurve.hpp"
36 35  
37 36 /* __ Forward declarations _______________________________________________ */
38 37  
... ... @@ -72,6 +71,8 @@ public:
72 71 void tstart(const GTime& tstart);
73 72 const GTime& tstop(void) const;
74 73 void tstop(const GTime& tstop);
  74 + const std::string& timesys(void) const;
  75 + void timesys(const std::string& timesys);
75 76 GTime t0(void) const;
76 77 void t0(const GTime& t0);
77 78 double phase(void) const;
... ... @@ -82,7 +83,8 @@ public:
82 83 void f1(const double& f1);
83 84 double f2(void) const;
84 85 void f2(const double& f2);
85   - double phase(const GTime& time) const;
  86 + double phase(const GTime& time,
  87 + const std::string& timesys) const;
86 88  
87 89 protected:
88 90 // Protected methods
... ... @@ -91,11 +93,16 @@ protected:
91 93 void free_members(void);
92 94  
93 95 // Protected members
94   - std::string m_name; //!< Pulsar name
95   - GTime m_tstart; //!< Validity start time
96   - GTime m_tstop; //!< Validity stop time
97   - GSkyDir m_dir; //!< Pulsar sky direction
98   - GModelTemporalPhaseCurve m_phase_curve; //!< Pulsar phase curve
  96 + std::string m_name; //!< Pulsar name
  97 + GTime m_tstart; //!< Validity start time
  98 + GTime m_tstop; //!< Validity stop time
  99 + GSkyDir m_dir; //!< Pulsar sky direction
  100 + std::string m_timesys; //!< Time system of pulsar ephemeris
  101 + GTime m_t0; //!< Reference epoch of pulsar ephemeris
  102 + double m_phase; //!< Pulse phase
  103 + double m_f0; //!< Pulsar frequency (Hz)
  104 + double m_f1; //!< Pulsar frequency derivative (s^-2)
  105 + double m_f2; //!< Pulsar second frequency derivative (s^-3)
99 106 };
100 107  
101 108  
... ... @@ -216,6 +223,48 @@ void GPulsarEphemeris::tstop(const GTime&amp; tstop)
216 223  
217 224  
218 225 /***********************************************************************//**
  226 + * @brief Return class name
  227 + *
  228 + * @return String containing the class name ("GPulsarEphemeris").
  229 + ***************************************************************************/
  230 +inline
  231 +std::string GPulsarEphemeris::classname(void) const
  232 +{
  233 + return ("GPulsarEphemeris");
  234 +}
  235 +
  236 +
  237 +/***********************************************************************//**
  238 + * @brief Returns pulsar ephemeris time system
  239 + *
  240 + * @return Pulsar ephemeris time system.
  241 + *
  242 + * The pulsar ephemeris time system is one of "TT", "TAI" or "UTC".
  243 + ***************************************************************************/
  244 +inline
  245 +const std::string& GPulsarEphemeris::timesys(void) const
  246 +{
  247 + // Return
  248 + return (m_timesys);
  249 +}
  250 +
  251 +
  252 +/***********************************************************************//**
  253 + * @brief Set pulsar ephemeris time system
  254 + *
  255 + * @param[in] timesys Pulsar ephemeris time system.
  256 + *
  257 + * The pulsar ephemeris time system is one of "TT", "TAI" or "UTC".
  258 + ***************************************************************************/
  259 +inline
  260 +void GPulsarEphemeris::timesys(const std::string& timesys)
  261 +{
  262 + m_timesys = timesys;
  263 + return;
  264 +}
  265 +
  266 +
  267 +/***********************************************************************//**
219 268 * @brief Returns reference epoch of pulsar ephemeris
220 269 *
221 270 * @return Reference epoch of pulsar ephemeris.
... ... @@ -223,7 +272,7 @@ void GPulsarEphemeris::tstop(const GTime&amp; tstop)
223 272 inline
224 273 GTime GPulsarEphemeris::t0(void) const
225 274 {
226   - return (m_phase_curve.mjd());
  275 + return (m_t0);
227 276 }
228 277  
229 278  
... ... @@ -235,7 +284,7 @@ GTime GPulsarEphemeris::t0(void) const
235 284 inline
236 285 void GPulsarEphemeris::t0(const GTime& t0)
237 286 {
238   - m_phase_curve.mjd(t0);
  287 + m_t0 = t0;
239 288 return;
240 289 }
241 290  
... ... @@ -248,7 +297,7 @@ void GPulsarEphemeris::t0(const GTime&amp; t0)
248 297 inline
249 298 double GPulsarEphemeris::phase(void) const
250 299 {
251   - return (m_phase_curve.phase());
  300 + return (m_phase);
252 301 }
253 302  
254 303  
... ... @@ -260,7 +309,7 @@ double GPulsarEphemeris::phase(void) const
260 309 inline
261 310 void GPulsarEphemeris::phase(const double& phase)
262 311 {
263   - m_phase_curve.phase(phase);
  312 + m_phase = phase;
264 313 return;
265 314 }
266 315  
... ... @@ -273,7 +322,7 @@ void GPulsarEphemeris::phase(const double&amp; phase)
273 322 inline
274 323 double GPulsarEphemeris::f0(void) const
275 324 {
276   - return (m_phase_curve.f0());
  325 + return (m_f0);
277 326 }
278 327  
279 328  
... ... @@ -285,7 +334,7 @@ double GPulsarEphemeris::f0(void) const
285 334 inline
286 335 void GPulsarEphemeris::f0(const double& f0)
287 336 {
288   - m_phase_curve.f0(f0);
  337 + m_f0 = f0;
289 338 return;
290 339 }
291 340  
... ... @@ -298,7 +347,7 @@ void GPulsarEphemeris::f0(const double&amp; f0)
298 347 inline
299 348 double GPulsarEphemeris::f1(void) const
300 349 {
301   - return (m_phase_curve.f1());
  350 + return (m_f1);
302 351 }
303 352  
304 353  
... ... @@ -310,7 +359,7 @@ double GPulsarEphemeris::f1(void) const
310 359 inline
311 360 void GPulsarEphemeris::f1(const double& f1)
312 361 {
313   - m_phase_curve.f1(f1);
  362 + m_f1 = f1;
314 363 return;
315 364 }
316 365  
... ... @@ -323,7 +372,7 @@ void GPulsarEphemeris::f1(const double&amp; f1)
323 372 inline
324 373 double GPulsarEphemeris::f2(void) const
325 374 {
326   - return (m_phase_curve.f2());
  375 + return (m_f2);
327 376 }
328 377  
329 378  
... ... @@ -335,20 +384,8 @@ double GPulsarEphemeris::f2(void) const
335 384 inline
336 385 void GPulsarEphemeris::f2(const double& f2)
337 386 {
338   - m_phase_curve.f2(f2);
  387 + m_f2 = f2;
339 388 return;
340 389 }
341 390  
342   -
343   -/***********************************************************************//**
344   - * @brief Returns pulsar phase
345   - *
346   - * @return Pulsar phase.
347   - ***************************************************************************/
348   -inline
349   -double GPulsarEphemeris::phase(const GTime& time) const
350   -{
351   - return (m_phase_curve.phase(time));
352   -}
353   -
354 391 #endif /* GPULSAREPHEMERIS_HPP */
... ...
pyext/GPulsarEphemeris.i
  View file @571bc6f
... ... @@ -56,6 +56,8 @@ public:
56 56 void tstart(const GTime& tstart);
57 57 const GTime& tstop(void) const;
58 58 void tstop(const GTime& tstop);
  59 + const std::string& timesys(void) const;
  60 + void timesys(const std::string& timesys);
59 61 GTime t0(void) const;
60 62 void t0(const GTime& t0);
61 63 double phase(void) const;
... ... @@ -66,7 +68,8 @@ public:
66 68 void f1(const double& f1);
67 69 double f2(void) const;
68 70 void f2(const double& f2);
69   - double phase(const GTime& time) const;
  71 + double phase(const GTime& time,
  72 + const std::string& timesys) const;
70 73 };
71 74  
72 75  
... ...
src/obs/GPulsarEphemeris.cpp
  View file @571bc6f
... ... @@ -29,6 +29,7 @@
29 29 #include <config.h>
30 30 #endif
31 31 #include "GPulsarEphemeris.hpp"
  32 +#include "GTools.hpp"
32 33  
33 34 /* __ Method name definitions ____________________________________________ */
34 35  
... ... @@ -157,6 +158,53 @@ GPulsarEphemeris* GPulsarEphemeris::clone(void) const
157 158  
158 159  
159 160 /***********************************************************************//**
  161 + * @brief Returns pulsar phase
  162 + *
  163 + * @param[in] time Time.
  164 + * @param[in] timesys Time system for evaluation of time.
  165 + * @return Pulsar phase.
  166 + *
  167 + * Computes
  168 + *
  169 + * \f[
  170 + * \Phi(t) = \Phi_0 + f(t-t_0) + \frac{1}{2}\dot{f} (t-t_0)^2 +
  171 + * \frac{1}{6}\ddot{f} (t-t_0)^3
  172 + * \f]
  173 + *
  174 + * where
  175 + * \f$t\f$ is the specified @p time in seconds,
  176 + * \f$t_0\f$ is a reference time in seconds,
  177 + * \f$\Phi_0\f$ is the phase at the reference time,
  178 + * \f$f\f$ is the variation frequency at the reference time,
  179 + * \f$\dot{f}\f$ is the first derivative of the variation frequency at the
  180 + * reference time, and
  181 + * \f$\dot{f}\f$ is the second derivative of the variation frequency at the
  182 + * reference time.
  183 + *
  184 + * The phase \f$\Phi(t)\f$ is in the interval [0,1].
  185 + ***************************************************************************/
  186 +double GPulsarEphemeris::phase(const GTime& time,
  187 + const std::string& timesys) const
  188 +{
  189 + // Set constants
  190 + const double c1 = 0.5;
  191 + const double c2 = 1.0 / 6.0;
  192 +
  193 + // Compute time since reference time in seconds
  194 + double dt = time.secs(timesys) - m_t0.secs(m_timesys);
  195 +
  196 + // Computes phase
  197 + double phase = m_phase + dt * (m_f0 + dt * (c1 * m_f1 + c2 * m_f2 * dt));
  198 +
  199 + // Put phase into interval [0,1]
  200 + phase -= floor(phase);
  201 +
  202 + // Return phase
  203 + return phase;
  204 +}
  205 +
  206 +
  207 +/***********************************************************************//**
160 208 * @brief Print Pulsar ephemeris
161 209 *
162 210 * @param[in] chatter Chattiness.
... ... @@ -187,9 +235,13 @@ std::string GPulsarEphemeris::print(const GChatter&amp; chatter) const
187 235 double f0 = this->f0();
188 236  
189 237 // Append phase information
190   - result.append("\n"+gammalib::parformat("Time of phase 0"));
  238 + result.append("\n"+gammalib::parformat("Time system of ephemeris"));
  239 + result.append(m_timesys);
  240 + result.append("\n"+gammalib::parformat("Reference epoch"));
191 241 result.append("MJD ");
192   - result.append(gammalib::str(t0().mjd()));
  242 + result.append(gammalib::str(t0().mjd(m_timesys)));
  243 + result.append("\n"+gammalib::parformat("Phase"));
  244 + result.append(gammalib::str(phase()));
193 245 result.append("\n"+gammalib::parformat("Frequency"));
194 246 result.append(gammalib::str(f0));
195 247 result.append(" Hz");
... ... @@ -217,9 +269,9 @@ std::string GPulsarEphemeris::print(const GChatter&amp; chatter) const
217 269  
218 270 // Append validity information
219 271 result.append("\n"+gammalib::parformat("Validity MJD range"));
220   - result.append(gammalib::str(tstart().mjd()));
  272 + result.append(gammalib::str(tstart().mjd(m_timesys)));
221 273 result.append(" - ");
222   - result.append(gammalib::str(tstop().mjd()));
  274 + result.append(gammalib::str(tstop().mjd(m_timesys)));
223 275 result.append("\n"+gammalib::parformat("Validity UTC range"));
224 276 result.append(tstart().utc());
225 277 result.append(" - ");
... ... @@ -248,16 +300,12 @@ void GPulsarEphemeris::init_members(void)
248 300 m_tstart.clear();
249 301 m_tstop.clear();
250 302 m_dir.clear();
251   - m_phase_curve.clear();
252   -
253   - // Remove range from phase curve parameters
254   - m_phase_curve["MJD"].remove_range();
255   - m_phase_curve["F0"].remove_range();
256   - m_phase_curve["F1"].remove_range();
257   - m_phase_curve["F2"].remove_range();
258   -
259   - // Make sure that reference phase is zero
260   - m_phase_curve["Phase"].value(0.0);
  303 + m_t0.clear();
  304 + m_timesys = "TT";
  305 + m_phase = 0.0;
  306 + m_f0 = 0.0;
  307 + m_f1 = 0.0;
  308 + m_f2 = 0.0;
261 309  
262 310 // Return
263 311 return;
... ... @@ -272,11 +320,16 @@ void GPulsarEphemeris::init_members(void)
272 320 void GPulsarEphemeris::copy_members(const GPulsarEphemeris& ephemeris)
273 321 {
274 322 // Copy members
275   - m_name = ephemeris.m_name;
276   - m_tstart = ephemeris.m_tstart;
277   - m_tstop = ephemeris.m_tstop;
278   - m_dir = ephemeris.m_dir;
279   - m_phase_curve = ephemeris.m_phase_curve;
  323 + m_name = ephemeris.m_name;
  324 + m_tstart = ephemeris.m_tstart;
  325 + m_tstop = ephemeris.m_tstop;
  326 + m_dir = ephemeris.m_dir;
  327 + m_timesys = ephemeris.m_timesys;
  328 + m_t0 = ephemeris.m_t0;
  329 + m_phase = ephemeris.m_phase;
  330 + m_f0 = ephemeris.m_f0;
  331 + m_f1 = ephemeris.m_f1;
  332 + m_f2 = ephemeris.m_f2;
280 333  
281 334 // Return
282 335 return;
... ...