# How to calculate planet positions using VSOP2013

I've tried to figure out how to calculate planet positions using the VSOP2013 files, but I'm completely lost how to do the calculations for each term. VSOP87 is straightforward, but that's mainly because of Meeus' explanations on how to calculate the terms, but I am still not able to find a source that clearly and simply explain step-by-step how to calculate planet positions using VSOP2013. The VSOP2013 readme and published article do give a terse explanation, but I have no clue how to actually implement it in practice mainly because I don't understand the mathematical lingo being used.

Thank you so much for everyone's input, especially for Greg Miller and Pierre Paquette. It must've taking a very long time for you guys to figure all of this out, and it is immensely appreciated.

I just want to clarify regarding the matrix to calculate $$X_Q$$, $$Y_Q$$, and $$Z_Q$$ as I am not proficient with matrices: are the calculations for each value to made as follows (just using the calculation of $$X_Q$$ as an example):

$$\begin{equation*} X_Q= x \cdot (\cos\phi) + y \cdot (-\sin \phi \cos \phi) + z \cdot (\sin \phi \cos \epsilon) \end{equation*}$$

And then I also just want to clarify regarding converting the values of $$k$$, $$h$$, $$q$$, and $$p$$ to standard orbital elements. Are the formulas as follows (as I understand from the PHP code):

\begin{align*} \tan\varpi &= \frac{h}{k} \\[5pt] \tan\Omega &=\frac{p}{q} \\[5pt] \omega &=\varpi-\Omega \\[5pt] M &= \lambda-\varpi\\[5pt] \sin \frac{1}{2}i &= \frac{q}{\cos\Omega}\\[5pt] e &= \frac{k}{\cos \varpi}\\ \end{align*}

My sincere apologies for these very basic questions.

• I have not looked into it myself much, but I believe it is quite similar to the JPL Development Ephemeris. Might not be exact, but might point you in the right direction: celestialprogramming.com/jpl-ephemeris-format/… Commented Sep 2 at 23:51
• After looking at the VSOP2013 format, it is nearly identical to the JPL Development Ephemeris. The only differences is the ASCII file has a short header (which can be ignored), and the number of coefficients are different (the number of coefficients also vary for DE versions, so no big deal). If you read the article I linked above, and get the coefficient counts from the VSOP2013 Readme, you should have no issue putting it together. But I'm kinda curious why you'd want to use VSOP2013 rather than the JPL DE, which is much more widely used. Commented Sep 3 at 2:44
• Thank you so much,Greg. The main reason for using VSOP is because my impression was that it was easier to do the calculations than using the DE files. I don't have much knowledge of programming, and I'm only an amateur backyard astronomer doing the calculations for my own interests. Because of these, I'm happy with a model that's less accurate than the JPL, but accurate enough to give a planet position to 0.1 arcseconds in ecliptic long. I easily used Excel to create solutions to the VSOP87 theory (really easy, and excel does very well!), so I was wondering if I can do the same with VSOP2013. Commented Sep 3 at 4:59
• It looks like there are two methods of using VSOP2013. The URL below has two directories, the "ephemerides" one works like the JPL DE, the "solution" one works more like their ELP/82 and other theories. It's possible that one might be smaller for a longer range than the JPL DE, so might have some use. Just note that VSOP2013 only has the Earth-Moon Barycenter, and doesn't actually have the location of the Earth or Moon. So, even though it's smaller, it will be less useful. The implementation is fairly straight forward though, I'll post an answer later tonight when I have time. Commented Sep 3 at 12:06
• ftp.imcce.fr/pub/ephem/planets/vsop2013 Commented Sep 3 at 12:16

As noted in the comments, VSOP2013 has two options for implementation, one resembles that of the JPL Development Ephemeris, the implementation of which is described in the article The JPL Ephemeris Format. Since such an implementation has little advantage over using the JPL DE directly, and the changes required should be quite obvious, won't re-explain implementing it here.

The other option is a periodic solution, and resembles the process used in other algorithms like the IAU 2006 Nutation series, or ELP82-2000. So there is some use in learning how to implement it even if you don't plan to use it.

The output of the VSOP2013 algorithm is the Keplerian orbital elements (a, $$\lambda$$, k, h, q, p). The user could then use these directly, or to compute the cartesian position and velocity. It includes data for the planets Mercury through Pluto, except it only included the Earth-Moon Barycenter, not the Earth or Moon itself.

The primary arguments to the algorithm, is T (the time), and the planet number 1-9. T can be computed from a Julian Date: $$T=\frac{jd-2451545.0}{365250}$$

Each planet is given its own file named VSOP2013pX.data (where X is the planet number 1-9). The files are divided into sections, based on the exponent of T which they will be multiplied by. The beginning of a section is denoted by a header like:

 VSOP2013  2  1  0  32412    VENUS      VARIABLE A   *T*00


