Skip to content

Commit

Permalink
CalcInterpolatedElements accepts user passed data
Browse files Browse the repository at this point in the history 
This allows to cleanup the code in l1.c, so that we don't have to rely
on a static variable to pass the moon index to the computation function.
  • Loading branch information
@guillaumechereau
guillaumechereau committed Dec 5, 2017
1 parent 2b5eeef commit 4cca649
Show file tree
Hide file tree
Showing 8 changed files with 41 additions and 34 deletions.
24 changes: 13 additions & 11 deletions src/core/planetsephems/calc_interpolated_elements.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,11 +28,13 @@ For documentation see the header file.

void CalcInterpolatedElements(const double t,double elem[],
const int dim,
void (*calc_func)(const double t,double elem[]),
void (*calc_func)(const double t,double elem[],
void *user),
const double delta_t,
double *t0,double e0[],
double *t1,double e1[],
double *t2,double e2[]) {
double *t2,double e2[],
void *user) {
/*
printf("CalcInterpolatedElements: %12.9f %12.9f %12.9f %12.9f\n",t,*t0,*t1,*t2);
*/
Expand All @@ -41,30 +43,30 @@ printf("CalcInterpolatedElements: %12.9f %12.9f %12.9f %12.9f\n",t,*t0,*t1,*t2);
*t0 = -1e100;
*t2 = -1e100;
*t1 = t;
(*calc_func)(*t1,e1);
(*calc_func)(*t1,e1,user);
for (i=0;i<dim;i++) elem[i] = e1[i];
return;
}
if (t <= *t1) {
if (*t1 - delta_t <= t) { /* interpolate */
if (*t0 < -1e99) {
*t0 = *t1 - delta_t;
(*calc_func)(*t0,e0);
(*calc_func)(*t0,e0,user);
}
} else if (*t1 - 2.0*delta_t <= t) { /* interpolate */
if (*t0 < -1e99) {
*t0 = *t1 - delta_t;
(*calc_func)(*t0,e0);
(*calc_func)(*t0,e0,user);
}
*t2 = *t1;*t1 = *t0;
for (i=0;i<dim;i++) {e2[i] = e1[i];e1[i] = e0[i];}
*t0 = *t1 - delta_t;
(*calc_func)(*t0,e0);
(*calc_func)(*t0,e0,user);
} else {
*t0 = -1e100;
*t2 = -1e100;
*t1 = t;
(*calc_func)(*t1,e1);
(*calc_func)(*t1,e1,user);
for (i=0;i<dim;i++) elem[i] = e1[i];
return;
}
Expand All @@ -78,22 +80,22 @@ printf("CalcInterpolatedElements: %12.9f %12.9f %12.9f %12.9f\n",t,*t0,*t1,*t2);
if (*t1 + delta_t >= t) { /* interpolate */
if (*t2 < -1e99) {
*t2 = *t1 + delta_t;
(*calc_func)(*t2,e2);
(*calc_func)(*t2,e2,user);
}
} else if (*t1 + 2.0*delta_t >= t) { /* interpolate */
if (*t2 < -1e99) {
*t2 = *t1 + delta_t;
(*calc_func)(*t2,e2);
(*calc_func)(*t2,e2,user);
}
*t0 = *t1;*t1 = *t2;
for (i=0;i<dim;i++) {e0[i] = e1[i];e1[i] = e2[i];}
*t2 = *t1 + delta_t;
(*calc_func)(*t2,e2);
(*calc_func)(*t2,e2,user);
} else {
*t0 = -1e100;
*t2 = -1e100;
*t1 = t;
(*calc_func)(*t1,e1);
(*calc_func)(*t1,e1,user);
for (i=0;i<dim;i++) elem[i] = e1[i];
return;
}
Expand Down
13 changes: 9 additions & 4 deletions src/core/planetsephems/calc_interpolated_elements.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,11 +25,13 @@ SOFTWARE.
extern
void CalcInterpolatedElements(const double t,double elem[],
const int dim,
void (*calc_func)(const double t,double elem[]),
void (*calc_func)(const double t,double elem[],
void *user),
const double delta_t,
double *t0,double e0[],
double *t1,double e1[],
double *t2,double e2[]);
double *t2,double e2[],
void *user);

