GammaLib  2.0.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelTemporalPhaseCurve.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelTemporalPhaseCurve.cpp - Temporal phase curve model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-2018 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GModelTemporalPhaseCurve.cpp
23  * @brief Temporal phase curve model class interface implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GTools.hpp"
33 #include "GException.hpp"
36 #include "GFits.hpp"
37 #include "GFitsTable.hpp"
38 #include "GFitsTableCol.hpp"
39 
40 /* __ Constants __________________________________________________________ */
41 
42 /* __ Globals ____________________________________________________________ */
44 const GModelTemporalRegistry g_temporal_phase_registry(&g_temporal_phase_seed);
45 
46 /* __ Method name definitions ____________________________________________ */
47 #define G_READ "GModelTemporalPhaseCurve::read(GXmlElement&)"
48 #define G_WRITE "GModelTemporalPhaseCurve::write(GXmlElement&)"
49 #define G_LOAD_NODES "GModelTemporalPhaseCurve::load_nodes(GFilename&)"
50 #define G_NORMALISE_NODES "GModelTemporalPhaseCurve::normalize_nodes()"
51 
52 /* __ Macros _____________________________________________________________ */
53 
54 /* __ Coding definitions _________________________________________________ */
55 
56 /* __ Debug definitions __________________________________________________ */
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void constructor
67  ***************************************************************************/
69 {
70  // Initialise members
71  init_members();
72 
73  // Return
74  return;
75 }
76 
77 
78 /***********************************************************************//**
79  * @brief File constructor
80  *
81  * @param[in] filename File name of phase curve nodes.
82  * @param[in] mjd Reference time.
83  * @param[in] phase Phase at reference time.
84  * @param[in] f0 Frequency at reference time (Hz).
85  * @param[in] f1 First frequency derivative at reference time.
86  * @param[in] f2 Second frequency derivative at reference time.
87  * @param[in] norm Normalization factor.
88  * @param[in] normalize Normalize phase curve?
89  *
90  * Constructs phase curve from a list of nodes that is found in the specified
91  * FITS file, a reference time and phase information at the reference time.
92  * See the load_nodes() method for more information about the expected
93  * structure of the file.
94  ***************************************************************************/
96  const GTime& mjd,
97  const double& phase,
98  const double& f0,
99  const double& f1,
100  const double& f2,
101  const double& norm,
102  const bool& normalize) :
104 {
105  // Initialise members
106  init_members();
107 
108  // Set parameters
109  this->mjd(mjd);
110  this->phase(phase);
111  this->f0(f0);
112  this->f1(f1);
113  this->f2(f2);
114  this->norm(norm);
115 
116  // Set normalization flag
118 
119  // Load nodes
120  load_nodes(filename);
121 
122  // Return
123  return;
124 }
125 
126 
127 /***********************************************************************//**
128  * @brief XML constructor
129  *
130  * @param[in] xml XML element.
131  *
132  * Constructs phase curve by extracting information from an XML element. See
133  * the read() method for more information about the expected structure of the
134  * XML element.
135  ***************************************************************************/
138 {
139  // Initialise members
140  init_members();
141 
142  // Read information from XML element
143  read(xml);
144 
145  // Return
146  return;
147 }
148 
149 
150 /***********************************************************************//**
151  * @brief Copy constructor
152  *
153  * @param[in] model Phase curve.
154  ***************************************************************************/
156  GModelTemporal(model)
157 {
158  // Initialise members
159  init_members();
160 
161  // Copy members
162  copy_members(model);
163 
164  // Return
165  return;
166 }
167 
168 
169 /***********************************************************************//**
170  * @brief Destructor
171  ***************************************************************************/
173 {
174  // Free members
175  free_members();
176 
177  // Return
178  return;
179 }
180 
181 
182 /*==========================================================================
183  = =
184  = Operators =
185  = =
186  ==========================================================================*/
187 
188 /***********************************************************************//**
189  * @brief Assignment operator
190  *
191  * @param[in] model Phase curve.
192  * @return Phase curve.
193  ***************************************************************************/
195 {
196  // Execute only if object is not identical
197  if (this != &model) {
198 
199  // Copy base class members
200  this->GModelTemporal::operator=(model);
201 
202  // Free members
203  free_members();
204 
205  // Initialise members
206  init_members();
207 
208  // Copy members
209  copy_members(model);
210 
211  } // endif: object was not identical
212 
213  // Return
214  return *this;
215 }
216 
217 
218 /*==========================================================================
219  = =
220  = Public methods =
221  = =
222  ==========================================================================*/
223 
224 /***********************************************************************//**
225  * @brief Clear phase curve
226  ***************************************************************************/
228 {
229  // Free class members (base and derived classes, derived class first)
230  free_members();
232 
233  // Initialise members
235  init_members();
236 
237  // Return
238  return;
239 }
240 
241 
242 /***********************************************************************//**
243  * @brief Clone phase curve
244  *
245  * @return Pointer to deep copy of phase curve.
246  ***************************************************************************/
248 {
249  // Clone constant temporal model
250  return new GModelTemporalPhaseCurve(*this);
251 }
252 
253 
254 /***********************************************************************//**
255  * @brief Evaluate function
256  *
257  * @param[in] srcTime True photon arrival time.
258  * @param[in] gradients Compute gradients?
259  * @return Value of phase curve.
260  *
261  * Computes
262  *
263  * \f[
264  * S_{\rm t}(t) = r(\Phi(t)) \times {\tt m\_norm}
265  * \f]
266  *
267  * where
268  *
269  * \f$r(\Phi(t))\f$ is the phase dependent rate, defined by linear
270  * interpolation between the nodes in a FITS file, \f$\Phi(t)\f$ is the time
271  * dependent phase, and \f${\tt m\_norm}\f$ is a normalisation constant.
272  ***************************************************************************/
273 double GModelTemporalPhaseCurve::eval(const GTime& srcTime,
274  const bool& gradients) const
275 {
276  // Get phase
277  double phase = this->phase(srcTime);
278 
279  // Interpolate phase curve nodes
280  double func = m_nodes.interpolate(phase, m_values);
281 
282  // Compute function value
283  double value = m_norm.value() * func;
284 
285  // Optionally compute gradients
286  if (gradients) {
287 
288  // Compute partial derivatives of the parameter values
289  double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0;
290 
291  // Set gradients
292  m_norm.factor_gradient(g_norm);
293 
294  } // endif: gradient computation was requested
295 
296  // Return
297  return value;
298 }
299 
300 
301 /***********************************************************************//**
302  * @brief Return vector of random event times
303  *
304  * @param[in] rate Mean event rate (events per second).
305  * @param[in] tmin Minimum event time.
306  * @param[in] tmax Maximum event time.
307  * @param[in,out] ran Random number generator.
308  *
309  * Returns a vector of random event times between @p tmin and @p tmax for
310  * the phase curve.
311  ***************************************************************************/
312 GTimes GModelTemporalPhaseCurve::mc(const double& rate, const GTime& tmin,
313  const GTime& tmax, GRan& ran) const
314 {
315  // Allocates empty vector of times
316  GTimes times;
317 
318  // Compute event rate (in events per seconds)
319  double lambda = rate * norm() * m_scale;
320 
321  // Initialise start and stop times in seconds
322  double time = tmin.secs();
323  double tstop = tmax.secs();
324 
325  // Generate events until maximum event time is exceeded
326  while (time <= tstop) {
327 
328  // Simulate next event time
329  time += ran.exp(lambda);
330 
331  // Break if time is beyond the stop time
332  if (time > tstop) {
333  break;
334  }
335 
336  // Set time object
337  GTime srcTime;
338  srcTime.secs(time);
339 
340  // Get value for current time
341  double value = eval(srcTime);
342 
343  // Get uniform random number
344  double uniform = ran.uniform() * m_scale;
345 
346  // If the random number is not larger than the phase curve value
347  // then accept the time
348  if (uniform <= value) {
349  times.append(srcTime);
350  }
351 
352  } // endwhile: loop until stop time is reached
353 
354  // Return vector of times
355  return times;
356 }
357 
358 
359 /***********************************************************************//**
360  * @brief Read model from XML element
361  *
362  * @param[in] xml XML element.
363  *
364  * Reads the temporal information from an XML element. The XML element should
365  * have the format
366  *
367  * <temporalModel type="PhaseCurve" file="phase.fits">
368  * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
369  * <parameter name="MJD" scale="1" value="51544.5" min="0.0" max="1e10" free="0"/>
370  * <parameter name="Phase" scale="1" value="0.0" min="0.0" max="1e10" free="1"/>
371  * <parameter name="F0" scale="1" value="1.0" min="0.0" max="1e10" free="1"/>
372  * <parameter name="F1" scale="1" value="0.1" min="0.0" max="1e10" free="1"/>
373  * <parameter name="F2" scale="1" value="0.01" min="0.0" max="1e10" free="1"/>
374  * </temporalModel>
375  ***************************************************************************/
377 {
378  // Get parameter pointers
385 
386  // Read parameters
387  m_norm.read(*norm);
388  m_mjd.read(*mjd);
389  m_phase.read(*phase);
390  m_f0.read(*f0);
391  m_f1.read(*f1);
392  m_f2.read(*f2);
393 
394  // Get optional normalization attribute
395  m_normalize = true;
396  m_has_normalize = false;
397  if (xml.has_attribute("normalize")) {
398  m_has_normalize = true;
399  std::string arg = xml.attribute("normalize");
400  if (arg == "0" || gammalib::tolower(arg) == "false") {
401  m_normalize = false;
402  }
403  }
404 
405  // Load nodes from file
407 
408  // Return
409  return;
410 }
411 
412 
413 /***********************************************************************//**
414  * @brief Write model into XML element
415  *
416  * @param[in] xml XML element.
417  *
418  * @exception GExpection::invalid_value
419  * Invalid XML format encountered.
420  *
421  * Writes the temporal information into an XML element. The XML element will
422  * have the format
423  *
424  * <temporalModel type="PhaseCurve" file="phase.fits">
425  * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
426  * <parameter name="MJD" scale="1" value="51544.5" min="0.0" max="1e10" free="0"/>
427  * <parameter name="Phase" scale="1" value="0.0" min="0.0" max="1e10" free="1"/>
428  * <parameter name="F0" scale="1" value="1.0" min="0.0" max="1e10" free="1"/>
429  * <parameter name="F1" scale="1" value="0.1" min="0.0" max="1e10" free="1"/>
430  * <parameter name="F2" scale="1" value="0.01" min="0.0" max="1e10" free="1"/>
431  * </temporalModel>
432  ***************************************************************************/
434 {
435  // Set model type
436  if (xml.attribute("type") == "") {
437  xml.attribute("type", type());
438  }
439 
440  // Verify model type
441  if (xml.attribute("type") != type()) {
442  std::string msg = "Temporal model of type "+xml.attribute("type")+
443  " encountered while \""+type()+"\" was expected.";
445  }
446 
447  // Get XML parameters
454 
455  // Write parameters
456  m_norm.write(*norm);
457  m_mjd.write(*mjd);
458  m_phase.write(*phase);
459  m_f0.write(*f0);
460  m_f1.write(*f1);
461  m_f2.write(*f2);
462 
463  // Set file attribute
465 
466  // Set optional normalization attribute
467  if (m_has_normalize || !m_normalize) {
468  if (m_normalize) {
469  xml.attribute("normalize", "1");
470  }
471  else {
472  xml.attribute("normalize", "0");
473  }
474  }
475 
476  // Return
477  return;
478 }
479 
480 
481 /***********************************************************************//**
482  * @brief Computes the phase for a given time
483  *
484  * @param[in] time Time.
485  *
486  * Computes
487  *
488  * \f[
489  * \Phi(t) = \Phi_0 + f(t-t_0) + \frac{1}{2}\dot{f} (t-t_0)^2 +
490  * \frac{1}{6}\ddot{f} (t-t_0)^3
491  * \f]
492  *
493  * where
494  * \f$t_0\f$ is a reference time,
495  * \f$\Phi_0\f$ is the phase at the reference time,
496  * \f$f\f$ is the variation frequency at the reference time,
497  * \f$\dot{f}\f$ is the first derivative of the variation frequency at the
498  * reference time, and
499  * \f$\dot{f}\f$ is the second derivative of the variation frequency at the
500  * reference time.
501  *
502  * The phase \f$\Phi(t)\f$ is in the interval [0,1].
503  ***************************************************************************/
504 double GModelTemporalPhaseCurve::phase(const GTime& time) const
505 {
506  // Set constants
507  const double c1 = 0.5;
508  const double c2 = 1.0 / 6.0;
509 
510  // Compute time since reference time in seconds
511  double t = time - mjd();
512 
513  // Computes phase
514  double phase = this->phase() + t * (f0() + t * (c1 * f1() + c2 * f2() * t));
515 
516  // Put phase into interval [0,1]
517  phase -= floor(phase);
518 
519  // Return phase
520  return phase;
521 }
522 
523 
524 /***********************************************************************//**
525  * @brief Evaluate phase curve value for a given phase
526  *
527  * @param[in] phase Phase.
528  * @return Value of phase curve.
529  *
530  * Computes
531  *
532  * \f[
533  * S_{\rm t}(\Phi) = r(\Phi) \times {\tt m\_norm}
534  * \f]
535  *
536  * where
537  *
538  * \f$r(\Phi)\f$ is the phase dependent rate, defined by linear interpolation
539  * between the nodes in a FITS file, \f$\Phi\f$ is the phase, and
540  * \f${\tt m\_norm}\f$ is a normalisation constant.
541  ***************************************************************************/
542 double GModelTemporalPhaseCurve::value(const double& phase) const
543 {
544  // Compute function value
545  double value = m_norm.value() * m_nodes.interpolate(phase, m_values);
546 
547  // Return
548  return value;
549 }
550 
551 
552 /***********************************************************************//**
553  * @brief Print phase curve information
554  *
555  * @param[in] chatter Chattiness.
556  * @return String containing phase curve information.
557  ***************************************************************************/
558 std::string GModelTemporalPhaseCurve::print(const GChatter& chatter) const
559 {
560  // Initialise result string
561  std::string result;
562 
563  // Continue only if chatter is not silent
564  if (chatter != SILENT) {
565 
566  // Append header
567  result.append("=== GModelTemporalPhaseCurve ===");
568 
569  // Append information
570  result.append("\n"+gammalib::parformat("Function file"));
571  result.append(m_filename.url());
572  if (normalize()) {
573  result.append(" [normalized]");
574  }
575  result.append("\n"+gammalib::parformat("Number of parameters"));
576  result.append(gammalib::str(size()));
577  for (int i = 0; i < size(); ++i) {
578  result.append("\n"+m_pars[i]->print(chatter));
579  }
580  result.append("\n"+gammalib::parformat("Number of phase nodes"));
581  result.append(gammalib::str(m_nodes.size()));
582 
583  // EXPLICIT: Append node values
584  if (chatter >= EXPLICIT) {
585  for (int i = 0; i < m_nodes.size(); ++i) {
586  result.append("\n");
587  result.append(gammalib::parformat(" Phase "+gammalib::str(m_nodes[i])));
588  result.append(gammalib::str(m_values[i]));
589  }
590  }
591 
592  } // endif: chatter was not silent
593 
594  // Return result
595  return result;
596 }
597 
598 
599 /*==========================================================================
600  = =
601  = Private methods =
602  = =
603  ==========================================================================*/
604 
605 /***********************************************************************//**
606  * @brief Initialise class members
607  ***************************************************************************/
609 {
610  // Initialise normalisation parameter
611  m_norm.clear();
612  m_norm.name("Normalization");
613  m_norm.unit("(relative value)");
614  m_norm.scale(1.0);
615  m_norm.value(1.0);
616  m_norm.range(0.0,1000.0);
617  m_norm.fix();
618  m_norm.gradient(0.0);
619  m_norm.has_grad(true);
620 
621  // Initialise reference MJD parameter
622  m_mjd.clear();
623  m_mjd.name("MJD");
624  m_mjd.unit("days");
625  m_mjd.scale(1.0);
626  m_mjd.value(51544.5);
627  m_mjd.range(1.0e-10,1.0e+10);
628  m_mjd.fix();
629  m_mjd.gradient(0.0);
630  m_mjd.has_grad(false);
631 
632  // Initialise phase at reference MJD parameter
633  m_phase.clear();
634  m_phase.name("Phase");
635  m_phase.unit("");
636  m_phase.scale(1.0);
637  m_phase.value(0.0);
638  m_phase.range(-1.0,1.0);
639  m_phase.free();
640  m_phase.gradient(0.0);
641  m_phase.has_grad(false);
642 
643  // Initialise frequency at reference MJD parameter
644  m_f0.clear();
645  m_f0.name("F0");
646  m_f0.unit("Hz");
647  m_f0.scale(1.0);
648  m_f0.value(1.0);
649  m_f0.range(0.0,1000.0);
650  m_f0.free();
651  m_f0.gradient(0.0);
652  m_f0.has_grad(false);
653 
654  // Initialise first frequency derivative at reference MJD parameter
655  m_f1.clear();
656  m_f1.name("F1");
657  m_f1.unit("s^-2");
658  m_f1.scale(1.0);
659  m_f1.value(0.1);
660  m_f1.range(-1000.0,1000.0);
661  m_f1.free();
662  m_f1.gradient(0.0);
663  m_f1.has_grad(false);
664 
665  // Initialise second frequency derivative at reference MJD parameter
666  m_f2.clear();
667  m_f2.name("F2");
668  m_f2.unit("s^-3");
669  m_f2.scale(1.0);
670  m_f2.value(0.01);
671  m_f2.range(-1000.0,1000.0);
672  m_f2.free();
673  m_f2.gradient(0.0);
674  m_f2.has_grad(false);
675 
676  // Set parameter pointer(s)
677  m_pars.clear();
678  m_pars.push_back(&m_norm);
679  m_pars.push_back(&m_mjd);
680  m_pars.push_back(&m_phase);
681  m_pars.push_back(&m_f0);
682  m_pars.push_back(&m_f1);
683  m_pars.push_back(&m_f2);
684 
685  // Initialise members
686  m_nodes.clear();
687  m_values.clear();
688  m_filename.clear();
689  m_scale = 1.0;
690  m_normalize = true;
691  m_has_normalize = false;
692 
693  // Return
694  return;
695 }
696 
697 
698 /***********************************************************************//**
699  * @brief Copy class members
700  *
701  * @param[in] model Phase curve.
702  ***************************************************************************/
704 {
705  // Copy members
706  m_norm = model.m_norm;
707  m_mjd = model.m_mjd;
708  m_phase = model.m_phase;
709  m_f0 = model.m_f0;
710  m_f1 = model.m_f1;
711  m_f2 = model.m_f2;
712  m_nodes = model.m_nodes;
713  m_values = model.m_values;
714  m_filename = model.m_filename;
715  m_scale = model.m_scale;
716  m_normalize = model.m_normalize;
718 
719  // Set parameter pointer(s)
720  m_pars.clear();
721  m_pars.push_back(&m_norm);
722  m_pars.push_back(&m_mjd);
723  m_pars.push_back(&m_phase);
724  m_pars.push_back(&m_f0);
725  m_pars.push_back(&m_f1);
726  m_pars.push_back(&m_f2);
727 
728  // Return
729  return;
730 }
731 
732 
733 /***********************************************************************//**
734  * @brief Delete class members
735  ***************************************************************************/
737 {
738  // Return
739  return;
740 }
741 
742 
743 /***********************************************************************//**
744  * @brief Load nodes from file
745  *
746  * @param[in] filename File name.
747  *
748  * @exception GException::invalid_value
749  * Phase curve FITS file is invalid
750  *
751  * Load the phase curve nodes from a FITS file. The information is expected
752  * in the first extension of the FITS file, containing two mandatory columns
753  * with names "PHASE" and "NORM".
754  ***************************************************************************/
756 {
757  // Set maximum phase curve normalization value, including a small margin
758  const double max_norm = 1.0 + 1.0e-8;
759 
760  // Clear nodes and values
761  m_nodes.clear();
762  m_values.clear();
763 
764  // Set filename
766 
767  // Continue only if filename is not empty
768  if (!filename.is_empty()) {
769 
770  // Load FITS file
771  GFits fits = GFits(filename);
772 
773  // Extract binary table (so far always load extension 1 as table)
774  GFitsTable* table = fits.table(1);
775 
776  // Extract columns
777  GFitsTableCol* phase_col = (*table)["PHASE"];
778  GFitsTableCol* norm_col = (*table)["NORM"];
779 
780  // Check that there are at least two nodes in table
781  if (phase_col->nrows() < 2) {
782  std::string msg = "\"PHASE\" column contains "+
783  gammalib::str(phase_col->nrows())+" rows but at "
784  "least two rows are required. Please specify a valid "
785  "phase curve file.";
787  }
788 
789  // Check that both columns are consistent
790  if (phase_col->nrows() != norm_col->nrows()) {
791  std::string msg = "\"PHASE\" and \"NORM\" columns have inconsistent "
792  "number of rows ("+
793  gammalib::str(phase_col->nrows())+", "+
794  gammalib::str(norm_col->nrows())+"). Please "
795  "specify a valid phase curve file.";
797  }
798 
799  // Set number of nodes
800  int nodes = phase_col->nrows();
801 
802  // Check that phase values are in ascending order and comprised between
803  // 0 and 1, and that no node is larger than 1
804  double last_phase = -1.0;
805  for (int i = 0; i < nodes; ++i) {
806 
807  // Check if phase has increased
808  if (last_phase >= 0.0 && phase_col->real(i) <= last_phase) {
809  std::string msg = "Phase value "+gammalib::str(phase_col->real(i))+
810  " in row "+gammalib::str(i+1)+" of \"PHASE\" "
811  "column is equal to or smaller than preceeding "
812  "value. Please provide a phase curve file with "
813  "monotonically increasing phase values.";
815  }
816 
817  // Check if phase is within [0,1]
818  if (phase_col->real(i) < 0.0 || phase_col->real(i) > 1.0) {
819  std::string msg = "Phase value "+gammalib::str(phase_col->real(i))+
820  " outside range [0,1]. Please provide a phase "
821  "curve file with phase values comprised in the "
822  "interval [0,1].";
824  }
825 
826  // Check if value is smaller than maximum allowed normalisation
827  if (norm_col->real(i) > max_norm) {
828  std::string msg = "Value "+gammalib::str(norm_col->real(i))+" at "
829  "phase "+gammalib::str(phase_col->real(i))+" is "
830  "larger than 1. Please provide a phase curve file "
831  "with normalizations not exceeding 1.";
833  }
834 
835  } // endfor: looped over nodes
836 
837  // If first phase is larger than zero then add last node with phase-1 as
838  // first node. This makes sure that a valid normalization exists for
839  // phase=0.
840  if (phase_col->real(0) > 0.0) {
841  m_nodes.append(phase_col->real(nodes-1)-1.0);
842  m_values.push_back(norm_col->real(nodes-1));
843  }
844 
845  // Extract nodes
846  for (int i = 0; i < nodes; ++i) {
847  m_nodes.append(phase_col->real(i));
848  m_values.push_back(norm_col->real(i));
849  }
850 
851  // If last node is smaller than one then add first node with phase+1 as
852  // last node. This makes sure that a valid normalization exists for
853  // phase=1.
854  if (phase_col->real(nodes-1) < 1.0) {
855  m_nodes.append(phase_col->real(0)+1.0);
856  m_values.push_back(norm_col->real(0));
857  }
858 
859  // Make sure that no node exceeds 1
860  for (int i = 0; i < m_values.size(); ++i) {
861  if (m_values[i] > 1.0) {
862  m_values[i] = 1.0;
863  }
864  }
865 
866  // Close FITS file
867  fits.close();
868 
869  // Optionally normalize nodes
870  if (m_normalize) {
871  normalize_nodes();
872  }
873 
874  } // endif: filename was not empty
875 
876  // Return
877  return;
878 }
879 
880 
881 /***********************************************************************//**
882  * @brief Normalise nodes
883  *
884  * @exception GException::invalid_value
885  * Phase curve FITS file is invalid
886  *
887  * Normalise the node values so that the integral over the phase interval
888  * [0,1] gives an average normalisation of 1.
889  ***************************************************************************/
891 {
892  // Initialise node integral
893  double sum = 0.0;
894 
895  // Loop over all nodes-1
896  for (int i = 0; i < m_nodes.size()-1; ++i) {
897 
898  // Compute phase values
899  double a = (m_nodes[i] < 0.0) ? 0.0 : m_nodes[i];
900  double b = (m_nodes[i+1] > 1.0) ? 1.0 : m_nodes[i+1];
901 
902  // Compute normalisation at phase values
903  double fa = m_nodes.interpolate(a, m_values);
904  double fb = m_nodes.interpolate(b, m_values);
905 
906  // Compute intergal for this interval using the Trapezoid rule
907  sum += 0.5 * (b - a) * (fa + fb);
908 
909  }
910 
911  // Throw an exception if the integral is not positive
912  if (sum <= 0.0) {
913  std::string msg = "Integral over phase curve is not positive ("+
914  gammalib::str(sum)+"). Please provide a valid phase "
915  "curve file.";
917  }
918 
919  // Set scale factor
920  m_scale = 1.0 / sum;
921 
922  // Normalise the node values
923  for (int i = 0; i < m_values.size(); ++i) {
924  m_values[i] *= m_scale;
925  }
926 
927  // Return
928  return;
929 }
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:192
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
bool m_has_normalize
XML has normalize attribute.
const double & factor_gradient(void) const
Return parameter factor gradient.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:828
std::vector< GModelPar * > m_pars
Parameter pointers.
bool normalize(void) const
Return normalization flag.
const std::string & name(void) const
Return parameter name.
GNodeArray m_nodes
Phase values of nodes.
double gradient(void) const
Return parameter gradient.
Abstract temporal model base class.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
GTime mjd(void) const
Return reference Modified Julian Day.
void init_members(void)
Initialise class members.
GModelPar m_f0
Frequency at reference MJD.
bool m_normalize
Normalize phase curve (default: true)
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
XML element node class.
Definition: GXmlElement.hpp:47
double m_scale
Scale factor due to normalization.
void nrows(const int &nrows)
Set number of rows in column.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:908
Random number generator class.
Definition: GRan.hpp:44
Time class.
Definition: GTime.hpp:54
std::vector< double > m_values
Function values at nodes.
Gammalib tools definition.
FITS file class.
Definition: GFits.hpp:63
#define G_READ
FITS file class interface definition.
double f2(void) const
Return second frequency derivative at reference Modified Julian Day.
virtual GModelTemporalPhaseCurve & operator=(const GModelTemporalPhaseCurve &model)
Assignment operator.
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1733
virtual GTimes mc(const double &rate, const GTime &tmin, const GTime &tmax, GRan &ran) const
Return vector of random event times.
bool is_free(void) const
Signal if parameter is free.
FITS table column abstract base class definition.
Time container class.
Definition: GTimes.hpp:45
const double & scale(void) const
Return parameter scale.
virtual std::string type(void) const
Return model type.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
void free(void)
Free a parameter.
Filename class.
Definition: GFilename.hpp:62
Temporal phase curve model class interface definition.
Abstract interface for FITS table column.
void fix(void)
Fix a parameter.
bool has_attribute(const std::string &name) const
Check if element has a given attribute.
void load_nodes(const GFilename &filename)
Load nodes from file.
int size(void) const
Return number of parameters.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1527
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
double f0(void) const
Return frequency at reference Modified Julian Day.
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
void copy_members(const GModelTemporalPhaseCurve &model)
Copy class members.
GChatter
Definition: GTypemaps.hpp:33
void free_members(void)
Delete class members.
virtual void read(const GXmlElement &xml)
Read model from XML element.
virtual ~GModelTemporalPhaseCurve(void)
Destructor.
GModelPar m_f1
First freq. derivative at reference MJD.
GModelPar m_norm
Normalization factor.
double phase(void) const
Return phase at reference Modified Julian Day.
#define G_NORMALISE_NODES
Interface definition for the temporal model registry class.
double value(const double &phase) const
Evaluate phase curve value for a given phase.
const GFilename & filename(void) const
Return file name.
const double & secs(void) const
Return time in seconds in native reference (TT)
Definition: GTime.hpp:153
#define G_WRITE
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
double norm(void) const
Return normalization factor.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
virtual std::string print(const GChatter &chatter=NORMAL) const
Print phase curve information.
GModelPar m_f2
Second freq. derivative at reference MJD.
Temporal model registry class definition.
virtual double real(const int &row, const int &inx=0) const =0
GModelTemporalPhaseCurve(void)
Void constructor.
void range(const double &min, const double &max)
Set minimum and maximum parameter boundaries.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
virtual GModelTemporal & operator=(const GModelTemporal &model)
Assignment operator.
const std::string & unit(void) const
Return parameter unit.
Exception handler interface definition.
double f1(void) const
Return first frequency derivative at reference Modified Julian Day.
virtual double eval(const GTime &srcTime, const bool &gradients=false) const
Evaluate function.
std::string tolower(const std::string &s)
Convert string to lower case.
Definition: GTools.cpp:845
virtual GModelTemporalPhaseCurve * clone(void) const
Clone phase curve.
void init_members(void)
Initialise class members.
GModelPar m_phase
Phase at reference MJD.
GModelPar m_mjd
Reference MJD.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1033
#define G_LOAD_NODES
void normalize_nodes(void)
Normalise nodes.
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
GFilename m_filename
Name of phase file function.
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
Temporal phase curve model class.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1579
void append(const GTime &time)
Append time to container.
Definition: GTimes.cpp:214
const GModelTemporalPhaseCurve g_temporal_phase_seed
virtual void write(GXmlElement &xml) const
Write model into XML element.
void free_members(void)
Delete class members.
virtual void clear(void)
Clear phase curve.
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1703
double exp(const double &lambda)
Returns exponential deviates.
Definition: GRan.cpp:291
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415
FITS table abstract base class interface definition.