Here the first field is just the VSOP2013 name, the next field (2) is the planet number (2=Venus), the next field is the variable number (1=A), and the next is the exponent for T, so every computation for this section gets multiplied by $$T^{exp}$$ after being summed together. The next field is the count of lines in this section.

Here is an example line from a section (I choose one far down to avoid one with mostly or all 0's).

  152   0  3  0  0   0  0  0  0  0   -1  -5   0   0      0   0  0  0 -0.3941567052167408 -09  0.2886541545155104 -08


Here, the first field is just a counter, the next 17 fields are coefficients to the lambda functions (explained below). The columns are an odd way to format floating point numbers, to parse them, add an "E" where there's a space. E.g. -0.3941567052167408E-09 and 0.2886541545155104E-08, these are S and C mentioned on page 2 of the readme, the coefficients to the sin and cos.

The README on page 2 gives 17 equations for $$\lambda$$l(). They are in the source code below, so I won't repeat them here. These 17 equations pair up with the 17 fields above, so the result of each equation is multiplied by the integer from its corresponding column. Since most are 0, they will have no effect, the line above should be computed as

$$\Phi = 3\lambda l(2) - \lambda l(10) - 5\lambda l(11)$$

And a running sum of all rows in a section is computed as $$sum = sum + -0.3941567052167408*10^9\sin(\Phi) + 0.2886541545155104*10^8 \cos(\Phi)$$

After summing all rows in a section, multiply the sum by $$T^{exp}$$, for the first section exp=0

$$a=sum * T^0$$

There will be multiple sections for different exponents of T per variable. For example, Venus has 9 sections for the exponents of T for the A variable. The counts for all variables for all planets is noted in the coeefcounts variable in the source code. The results from each section are just added together:

\begin{align*} a&=0 \\ a&=a + sum0 * T^0 \\ a&=a + sum1 * T^1 \\ a&=a + sum2 * T^2 \\ &... \\ a&=a + sum9 * T^9 \end{align*}

The process is continued for each of the 6 variables in the file. The code at the end of this post implements the full algorithm. An example of how to call it (assuming the code is in a file named vsop2013.py):

from vsop2013 import VSOP2013ComputePlanetElements
jd=2411545.0
elements=VSOP2013ComputePlanetElements(jd,1)
print(elements["a"])
print(elements["l"])
print(elements["k"])
print(elements["h"])
print(elements["q"])
print(elements["p"])


The full code in Python:

#Greg Miller [email protected] http://celestialprogramming.com
#released as public domain 2024
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.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):
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)

return {"a": a, "l": l%(2*math.pi), "k": k, "h": h, "q": q, "p": p}

path="./" #set to folder containing VSOP2013p*.dat files (always end with /)


An example returning the cartesian (XYZ) position and velocity in the J2000 reference frame is below. I was unfamiliar with the algorithm they use to solve Kepler's equation, so I converted their Fortran routine for an exact match rather than using a generic routine.

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])]

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))
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

jd=2411545.0
elements=VSOP2013ComputePlanetElements(jd,1)
xyz=ELLXYZ(1,elements)
print(rotateEclipticToICRS(xyz))

#result
#[0.3493878714121343, -0.13020772657955196, -0.10587303630039294, 0.006318722188132302, 0.023978752991141696, 0.012146771237284295]


