This is an example implemenation of VSOP2013, it reads the ASCII Files provided by IMCCE for the computation. For practical use, it would be highly recommended to convert the ASCII files to some binary form to speed up computation. You will need the files named solution/VSOP2013pX.dat (where X=1 to 9) from the link above to run the programs. Just for reference, a binary file using one byte for each of the 17 arguments, and 8 for the two coefficients would be about 66Mb. That makes it a bit impractical for use in a webpage, so the code below is in Python.
The code is divded into to programs, one for the main VSOP2013 computation wich produces Kelperian elements. The other is a test program which calls the main VSOP2013 routine and converts the Kelperian elements into ICRS/J2000 cartesian coordinates for both position and velocity. You may also find A more detailed breakdown of the implementation useful.
Note, VSOP2013 does not provide the actual position of the Earth, only the Earth-Moon Barycenter. To compute the position of the Earth, you will need an ephemeris algorithm like ELP/MPP02 which provides the position of the Moon, and use that to compute the Earth's position.
Test.py
from vsop2013 import VSOP2013ComputePlanetElements
import math
def dcmplx(r,i):
return [r,i]
def cdabs(c):
return math.sqrt(c[0]*c[0] + c[1]*c[1])
def dimag(c):
return c[1]
def dreal(c):
return c[0]
def dconjg(c):
return [c[0],-c[1]]
def mul(c1,c2):
r=c1[0]*c2[0] - c1[1]*c2[1]
i=c1[1]*c2[0] + c1[0]*c2[1]
return [r,i]
def scalar(c,s):
return [c[0]*s,c[1]*s]
def cdexp(c):
return [math.exp(c[0])*math.cos(c[1]), math.exp(c[0])*math.sin(c[1])]
def add(c1,c2):
return [c1[0]+c2[0],c1[1]+c2[1]]
def ELLXYZ (ibody,v):
dpi=6.283185307179586e0
gmp=[
4.9125474514508118699e-11,
7.2434524861627027000e-10,
8.9970116036316091182e-10,
9.5495351057792580598e-11,
2.8253458420837780000e-07,
8.4597151856806587398e-08,
1.2920249167819693900e-08,
1.5243589007842762800e-08,
2.1886997654259696800e-12
]
gmsol=2.9591220836841438269e-04
rgm=math.sqrt(gmp[ibody-1]+gmsol)
xa=v["a"]
xl=v["l"]
xk=v["k"]
xh=v["h"]
xq=v["q"]
xp=v["p"]
xfi=math.sqrt(1.e0-xk*xk-xh*xh)
xki=math.sqrt(1.e0-xq*xq-xp*xp)
u=1.e0/(1.e0+xfi)
z=dcmplx(xk,xh)
ex=cdabs(z)
ex2=ex*ex
ex3=ex2*ex
z1=dconjg(z)
gl=xl%dpi
gm=gl-math.atan2(xh,xk)
e=gl+(ex-0.125e0*ex3)*math.sin(gm) + 0.5e0*ex2*math.sin(2.e0*gm) + 0.375e0*ex3*math.sin(3.e0*gm)
while True:
z2=dcmplx(0.e0,e)
zteta=cdexp(z2)
z3=mul(z1,zteta)
dl=gl-e+dimag(z3)
rsa=1.e0-dreal(z3)
e=e+dl/rsa
if (abs(dl)<1.e-15): break
z1=scalar(z,u*dimag(z3))
z2=dcmplx(dimag(z1),-dreal(z1))
zto=scalar(add(add(scalar(z,-1),zteta),z2),1/rsa)
xcw=dreal(zto)
xsw=dimag(zto)
xm=xp*xcw-xq*xsw
xr=xa*rsa
w=[0,0,0,0,0,0]
w[0]=xr*(xcw-2.e0*xp*xm)
w[1]=xr*(xsw+2.e0*xq*xm)
w[2]=-2.e0*xr*xki*xm
xms=xa*(xh+xsw)/xfi
xmc=xa*(xk+xcw)/xfi
xn=rgm/pow(xa,1.5e0)
w[3]=xn*((2.e0*xp*xp-1.e0)*xms+2.e0*xp*xq*xmc)
w[4]=xn*((1.e0-2.e0*xq*xq)*xmc-2.e0*xp*xq*xms)
w[5]=2.e0*xn*xki*(xp*xms+xq*xmc)
return w
def rotateEclipticToICRS(v):
eps=(23 + 26/60.0 + 21.41136/60.0/60.0)*math.pi/180
l=(-0.05188/60.0/60.0)*math.pi/180
coseps=math.cos(eps)
sineps=math.sin(eps)
cosl=math.cos(l)
sinl=math.sin(l)
r=[0,0,0,0,0,0]
r[0]=cosl*v[0] - sinl*coseps*v[1] + sinl*sineps*v[2]
r[3]=cosl*v[3] - sinl*coseps*v[4] + sinl*sineps*v[5]
r[1]=sinl*v[0] + cosl*coseps*v[1] - cosl*sineps*v[2]
r[4]=sinl*v[3] + cosl*coseps*v[4] - cosl*sineps*v[5]
r[2]=sineps*v[1] + coseps*v[2]
r[5]=sineps*v[4] + coseps*v[5]
return r
aukm=1.495979e+8
jd=2411545.0
elements=VSOP2013ComputePlanetElements(jd,1)
xyz=ELLXYZ(1,elements)
icrs=rotateEclipticToICRS(xyz)
print(icrs)
#result
#[0.3493878714121343, -0.13020772657955196, -0.10587303630039294, 0.006318722188132302, 0.023978752991141696, 0.012146771237284295]
vsop2013.py
import math
def parseFields(line):
fields=[]
fields.append(line[ 5:9])
fields.append(line[ 9:12])
fields.append(line[12:15])
fields.append(line[15:18])
fields.append(line[18:22])
fields.append(line[22:25])
fields.append(line[25:28])
fields.append(line[28:31])
fields.append(line[31:34])
fields.append(line[34:39])
fields.append(line[39:43])
fields.append(line[43:47])
fields.append(line[47:51])
fields.append(line[51:58])
fields.append(line[58:62])
fields.append(line[62:65])
fields.append(line[65:68])
fields.append(line[68:88]+"E"+line[89:92])
fields.append(line[92:112]+"E"+line[113:116])
return fields
def computeExponentSum(file,lambdas,count):
sum=0
lines=[]
for i in range(count):
lines.append(file.readline())
lines.reverse()
for line in lines:
fields=parseFields(line)
phi=0
for i in range(17):
phi+=lambdas[i]*int(fields[i])
S=float(fields[17])
C=float(fields[18])
sum+=(S*math.sin(phi) + C*math.cos(phi))
return sum
def computeVariableSum(file,T,exps,lambdas):
sum=0
for x in range(exps):
l=file.readline()
fields=l.split()
exponent=fields[3]
count=fields[4]
sum+=pow(T,int(exponent))*computeExponentSum(file,lambdas,int(count))
return sum
def VSOP2013ComputePlanetElements(jd,planet):
T=(jd-2451545.0)/365250
coeefcounts=[
{},
{"Name": "MERCURY", "A":10, "L":10, "K":10, "H":10, "Q": 8, "P": 8},
{"Name": "VENUS", "A": 9, "L":10, "K":10, "H": 9, "Q": 8, "P": 8},
{"Name": "EARTH-MOON","A": 8, "L":10, "K": 9, "H": 9, "Q": 7, "P": 7},
{"Name": "MARS", "A": 8, "L": 9, "K": 8, "H": 8, "Q": 7, "P": 6},
{"Name": "JUPITER", "A":20, "L":18, "K":10, "H":10, "Q": 8, "P": 8},
{"Name": "SATURN", "A":20, "L":18, "K":10, "H":10, "Q": 8, "P": 8},
{"Name": "URANUS", "A":10, "L":10, "K":10, "H":10, "Q": 6, "P": 7},
{"Name": "NEPTUNE", "A":11, "L":10, "K":10, "H":10, "Q": 8, "P":10},
{"Name": "PLUTO", "A":12, "L":12, "K":12, "H":12, "Q":12, "P":12}
]
lambdas=[
4.402608631669 + 26087.90314068555 * T,
3.176134461576 + 10213.28554743445 * T,
1.753470369433 + 6283.075850353215 * T,
6.203500014141 + 3340.612434145457 * T,
4.091360003050 + 1731.170452721855 * T,
1.713740719173 + 1704.450855027201 * T,
5.598641292287 + 1428.948917844273 * T,
2.805136360408 + 1364.756513629990 * T,
2.326989734620 + 1361.923207632842 * T,
0.599546107035 + 529.6909615623250 * T,
0.874018510107 + 213.2990861084880 * T,
5.481225395663 + 74.78165903077800 * T,
5.311897933164 + 38.13297222612500 * T,
0.3595362285049309 * T,
5.198466400630 + 77713.7714481804 * T,
1.627905136020 + 84334.6615717837 * T,
2.355555638750 + 83286.9142477147 * T,
]
f=open(path+"VSOP2013p"+str(planet)+".dat")
a=computeVariableSum(f,T,coeefcounts[planet]["A"],lambdas)
l=computeVariableSum(f,T,coeefcounts[planet]["L"],lambdas)
k=computeVariableSum(f,T,coeefcounts[planet]["K"],lambdas)
h=computeVariableSum(f,T,coeefcounts[planet]["H"],lambdas)
q=computeVariableSum(f,T,coeefcounts[planet]["Q"],lambdas)
p=computeVariableSum(f,T,coeefcounts[planet]["P"],lambdas)
f.close()
return {"a": a, "l": l%(2*math.pi), "k": k, "h": h, "q": q, "p": p}
path="orig/" #set to folder containing VSOP2013p*.dat files (always end with /)