Planetary Positions
05-17-2018, 02:46 PM (This post was last modified: 05-17-2018 02:49 PM by toml_12953.)
Post: #1
 toml_12953 Senior Member Posts: 1,138 Joined: Dec 2013
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

French phrases you rarely hear: "trop de beurre" (too much butter)
 « Next Oldest | Next Newest »

 Messages In This Thread Planetary Positions - toml_12953 - 05-17-2018 02:46 PM

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