<p>As noted in the comments, <a href="https://ftp.imcce.fr/pub/ephem/planets/vsop2013" rel="nofollow noreferrer">VSOP2013</a> has two options for implementation, one resembles that of the JPL Development Ephemeris, the implementation of which is described in the article The <a href="https://www.celestialprogramming.com/jpl-ephemeris-format/jpl-ephemeris-format.html" rel="nofollow noreferrer">JPL Ephemeris Format</a>.  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.</p>
<p>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.</p>
<p>The output of the VSOP2013 algorithm is the Keplerian orbital elements (a, <span class="math-container">$ \lambda $</span>, 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.</p>
<p>The primary arguments to the algorithm, is T (the time), and the planet number 1-9.  T can be computed from a Julian Date: <span class="math-container">$T=\frac{jd-2451545.0}{365250} $</span></p>
<p>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:</p>
<pre><code> VSOP2013  2  1  0  32412    VENUS      VARIABLE A   *T*00
</code></pre>
<p>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 <span class="math-container">$T^{exp}$</span> after being summed together.  The next field is the count of lines in this section.</p>
<p>Here is an example line from a section (I choose one far down to avoid one with mostly or all 0's).</p>
<pre><code>  152   0  3  0  0   0  0  0  0  0   -1  -5   0   0      0   0  0  0 -0.3941567052167408 -09  0.2886541545155104 -08
</code></pre>
<p>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 &quot;E&quot; 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.</p>
<p>The <a href="https://ftp.imcce.fr/pub/ephem/planets/vsop2013/solution/README.pdf" rel="nofollow noreferrer">README</a> on page 2 gives 17 equations for <span class="math-container">$\lambda$</span>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</p>
<p><span class="math-container">$$
\Phi = 3\lambda l(2) - \lambda l(10) - 5\lambda l(11)
$$</span></p>
<p>And a running sum of all rows in a section is computed as
<span class="math-container">$$
sum = sum + -0.3941567052167408*10^9\sin(\Phi) + 0.2886541545155104*10^8 \cos(\Phi)
$$</span></p>
<p>After summing all rows in a section, multiply the sum by <span class="math-container">$T^{exp}$</span>, for the first section exp=0</p>
<p><span class="math-container">$$
a=sum * T^0
$$</span></p>
<p>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:</p>
<p><span class="math-container">$$
\begin{align*}
a&amp;=0 \\
a&amp;=a + sum0 * T^0 \\
a&amp;=a + sum1 * T^1 \\
a&amp;=a + sum2 * T^2 \\
 &amp;... \\
a&amp;=a + sum9 * T^9
\end{align*}
$$</span></p>
<p>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):</p>
<pre><code>from vsop2013 import VSOP2013ComputePlanetElements
jd=2411545.0
elements=VSOP2013ComputePlanetElements(jd,1)
print(elements[&quot;a&quot;])
print(elements[&quot;l&quot;])
print(elements[&quot;k&quot;])
print(elements[&quot;h&quot;])
print(elements[&quot;q&quot;])
print(elements[&quot;p&quot;])
</code></pre>
<p>The full code in Python:</p>
<pre><code>#Greg Miller gmiller@gregmiller.net 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]+&quot;E&quot;+line[89:92])
    fields.append(line[92:112]+&quot;E&quot;+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=[
        {},
        {&quot;Name&quot;: &quot;MERCURY&quot;,   &quot;A&quot;:10, &quot;L&quot;:10, &quot;K&quot;:10, &quot;H&quot;:10, &quot;Q&quot;: 8, &quot;P&quot;: 8},
        {&quot;Name&quot;: &quot;VENUS&quot;,     &quot;A&quot;: 9, &quot;L&quot;:10, &quot;K&quot;:10, &quot;H&quot;: 9, &quot;Q&quot;: 8, &quot;P&quot;: 8},
        {&quot;Name&quot;: &quot;EARTH-MOON&quot;,&quot;A&quot;: 8, &quot;L&quot;:10, &quot;K&quot;: 9, &quot;H&quot;: 9, &quot;Q&quot;: 7, &quot;P&quot;: 7},
        {&quot;Name&quot;: &quot;MARS&quot;,      &quot;A&quot;: 8, &quot;L&quot;: 9, &quot;K&quot;: 8, &quot;H&quot;: 8, &quot;Q&quot;: 7, &quot;P&quot;: 6},
        {&quot;Name&quot;: &quot;JUPITER&quot;,   &quot;A&quot;:20, &quot;L&quot;:18, &quot;K&quot;:10, &quot;H&quot;:10, &quot;Q&quot;: 8, &quot;P&quot;: 8},
        {&quot;Name&quot;: &quot;SATURN&quot;,    &quot;A&quot;:20, &quot;L&quot;:18, &quot;K&quot;:10, &quot;H&quot;:10, &quot;Q&quot;: 8, &quot;P&quot;: 8},
        {&quot;Name&quot;: &quot;URANUS&quot;,    &quot;A&quot;:10, &quot;L&quot;:10, &quot;K&quot;:10, &quot;H&quot;:10, &quot;Q&quot;: 6, &quot;P&quot;: 7},
        {&quot;Name&quot;: &quot;NEPTUNE&quot;,   &quot;A&quot;:11, &quot;L&quot;:10, &quot;K&quot;:10, &quot;H&quot;:10, &quot;Q&quot;: 8, &quot;P&quot;:10},
        {&quot;Name&quot;: &quot;PLUTO&quot;,     &quot;A&quot;:12, &quot;L&quot;:12, &quot;K&quot;:12, &quot;H&quot;:12, &quot;Q&quot;:12, &quot;P&quot;: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+&quot;VSOP2013p&quot;+str(planet)+&quot;.dat&quot;)
    a=computeVariableSum(f,T,coeefcounts[planet][&quot;A&quot;],lambdas)
    l=computeVariableSum(f,T,coeefcounts[planet][&quot;L&quot;],lambdas)
    k=computeVariableSum(f,T,coeefcounts[planet][&quot;K&quot;],lambdas)
    h=computeVariableSum(f,T,coeefcounts[planet][&quot;H&quot;],lambdas)
    q=computeVariableSum(f,T,coeefcounts[planet][&quot;Q&quot;],lambdas)
    p=computeVariableSum(f,T,coeefcounts[planet][&quot;P&quot;],lambdas)

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

path=&quot;./&quot; #set to folder containing VSOP2013p*.dat files (always end with /)
</code></pre>
<p>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.</p>
<pre><code>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[&quot;a&quot;]
    xl=v[&quot;l&quot;]
    xk=v[&quot;k&quot;]
    xh=v[&quot;h&quot;]
    xq=v[&quot;q&quot;]
    xp=v[&quot;p&quot;]

    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)&lt;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]
</code></pre>

CC BY-SA 4.0
2024-12-11T05:47:28.817Z
Nice. I was intending to come around to this but you have already done it simplifying my work. Thanks!  BTW in place of `dcmplx(r,i)` you can use Python's builtin `complex(r, i)` and its methods.
jamadagni
CC BY-SA 4.0
2024-12-11T16:41:17.353Z
I avoided any uncommon language features to make it easier to port to any other language.
Greg Miller
CC BY-SA 4.0
2024-09-04T02:20:09.553Z
2024-09-04T03:34:56.440Z
2024-09-04T03:34:56.440Z
Greg Miller
<p>VSOP 2013 does not provide <strong>positions</strong> but <strong>orbital elements</strong>, 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.</p>
<pre><code>    function VSOP2013(<span class="math-container">$planète, $</span>JD) { // EARTH-MOON BARYCENTRE ONLY! NOT EARTH BY ITSELF
    // Provides orbital elements… from which are derived positions
    <span class="math-container">$T = ($</span>JD - 2451545) / 365250; <span class="math-container">$noLigne = 0; $</span>A = 0; <span class="math-container">$L = 0; $</span>K = 0; <span class="math-container">$H = 0; $</span>Q = 0; <span class="math-container">$P = 0;
    $</span>λ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(<span class="math-container">$i = 1; $</span>i &lt;= 17; <span class="math-container">$i ++) { $λl[$i] = $λls[$i][0] + $λls[$i][1] * $T; }
    $</span>file = fopen(&quot;VSOP2013/VSOP2013p&quot; . <span class="math-container">$planète . ".dat", "r");
    $</span>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(<span class="math-container">$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 &lt;= 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($</span>file);
    while(<span class="math-container">$L &gt; 2 * pi()) $</span>L -= 2 * pi(); while(<span class="math-container">$L &lt; 0) $</span>L += 2 * pi();
    <span class="math-container">$π = atan2($</span>H, <span class="math-container">$K); $</span>Ω = atan2(<span class="math-container">$P, $</span>Q); <span class="math-container">$ω = $</span>π - <span class="math-container">$Ω; $</span>M = <span class="math-container">$L - $</span>π; <span class="math-container">$in = 2 * asin($</span>Q / cos(<span class="math-container">$Ω)); $</span>e = <span class="math-container">$K / cos($</span>π);
    <span class="math-container">$ν = deg2rad(Kepler($</span>e, rad2deg(<span class="math-container">$M))[1]); $</span>r = (<span class="math-container">$A * (1 - ²($</span>e))) / (1 + <span class="math-container">$e * cos($</span>ν)); <span class="math-container">$u = $</span>ω + <span class="math-container">$ν;
    $</span>x = <span class="math-container">$r * (cos($</span>Ω) * cos(<span class="math-container">$u) - sin($</span>Ω) * sin(<span class="math-container">$u) * cos($</span>in));
    <span class="math-container">$y = $</span>r * (sin(<span class="math-container">$Ω) * cos($</span>u) + cos(<span class="math-container">$Ω) * sin($</span>u) * cos(<span class="math-container">$in));
    $</span>z = <span class="math-container">$r * sin($</span>in) * sin(<span class="math-container">$u);
    $</span>ε = deg2rad(23 + 26 / 60 + 21.41136 / 3600); <span class="math-container">$φ = deg2rad(-.05188 / 3600);
    $</span>matrice = [[cos(<span class="math-container">$φ), -sin($</span>φ) * cos(<span class="math-container">$ε), sin($</span>φ) * cos(<span class="math-container">$ε)], [sin($</span>φ), cos(<span class="math-container">$φ) * cos($</span>ε), -cos(<span class="math-container">$φ) * sin($</span>ε)], [0, sin(<span class="math-container">$ε), cos($</span>ε)]];
    [<span class="math-container">$XQ, $</span>YQ, <span class="math-container">$ZQ] = m33m31($</span>matrice, [<span class="math-container">$x, $</span>y, <span class="math-container">$z]);
    $</span>longitude = atan2(<span class="math-container">$y, $</span>x); <span class="math-container">$latitude = asin($</span>z / <span class="math-container">$r);
    while($</span>longitude &gt; 2 * pi()) <span class="math-container">$longitude -= 2 * pi(); while($</span>longitude &lt; 0) <span class="math-container">$longitude += 2 * pi();
    return [$</span>A, <span class="math-container">$L, $</span>K, <span class="math-container">$H, $</span>Q, <span class="math-container">$P, $</span>M, <span class="math-container">$Ω, $</span>ω, <span class="math-container">$in, $</span>e, <span class="math-container">$r, $</span>ν, <span class="math-container">$x, $</span>y, <span class="math-container">$z, $</span>XQ, <span class="math-container">$YQ, $</span>ZQ, <span class="math-container">$longitude, $</span>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!!!
}
</code></pre>
<p>Here’s Kepler:</p>
<pre><code>    function Kepler(<span class="math-container">$e, $</span>M, <span class="math-container">$dr = "", $</span>max = 57) {
    <span class="math-container">$M = deg2rad($</span>M);
    <span class="math-container">$F = sign($</span>M);
    <span class="math-container">$M = abs($</span>M) / (2 * PI());
    <span class="math-container">$M = ($</span>M - floor(<span class="math-container">$M)) * 2 * PI() * $</span>F;
    if(<span class="math-container">$M &lt; 0) $</span>M = <span class="math-container">$M + 2 * PI();
    $</span>F = 1;
    if(<span class="math-container">$M &gt; PI()) { $F = -1; $M = 2 * PI() - $M; }
    $</span>E0 = PI() / 2;
    <span class="math-container">$D = PI() / 4;
    for($</span>J = 1; <span class="math-container">$J &lt;= $</span>max; <span class="math-container">$J ++) {
        $M1 = $E0 - $e * sin($E0);
        $E0 = $E0 + $D * sign($M - $M1);
        $D = $D / 2;
    }
    $</span>E0 = <span class="math-container">$E0 * $</span>F;
    if(<span class="math-container">$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 [$</span>E0, $ν];
}
</code></pre>
<p>Here’s m33m31:</p>
<pre><code>PHP. I wrote everything myself, from various sources. Sorry for the missing dependencies. Here’s m33m31:
function m33m31(<span class="math-container">$a, $</span>b) {
    <span class="math-container">$r1 = $</span>a[0][0] * <span class="math-container">$b[0] + $</span>a[0][1] * <span class="math-container">$b[1] + $</span>a[0][2] * <span class="math-container">$b[2];
    $</span>r2 = <span class="math-container">$a[1][0] * $</span>b[0] + <span class="math-container">$a[1][1] * $</span>b[1] + <span class="math-container">$a[1][2] * $</span>b[2];
    <span class="math-container">$r3 = $</span>a[2][0] * <span class="math-container">$b[0] + $</span>a[2][1] * <span class="math-container">$b[1] + $</span>a[2][2] * <span class="math-container">$b[2];
    return [$</span>r1, <span class="math-container">$r2, $</span>r3];
}
</code></pre>

CC BY-SA 4.0
2024-09-04T01:51:47.213Z
(This code was written by myself from what can be found in the IMCCE’s readme file.)
Pierre Paquette
CC BY-SA 4.0
2024-09-04T12:11:14.207Z
Is that PHP or Perl?  Also, are the functions Kepler and m33m31 available?
Greg Miller
CC BY-SA 4.0
2024-09-04T19:37:25.820Z
PHP. I wrote everything myself, from various sources.
Pierre Paquette
CC BY-SA 4.0
2024-09-04T01:51:02.683Z
2024-09-04T19:36:33.550Z
2024-09-04T19:36:33.550Z
Pierre Paquette
<p>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.</p>
<p>Would someone be able to help me out explaining how to go about this?</p>
<hr />
<p>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.</p>
<p>I just want to clarify regarding the matrix to calculate <span class="math-container">$X_Q$</span>, <span class="math-container">$Y_Q$</span>, and <span class="math-container">$Z_Q$</span> as I am not proficient with matrices: are the calculations for each value to made as follows (just using the calculation of <span class="math-container">$X_Q$</span> as an example):</p>
<p><span class="math-container">\begin{equation*}
X_Q= x \cdot (\cos\phi) + y \cdot (-\sin \phi \cos \phi) + z \cdot (\sin \phi \cos \epsilon)
\end{equation*}</span></p>
<p>And then I also just want to clarify regarding converting the values of <span class="math-container">$k$</span>, <span class="math-container">$h$</span>, <span class="math-container">$q$</span>, and <span class="math-container">$p$</span> to standard orbital elements. Are the formulas as follows (as I understand from the PHP code):</p>
<p><span class="math-container">\begin{align*}
\tan\varpi &amp;= \frac{h}{k} \\[5pt]
\tan\Omega &amp;=\frac{p}{q} \\[5pt]
\omega &amp;=\varpi-\Omega \\[5pt]
M &amp;= \lambda-\varpi\\[5pt]
\sin \frac{1}{2}i &amp;= \frac{q}{\cos\Omega}\\[5pt]
e &amp;= \frac{k}{\cos \varpi}\\
\end{align*}</span></p>
<p>My sincere apologies for these very basic questions.</p>

CC BY-SA 4.0
2024-09-02T23:51:49.083Z
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: https://www.celestialprogramming.com/jpl-ephemeris-format/jpl-ephemeris-format.html
Greg Miller
CC BY-SA 4.0
2024-09-03T02:44:20.737Z
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.
Greg Miller
CC BY-SA 4.0
2024-09-03T04:59:08.163Z
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.
Coenraad van Schoor
CC BY-SA 4.0
2024-09-03T12:06:12.080Z
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.
Greg Miller
CC BY-SA 4.0
2024-09-03T12:16:08.310Z
https://ftp.imcce.fr/pub/ephem/planets/vsop2013/
Greg Miller
CC BY-SA 4.0
2024-09-03T14:25:25.123Z
Thanks Greg. The imcce ftp site is where I got the VSOP2013 files from, I'm just not sure how to do the calculations though. Each file for each planet has the terms for the different parameters, and there is a README as well with instructions, but I'm struggling to make out heads or tails on how to implement it.
Coenraad van Schoor
<planet><positional-astronomy><celestial-mechanics>
How to calculate planet positions using VSOP2013
Coenraad van Schoor
# How to calculate planet positions using VSOP2013

(Asked by: username_0 on 2024-09-02)

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 username_1 username_1 and username_3. 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.

**Tags:** planet, positional-astronomy, celestial-mechanics

> **Comments:**
> **username_1**: 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: https://www.celestialprogramming.com/jpl-ephemeris-format/jpl-ephemeris-format.html  
> **username_1**: 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.  
> **username_0**: Thank you so much,username_1. 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.  
> **username_1**: 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.  
> **username_1**: https://ftp.imcce.fr/pub/ephem/planets/vsop2013/  
> **username_0**: Thanks username_1. The imcce ftp site is where I got the VSOP2013 files from, I'm just not sure how to do the calculations though. Each file for each planet has the terms for the different parameters, and there is a README as well with instructions, but I'm struggling to make out heads or tails on how to implement it.  

---

## 2 answers

### Answer by username_1 on 2024-09-04. Score: 2

As noted in the comments, [VSOP2013](https://ftp.imcce.fr/pub/ephem/planets/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](https://www.celestialprogramming.com/jpl-ephemeris-format/jpl-ephemeris-format.html). 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](https://ftp.imcce.fr/pub/ephem/planets/vsop2013/solution/README.pdf) 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:

```
#username_1 username_1 gmiller@username_1.net 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]

```

> **Comments:**
> **username_2**: Nice. I was intending to come around to this but you have already done it simplifying my work. Thanks!  BTW in place of `dcmplx(r,i)` you can use Python's builtin `complex(r, i)` and its methods.  
> **username_1**: I avoided any uncommon language features to make it easier to port to any other language.  

---

### Answer by username_3 on 2024-09-04. Score: 1

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

```

> **Comments:**
> **username_3**: (This code was written by myself from what can be found in the IMCCE’s readme file.)  
> **username_1**: Is that PHP or Perl?  Also, are the functions Kepler and m33m31 available?  
> **username_3**: PHP. I wrote everything myself, from various sources.  

---

