📄 readjpl_eph.cpp
字号:
// ReadJPL_EPH.cpp: implementation of the CReadJPL_EPH class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
//#include "spdg.h"
#include "ReadJPL_EPH.h"
#include "stdio.h"
#include "math.h"
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////////
// NTARG = INTEGER NUMBER OF 'TARGET' POINT.
//
// NCENT = INTEGER NUMBER OF CENTER POINT.
//
// THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS:
//
// 1 = MERCURY 8 = NEPTUNE
// 2 = VENUS 9 = PLUTO
// 3 = EARTH 10 = MOON
// 4 = MARS 11 = SUN
// 5 = JUPITER 12 = SOLAR-SYSTEM BARYCENTER
// 6 = SATURN 13 = EARTH-MOON BARYCENTER
// 7 = URANUS 14 = NUTATIONS (LONGITUDE AND OBLIQ)
// 15 = LIBRATIONS, IF ON EPH FILE
//
// (IF NUTATIONS ARE WANTED, SET NTARG = 14. FOR LIBRATIONS,
// SET NTARG = 15. SET NCENT=0.)
//
// RRD = OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY
// OF POINT 'NTARG' RELATIVE TO 'NCENT'. THE UNITS ARE AU AND
// AU/DAY. FOR LIBRATIONS THE UNITS ARE RADIANS AND RADIANS
// PER DAY. IN THE CASE OF NUTATIONS THE FIRST FOUR WORDS OF
// RRD WILL BE SET TO NUTATIONS AND RATES, HAVING UNITS OF
// RADIANS AND RADIANS/DAY.
//
// The option is available to have the units in km and km/sec.
// For this, set km=.true. in the STCOMX common block.
//////////////////////////////////////////////////////////////////////////
#define DENUM 405
/*#define DENUM 406*/
//#if DENUM==200
//#define KSIZE 1652
//#elif DENUM==403 || DENUM==405
#define KSIZE 2036
//#elif DENUM==404 || DENUM==406
//#define KSIZE 1456
//#endif
#define NRECL 4
#define RECSIZE (NRECL*KSIZE)
#define NCOEFF (KSIZE/2)
#define TRUE 1
#define FALSE 0
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CReadJPL_EPH::CReadJPL_EPH()
{
}
CReadJPL_EPH::~CReadJPL_EPH()
{
}
void CReadJPL_EPH::constan(char NAM[][6], double VAL[], double SSS[], int *N)
{
int i,j;
char spare1[RECSIZE-sizeof(TTL)+sizeof(CNAM)+sizeof(double)*5+sizeof(long)*41];
char spare2[RECSIZE-sizeof(double)*400];
fread(TTL,sizeof(TTL),1,Fp);
fread(CNAM,sizeof(CNAM),1,Fp);
fread(SS,sizeof(SS),1,Fp);
fread(&NCON,sizeof(long int),1,Fp);
fread(&AU,sizeof(double),1,Fp);
fread(&EMRAT,sizeof(double),1,Fp);
fread(IPT,sizeof(IPT),1,Fp);
fread(&NUMDE,sizeof(long),1,Fp);
fread(LPT,sizeof(LPT),1,Fp);
fread(spare1,sizeof(spare1),1,Fp);
fread(CVAL,sizeof(CVAL),1,Fp);
fread(spare2,sizeof(spare2),1,Fp);
*N=NCON;
for(i=0;i<3;i++) SSS[i]=SS[i];
for(i=0;i<=*N;i++)
{
VAL[i] = CVAL[i];
for(j=0;j<6;j++) NAM[i][j] = CNAM[i][j];
}
}
void CReadJPL_EPH::interp(double coef[], double t[], int ncf, int ncm, int na, int ifl, double posvel[])
{
static double pc[18],vc[18];
static int np=2, nv=3, first=1;
static double twot=0.0;
double dna,dt1,temp,tc,vfac,temp1;
int l,i,j;
if(first)
{ /* initialize static vectors when called first time */
pc[0]=1.0;
pc[1]=0.0;
vc[1]=1.0;
first=0;
}
/* entry point. get correct sub-interval number for this set
of coefficients and then get normalized chebyshev time
within that subinterval. */
dna=(double)na;
modf(t[0],&dt1);
temp=dna*t[0];
l=(int)(temp-dt1);
/* tc is the normalized chebyshev time (-1 <= tc <= 1) */
tc=2.0*(modf(temp,&temp1)+dt1)-1.0;
/* check to see whether chebyshev time has changed,
and compute new polynomial values if it has.
(the element pc[1] is the value of t1[tc] and hence
contains the value of tc on the previous call.) */
if(tc != pc[1])
{
np=2;
nv=3;
pc[1]=tc;
twot=tc+tc;
}
/* be sure that at least 'ncf' polynomials have been evaluated
and are stored in the array 'pc'. */
if(np < ncf)
{
for(i=np;i<ncf;i++) pc[i]=twot*pc[i-1]-pc[i-2];
np=ncf;
}
/* interpolate to get position for each component */
for(i=0;i<ncm;i++) /* ncm is a number of coordinates */
{
posvel[i]=0.0;
for(j=ncf-1;j>=0;j--)
posvel[i]=posvel[i]+pc[j]*coef[j+i*ncf+l*ncf*ncm];
}
if(ifl <= 1) return;
/* if velocity interpolation is wanted, be sure enough
derivative polynomials have been generated and stored. */
vfac=(dna+dna)/t[1];
vc[2]=twot+twot;
if(nv < ncf)
{
for(i=nv;i<ncf;i++) vc[i]=twot*vc[i-1]+pc[i-1]+pc[i-1]-vc[i-2];
nv=ncf;
}
/* interpolate to get velocity for each component */
for(i=0;i<ncm;i++)
{
posvel[i+ncm]=0.0;
for(j=ncf-1;j>0;j--)
posvel[i+ncm]=posvel[i+ncm]+vc[j]*coef[j+i*ncf+l*ncf*ncm];
posvel[i+ncm]=posvel[i+ncm]*vfac;
}
return;
}
void CReadJPL_EPH::split(double tt, double fr[])
{
/* main entry -- get integer and fractional parts */
fr[1]=modf(tt,&fr[0]);
if(tt >= 0.0 || fr[1] == 0.0) return;
/* make adjustments for negative input number */
fr[0]=fr[0]-1.0;
fr[1]=fr[1]+1.0;
return;
}
void CReadJPL_EPH::state(double et2[], int list[], double pv[][6], double nut[])
{
int i,j;
static int ipt[13][3], first=TRUE;
long int nr;
static long int nrl=0;
double pjd[4];
static double buf[NCOEFF];
double s,t[2],aufac;
double pefau[6];
if(first)
{
first=FALSE;
for(i=0;i<3;i++)
{
for(j=0;j<12;j++) ipt[j][i]=(int)IPT[j][i];
ipt[12][i] = (int)LPT[i];
}
}
// ********** main entry point ********** //
s=et2[0] - 0.5;
split(s,&pjd[0]);
split(et2[1],&pjd[2]);
pjd[0]=pjd[0]+pjd[2]+0.5;
pjd[1]=pjd[1]+pjd[3];
split(pjd[1],&pjd[2]);
pjd[0]=pjd[0]+pjd[2];
// here pjd[0] contains last midnight before epoch desired (in JED: *.5) //
// and pjd[3] contains the remaining, fractional part of the epoch //
if( (pjd[0]+pjd[3]) < SS[0] || (pjd[0]+pjd[3]) > SS[1] )
{
AfxMessageBox("Requested JED not within ephemeris limits!");
return;
}
nr=(long)((pjd[0]-SS[0])/SS[2])+2; // add 2 to adjus for the first two records containing header data
if(pjd[0] == SS[1]) nr=nr-1;
t[0]=( pjd[0]-( (1.0*nr-2.0)*SS[2]+SS[0] ) +
pjd[3] )/SS[2];
// read correct record if not in core (static vector buf[]) //
if(nr!=nrl)
{
nrl=nr;
fseek(Fp,nr*RECSIZE,SEEK_SET);
fread(buf,sizeof(buf),1,Fp);
}
if(KM)
{
t[1]=SS[2]*86400.0;
aufac=1.0;
}
else
{
t[1]=SS[2];
aufac=1.0/AU;
}
// every time interpolate Solar System barycentric sun state //
interp(&buf[ipt[10][0]-1],t,ipt[10][1],3,ipt[10][2],2,pefau);
for(i=0;i<6;i++) PVSUN[i]=pefau[i]*aufac;
// check and interpolate whichever bodies are requested //
for(i=0;i<10;i++)
{
if(list[i] == 0) continue;
interp(&buf[ipt[i][0]-1],t,ipt[i][1],3,ipt[i][2],list[i],pefau);
for(j=0;j<6;j++)
{
if(i < 9 && !BARY) pv[i][j]=pefau[j]*aufac-PVSUN[j];
else pv[i][j]=pefau[j]*aufac;
}
}
// do nutations if requested (and if on file) //
if(list[10] > 0 && ipt[11][1] > 0)
interp(&buf[ipt[11][0]-1],t,ipt[11][1],2,ipt[11][2],list[10],nut);
// get librations if requested (and if on file) //
if(list[11] > 0 && ipt[12][1] > 0)
{
interp(&buf[ipt[12][0]-1],t,ipt[12][1],3,ipt[12][2],list[11],pefau);
for(j=0;j<6;j++) pv[10][j]=pefau[j];
}
return;
}
void CReadJPL_EPH::pleph(double et, int ntarg, int ncent, double rrd[])
{
double et2[2],pv[13][6];
/* pv is the position/velocity array
NUMBERED FROM ZERO: 0=Mercury,1=Venus,...
8=Pluto,9=Moon,10=Sun,11=SSBary,12=EMBary
First 10 elements (0-9) are affected by state(),
all are adjusted here.
*/
int bsave,i,k;
int list[12];
/* list is a vector denoting, for which "body"
ephemeris values should be calculated by state():
0=Mercury,1=Venus,2=EMBary,...,8=Pluto,
9=geocentric Moon, 10=nutations in long. & obliq.
11= lunar librations
*/
/* initialize et2 for 'state' and set up component count */
et2[0]=et;
et2[1]=0.0;
for(i=0;i<6;i++) rrd[i]=0.0;
if(ntarg == ncent) return;
for(i=0;i<12;i++) list[i]=0;
/* check for nutation call */
if(ntarg == 14)
{
if(IPT[11][1] > 0) /* there is nutation on ephemeris */
{
list[10]=2;
state(et2,list,pv,rrd);
}
else AfxMessageBox("No nutations on the ephemeris file!");
return;
}
/* check for librations */
if(ntarg == 15)
{
if(LPT[1] > 0) /* there are librations on ephemeris file */
{
list[11]=2;
state(et2,list,pv,rrd);
for(i=0;i<6;i++) rrd[i]=pv[10][i]; /* librations */
}
else AfxMessageBox("No librations on the ephemeris file");
return;
}
/* force barycentric output by 'state' */
bsave=BARY;
BARY= TRUE;
/* set up proper entries in 'list' array for state call */
for(i=0;i<2;i++) /* list[] IS NUMBERED FROM ZERO ! */
{
k=ntarg-1;
if(i == 1) k=ncent-1; /* same for ntarg & ncent */
if(k <= 9) list[k]=2; /* Major planets */
if(k == 9) list[2]=2; /* for moon state earth state is necessary*/
if(k == 2) list[9]=2; /* for earth state moon state is necessary*/
if(k == 12) list[2]=2; /* EMBary state additionaly */
}
/* make call to state */
state(et2,list,pv,rrd);
/* Solar System barycentric Sun state goes to pv[10][] */
if(ntarg == 11 || ncent == 11) for(i=0;i<6;i++) pv[10][i]=PVSUN[i];
/* Solar System Barycenter coordinates & velocities equal to zero */
if(ntarg == 12 || ncent == 12) for(i=0;i<6;i++) pv[11][i]=0.0;
/* Solar System barycentric EMBary state: */
if(ntarg == 13 || ncent == 13) for(i=0;i<6;i++) pv[12][i]=pv[2][i];
/* if moon from earth or earth from moon ..... */
if( (ntarg*ncent) == 30 && (ntarg+ncent) == 13)
for(i=0;i<6;i++) pv[2][i]=0.0;
else
{
if(list[2] == 2) /* calculate earth state from EMBary */
for(i=0;i<6;i++) pv[2][i] -= pv[9][i]/(1.0+EMRAT);
if(list[9] == 2) /* calculate Solar System barycentric moon state */
for(i=0;i<6;i++) pv[9][i] += pv[2][i];
}
for(i=0;i<6;i++) rrd[i]=pv[ntarg-1][i]-pv[ncent-1][i];
BARY=bsave;
return;
}
void CReadJPL_EPH::InitJPL_EPH()
{
KM=1;
BARY=0;
char Path[MAX_PATH];
CString FullPath;
GetCurrentDirectory(MAX_PATH,Path);
FullPath=Path;
FullPath+="\\Model\\JPLEPH";
Fp=fopen(FullPath,"rb");
if(Fp==NULL) {AfxMessageBox("don't find the file:JPLEPH");return;}
char NAM[400][6];
double VAL[400];
double SSS[3];
int N;
constan(NAM,VAL,SSS,&N);
}
void CReadJPL_EPH::CloseJPL_EPH()
{
fclose(Fp);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -