Add GModelTemporalPhaseCurve::value() method and add unit tests (#2154)

Jürgen Knödlseder
1 parent 4c6f4ae3
Showing 5 changed files with 93 additions and 36 deletions
NEWS
include/GModelTemporalPhaseCurve.hpp
pyext/GModelTemporalPhaseCurve.i
src/model/GModelTemporalPhaseCurve.cpp
test/test_GModel.cpp
@@ -23,7 +23,8 @@ The following classes have been renamed:
 The following methods have been added:
 - GSkyRegion::overlaps(GSkyRegions&)
 - GSkyMap::overlaps(GSkyRegion&)
-- GModelTemporalPhaseCurve::normalise()
+- GModelTemporalPhaseCurve::normalize()
+- GModelTemporalPhaseCurve::value(double&)
 The following methods have been removed:
 - n.a.
@@ -116,6 +116,7 @@ public:
     double           phase(void) const;
     double           phase(const GTime& time) const;
     void             phase(const double& phase);
+    double           value(const double& phase) const;
     double           f0(void) const;
     void             f0(const double& f0);
     double           f1(void) const;
@@ -130,7 +131,7 @@ protected:
     void copy_members(const GModelTemporalPhaseCurve& model);
     void free_members(void);
     void load_nodes(const GFilename& filename);
-    void normalise_nodes(void);
+    void normalize_nodes(void);
     // Protected parameter members
     GModelPar           m_norm;          //!< Normalization factor
@@ -73,6 +73,7 @@ public:
     double           phase(void) const;
     double           phase(const GTime& time) const;
     void             phase(const double& phase);
+    double           value(const double& phase) const;
     double           f0(void) const;
     void             f0(const double& f0);
     double           f1(void) const;
@@ -47,7 +47,7 @@ const GModelTemporalRegistry   g_temporal_phase_registry(&amp;g_temporal_phase_seed)
 #define G_READ                 "GModelTemporalPhaseCurve::read(GXmlElement&)"
 #define G_WRITE               "GModelTemporalPhaseCurve::write(GXmlElement&)"
 #define G_LOAD_NODES       "GModelTemporalPhaseCurve::load_nodes(GFilename&)"
-#define G_NORMALISE_NODES       "GModelTemporalPhaseCurve::normalise_nodes()"
+#define G_NORMALISE_NODES       "GModelTemporalPhaseCurve::normalize_nodes()"
 /* __ Macros _____________________________________________________________ */
@@ -522,6 +522,34 @@ double GModelTemporalPhaseCurve::phase(const GTime&amp; time) const
 /***********************************************************************//**
+ * @brief Evaluate phase curve value for a given phase
+ *
+ * @param[in] phase Phase.
+ * @return Value of phase curve.
+ *
+ * Computes
+ *
+ * \f[
+ *    S_{\rm t}(\Phi) = r(\Phi) \times {\tt m\_norm}
+ * \f]
+ *
+ * where
+ *
+ * \f$r(\Phi)\f$ is the phase dependent rate, defined by linear interpolation
+ * between the nodes in a FITS file, \f$\Phi\f$ is the phase, and
+ * \f${\tt m\_norm}\f$ is a normalisation constant.
+ ***************************************************************************/
+double GModelTemporalPhaseCurve::value(const double& phase) const
+{
+    // Compute function value
+    double value  = m_norm.value() * m_nodes.interpolate(phase, m_values);
+
+    // Return
+    return value;
+}
+
+
+/***********************************************************************//**
  * @brief Print phase curve information
  *
  * @param[in] chatter Chattiness.
@@ -833,9 +861,9 @@ void GModelTemporalPhaseCurve::load_nodes(const GFilename&amp; filename)
     // Close FITS file
     fits.close();
-    // Optionally normalise nodes
+    // Optionally normalize nodes
     if (m_normalize) {
-        normalise_nodes();
+        normalize_nodes();
     }
     // Return
@@ -852,7 +880,7 @@ void GModelTemporalPhaseCurve::load_nodes(const GFilename&amp; filename)
  * Normalise the node values so that the integral over the phase interval
  * [0,1] gives an average normalisation of 1.
  ***************************************************************************/
-void GModelTemporalPhaseCurve::normalise_nodes(void)
+void GModelTemporalPhaseCurve::normalize_nodes(void)
 {
     // Initialise node integral
     double sum = 0.0;
@@ -2322,52 +2322,78 @@ void TestGModel::test_temp_phasecurve(void)
     test_value(model2.f1(), 0.1);
     test_value(model2.f2(), 0.01);
     test_value(model2.norm(), 2.0);
-   
-    // Test XML constructor
-    GXml         xml(m_xml_model_point_temp_phasecurve);
-    GXmlElement* element = xml.element(0)->element(0)->element("temporal", 0);
-    GModelTemporalPhaseCurve model3(*element);
-    test_value(model3.size(), 6);
-    test_value(model3.type(), "PhaseCurve", "Check model type");
+    test_value(model2.value(0.1), 0.0);
+    test_value(model2.value(0.2), 9.09090909090909);
+    test_value(model2.value(0.6), 8.18181818181818);
+    test_assert(model2.normalize(), "Check that phase curve is normalised");
+
+    // Test value constructor without phase curve normalization
+    GModelTemporalPhaseCurve model3(m_temp_phasecurve,
+                                    GTime(3.0, "secs"),
+                                    0.0,
+                                    1.0,
+                                    0.1,
+                                    0.01,
+                                    2.0,
+                                    false);
     test_value(model3.filename().url(), m_temp_phasecurve,
                "Check phase curve file name");
-    test_value(model3.norm(), 1.0);
-    test_value(model3.mjd().mjd(), 51544.5);
+    test_value(model3.mjd().secs(), 3.0);
     test_value(model3.phase(), 0.0);
     test_value(model3.f0(), 1.0);
     test_value(model3.f1(), 0.1);
     test_value(model3.f2(), 0.01);
+    test_value(model3.norm(), 2.0);
+    test_value(model3.value(0.1), 0.0);
+    test_value(model3.value(0.2), 2.0);
+    test_value(model3.value(0.6), 1.8);
+    test_assert(!model3.normalize(), "Check that phase curve is not normalised");
+
+    // Test XML constructor
+    GXml         xml(m_xml_model_point_temp_phasecurve);
+    GXmlElement* element = xml.element(0)->element(0)->element("temporal", 0);
+    GModelTemporalPhaseCurve model4(*element);
+    test_value(model4.size(), 6);
+    test_value(model4.type(), "PhaseCurve", "Check model type");
+    test_value(model4.filename().url(), m_temp_phasecurve,
+               "Check phase curve file name");
+    test_value(model4.norm(), 1.0);
+    test_value(model4.mjd().mjd(), 51544.5);
+    test_value(model4.phase(), 0.0);
+    test_value(model4.f0(), 1.0);
+    test_value(model4.f1(), 0.1);
+    test_value(model4.f2(), 0.01);
     // Test filename method
-    model3.filename(m_temp_phasecurve);
-    test_value(model3.filename().url(), m_temp_phasecurve,
+    model4.filename(m_temp_phasecurve);
+    test_value(model4.filename().url(), m_temp_phasecurve,
                "Check phase curve file name");
     // Test parameter methods
-    model3.norm(3.0);
-    test_value(model3.norm(), 3.0);
-    model3.mjd(GTime(3.0, "s"));
-    test_value(model3.mjd().secs(), 3.0);
-    model3.phase(0.3);
-    test_value(model3.phase(), 0.3);
-    model3.f0(3.0);
-    test_value(model3.f0(), 3.0);
-    model3.f1(3.0);
-    test_value(model3.f1(), 3.0);
-    model3.f2(3.0);
-    test_value(model3.f2(), 3.0);
+    model4.norm(3.0);
+    test_value(model4.norm(), 3.0);
+    model4.mjd(GTime(3.0, "s"));
+    test_value(model4.mjd().secs(), 3.0);
+    model4.phase(0.3);
+    test_value(model4.phase(), 0.3);
+    model4.f0(3.0);
+    test_value(model4.f0(), 3.0);
+    model4.f1(3.0);
+    test_value(model4.f1(), 3.0);
+    model4.f2(3.0);
+    test_value(model4.f2(), 3.0);
     // Test operator access
     const char* strarray[] = {"Normalization", "MJD", "Phase", "F0", "F1", "F2"};
     for (int i = 0; i < 6; ++i) {
         std::string keyname(strarray[i]);
-        model3[keyname].remove_range(); // To allow setting of any value
-        model3[keyname].value(2.1);
-        model3[keyname].error(1.9);
-        model3[keyname].gradient(0.8);
-        test_value(model3[keyname].value(), 2.1);
-        test_value(model3[keyname].error(), 1.9);
-        test_value(model3[keyname].gradient(), 0.8);
+        model4[keyname].remove_range(); // To allow setting of any value
+        model4[keyname].value(2.1);
+        model4[keyname].error(1.9);
+        model4[keyname].gradient(0.8);
+        test_value(model4[keyname].value(), 2.1);
+        test_value(model4[keyname].error(), 1.9);
+        test_value(model4[keyname].gradient(), 0.8);
     }
     // Return