Post Reply 
Planetary Positions
05-17-2018, 02:46 PM (This post was last modified: 05-17-2018 02:49 PM by toml_12953.)
Post: #1
Planetary Positions
Here's a program by Stephen R. Schmitt (converted to HPPL by me) to compute the positions of the planets (and one dwarf planet!) Why you'd want to do this IDK but obviously someone does, so enjoy!


Code:
day_number();
get_coord();
mean_elements();
true_anomaly();
ha2str();
dec2str();
pname();
get_int();
abs_floor();
mod2pi();

degs := 180 / pi;
rads := pi / 180;
eps  := 1.0e-10;

ai :=1;  // semi-major axis [AU]
ei :=2;  // eccentricity of orbit
ii :=3;  // inclination of orbit [deg]
Oi :=4;  // longitude of the ascending node [deg]
wi :=5;  // longitude of perihelion [deg]
Li :=6;  // mean longitude [deg]

ra:=1; dec:=2; rvec:=3;

export PlanetPostions()             
begin

    local year, month, day, hour, mins, p;
    
    // get the date and universal time
    year  := get_int( "year" );
    month := get_int( "month" );
    day   := get_int( "day" );
    hour  := get_int( "hour" );
    mins  := get_int( "minute" );

    // compute day number for date/time
    local d := day_number( year, month, day, hour, mins );
    print();
    print("Date: "+ month+ "/"+ day+ "/"+ year+ " UT: "+ hour+ ":"+ mins);
    print("days since J2000: "+ d);
    print(" ");    
    print("Object      RA        DEC       Distance");
    print("---------------------------------------");

    // compute location of objects
    for p := 1 to 9 do  

        local obj := makemat(0,3);
        obj := get_coord( obj, p, d );

        print(pname(p)+" "+ ha2str(obj(ra)*degs)+" "+ dec2str(obj(dec)*degs)+" "+round(obj(rvec),6));

    end;

    print(" ");

end;

// day number to/from J2000 (Jan 1.5, 2000)
day_number( y, m, d, hour, mins )
begin

    local rv;
    local h := hour + mins / 60;
    
    rv := 367*y - ip(7*ip(y + ip(m + 9) / 12) / 4) 
        + ip(275*m / 9) + d - 730531.5 + h / 24;

    return rv;
    
end;

// compute RA, DEC, and distance of planet-p for day number-d
// result returned in structure obj
get_coord( obj, p, d )
begin

    local planet := makemat(0,6);
    planet := mean_elements( planet, p, d );
    local ap := planet(ai);
    local ep := planet(ei);
    local ipp := planet(ii);
    local op := planet(Oi);
    local pp := planet(wi);
    local lp := planet(Li);

    local earth := makemat(0,6);
    earth := mean_elements( earth, 3, d ); 
    local ae := earth(ai);
    local ee := earth(ei);
    local ie := earth(ii);
    local oe := earth(Oi);
    local pe := earth(wi);
    local le := earth(Li);
    
    // position of Earth in its orbit
    local me := mod2pi(le - pe);
    local ve := true_anomaly(me, ee);
    local re := ae*(1 - ee*ee) / (1 + ee*cos(ve));
    
    // heliocentric rectangular coordinates of Earth
    local xe := re*cos(ve + pe);
    local ye := re*sin(ve + pe);
    local ze := 0;
    
    // position of planet in its orbit
    local mp := mod2pi(lp - pp);
    local vp := true_anomaly(mp, planet(ei));
    local rp := ap*(1 - ep*ep) / (1 + ep*cos(vp));
    
    // heliocentric rectangular coordinates of planet
    local xh := rp*(cos(op)*cos(vp + pp - op) - sin(op)*sin(vp + pp - op)*cos(ipp));
    local yh := rp*(sin(op)*cos(vp + pp - op) + cos(op)*sin(vp + pp - op)*cos(ipp));
    local zh := rp*(sin(vp + pp - op)*sin(ipp));

    if p = 3 then                       // earth --> compute sun
        xh := 0;
        yh := 0;
        zh := 0;
    end;
    
    // convert to geocentric rectangular coordinates
    local xg := xh - xe;
    local yg := yh - ye;
    local zg := zh - ze;
    
    // rotate around x axis from ecliptic to equatorial coords
    local ecl := 23.439281*rads;    //value for J2000.0 frame
    local xeq := xg;
    local yeq := yg*cos(ecl) - zg*sin(ecl);
    local zeq := yg*sin(ecl) + zg*cos(ecl);
    
    // find the RA and DEC from the rectangular equatorial coords
    obj(ra)   := mod2pi( arg((xeq, yeq)) );
    obj(dec)  := atan(zeq / sqrt(xeq*xeq + yeq*yeq));
    obj(rvec) := sqrt(xeq*xeq + yeq*yeq + zeq*zeq);
    return obj;
