fixed bug in pointing interpolation (now gives reasonable answers in test case)

Karl Kosack
1 parent b95897e9
Showing 5 changed files with 55 additions and 8 deletions
include/GHorizDir.hpp
inst/cta/include/GCTAPointing.hpp
inst/cta/src/GCTAPointing.cpp
inst/cta/test/test_CTA.cpp
inst/cta/test/test_CTA.hpp
@@ -31,6 +31,7 @@
 #include <string>
 #include "GBase.hpp"
 #include "GVector.hpp"
+#include "GMath.hpp"
  
 /* __ Compile options ____________________________________________________ */
 #define G_SINCOS_CACHE
@@ -33,6 +33,7 @@
 #include "GTime.hpp"
 #include "GMatrix.hpp"
 #include "GNodeArray.hpp"
+#include "GHorizDir.hpp"
  
 /***********************************************************************//**
  * @class GCTAPointing
@@ -91,8 +92,10 @@ protected:
  
     bool                m_has_table; //!< table is loaded
     GNodeArray          m_table_nodes;
-    std::vector<double> m_table_az;
-    std::vector<double> m_table_alt;
+    std::vector<double> m_table_az; //!< table of azimuths (rad)
+    std::vector<double> m_table_alt; //!< table of altitudes (rad)
+    GTime               m_table_tmin; //!< min time bound in table
+    GTime               m_table_tmax; //!<max time bound in table
  
     // Cached members
     mutable bool    m_has_cache;  //!< Has transformation cache
@@ -68,8 +68,11 @@ GCTAPointing::load_pointing_table(std::string filename)
     GFitsTableCol *stop  = (*table)["STOP"];
     GFitsTableCol *alt    = (*table)["ALT_PNT"];
     GFitsTableCol *az     = (*table)["AZ_PNT"];
-    //  GFitsTableCol *ra     = (*table)["RA_PNT"];
-    //GFitsTableCol *dec    = (*table)["DEC_PNT"];
+
+    // later on, may also interpolate RA/Dec, in the case of a drift
+    //  scan or other tracking mode:
+    //     GFitsTableCol *ra     = (*table)["RA_PNT"];
+    //     GFitsTableCol *dec    = (*table)["DEC_PNT"];
  
  
     // Loop over the elements and build a lookup table:
@@ -79,7 +82,7 @@ GCTAPointing::load_pointing_table(std::string filename)
  
     for (size_t ii=0; ii < table->nrows(); ii++) {
  
-      double midtime = 0.5*(stop->real(ii) - start->real(ii));
+      double midtime = 0.5*(stop->real(ii) + start->real(ii));
       GTime tt(midtime, "sec" );
       m_table_nodes.append( tt.secs() );
       m_table_az[ii] = az->real(ii) * gammalib::deg2rad;
@@ -187,7 +187,11 @@ void TestGCTAPointing::set(void)
     name("GCTAPointing");
  
     // Append tests to test suite
-    append(static_cast<pfunction>(&TestGCTAPointing::test_load_table), "Test load pointing from table");
+    append(static_cast<pfunction>(&TestGCTAPointing::test_load_table), 
+           "Test load pointing from table");
+
+    append(static_cast<pfunction>(&TestGCTAPointing::test_interpolate_altaz),
+           "Test alt/az interpolation given a time");
  
     return;
 }
@@ -206,9 +210,8 @@ TestGCTAPointing* TestGCTAPointing::clone(void) const
  
  
 /***********************************************************************//**
- * @brief Test CTA Npred computation
+ * @brief Test ability to load a CTA pointing table
  *
-
  ***************************************************************************/
 void TestGCTAPointing::test_load_table(void)
 {
@@ -221,6 +224,42 @@ void TestGCTAPointing::test_load_table(void)
  
  
 /***********************************************************************//**
+ * @brief Test interpolation of alt/az pointing dir as a function of
+ * time
+ *
+ ***************************************************************************/
+void TestGCTAPointing::test_interpolate_altaz(void)
+{
+
+    GCTAObservation run;
+
+    GCTAPointing pnt;
+    pnt.load_pointing_table( cta_point_table );
+
+    // test an out-of bounds time
+    
+
+    try {
+      GTime time(155470378.7, "sec"); 
+      GHorizDir dir = pnt.dir_horiz( time );
+    }
+    catch (GException::out_of_range &exc) {
+      std::cout << exc.what() << std::endl;
+    }
+
+    // Get a time that is somewhere in the run:
+
+    GTime time(155470378.7, "sec"); 
+    GHorizDir dir = pnt.dir_horiz( time );
+    
+    std::cout << "time="<<time<<" dir="<<dir<<std::endl;
+
+    return;
+}
+
+
+
+/***********************************************************************//**
  * @brief Test CTA Aeff computation
  ***************************************************************************/
 void TestGCTAResponse::test_response_aeff(void)
@@ -137,6 +137,7 @@ public:
     virtual void              set(void);
     virtual TestGCTAPointing* clone(void) const;
     void                      test_load_table(void);
+    void                      test_interpolate_altaz(void);
 };