Eirik Hatlen / gammalib

Commit b95897e93bc71944fcdcfb22f402fe89da0fad10

Authored by Karl Kosack
1 parent 0c9205e3

use GTime in interpolation

Inline Side-by-side
Showing 1 changed file with 33 additions and 7 deletions
inst/cta/src/GCTAPointing.cpp
  View file @b95897e
... ... @@ -31,6 +31,7 @@
31 31 #include "GCTAPointing.hpp"
32 32 #include "GTools.hpp"
33 33 #include "GFits.hpp"
  34 +#include "GHorizDir.hpp"
34 35  
35 36 /* __ Method name definitions ____________________________________________ */
36 37  
... ... @@ -45,11 +46,15 @@
45 46 * The advantage of this method is that ctools does not need to
46 47 * implement ra/dec to alt/az coordinate conversions, since the values
47 48 * are pre-calculated.
  49 + *
  50 + * \bug: GTimeReference should be used to convert the time from the
  51 + * FITS file to a GTime in the correct system
48 52 ************************************************************************/
49 53 void
50 54 GCTAPointing::load_pointing_table(std::string filename)
51 55 {
52 56  
  57 +
53 58 std::cout << "load pointing from: " << filename << std::endl;
54 59  
55 60 // Open the pointing file and find the pointing table:
... ... @@ -74,12 +79,23 @@ GCTAPointing::load_pointing_table(std::string filename)
74 79  
75 80 for (size_t ii=0; ii < table->nrows(); ii++) {
76 81  
77   - m_table_nodes.append( 0.5*(stop->real(ii) - start->real(ii)));
78   - m_table_az[ii] = az->real(ii);
79   - m_table_alt[ii] = alt->real(ii);
  82 + double midtime = 0.5*(stop->real(ii) - start->real(ii));
  83 + GTime tt(midtime, "sec" );
  84 + m_table_nodes.append( tt.secs() );
  85 + m_table_az[ii] = az->real(ii) * gammalib::deg2rad;
  86 + m_table_alt[ii] = alt->real(ii) * gammalib::deg2rad;
80 87  
  88 + if (ii==0) {
  89 + m_table_tmin = tt;
  90 + }
  91 +
  92 + if (ii==table->nrows()-1){
  93 + m_table_tmax = tt;
  94 + }
  95 +
81 96 }
82 97  
  98 +
83 99 m_has_table = true;
84 100  
85 101 }
... ... @@ -100,17 +116,27 @@ GCTAPointing::dir_horiz( const GTime &amp;time ) const {
100 116 // average direction...
101 117 }
102 118  
  119 + // first check if time is inside table bounds
  120 +
  121 + if (time<m_table_tmin || time>m_table_tmax) {
  122 + throw GException::out_of_range( __func__, time.secs(),
  123 + m_table_tmin.secs(),
  124 + m_table_tmax.secs());
  125 + }
  126 +
103 127 // get interpolated alt and az for the given time using the pointing
104 128 // table:
105 129  
106   - alt = m_table_nodes.interpolate( time.sec(), m_table_alt );
107   - az = m_table_nodes.interpolate( time.sec(), m_table_az );
  130 + double alt = m_table_nodes.interpolate( time.secs(), m_table_alt );
  131 + double az = m_table_nodes.interpolate( time.secs(), m_table_az );
108 132  
109 133 // construct a GHorizDir and return it:
110 134  
111   - return GHorizDir(alt,az);
  135 + GHorizDir dir;
  136 + dir.altaz(alt,az);
  137 + return dir;
112 138  
113   -}
  139 + }
114 140  
115 141  
116 142 /* __ Macros _____________________________________________________________ */
... ...