ephcom_wrapper.cpp

#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);
}