end;

// Compute the elements of the orbit for planet-i at day number-d
// result is returned in structure p
mean_elements( p, i1, d )
begin

    local cy := d / 36525;          // centuries since J2000

    case
    if i1 == 1 then // Mercury
        p(ai) := 0.38709893 + 0.00000066*cy;
        p(ei) := 0.20563069 + 0.00002527*cy;
        p(ii) := ( 7.00487  -  23.51*cy/3600)*rads;
        p(Oi) := (48.33167  - 446.30*cy/3600)*rads;
        p(wi) := (77.45645  + 573.57*cy/3600)*rads;
        p(Li) := mod2pi((252.25084 + 538101628.29*cy/3600)*rads);
    end;
    if i1 == 2 then // Venus
        p(ai) := 0.72333199 + 0.00000092*cy;
        p(ei) := 0.00677323 - 0.00004938*cy;
        p(ii) := (  3.39471 -   2.86*cy/3600)*rads;
        p(Oi) := ( 76.68069 - 996.89*cy/3600)*rads;
        p(wi) := (131.53298 - 108.80*cy/3600)*rads;
        p(Li) := mod2pi((181.97973 + 210664136.06*cy/3600)*rads);
    end;
    if i1 ==  3 then // Earth/Sun
        p(ai) := 1.00000011 - 0.00000005*cy;
        p(ei) := 0.01671022 - 0.00003804*cy;
        p(ii) := (  0.00005 -    46.94*cy/3600)*rads;
        p(Oi) := (-11.26064 - 18228.25*cy/3600)*rads;
        p(wi) := (102.94719 +  1198.28*cy/3600)*rads;
        p(Li) := mod2pi((100.46435 + 129597740.63*cy/3600)*rads);
    end;
    if i1 == 4 then // Mars
        p(ai) := 1.52366231 - 0.00007221*cy;
        p(ei) := 0.09341233 + 0.00011902*cy;
        p(ii) := (  1.85061 -   25.47*cy/3600)*rads;
        p(Oi) := ( 49.57854 - 1020.19*cy/3600)*rads;
        p(wi) := (336.04084 + 1560.78*cy/3600)*rads;
        p(Li) := mod2pi((355.45332 + 68905103.78*cy/3600)*rads);
    end;
    if i1 == 5 then // Jupiter
        p(ai) := 5.20336301 + 0.00060737*cy;
        p(ei) := 0.04839266 - 0.00012880*cy;
        p(ii) := (  1.30530 -    4.15*cy/3600)*rads;
        p(Oi) := (100.55615 + 1217.17*cy/3600)*rads;
        p(wi) := ( 14.75385 +  839.93*cy/3600)*rads;
        p(Li) := mod2pi((34.40438 + 10925078.35*cy/3600)*rads);
    end;
    if i1 == 6 then // Saturn
        p(ai) := 9.53707032 - 0.00301530*cy;
        p(ei) := 0.05415060 - 0.00036762*cy;
        p(ii) := (  2.48446 +    6.11*cy/3600)*rads;
        p(Oi) := (113.71504 - 1591.05*cy/3600)*rads;
        p(wi) := ( 92.43194 - 1948.89*cy/3600)*rads;
        p(Li) := mod2pi((49.94432 + 4401052.95*cy/3600)*rads);
    end;
    if i1 == 7 then // Uranus
        p(ai) := 19.19126393 + 0.00152025*cy;
        p(ei) :=  0.04716771 - 0.00019150*cy;
        p(ii) := (  0.76986  -    2.09*cy/3600)*rads;
        p(Oi) := ( 74.22988  - 1681.40*cy/3600)*rads;
        p(wi) := (170.96424  + 1312.56*cy/3600)*rads;
        p(Li) := mod2pi((313.23218 + 1542547.79*cy/3600)*rads);
    end;
    if i1 == 8 then // Neptune
        p(ai) := 30.06896348 - 0.00125196*cy;
        p(ei) :=  0.00858587 + 0.00002510*cy;
        p(ii) := (  1.76917  -   3.64*cy/3600)*rads;
        p(Oi) := (131.72169  - 151.25*cy/3600)*rads;
        p(wi) := ( 44.97135  - 844.43*cy/3600)*rads;
        p(Li) := mod2pi((304.88003 + 786449.21*cy/3600)*rads);
    end;
    if i1 == 9 then // Pluto
        p(ai) := 39.48168677 - 0.00076912*cy;
        p(ei) :=  0.24880766 + 0.00006465*cy;
        p(ii) := ( 17.14175  +  11.07*cy/3600)*rads;
        p(Oi) := (110.30347  -  37.33*cy/3600)*rads;
        p(wi) := (224.06676  - 132.25*cy/3600)*rads;
        p(Li) := mod2pi((238.92881 + 522747.90*cy/3600)*rads);
    end;
    default
        //assert false
    end;
    return p;