VSOP 2013 does not provide positions but orbital elements, from which may be derived positions. Using the files found in the IMCCE’s FTP site linked above by Greg Miller, the following PHP function performs the required calculations. You may adapt it to your needs.

    function VSOP2013($$planète,$$JD) { // EARTH-MOON BARYCENTRE ONLY! NOT EARTH BY ITSELF
// Provides orbital elements… from which are derived positions
$$T = ($$JD - 2451545) / 365250; $$noLigne = 0;$$A = 0; $$L = 0;$$K = 0; $$H = 0;$$Q = 0; $$P = 0;$$λls = [[], [4.402608631669, 26087.90314068555], [3.176134461576, 10213.28554743445], [1.753470369433, 6283.075850353215], [6.203500014141, 3340.612434145457], [4.09136000305, 1731.170452721855], [1.713740719173, 1704.450855027201], [5.598641292287, 1428.948917844273], [2.805136360408, 1364.756513629990], [2.32698973462, 1361.923207632842], [.599546107035, 529.690961562325], [.874018510107, 213.299086108488], [5.481225395663, 74.781659030778], [5.311897933164, 38.132972226125], [0, .3595362285049309], [5.19846640063, 77713.7714481804], [1.62790513602, 84334.6615717837], [2.35555563875, 83286.9142477147]];
// Mercury 1, Venus 2, Earth-Moon 3, Mars 4, Vesta, Iris, Bamberga, Ceres, Pallas, Jupiter 5, Saturn 6, Uranus 7, Neptune 8, Pluto 9, D Moon, F Moon, ℓ Moon
for($$i = 1;$$i <= 17; $$i ++) { λl[i] = λls[i][0] + λls[i][1] * T; }$$file = fopen("VSOP2013/VSOP2013p" . $$planète . ".dat", "r");$$foras = [[0, 0], [6, 3], [10, 2], [13, 2], [16, 2], [19, 3], [23, 2], [26, 2], [29, 2], [32, 2], [35, 4], [40, 3], [44, 3], [48, 3], [52, 6], [59, 2], [63, 2], [66, 2]];
while(!feof($$file)) { ligne = fgets(file); if(substr(ligne, 1, 4) == "VSOP") { variable = substr(ligne, 49, 1); exposant = floatval(substr(ligne, 56, 2)); } if(substr(ligne, 1, 4) != "VSOP") { Φ = 0; for(lea = 1; lea <= 17; lea ++) { a[lea] = floatval(substr(ligne, foras[lea][0], foras[lea][1])); Φ += a[lea] * λl[lea]; } S = floatval(substr(ligne, 69, 19)); Sexp = floatval(substr(ligne, 89, 3)); S = S * pow(10, Sexp); C = floatval(substr(ligne, 93, 19)); Cexp = floatval(substr(ligne, 113, 3)); C = C * pow(10, Cexp); variable += pow(T, exposant) * (S * sin(Φ) + C * cos(Φ)); } } fclose($$file);
while($$L > 2 * pi())$$L -= 2 * pi(); while($$L < 0)$$L += 2 * pi();
$$π = atan2($$H, $$K);$$Ω = atan2($$P,$$Q); $$ω =$$π - $$Ω;$$M = $$L -$$π; $$in = 2 * asin($$Q / cos($$Ω));$$e = $$K / cos($$π);
$$ν = deg2rad(Kepler($$e, rad2deg($$M))[1]);$$r = ($$A * (1 - ²($$e))) / (1 + $$e * cos($$ν)); $$u =$$ω + $$ν;$$x = $$r * (cos($$Ω) * cos($$u) - sin($$Ω) * sin($$u) * cos($$in));
$$y =$$r * (sin($$Ω) * cos($$u) + cos($$Ω) * sin($$u) * cos($$in));$$z = $$r * sin($$in) * sin($$u);$$ε = deg2rad(23 + 26 / 60 + 21.41136 / 3600); $$φ = deg2rad(-.05188 / 3600);$$matrice = [[cos($$φ), -sin($$φ) * cos($$ε), sin($$φ) * cos($$ε)], [sin($$φ), cos($$φ) * cos($$ε), -cos($$φ) * sin($$ε)], [0, sin($$ε), cos($$ε)]];
[$$XQ,$$YQ, $$ZQ] = m33m31($$matrice, [$$x,$$y, $$z]);$$longitude = atan2($$y,$$x); $$latitude = asin($$z / $$r); while($$longitude > 2 * pi()) $$longitude -= 2 * pi(); while($$longitude < 0) $$longitude += 2 * pi(); return [$$A, $$L,$$K, $$H,$$Q, $$P,$$M, $$Ω,$$ω, $$in,$$e, $$r,$$ν, $$x,$$y, $$z,$$XQ, $$YQ,$$ZQ, $$longitude,$$latitude, $π]; // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 // VALUES ARE GIVEN IN RADIANS WHEN APPLICABLE!!! }  Here’s Kepler:  function Kepler($$e,$$M, $$dr = "",$$max = 57) { $$M = deg2rad($$M); $$F = sign($$M); $$M = abs($$M) / (2 * PI()); $$M = ($$M - floor($$M)) * 2 * PI() *$$F; if($$M < 0)$$M = $$M + 2 * PI();$$F = 1; if($$M > PI()) { F = -1; M = 2 * PI() - M; }$$E0 = PI() / 2; $$D = PI() / 4; for($$J = 1; $$J <=$$max; $$J ++) { M1 = E0 - e * sin(E0); E0 = E0 + D * sign(M - M1); D = D / 2; }$$E0 = $$E0 *$$F; if($$dr == "") { E0 = e360(rad2deg(E0)); ν = e360(2 * atand(sqrt((1 + e) / (1 - e)) * tand(E0 / 2))); } else { ν = 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E0 / 2)); } return [$$E0,$ν];
}


Here’s m33m31:

PHP. I wrote everything myself, from various sources. Sorry for the missing dependencies. Here’s m33m31:
function m33m31($$a,$$b) {
$$r1 =$$a[0][0] * $$b[0] +$$a[0][1] * $$b[1] +$$a[0][2] * $$b[2];$$r2 = $$a[1][0] *$$b[0] + $$a[1][1] *$$b[1] + $$a[1][2] *$$b[2];
$$r3 =$$a[2][0] * $$b[0] +$$a[2][1] * $$b[1] +$$a[2][2] * $$b[2]; return [$$r1, $$r2,$$r3];
}

• (This code was written by myself from what can be found in the IMCCE’s readme file.) Commented Sep 4 at 1:51
• Is that PHP or Perl? Also, are the functions Kepler and m33m31 available? Commented Sep 4 at 12:11
• PHP. I wrote everything myself, from various sources. Commented Sep 4 at 19:37