How to calculate planet positions using VSOP2013

Question

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.

Would someone be able to help me out explaining how to go about this?

---

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.

Answer 1

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

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

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

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]

Answer 2

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];
}