end;

// compute the true anomaly from mean anomaly using iteration
//  M - mean anomaly in radians
//  e - orbit eccentricity
true_anomaly( M, e1 )
begin    

    local V, E, E1;

    // initial approximation of eccentric anomaly
    E := M + e1*sin(M)*(1.0 + e1*cos(M));

    // iterate to improve accuracy
    repeat
        E1 := E;
        E := E1 - (E1 - e1*sin(E1) - M)/(1 - e1*cos(E1));
    until abs( E - E1 ) < eps;

    // convert eccentric anomaly to true anomaly
    V := 2*atan(sqrt((1 + e1)/(1 - e1))*tan(0.5*E));

    // modulo 2pi
    if (V < 0) then
        V := V + (2*pi); 
    end;
    
    return V;

end;

// converts hour angle in degrees into hour angle string
ha2str( x )
begin

    // assert (0 <= x) and (x < 360)       // assure valid range
    local ra := x / 15;             // degrees to hours
    local h := floor( ra );
    local m := 60 * ( ra - h );
    return string( h ) + "h " + round( m,1 ) + "m";

end;

// converts declination angle in degrees into string
dec2str( x )
begin

    // assert (-90 <= x) and (x <= +90)     // assure valid range
    local dec := abs( x );           // work with absolute value
    local d := floor( dec );
    local m := 60 * ( dec - d );
    return string( sign(x)*d ) + "° " + round( m,1 ) + "'";

end;

// get object name
pname( i1 )
begin

    case
    if i1 == 1 then return "Mercury";end;
    if i1 == 2 then return "Venus  ";end;
    if i1 == 3 then return "Sun    ";end;
    if i1 == 4 then return "Mars   ";end;
    if i1 == 5 then return "Jupiter";end;
    if i1 == 6 then return "Saturn ";end;
    if i1 == 7 then return "Uranus ";end;
    if i1 == 8 then return "Neptune";end;
    if i1 == 9 then return "Pluto  ";end;
    end;

    return "not in range";

end;

// prompts for input of an int
get_int( prompt  )
begin

    local n;
    input(n, prompt, prompt);
    return n;

end;

// return the integer part of a number
abs_floor( x )
begin

    local r;

    if x >= 0.0 then
        r := floor( x );
    else
        r := ceiling( x );
    end;

    return r;

end;

// return an angle in the range 0 to 2pi
mod2pi( x )     // radians
begin    

    local b := x / (2*pi);
    local a := (2*pi)*(b - abs_floor(b));
    
    if a < 0 then
        a := (2*pi) + a; 
    end;
    
    return a;
    
end;

Tom L

DM42 SN: 00025 (Beta)
SN: 00221 (Production)
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)