#include "ephcom_wrapper.hpp" Ephcom::Ephcom(std::string& DE_filename){ this->infp = NULL; this->ORIG = EPHCOM_EARTH; this->init(DE_filename); } Ephcom::~Ephcom(){ if( infp != NULL ){ fclose(infp); infp = NULL; } } void Ephcom::init(std::string& DE_filename){ infp = NULL; infp = fopen(DE_filename.c_str(),"r"); if( infp == NULL ){ std::cout << __FILE__ << ":" << __func__ << ":L:" << __LINE__ << ":" << " failed to open: " << DE_filename << std::endl; } else { std::cout << "Successfully loaded ephcom data from: " << DE_filename << std::endl; } ephcom_readbinary_header(infp, &(this->header)); this->coords.km = 1; // not AU, use kilometers this->coords.seconds = 0; // Timescale is days, not seconds this->coords.bary = 1; // Center is Solar System Barycenter this->coords.coordtype = 0; // No correction for light travel time or relativistic effects from Sun. } /* targ: one of the object, SUN, ..., MERCURY (see ephcom_wrapper.hpp) jd : Julian date pos : position of the object vel : velocity of the object */ void Ephcom::getPosVel(int targ, double jd, double pos[3], double vel[3]){ if( (targ < EPHCOM_MERCURY) || (targ > EPHCOM_SUN) ){ std::cout << __FILE__ << ":" << __func__ << ":L:" << __LINE__ << ":" << "targ must >= EPHCOM_MERCURY and <= EPHCOM_SUN" << std::endl; } double res[6]; this->coords.et2[0] = jd; this->coords.et2[1] = 0.0; ephcom_readbinary_header(this->infp, &this->header); double *datablock = (double *) malloc(this->header.ncoeff * sizeof(double)); if ( ephcom_get_coords(this->infp, &this->header, &this->coords, datablock) == 0 ) { ephcom_pleph(&this->coords, targ, this->ORIG, res); pos[0] = res[0]; pos[1] = res[1]; pos[2] = res[2]; vel[0] = res[3] / 86400; // convert to km/sec vel[1] = res[4] / 86400; vel[2] = res[5] / 86400; } else { std::cout << __FILE__ << ":" << __func__ << ":L:" << __LINE__ << ":" << "Coordinates not found for Julian date:" << jd << std::endl; } free(datablock); }