Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#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);
}