06-20-2015, 05:53 PM
hi everybody,
waiting for the evolution of the very interesting app by Marcel (I'm already using!), I'm developing also this program to get astronomical results for Sun, Moon and Planets...
Effemeridi (Astrolabio) is a program that I'm adapted from an old version I made many years ago for a TI pocket calculator.
The program has still some incorrect output (very few), so I'm glad to have your advice and hints to improve it, thank you.
The program calculates Sidereal Time, Longitude, Latitude (various coordinates), Rect Ascension and Declination, Azimuth and Height, phases of the Moon, Retrogradation of planets, Ascendent, after a book of S. Bouiges...
I'm working to add eclipses, and improving program also in "astrological" sense (I was asked for this) and Satellites of Jupiter and Saturn.
Also data for Pluto are still not so precise as those in the app of Marcel...
Some routines and data are still in Italian, sorry.
I need your help, please, test, and thank you in advance!
Salvo
The code (I'll post here the new versions for now):
waiting for the evolution of the very interesting app by Marcel (I'm already using!), I'm developing also this program to get astronomical results for Sun, Moon and Planets...
Effemeridi (Astrolabio) is a program that I'm adapted from an old version I made many years ago for a TI pocket calculator.
The program has still some incorrect output (very few), so I'm glad to have your advice and hints to improve it, thank you.
The program calculates Sidereal Time, Longitude, Latitude (various coordinates), Rect Ascension and Declination, Azimuth and Height, phases of the Moon, Retrogradation of planets, Ascendent, after a book of S. Bouiges...
I'm working to add eclipses, and improving program also in "astrological" sense (I was asked for this) and Satellites of Jupiter and Saturn.
Also data for Pluto are still not so precise as those in the app of Marcel...
Some routines and data are still in Italian, sorry.
I need your help, please, test, and thank you in advance!
Salvo
The code (I'll post here the new versions for now):
Code:
data();
calcN();
calcTS();
sunCalc();
planetCalc();
planets();
sun();
moon();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendent();
zodiaco();
trunc_angle();
atan2();
thisday:= 2000.0101;
lstd:=0; gmt:= 12;
long:= 0; lat:= 0;
m:= 1; dst:= 0; utc:=0;
sunlist:= MAKELIST(X,X,1,25);
planet_const:= MAKELIST(X,X,1,10);
ε:= →HMS(23.4392911); // approx 23.45
EXPORT Effemeridi()
// from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015
// Credits: from a book, Calcule astronomique pour amateurs, by Serge Bouiges
// Thanks Didier Lachieze, Marcel, Eried
BEGIN
LOCAL ch;
HAngle:=1;
data(); // input data and calc Sun
CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc N", "Calc TS", "Ascendent", "Quit");
CASE
IF ch==1 THEN sun(); END;
IF ch==2 THEN moon(); END;
IF ch==3 THEN planets(); END;
IF ch==4 THEN calcN(); END;
IF ch==5 THEN calcTS(); END;
IF ch==6 THEN ascendent(); END;
IF ch==7 THEN HAngle:=0; RETURN; END;
DEFAULT
END; // case
END;
data()
BEGIN
LOCAL dd, mm, yy, hh, mn, ss, loct;
yy:= IP(Date);
mm:= IP(FP(Date)*100);
dd:= FP(Date*100)*100;
hh:= IP(HMS→(Time));
mn:= IP(FP(HMS→(Time))*60);
ss:= IP(FP(FP(FP(HMS→(Time))*60))*60);
loct:= →HMS(hh+mn/60+ss/3600);
INPUT({ {m,[0],{15,15,1}}, {D,[0],{50,15,1}}, {Y,[0],{80,15,1}},
{lstd,[0],{65,30,2}},
{dst,0,{40,15,3}}, {utc,[0],{80,15,3}},
{long,[0],{20,25,4}},
{lat,[0],{70,25,4}} },
"Data: Use Shift+9 for H°M′S″",
{"Month:","Day :","Year :","Local Time (24):", "DST", "GMT",
"Long (-E):","Lat (+N)"},
{"Enter month", "Enter day", "Enter year (1901-2099)",
"Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC lag",
"Longitude", "Latitude"},
{1,1,1901,0, 0,0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″});
// Ragusa (-14°43′41″, 36°55′15″) -
// Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)
D:=EVAL(D); m:=EVAL(m); Y:=EVAL(Y);
lstd:= lstd;
// Greenwich Mean Time
gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)
thisday:= Y+m/100+D/10000;
sunCalc(); // calculate Sun data
END;
planets()
BEGIN
LOCAL ch, nameplanet;
INPUT({{ch,{"Mercurius","Venus","Mars", "Jupiter", "Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},
"Choose planet",
"planet: ", "Choose the planet to calculate an press OK", 0,5 );
CASE
IF ch==1 THEN nameplanet:="Mercurius"; planetCalc(nameplanet); END;
IF ch==2 THEN nameplanet:="Venus"; planetCalc(nameplanet); END;
IF ch==3 THEN nameplanet:="Mars"; planetCalc(nameplanet); END;
IF ch==4 THEN nameplanet:="Jupiter"; planetCalc(nameplanet); END;
IF ch==5 THEN nameplanet:="Saturn"; planetCalc(nameplanet); END;
IF ch==6 THEN nameplanet:="Uranus"; planetCalc(nameplanet); END;
IF ch==7 THEN nameplanet:="Neptunus"; planetCalc(nameplanet); END;
IF ch==8 THEN nameplanet:="Pluto"; planetCalc(nameplanet); END;
DEFAULT
END; // case
END;
calcN()
BEGIN
N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT
RETURN N;
END;
calcTS()
BEGIN
LOCAL TS, TSd, N0, TSL;
calcN();
N0:= DDAYS(1901.0101, thisday)+1+0; // N at 0 GMT (for TS)
TS:= ((23750.3 + 236.555362*(N0+gmt/24)) / 3600) MOD 24; // Sideral Time (GMT) in seconds
TS:= →HMS(TS); // in HMS
TSL:= (((23750.3 + 236.555362*N+lstd/24) / 3600) MOD 24); // local ST
TSL:= (TSL + gmt - 4*(HMS→(long))/60) MOD 24;
TSL:= →HMS(TSL); // in HMS
TSd:= →HMS((98.965 + 0.985647342*N0) MOD 360); // ST at 0 GMT in degrees
RETURN {TS, TSL, TSd, gmt};
END;
ascendent()
BEGIN
LOCAL ASC, ts, segno;
ts:= calcTS();
ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));
IF ASC<0 THEN ASC:= ASC+360; END; IF ASC>360 THEN ASC:=ASC-360; END;
segno:= zodiaco(ASC);
RETURN {→HMS(ASC), segno};
END;
zodiaco(astro)
BEGIN
LOCAL segno;
CASE // Segni e case
IF astro>=0 AND astro <30 THEN segno:="Ariete"; END;
IF astro>=30 AND astro <60 THEN segno:="Toro"; END;
IF astro>=60 AND astro <90 THEN segno:="Gemelli"; END;
IF astro>=90 AND astro <120 THEN segno:="Cancro"; END;
IF astro>=120 AND astro <150 THEN segno:="Leone"; END;
IF astro>=150 AND astro <180 THEN segno:="Vergine"; END;
IF astro>=180 AND astro <210 THEN segno:="Bilancia"; END;
IF astro>=210 AND astro <240 THEN segno:="Scorpione"; END;
IF astro>=240 AND astro <270 THEN segno:="Sagittario"; END;
IF astro>=270 AND astro <300 THEN segno:="Capricorno"; END;
IF astro>=300 AND astro <330 THEN segno:="Acquario"; END;
IF astro>=330 AND astro <360 THEN segno:="Pesci"; END;
END; // case
RETURN segno;
END;
atan2(y,x)
BEGIN
// atan2 returns coordinates in correct angle for all four quadrants (DEG)
CASE
IF x==0 AND y==0 then return "undefined"; end; // x=0, y=0
IF x>0 then return atan(y/x); end; // x>0
IF x<0 AND y>=0 then return atan(y/x)+180; end; // x<0, y>=0
IF x==0 AND y>0 then return 90; end; // x=0, y>0
IF x<0 AND y<=0 then return atan(y/x)-180; end; // x<0, y<0
IF x==0 AND y<0 then return -90; end; // x=0, y<0
END;
RETURN;
END;
trunc_angle(ang)
// trunc decimal from a DEG form angle
BEGIN
→HMS((IP(HMS→(ang)*3600)/3600));
END;
sunCalc()
BEGIN
LOCAL ω, l, v, r, x, y, dz;
LOCAL α, δ, dayY, H, ts, az, z, h;
LOCAL SAD, DA, sorg, tram, culm, Hd, lt;
LOCAL as, ac, zc;
calcN();
L:= →HMS((278.965+0.985647342*N)) MOD 360; // long media
ω:= →HMS((281.235 + 0.0000469*N)) MOD 360; // long del perielio
M:= (L - ω); // anomalia vera
v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio
l:= v + ω; // Longitudine Sole
lt:= l + 180; // Longitudine eliocentrica della Terra
r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)
x:= r * cos(l); y:= r * sin(l);
dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number
// δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper
δ:= →HMS((asin(sin(ε)*sin(l))) );
// α:= →HMS((acos(cos(l)/cos(δ))/15) MOD 24); // ascensione retta
α:= (acos(cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta
IF l>180 THEN α:= α +180; END;
α:= →HMS(α/15) MOD 24; // α in hms
ts:= calcTS(); // calcola tempo siderale
H:= (ts(2) - α) MOD 24; // Angolo orario
Hd:= →HMS(15*H); // Angolo orario in gradi
zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt
dz:= →HMS(acos(zc)) ; // distanza zenitale
// alt:= →HMS(90 - z); // Altezza
h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza
// as:= cos(δ)*sin(Hd)/sin(dz); // calc az
// ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);
// az:= →HMS(ATAN(as/ac) ); // Azimuth
// IF (ac<0) THEN az:= (360 + az) MOD 360; END;
az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));
IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;
IF (az<0) THEN az:= 360+az; END;
E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)
DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
SAD:= (90 + DA); // Semiarco Diurno
// SAD:= →HMS(acos(-tan(lat)*tan(δ))); // Semiarco Diurno
culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione
sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)
sunlist:= {L,ω,l, M, α,δ, az, h, x,y,sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA, SAD, E, lt};
RETURN sunlist;
END;
sun()
BEGIN
LOCAL ω, l, v, r, x, y, dz;
LOCAL α, δ, dayY, H, ts, az, z, h;
LOCAL SAD, DA, sorg, tram, culm, Hd, lt;
LOCAL as, ac, zc, segno;
L:= sunlist(1); ω:= sunlist(2);
l:= sunlist(3); M:= sunlist(4);
α:= sunlist(5); δ:= sunlist(6);
az:= sunlist(7); h:= sunlist(8);
x:= sunlist(9); y:= sunlist(10);
sorg:= sunlist(11); tram:= sunlist(12); culm:= sunlist(13);
r:= sunlist(14); H:= sunlist(15); Hd:= sunlist(16);
zc:= sunlist(17); dz:= sunlist(18);
ts:= calcTS(); // no 19
DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);
segno:= zodiaco(l);
PRINT;
PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("");
PRINT("L longitudine media " + trunc_angle(L));
PRINT("ω long del perielio " + trunc_angle(ω));
PRINT("M anomalia vera " + trunc_angle(M));
PRINT("λ longitudine Sole " + trunc_angle(l) + " " + segno);
PRINT("r Distanza (UA) " + r);
PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
PRINT(" ");
PRINT("Ang. orario "+ H + " " + trunc_angle(Hd));
PRINT("TS (UTC) " + trunc_angle(ts(1)));
PRINT("TS locale " + trunc_angle(ts(2)));
PRINT("Equazione del tempo (min) " + trunc_angle(E));
PRINT("Dif asc." + DA + " SAD " + trunc_angle(SAD));
PRINT("");
PRINT("Culmina a: " + trunc_angle(culm) + " UTC");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");
PRINT("Long elioc. Terra "+ trunc_angle(lt));
// RETURN {L, ω, M, v, l, r, x, y};
RETURN [[L, ω], [l, M], [α,δ], [az, h], [x,y], [sorg, tram], [culm, r], [H,ts(1)]];
END;
moon()
BEGIN
LOCAL ω,Ω, l, v, r, x, y, i;
LOCAL Ls, Ms, LongM, LatM, P, d, s;
LOCAL α,δ, ts, Hd, zc, dz, h, az;
LOCAL as, ac, DA, SAD, DD, culm, sorg, tram;
LOCAL b, dayY, xs, ys, ph, fase, eta;
LOCAL dd, hh, mm, segno;
i:= 5.13333333; // inclinazione 5°8' sull'eclittica
calcN();
// Sun
Ls:= sunlist(1); // Longitudine media
Ms:= sunlist(4); // Anomalia
// end Sun
L:= (33.231 + 13.17639653*N) MOD 360; // long. media
Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente
M:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera
ω:= (L - M) MOD 360; // perigeo medio
DD:= L - Ls; F:= L-Ω; // for the variations
LongM:= L + 6.28875 * sin(M); // calculus (h/ and above) of Moon longitude
LongM:= LongM + 0.2136 * sin(2*M); // equazione del centro
LongM:= LongM + 0.6583 * sin(2*DD); // variazione
LongM:= LongM - 0.1856 * sin(Ms); // equazione annuale
LongM:= LongM + 1.2740 * sin(2*DD-M); // evezione
LongM:= LongM - 0.1143 * sin(2*F); // riduzione all'eclittica
LongM:= LongM + 0.0588 * sin(2*DD-2*M);
LongM:= LongM + 0.0572 * sin(2*DD-M-Ms);
LongM:= LongM + 0.0533 * sin(2*DD+M);
LongM:= LongM + 0.0459 * sin(2*DD-Ms);
LongM:= LongM + 0.0410 * sin(M-Ms);
LongM:= LongM - 0.0305 * sin(M+Ms);
LongM:= LongM - 0.0348 * sin(DD);
LongM:= LongM MOD 360;
LatM := 5.1280 * sin(F); // termine principale
LatM := LatM + 0.2806 * sin(M+F);
LatM := LatM + 0.2777 * sin(M-F); // equazione del centro
LatM := LatM + 0.1732 * sin(2*DD-F); // grande ineguaglianza
LatM := LatM + 0.0554 * sin(2*DD-M+F); // evezione
LatM := LatM + 0.0462 * sin(2*DD-M-F); // evezione (2)
LatM := LatM + 0.0326 * sin(2*DD+F);
LatM := LatM MOD 360;
b:= LatM;
dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number
δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(LongM))) ); // declinazione (coord locali)
// α:= →HMS((acos(cos(b)*cos(LongM)/cos(δ))/15) MOD 24); // ascensione retta
α:= (acos(cos(b)*cos(IFTE(LongM>180, LongM-180, LongM))/cos(δ))); // ascensione retta
IF LongM>180 THEN α:= α +180; END;
α:= →HMS(α/15) MOD 24; // α in hms
ts:= calcTS(); // calcola tempo siderale
H:= (ts(2) - α) MOD 24; // Angolo orario
Hd:= →HMS(15*H); // Angolo orario in gradi
zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt
dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza
h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)
as:= cos(δ)*sin(Hd)/sin(dz); // calc az
ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);
// az:= →HMS(ATAN(as/ac) ); // Azimuth
// IF (ac<0) THEN az:= ((360 + az) MOD 360); END;
az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));
IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;
DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
SAD:= (90 + DA); // Semiarco Diurno
// SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione
sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)
P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse
// Raggio quadratico medio della Terra: 6373,044737
d:= 1/sin(P); // distanza Terra-Luna in raggi terrestri
s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)
ph:= abs(LongM - Ls); // Phase
CASE
// IF ph == 0 THEN fase:="New Moon"; END;
IF ph>280 AND ph<350 THEN fase:="Decrescente"; END;
IF ph>=260 AND ph<=280 THEN fase:="Last Quart"; END;
IF ph>190 AND ph<260 THEN fase:="Gibbosa Cal."; END;
IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
IF ph>100 AND ph<170 THEN fase:="Gibbosa Cresc."; END;
IF ph>=80 AND ph<=100 THEN fase:="1st Quart"; END;
IF ph>10 AND ph<80 THEN fase:="Crescente"; END;
DEFAULT fase:="New Moon";
END; // case
// età della Luna
dd:= IP(HMS→(ph/15));
hh:= IP(FP(HMS→(ph/15))*60);
IF (hh>24) THEN dd:= dd+1; hh:= hh MOD 24; END;
mm:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);
eta:= →HMS(dd+hh/60+mm/3600); // Age of the Moon
segno:= zodiaco(LongM);
PRINT;
PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("L Long media " + trunc_angle(L));
PRINT("M anomalia vera " + trunc_angle(M));
PRINT("λ Longitudine " + trunc_angle(LongM) + " " + segno);
PRINT("β Latitudine " + trunc_angle(LatM));
PRINT("ω Perigeo medio " + trunc_angle(ω));
PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
PRINT("α Ascensione retta (gradi) " + trunc_angle((α*15) MOD 360));
PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
PRINT("");
PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));
PRINT("TS (UTC) " + trunc_angle(ts(1)));
PRINT("TS (locale) " + trunc_angle(ts(2)));
PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT("P Parallasse " + trunc_angle(P));
PRINT("d Distanza Terra - Luna "+ ROUND(d*6373.044737,3) + " km");
PRINT("Semidiametro "+ trunc_angle(→HMS(s)) + " :: " + ROUND(s*6373.044737,3) + " km");
PRINT(" ");
PRINT("Culmina a: " + trunc_angle(culm) + " UTC");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");
PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");
PRINT("Età " + eta);
RETURN [[L, M],[LongM, LatM],[ω,Ω],[α,δ], [az, h], [P,0],[d,s], [sorg, tram], [culm, ts(1)],[ph, eta] ];
END;
planetCalc(nameplanet)
BEGIN
LOCAL ω,Ω, l, v, r, x, y, z;
LOCAL Ls, Ms, b, xs, ys, b_ec;
LOCAL x1, y1, z1, R, l_ec, H, Hd;
LOCAL ec, in, a, δ, α, ts, az, h;
LOCAL DA, SAD, tram, sorg, culm, dayY;
LOCAL as, ac, zc, dz, θ, mov, segno;
// calc of Sun (as basis)
xs:= sunlist(9); ys:= sunlist(10);
Ls:= sunlist(1); // Longitudine media
Ms:= sunlist(4); // Anomalia
// end Sun
EVAL(EXPR(nameplanet+"()")); // call planet-name function
L:= planet_const(1); // Longitudine media
ω:= planet_const(2); // perielio
Ω:= planet_const(3); // nodo
M:= planet_const(4); // anomalia
v:= planet_const(5);
ec:= planet_const(6); // eccentricità
in:= planet_const(7); // inclinazione
a:= planet_const(8); // semiasse maggiore
θ:= planet_const(9); // θs to calc if retrogade
b_ec:= asin(sin(in)*sin(v+ω-Ω)); // argomento di latitudine del perielio
l_ec:= acos(cos(v+ω-Ω)/cos(b_ec)); // long from eclittic
IF b_ec<0 THEN l_ec:= 360 - l_ec; END; // argomento secondario
l_ec:= →HMS((l_ec + Ω)) MOD 360; // aggiunge Nodo e riporta a 2PI
r:= a*(1-ec^2)/(1+ec*cos(v)); // 9.524899/(1+ 0.0559*cos(v)); raggio vettore
x1:= r*cos(b_ec)*cos(l_ec); y1:= r*cos(b_ec)*sin(l_ec); z1:=r*sin(b_ec); // cartesian
x:= x1+xs; y:= y1+ys; z:=z1; // cohordinates cartesian->spheric
R:= sqrt(x^2+y^2+z^2);
b:=asin(z/R); // latitudine β
// l:=→HMS(atan(y/x)); // cart->spheric longtitudine λ
l:= atan2(y,x);
IF l<0 THEN l:= l+360; END; IF l>360 THEN l:=l-360; END;
dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number
// δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper
δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(l))) ); // declinazione (coord locali)
α:= (acos(cos(b)*cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta
IF l>180 THEN α:= α +180; END;
α:= →HMS(α/15) MOD 24; // α in hms
ts:= calcTS(); // calcola tempo siderale
H:= (ts(2) - α) MOD 24; // Angolo orario
Hd:= →HMS(15*H); // Angolo orario in gradi
zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt
dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza
h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)
as:= cos(δ)*sin(Hd)/sin(dz); // calc az
ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);
// az:= →HMS(ATAN(as/ac) ); // Azimuth
az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));
IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;
DA:= →HMS(asin(tan(lat)*tan(δ))); // Differenza ascensionale
SAD:= →HMS(90 + DA); // Semiarco Diurno
// SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) MOD 24 ; // Culminazione
sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)
IF abs(l-220.4)>θ THEN mov:="diretto"; ELSE mov:="retrogrado"; END;
segno:= zodiaco(l);
PRINT;
PRINT("Effemeridi di " + nameplanet + " a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("L Long media " + trunc_angle(L));
PRINT("ω Perielio " + trunc_angle(ω));
PRINT("Ω Nodo " + trunc_angle(Ω));
PRINT("M Anomalia vera " + trunc_angle(M));
PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));
PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(v));
PRINT("r Raggio vettore " + r);
PRINT("Cartes. elioc.: x' " + ROUND(x1,3) + " y' " + ROUND(y1,3) + " z' " + ROUND(z1,3));
PRINT("Cartes. geoc.: x " + ROUND(x,3) + " y " + ROUND(y,3) + " z " + ROUND(z,3));
PRINT("Spher.: λ Longitudine " + trunc_angle(l) + " " + segno); // spheric
PRINT("Spher.: β Latitudine " + trunc_angle(b));
PRINT("Spher.: R Distanza " + ROUND(R,3));
PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
PRINT("");
PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));
PRINT("TS (UTC) " + trunc_angle(ts(1)));
PRINT("TS (locale) " + trunc_angle(ts(2)));
PRINT("Diff. ascensionale " + trunc_angle(DA));
PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));
PRINT("Culmina a: " + trunc_angle(culm) + " UTC");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");
PRINT("Movimento " + mov);
RETURN [[L,ω,Ω ],[M,l_ec,b_ec], [v,r,0], [x1,y1,z1],[x,y,z],[l,b,R], [α,δ,H], [az,h,ts(1)],[sorg, tram, culm]];
END;
Mercurius()
BEGIN
// Costanti per Mercurio
LOCAL ω,Ω, L, M, v;
LOCAL ec, in, a, θ;
ec:= 0.205615; // eccentricità
in:= 7.003; // inclinazione
a:= 0.387098; // semiasse maggiore
θ:= 35.6; // θs to calc if retrogade
calcN();
// costanti planetarie
L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media
ω:= →HMS((75.913 + 0.00004253*N) MOD 360); // perielio
Ω:= →HMS((47.157 + 0.00003244*N) MOD 360); // nodo
M:= (L - ω); // Anomalia vera
v:= M + 23.4372*sin(M)+2.9810*sin(2*M)+0.5396*sin(3*M)+0.1098*sin(4*M);
planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
RETURN planet_const;
END;
Venus()
BEGIN
// Costanti per Venere
LOCAL ω,Ω, L, M, v;
LOCAL ec, in, a, θ;
ec:= 0.006816; // eccentricità
in:= 3.393636; // inclinazione
a:= 0.72333; // semiasse maggiore
θ:= 13.0; // θs to calc if retrogade
calcN();
// costanti planetarie
L:= →HMS((206.758 + 1.60216873*N) MOD 360); // Longitudine media
ω:= →HMS((130.154 + 0.00003757*N) MOD 360); // perielio
Ω:= →HMS((75.797 + 0.00002503*N) MOD 360); // nodo
M:= (L - ω); // Anomalia vera
v:= M + 0.7811*sin(M);
planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
RETURN planet_const;
END;
Mars()
BEGIN
// Costanti per Marte
LOCAL ω,Ω, L, M, v;
LOCAL ec, in, a, θ;
ec:= 0.093309; // eccentricità
in:= 1.850303; // inclinazione
a:= 1.523678; // semiasse maggiore
θ:= 16.8; // θs to calc if retrogade
calcN();
// costanti planetarie
L:= →HMS((124.765 + 0.52407108*N) MOD 360); // Longitudine media
ω:= →HMS((334.251 + 0.00005038*N) MOD 360); // perielio
Ω:= →HMS((48.794 + 0.00002127*N) MOD 360); // nodo
M:= (L - ω); // Anomalia vera
v:= M + 10.608*sin(M)+0.6216*sin(2*M)+0.0504*sin(3*M)+0.0047*sin(4*M);
planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
RETURN planet_const;
END;
Jupiter()
BEGIN
// Costanti per Giove
LOCAL ω,Ω, L, M, v;
LOCAL ec, in, a, θ;
ec:= 0.048335; // eccentricità
in:= 1.308736; // inclinazione
a:= 5.202561; // semiasse maggiore
θ:= 54.43; // θs to calc if retrogade
calcN();
// costanti planetarie
L:= →HMS((268.350 + 0.08312942*N) MOD 360); // Longitudine media
ω:= →HMS((12.737 + 0.00004408*N) MOD 360); // perielio
Ω:= →HMS((99.453 + 0.00002767*N) MOD 360); // nodo
M:= (L - ω); // Anomalia vera
v:= M + 5372*sin(M)+0.1662*sin(2*M)+0.007*sin(3*M);
planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
RETURN planet_const;
END;
Saturn()
BEGIN
// Costanti per Saturno
LOCAL ω,Ω, L, M, v;
LOCAL ec, in, a, θ;
ec:= 0.055892; // eccentricità
in:= 2.492519; // inclinazione
a:= 9.554747; // semiasse maggiore
θ:= 65.53; // θs to calc if retrogade
calcN();
// costanti planetarie
L:= →HMS((278.774 + 0.03349788*N) MOD 360); // Longitudine media
ω:= →HMS((91.117 + 0.00005362*N) MOD 360); // perielio
Ω:= →HMS((112.79 + 0.00002391*N) MOD 360); // nodo
M:= (L - ω); // Anomalia vera
v:= M + 6.4023*sin(M)+0.2235*sin(2*M)+0.011*sin(3*M);
planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
RETURN planet_const;
END;
Uranus()
BEGIN
// Costanti per Urano
LOCAL ω,Ω, L, M, v;
LOCAL ec, in, a, θ;
ec:= 0.046344; // eccentricità
in:= 0.772464; // inclinazione
a:= 19.21814; // semiasse maggiore
θ:= 73.92; // θs to calc if retrogade
calcN();
// costanti planetarie
L:= →HMS((248.487 + 0.01176902*N) MOD 360); // Longitudine media
ω:= →HMS((171.563 + 0.00004064*N) MOD 360); // perielio
Ω:= →HMS((73.482 + 0.00001365*N) MOD 360); // nodo
M:= (L - ω); // Anomalia vera
v:= M + 5.3092*sin(M)+0.1537*sin(2*M)+0.0062*sin(3*M);
planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
RETURN planet_const;
END;
Neptunus()
BEGIN
// Costanti per Nettuno
LOCAL ω,Ω, L, M, v;
LOCAL ec, in, a, θ;
ec:= 0.008997; // eccentricità
in:= 1.779242; // inclinazione
a:= 30.10957; // semiasse maggiore
θ:= 77.63; // θs to calc if retrogade
calcN();
// costanti planetarie
L:= →HMS((86.652 + 0.00602015*N) MOD 360); // Longitudine media
ω:= →HMS((46.742 + 0.000039*N) MOD 360); // perielio
Ω:= →HMS((130.693 + 0.00003009*N) MOD 360); // nodo
M:= (L - ω); // Anomalia vera
v:= M + 1.031*sin(M)+0.0058*sin(2*M);
planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
RETURN planet_const;
END;
Pluto()
BEGIN
// Costanti per Plutone
LOCAL ω,Ω, L, M, v;
LOCAL ec, in, a, θ;
ec:= 0.250236; // eccentricità
in:= 17.17047; // inclinazione
a:= 39.438712; // semiasse maggiore
θ:= 79.41; // θs to calc if retrogade
calcN();
// costanti planetarie
// Long, perielio, nodo, non precisi, considerati costanti
// L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media
M:= →HMS((230.67 + 0.00397943*N) MOD 360); // Anomalia vera (calcolata direttamente)
Ω:= →HMS((109.06 + 0.00003823*N) MOD 360); // nodo
ω:= →HMS((114.27 + Ω) MOD 360); // perielio
L:= →HMS((M + ω) MOD 360); // L (non sicuro)
v:= M + 28.45*sin(M)+4.38*sin(2*M)+0.97*sin(3*M)+0.24*sin(4*M);
planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
RETURN planet_const;
END;