24 #include "Tempo2PhaseModel.h" 25 #include "panda/Error.h" 26 #include "panda/CaseInsensitiveString.h" 27 #include <boost/algorithm/string/trim.hpp> 28 #include <boost/filesystem.hpp> 29 #include <boost/filesystem/fstream.hpp> 59 void resize(
unsigned x,
unsigned y);
63 long double& coeff(
unsigned x,
unsigned y) {
64 return _coeff[y*_nx+x];
68 long double evaluate(
long double x,
long double y)
const;
71 void construct_x_derivative(Cheby2D& result)
const;
74 long double evaluate_1d(
long double x,
int index)
const;
82 void Cheby2D::resize(
unsigned x,
unsigned y)
91 _coeff =
static_cast<long double*
>(malloc(_nx*_ny*
sizeof(
long double)));
94 long double Cheby2D::evaluate_1d(
long double x,
int index)
const 97 long double di = 0.0L;
98 long double di1 = 0.0L;
100 long double* coeffcurr = _coeff + (index + 1) * _nx - 1;
101 for(std::size_t i = _nx - 1; i > 0; --i) {
104 di = 2.0L * x * di1 - di2 + *(coeffcurr--);
107 return x * di - di1 + 0.5L * (*coeffcurr--);
110 long double Cheby2D::evaluate(
long double x,
long double y)
const 116 long double y_coeff, dj, dj1, dj2;
119 for (
int j = _ny - 1; j > 0; --j) {
120 y_coeff = evaluate_1d(x, j);
123 dj = 2.0L * y * dj1 - dj2 + y_coeff;
125 return y * dj - dj1 + 0.5L * evaluate_1d(x, 0);
128 void Cheby2D::construct_x_derivative(Cheby2D& derived_cheby)
const 131 assert(derived_cheby._nx == _nx && derived_cheby._ny == _ny);
134 for (
unsigned iy=0; iy < _ny; ++iy)
136 unsigned ioff = i+_nx;
137 derived_cheby._coeff[ioff-1] = 0.0;
138 derived_cheby._coeff[ioff-2] = 2 * (_nx-1) * _coeff[ioff-1];
142 derived_cheby._coeff[i+ix]= derived_cheby._coeff[i+ix+2] + 2 * (ix+1) * _coeff[i+ix+1];
152 void load(boost::filesystem::ifstream& is,
unsigned& line_count)
155 std::istringstream iss;
158 auto next_line = [&]() ->
bool {
159 while( std::getline(is, line) ) {
161 boost::trim_left(line);
162 if(line.empty() || line.front() ==
'#') {
171 auto parse_error = [&](std::string
const& msg)
173 ska::panda::Error e(msg);
174 e <<
" on line " << line_count;
179 ska::panda::CaseInsensitiveString keyword;
180 if(!next_line())
return;
182 if(keyword ==
"ChebyModel") {
184 if(keyword !=
"BEGIN") {
185 parse_error(
"expecting ChebyModel BEGIN");
188 while( next_line() ) {
190 if(keyword ==
"ChebyModel") {
192 if(keyword ==
"BEGIN") {
193 parse_error(
"unexpected ChebyModel BEGIN");
196 if(keyword ==
"END") {
197 cheby.construct_x_derivative(frequency_cheby);
201 if(keyword ==
"PSRNAME") {
205 if(keyword ==
"SITENAME") {
209 if(keyword ==
"TIME_RANGE")
215 if(keyword ==
"FREQ_RANGE")
221 if(keyword ==
"DISPERSION_CONSTANT")
223 iss >> dispersion_constant;
226 if(keyword ==
"NCOEFF_TIME")
231 if (keyword ==
"NCOEFF_FREQ")
236 if (keyword ==
"COEFFS")
240 parse_error(
"unexpected COEFFS keyword before NCOEFF_TIME or NCOEFF_FREQ");
242 frequency_cheby.resize(nx,ny);
249 iss >> cheby.coeff(ix, iy);
255 parse_error(std::string(
"unrecognised keyword :") + keyword.c_str());
261 std::string sitename;
262 long double mjd_start, mjd_end;
263 long double freq_start, freq_end;
264 long double dispersion_constant;
266 Cheby2D frequency_cheby;
270 class ChebyModelSet {
273 for(
auto ptr : _models) {
277 std::vector<ChebyModel*>& models() {
return _models; };
280 std::size_t size()
const {
return _models.size(); }
284 long double phase(
long double mjd,
long double freq)
const 286 ChebyModel
const* model = get_nearest_model(mjd);
289 return model->cheby.evaluate(
290 -1.0L + 2.0L * (mjd - model->mjd_start)/(model->mjd_end - model->mjd_start)
292 -1.0L + 2.0L * (freq - model->freq_start)/(model->freq_end - model->freq_start));
298 void load(boost::filesystem::ifstream& is)
301 unsigned line_number=0;
302 ska::panda::CaseInsensitiveString keyword;
303 unsigned model_number;
304 while( std::getline(is, line) ) {
306 boost::trim_left(line);
307 if(line.empty() || line.front() ==
'#') {
310 std::istringstream iss(line);
311 iss >> keyword >> model_number;
312 if(keyword !=
"ChebyModelSet") {
313 throw ska::panda::Error(
"expecting a ChebyModelSet format");
315 _models.reserve(model_number);
316 for(
unsigned num =0; num < model_number; ++num)
318 ChebyModel* model =
new ChebyModel();
319 _models.push_back(model);
320 model->load(is, line_number);
327 ChebyModel
const* get_nearest_model(
long double mjd)
const 329 ChebyModel
const* rv_model =
nullptr;
330 long double min_offset=std::numeric_limits<long double>::max(), offset;
331 for(ChebyModel
const* model : _models) {
333 if(mjd < model->mjd_start || mjd > model->mjd_end) {
336 offset = std::abs((model->mjd_start+model->mjd_end)/2 - mjd);
337 if( offset < min_offset ) {
346 std::vector<ChebyModel*> _models;
353 namespace generators {
365 long double phase(
long double mjd,
long double freq)
367 return _modelset->phase(mjd, freq);
371 void load(boost::filesystem::path
const& filename) {
372 boost::filesystem::ifstream is(filename);
374 panda::Error e(
"unable to open file: ");
378 if(_modelset ==
nullptr) _modelset =
new ChebyModelSet;
383 ChebyModelSet* _modelset;
393 void Tempo2PhaseModel::load(boost::filesystem::path
const& tempo2_predictor_file)
396 if(!boost::filesystem::is_regular_file(tempo2_predictor_file))
398 panda::Error e(
"file does not exist:");
399 e << tempo2_predictor_file;
402 PANDA_LOG <<
"tempo2 predictor model loading from file " << tempo2_predictor_file;
404 _pred->load(tempo2_predictor_file);
408 Tempo2PhaseModel::~Tempo2PhaseModel()
412 double Tempo2PhaseModel::operator()(utils::ModifiedJulianClock::time_point
const& mjd, boost::units::quantity<data::MegaHertz, double>
const mhz)
const 414 long double const mjd_dbl = mjd.time_since_epoch().count();
415 return _pred->phase(mjd_dbl, mhz.value());
configuration parameters for a tempo2 predictor model
Some limits and constants for FLDO.
boost::filesystem::path const & filename() const
return the tempo2 file containing the model parameters