/*
Simple interpolation routine with external cache.
Expand All @@ -38,8 +40,9 @@ The cache consists of 3 sets of values:
e1[0..dim-1] are the cached values at time *t1,
e2[0..dim-1] are the cached values at time *t2
delta_t is the time step: *t2-*t1 = *t1-*t0 = delta_t,
(*calc_func)(t,elem) calculates the values elem[0..dim-1] at time t,
t is the input parameter, elem[0..dim-1] are the output values.
(*calc_func)(t,elem,user) calculates the values elem[0..dim-1] at time t,
t is the input parameter, elem[0..dim-1] are the output values, user is the
user supplied data.
The user must supply *t0,*t1,*t2,e0,e1,e2.
The initial values must be *t0 = *t1 = *t2 = -1e100,
Expand All @@ -50,4 +53,6 @@ the user must never change them.
The user must always supply the same delta_t
for one set of (*t0,*t1,*t2,e0,e1,e2),
and of course the same dim and calc_func.
The user argument is passed to the calc_func callback.
*/
4 changes: 2 additions & 2 deletions src/core/planetsephems/elp82b.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162675,7 +162675,7 @@ static const double a0_div_ath_times_au =
384747.9806448954 / (384747.9806743165 * 149597870.691);

static
void GetElp82bSphericalCoor(const double t,double r[3]) {
void GetElp82bSphericalCoor(const double t,double r[3], void *user) {
int i,k;
double lambda[17];
double cos_sin_lambda[303*4];
Expand Down Expand Up @@ -162735,7 +162735,7 @@ void GetElp82bCoor(const double jd,double xyz[3]) {
const double t = (jd - 2451545.0) / 36525.0;
double r[3];
CalcInterpolatedElements(t,r,3,&GetElp82bSphericalCoor,DELTA_T,
&t_0,r_0,&t_1,r_1,&t_2,r_2);
&t_0,r_0,&t_1,r_1,&t_2,r_2,0);
{
const double rh = r[2] * cos(r[1]);
const double x3 = r[2] * sin(r[1]);
Expand Down
5 changes: 3 additions & 2 deletions src/core/planetsephems/gust86.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ const double phi[5] = {5.702313,
1.746237,
4.206896};

void CalcGust86Elem(double t,double elem[5*6]) {
void CalcGust86Elem(double t,double elem[5*6],void *user) {
double an[5],ae[5],ai[5];
int i;
for (i=0;i<5;i++) {
Expand Down Expand Up @@ -456,7 +456,8 @@ void GetGust86OsculatingCoor(const double jd0,const double jd,
&CalcGust86Elem,DELTA_T,
&t_0,gust86_elem_0,
&t_1,gust86_elem_1,
&t_2,gust86_elem_2);
&t_2,gust86_elem_2,
0);
/*
printf("GetGust86Coor(%d): %f %f %f %f %f %f\n",
body,
Expand Down
14 changes: 5 additions & 9 deletions src/core/planetsephems/l1.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -902,7 +902,8 @@ static const struct L1Body l1_bodies[4] = {
},
};

static void CalcL1Elem(const double t,const int body,double elem[6]) {
static void CalcL1Elem(const double t, double elem[6], void *user) {
const int body = *((int*)user);
int j;
const struct L1Body *const bp = l1_bodies + body;
const double *cheb = bp->cheb_coef;
Expand Down Expand Up @@ -992,11 +993,6 @@ static double l1_elem_2[4*6];
static double l1_jd0[4] = {-1e100,-1e100,-1e100,-1e100};
static double l1_elem[4*6];

static int ugly_static_parameter_body = -1;
static void CalcUglyStaticL1Elem(double t,double elem[6]) {
CalcL1Elem(t,ugly_static_parameter_body,elem);
}

void GetL1Coor(double jd, int body, double *xyz, double *xyzdot) {
double xyz6[6];
GetL1OsculatingCoor(jd,jd,body,xyz6);
Expand All @@ -1010,12 +1006,12 @@ void GetL1OsculatingCoor(const double jd0,const double jd,
if (jd0 != l1_jd0[body]) {
const double t0 = jd0 - 2433282.5;
l1_jd0[body] = jd0;
ugly_static_parameter_body = body;
CalcInterpolatedElements(t0,l1_elem+(body*6),6,
&CalcUglyStaticL1Elem,DELTA_T,
&CalcL1Elem, DELTA_T,
t_0+body,l1_elem_0+(body*6),
t_1+body,l1_elem_1+(body*6),
t_2+body,l1_elem_2+(body*6));
t_2+body,l1_elem_2+(body*6),
(void*)&body);
}
EllipticToRectangularA(l1_bodies[body].mu,l1_elem+(body*6),jd-jd0,x);
xyz[0] = L1toVsop87[0]*x[0]+L1toVsop87[1]*x[1]+L1toVsop87[2]*x[2];
Expand Down
5 changes: 3 additions & 2 deletions src/core/planetsephems/marssat.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -421,7 +421,7 @@ static double marssat_elem_2[2*6];
static double marssat_jd0 = -1e100;
static double marssat_elem[2*6];

static void CalcAllMarsSatElem(double t,double elem[12]) {
static void CalcAllMarsSatElem(double t,double elem[12], void *user) {
CalcMarsSatElem(t,0,elem+(0*6));
CalcMarsSatElem(t,1,elem+(1*6));
}
Expand All @@ -445,7 +445,8 @@ void GetMarsSatOsculatingCoor(const double jd0,const double jd,
&CalcAllMarsSatElem,DELTA_T,
&t_0,marssat_elem_0,
&t_1,marssat_elem_1,
&t_2,marssat_elem_2);
&t_2,marssat_elem_2,
NULL);
GenerateMarsSatToVSOP87(t0,mars_sat_to_vsop87);
}
EllipticToRectangularA(mars_sat_bodies[body].mu,marssat_elem+(body*6),jd-jd0,x);
Expand Down
5 changes: 3 additions & 2 deletions src/core/planetsephems/tass17.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3155,7 +3155,7 @@ static double tass17_elem_2[TASS17_DIM];
static double tass17_jd0 = -1e100;
static double tass17_elem[TASS17_DIM];

void CalcAllTass17Elem(const double t,double elem[TASS17_DIM])
void CalcAllTass17Elem(const double t,double elem[TASS17_DIM], void *user)
{
int body;
double lon[8];
Expand Down Expand Up @@ -3183,7 +3183,8 @@ void GetTass17OsculatingCoor(const double jd0,const double jd, const int body,do
&CalcAllTass17Elem,DELTA_T,
&t_0,tass17_elem_0,
&t_1,tass17_elem_1,
&t_2,tass17_elem_2);
&t_2,tass17_elem_2,
0);
/*
printf("GetTass17Coor(%d): %f %f %f %f %f %f\n",
body,
Expand Down
5 changes: 3 additions & 2 deletions src/core/planetsephems/vsop87.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137263,7 +137263,7 @@ void AccumulateVsop87Terms(const unsigned char *instructions,
}

static
void CalcVsop87Elem(const double t,double elem[8*6]) {
void CalcVsop87Elem(const double t,double elem[8*6], void *user) {
unsigned int i;
double lambda[12];
double cos_sin_lambda[203*4];
Expand Down Expand Up @@ -137347,7 +137347,8 @@ void GetVsop87OsculatingCoor(const double jd0,const double jd,const int body,dou
&CalcVsop87Elem,DELTA_T,
&t_0,vsop87_elem_0,
&t_1,vsop87_elem_1,
&t_2,vsop87_elem_2);
&t_2,vsop87_elem_2,
0);
}
EllipticToRectangularA(vsop87_mu[body],vsop87_elem+(body*6),jd-jd0,xyz);
}

0 comments on commit 4cca649

Please sign in to comment.