GammaLib  2.1.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-2021 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  // Verify number of model parameters
380 
381  // Get parameter pointers
388 
389  // Read parameters
390  m_norm.read(*norm);
391  m_mjd.read(*mjd);
392  m_phase.read(*phase);
393  m_f0.read(*f0);
394  m_f1.read(*f1);
395  m_f2.read(*f2);
396 
397  // Get optional normalization attribute
398  m_normalize = true;
399  m_has_normalize = false;
400  if (xml.has_attribute("normalize")) {
401  m_has_normalize = true;
402  std::string arg = xml.attribute("normalize");
403  if (arg == "0" || gammalib::tolower(arg) == "false") {
404  m_normalize = false;
405  }
406  }
407 
408  // Load nodes from file
410 
411  // Return
412  return;
413 }
414 
415 
416 /***********************************************************************//**
417  * @brief Write model into XML element
418  *
419  * @param[in] xml XML element.
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  // Verify model type
437 
438  // Get XML parameters
445 
446  // Write parameters
447  m_norm.write(*norm);
448  m_mjd.write(*mjd);
449  m_phase.write(*phase);
450  m_f0.write(*f0);
451  m_f1.write(*f1);
452  m_f2.write(*f2);
453 
454  // Set file attribute
456 
457  // Set optional normalization attribute
458  if (m_has_normalize || !m_normalize) {
459  if (m_normalize) {
460  xml.attribute("normalize", "1");
461  }
462  else {
463  xml.attribute("normalize", "0");
464  }
465  }
466 
467  // Return
468  return;
469 }
470 
471 
472 /***********************************************************************//**
473  * @brief Computes the phase for a given time
474  *
475  * @param[in] time Time.
476  *
477  * Computes
478  *
479  * \f[
480  * \Phi(t) = \Phi_0 + f(t-t_0) + \frac{1}{2}\dot{f} (t-t_0)^2 +
481  * \frac{1}{6}\ddot{f} (t-t_0)^3
482  * \f]
483  *
484  * where
485  * \f$t\f$ is the specified @p time in seconds,
486  * \f$t_0\f$ is a reference time in seconds,
487  * \f$\Phi_0\f$ is the phase at the reference time,
488  * \f$f\f$ is the variation frequency at the reference time,
489  * \f$\dot{f}\f$ is the first derivative of the variation frequency at the
490  * reference time, and
491  * \f$\dot{f}\f$ is the second derivative of the variation frequency at the
492  * reference time.
493  *
494  * The phase \f$\Phi(t)\f$ is in the interval [0,1].
495  ***************************************************************************/
496 double GModelTemporalPhaseCurve::phase(const GTime& time) const
497 {
498  // Set constants
499  const double c1 = 0.5;
500  const double c2 = 1.0 / 6.0;
501 
502  // Compute time since reference time in seconds
503  double dt = time - mjd();
504 
505  // Computes phase
506  double phase = this->phase() + dt * (f0() + dt * (c1 * f1() + c2 * f2() * dt));
507 
508  // Put phase into interval [0,1]
509  phase -= floor(phase);
510 
511  // Return phase
512  return phase;
513 }
514 
515 
516 /***********************************************************************//**
517  * @brief Evaluate phase curve value for a given phase
518  *
519  * @param[in] phase Phase.
520  * @return Value of phase curve.
521  *
522  * Computes
523  *
524  * \f[
525  * S_{\rm t}(\Phi) = r(\Phi) \times {\tt m\_norm}
526  * \f]
527  *
528  * where
529  *
530  * \f$r(\Phi)\f$ is the phase dependent rate, defined by linear interpolation
531  * between the nodes in a FITS file, \f$\Phi\f$ is the phase, and
532  * \f${\tt m\_norm}\f$ is a normalisation constant.
533  ***************************************************************************/
534 double GModelTemporalPhaseCurve::value(const double& phase) const
535 {
536  // Compute function value
537  double value = m_norm.value() * m_nodes.interpolate(phase, m_values);
538 
539  // Return
540  return value;
541 }
542 
543 
544 /***********************************************************************//**
545  * @brief Print phase curve information
546  *
547  * @param[in] chatter Chattiness.
548  * @return String containing phase curve information.
549  ***************************************************************************/
550 std::string GModelTemporalPhaseCurve::print(const GChatter& chatter) const
551 {
552  // Initialise result string
553  std::string result;
554 
555  // Continue only if chatter is not silent
556  if (chatter != SILENT) {
557 
558  // Append header
559  result.append("=== GModelTemporalPhaseCurve ===");
560 
561  // Append information
562  result.append("\n"+gammalib::parformat("Function file"));
563  result.append(m_filename.url());
564  if (normalize()) {
565  result.append(" [normalized]");
566  }
567  result.append("\n"+gammalib::parformat("Number of parameters"));
568  result.append(gammalib::str(size()));
569  for (int i = 0; i < size(); ++i) {
570  result.append("\n"+m_pars[i]->print(chatter));
571  }
572  result.append("\n"+gammalib::parformat("Number of phase nodes"));
573  result.append(gammalib::str(m_nodes.size()));
574 
575  // EXPLICIT: Append node values
576  if (chatter >= EXPLICIT) {
577  for (int i = 0; i < m_nodes.size(); ++i) {
578  result.append("\n");
579  result.append(gammalib::parformat(" Phase "+gammalib::str(m_nodes[i])));
580  result.append(gammalib::str(m_values[i]));
581  }
582  }
583 
584  } // endif: chatter was not silent
585 
586  // Return result
587  return result;
588 }
589 
590 
591 /*==========================================================================
592  = =
593  = Private methods =
594  = =
595  ==========================================================================*/
596 
597 /***********************************************************************//**
598  * @brief Initialise class members
599  ***************************************************************************/
601 {
602  // Initialise normalisation parameter
603  m_norm.clear();
604  m_norm.name("Normalization");
605  m_norm.unit("(relative value)");
606  m_norm.scale(1.0);
607  m_norm.value(1.0);
608  m_norm.range(0.0,1000.0);
609  m_norm.fix();
610  m_norm.gradient(0.0);
611  m_norm.has_grad(true);
612 
613  // Initialise reference MJD parameter
614  m_mjd.clear();
615  m_mjd.name("MJD");
616  m_mjd.unit("days");
617  m_mjd.scale(1.0);
618  m_mjd.value(51544.5);
619  m_mjd.range(1.0e-10,1.0e+10);
620  m_mjd.fix();
621  m_mjd.gradient(0.0);
622  m_mjd.has_grad(false);
623 
624  // Initialise phase at reference MJD parameter
625  m_phase.clear();
626  m_phase.name("Phase");
627  m_phase.unit("");
628  m_phase.scale(1.0);
629  m_phase.value(0.0);
630  m_phase.range(-1.0,1.0);
631  m_phase.free();
632  m_phase.gradient(0.0);
633  m_phase.has_grad(false);
634 
635  // Initialise frequency at reference MJD parameter
636  m_f0.clear();
637  m_f0.name("F0");
638  m_f0.unit("Hz");
639  m_f0.scale(1.0);
640  m_f0.value(0.0);
641  m_f0.range(0.0,1000.0);
642  m_f0.free();
643  m_f0.gradient(0.0);
644  m_f0.has_grad(false);
645 
646  // Initialise first frequency derivative at reference MJD parameter
647  m_f1.clear();
648  m_f1.name("F1");
649  m_f1.unit("s^-2");
650  m_f1.scale(1.0);
651  m_f1.value(0.0);
652  m_f1.range(-1000.0,1000.0);
653  m_f1.free();
654  m_f1.gradient(0.0);
655  m_f1.has_grad(false);
656 
657  // Initialise second frequency derivative at reference MJD parameter
658  m_f2.clear();
659  m_f2.name("F2");
660  m_f2.unit("s^-3");
661  m_f2.scale(1.0);
662  m_f2.value(0.0);
663  m_f2.range(-1000.0,1000.0);
664  m_f2.free();
665  m_f2.gradient(0.0);
666  m_f2.has_grad(false);
667 
668  // Set parameter pointer(s)
669  m_pars.clear();
670  m_pars.push_back(&m_norm);
671  m_pars.push_back(&m_mjd);
672  m_pars.push_back(&m_phase);
673  m_pars.push_back(&m_f0);
674  m_pars.push_back(&m_f1);
675  m_pars.push_back(&m_f2);
676 
677  // Initialise members
678  m_nodes.clear();
679  m_values.clear();
680  m_filename.clear();
681  m_scale = 1.0;
682  m_normalize = true;
683  m_has_normalize = false;
684 
685  // Return
686  return;
687 }
688 
689 
690 /***********************************************************************//**
691  * @brief Copy class members
692  *
693  * @param[in] model Phase curve.
694  ***************************************************************************/
696 {
697  // Copy members
698  m_norm = model.m_norm;
699  m_mjd = model.m_mjd;
700  m_phase = model.m_phase;
701  m_f0 = model.m_f0;
702  m_f1 = model.m_f1;
703  m_f2 = model.m_f2;
704  m_nodes = model.m_nodes;
705  m_values = model.m_values;
706  m_filename = model.m_filename;
707  m_scale = model.m_scale;
708  m_normalize = model.m_normalize;
710 
711  // Set parameter pointer(s)
712  m_pars.clear();
713  m_pars.push_back(&m_norm);
714  m_pars.push_back(&m_mjd);
715  m_pars.push_back(&m_phase);
716  m_pars.push_back(&m_f0);
717  m_pars.push_back(&m_f1);
718  m_pars.push_back(&m_f2);
719 
720  // Return
721  return;
722 }
723 
724 
725 /***********************************************************************//**
726  * @brief Delete class members
727  ***************************************************************************/
729 {
730  // Return
731  return;
732 }
733 
734 
735 /***********************************************************************//**
736  * @brief Load nodes from file
737  *
738  * @param[in] filename File name.
739  *
740  * @exception GException::invalid_value
741  * Phase curve FITS file is invalid
742  *
743  * Load the phase curve nodes from a FITS file. The information is expected
744  * in the first extension of the FITS file, containing two mandatory columns
745  * with names "PHASE" and "NORM".
746  ***************************************************************************/
748 {
749  // Set maximum phase curve normalization value, including a small margin
750  const double max_norm = 1.0 + 1.0e-8;
751 
752  // Clear nodes and values
753  m_nodes.clear();
754  m_values.clear();
755 
756  // Set filename
758 
759  // Continue only if filename is not empty
760  if (!filename.is_empty()) {
761 
762  // Load FITS file
763  GFits fits = GFits(filename);
764 
765  // Extract binary table (so far always load extension 1 as table)
766  GFitsTable* table = fits.table(1);
767 
768  // Extract columns
769  GFitsTableCol* phase_col = (*table)["PHASE"];
770  GFitsTableCol* norm_col = (*table)["NORM"];
771 
772  // Check that there are at least two nodes in table
773  if (phase_col->nrows() < 2) {
774  std::string msg = "\"PHASE\" column contains "+
775  gammalib::str(phase_col->nrows())+" rows but at "
776  "least two rows are required. Please specify a valid "
777  "phase curve file.";
779  }
780 
781  // Check that both columns are consistent
782  if (phase_col->nrows() != norm_col->nrows()) {
783  std::string msg = "\"PHASE\" and \"NORM\" columns have inconsistent "
784  "number of rows ("+
785  gammalib::str(phase_col->nrows())+", "+
786  gammalib::str(norm_col->nrows())+"). Please "
787  "specify a valid phase curve file.";
789  }
790 
791  // Set number of nodes
792  int nodes = phase_col->nrows();
793 
794  // Check that phase values are in ascending order and comprised between
795  // 0 and 1, and that no node is larger than 1
796  double last_phase = -1.0;
797  for (int i = 0; i < nodes; ++i) {
798 
799  // Check if phase has increased
800  if (last_phase >= 0.0 && phase_col->real(i) <= last_phase) {
801  std::string msg = "Phase value "+gammalib::str(phase_col->real(i))+
802  " in row "+gammalib::str(i+1)+" of \"PHASE\" "
803  "column is equal to or smaller than preceeding "
804  "value. Please provide a phase curve file with "
805  "monotonically increasing phase values.";
807  }
808 
809  // Check if phase is within [0,1]
810  if (phase_col->real(i) < 0.0 || phase_col->real(i) > 1.0) {
811  std::string msg = "Phase value "+gammalib::str(phase_col->real(i))+
812  " outside range [0,1]. Please provide a phase "
813  "curve file with phase values comprised in the "
814  "interval [0,1].";
816  }
817 
818  // Check if value is smaller than maximum allowed normalisation
819  if (norm_col->real(i) > max_norm) {
820  std::string msg = "Value "+gammalib::str(norm_col->real(i))+" at "
821  "phase "+gammalib::str(phase_col->real(i))+" is "
822  "larger than 1. Please provide a phase curve file "
823  "with normalizations not exceeding 1.";
825  }
826 
827  } // endfor: looped over nodes
828 
829  // If first phase is larger than zero then add last node with phase-1 as
830  // first node. This makes sure that a valid normalization exists for
831  // phase=0.
832  if (phase_col->real(0) > 0.0) {
833  m_nodes.append(phase_col->real(nodes-1)-1.0);
834  m_values.push_back(norm_col->real(nodes-1));
835  }
836 
837  // Extract nodes
838  for (int i = 0; i < nodes; ++i) {
839  m_nodes.append(phase_col->real(i));
840  m_values.push_back(norm_col->real(i));
841  }
842 
843  // If last node is smaller than one then add first node with phase+1 as
844  // last node. This makes sure that a valid normalization exists for
845  // phase=1.
846  if (phase_col->real(nodes-1) < 1.0) {
847  m_nodes.append(phase_col->real(0)+1.0);
848  m_values.push_back(norm_col->real(0));
849  }
850 
851  // Make sure that no node exceeds 1
852  for (int i = 0; i < m_values.size(); ++i) {
853  if (m_values[i] > 1.0) {
854  m_values[i] = 1.0;
855  }
856  }
857 
858  // Close FITS file
859  fits.close();
860 
861  // Optionally normalize nodes
862  if (m_normalize) {
863  normalize_nodes();
864  }
865 
866  } // endif: filename was not empty
867 
868  // Return
869  return;
870 }
871 
872 
873 /***********************************************************************//**
874  * @brief Normalise nodes
875  *
876  * @exception GException::invalid_value
877  * Phase curve FITS file is invalid
878  *
879  * Normalise the node values so that the integral over the phase interval
880  * [0,1] gives an average normalisation of 1.
881  ***************************************************************************/
883 {
884  // Initialise node integral
885  double sum = 0.0;
886 
887  // Loop over all nodes-1
888  for (int i = 0; i < m_nodes.size()-1; ++i) {
889 
890  // Compute phase values
891  double a = (m_nodes[i] < 0.0) ? 0.0 : m_nodes[i];
892  double b = (m_nodes[i+1] > 1.0) ? 1.0 : m_nodes[i+1];
893 
894  // Compute normalisation at phase values
895  double fa = m_nodes.interpolate(a, m_values);
896  double fb = m_nodes.interpolate(b, m_values);
897 
898  // Compute intergal for this interval using the Trapezoid rule
899  sum += 0.5 * (b - a) * (fa + fb);
900 
901  }
902 
903  // Throw an exception if the integral is not positive
904  if (sum <= 0.0) {
905  std::string msg = "Integral over phase curve is not positive ("+
906  gammalib::str(sum)+"). Please provide a valid phase "
907  "curve file.";
909  }
910 
911  // Set scale factor
912  m_scale = 1.0 / sum;
913 
914  // Normalise the node values
915  for (int i = 0; i < m_values.size(); ++i) {
916  m_values[i] *= m_scale;
917  }
918 
919  // Return
920  return;
921 }
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:482
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:864
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:48
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:944
Random number generator class.
Definition: GRan.hpp:44
Time class.
Definition: GTime.hpp:55
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:1946
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:1637
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.
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
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:156
#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:955
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:1143
#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:1342
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:1689
void append(const GTime &time)
Append time to container.
Definition: GTimes.cpp:216
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:1889
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:489
FITS table abstract base class interface definition.
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819