+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: Planetary Positions (/thread-10757.html)
Planetary Positions - toml_12953 - 05-17-2018 02:46 PM
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;