+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: HP Prime (/forum-5.html)
+--- Thread: Astronomy... (/thread-4195.html)
Astronomy... - salvomic - 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):
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;RE: Astronomy... - Marcel - 06-20-2015 07:18 PM
Hi Salvomic
We don't use the same technique for calculating the planetary position. I use the VSOP technique (Variation séculaire des orbites planétaires). In this technique, we use polynomials and series of DATA to obtain the heliocentric position of a planet...
I complete today two new functions for the location of Pluto. This is the year of Pluto because in 25 days, the probe New Horizon will arrive at Pluto. For this reason, I will post my astronomy app's just before that moment.
Good lock in your project, this is fun and we can use both programs.
Marcel
RE: Astronomy... - salvomic - 06-20-2015 08:34 PM
(06-20-2015 07:18 PM)Marcel Wrote: Hi Salvomic
We don't use the same technique for calculating the planetary position. I use the VSOP technique (Variation séculaire des orbites planétaires). In this technique, we use polynomials and series of DATA to obtain the heliocentric position of a planet...
I complete today two new functions for the location of Pluto. This is the year of Pluto because in 25 days, the probe New Horizon will arrive at Pluto. For this reason, I will post my astronomy app's just before that moment.
Good lock in your project, this is fun and we can use both programs.
Marcel
hi Marcel,
yes, I'm using a different approach, following Bouiges: using the constants for planets, Sun and Moon...
About Pluto, my approach is very few, for now, as I'm using the old data (they were written in 1973...) and then Pluto was not completely known about Longitude, for example, so my program for now use some values as constant, like perigeo or longitude and anomaly...
You are right: we could now give more attention to Pluto.
However we could perhaps to search if some functions are common to the two programs and to find a common point of view...
Beside Pluto, in your program, I'm looking forward to see also the Moon functions, phases and so on...
Salvo
RE: Astronomy... - Marcel - 06-20-2015 10:29 PM
Hi,
I send you my app's with Pluto, so you will be able to compare with your program. If you have Mathematica or WolframAlpha, check for precision. It is very good.
Marcel
RE: Astronomy... - salvomic - 06-21-2015 06:38 AM
(06-20-2015 10:29 PM)Marcel Wrote: Hi,
I send you my app's with Pluto, so you will be able to compare with your program. If you have Mathematica or WolframAlpha, check for precision. It is very good.
Marcel
thank Marcel,
I'm trying it!
for now with my program I get the same declination, but not the same RA, but my program has a problem to get right longitude of Pluto for now, your does it better...
I've both Mathematica and Wolfram, I'll check there also.
Salvo
EDIT: I reduced difference for Pluto between the to programs about (1' - 3'), for now: there was an error in a formula of mine...
RE: Astronomy... - salvomic - 06-22-2015 12:42 PM
some improvements (calculation of times, sidereal times, julian day, epsilon, data for Pluto...)
Code:
data();
calcN();
calcTS();
julianDay();
tempo();
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;
deltaT:= 68;
datetimelist:= MAKELIST(X,X,1,6);
sunlist:= MAKELIST(X,X,1,25);
planet_const:= MAKELIST(X,X,1,10);
ε:= →HMS(23.4392911); // approx 23.45 (23°26'21")
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 TS JD N", "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 tempo(); END;
IF ch==5 THEN ascendent(); END;
IF ch==6 THEN HAngle:=0; RETURN; END;
DEFAULT
END; // case
END;
data()
BEGIN
LOCAL dd, mm, yy, hh, mn, ss, loct;
LOCAL t, jd;
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}},
{deltaT, [0], {20,25,5}} },
"Data: Use Shift+9 for H°M′S″",
{"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",
"Long (-E):","Lat (+N)", "delta T"},
{"Enter month", "Enter day", "Enter year (1901-2099)",
"Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",
"Longitude", "Latitude", "Delta T (sec.)"},
{1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});
// 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);
datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };
// Greenwich Mean Time
gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)
thisday:= Y+m/100+D/10000;
jd:= julianDay();
t:=(jd(1)-2451545)/36525;
ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3;
ε:= →HMS(ε);
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, TSapp;
calcN();
N0:= DDAYS(1901.0101, thisday)+1+0; // N at 0 GMT (for TS)
TS:= ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds
TS:= →HMS(TS); // in HMS
TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST
TSL:= →HMS(TSL); // in HMS
TSapp:= (TS + gmt) MOD 24; // apparent ST (ST at Greenwitch at local our)
TSapp:= →HMS(TSapp); // in HMS
TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees
RETURN {trunc_angle(TS), trunc_angle(TSL), trunc_angle(TSapp), trunc_angle(TSd), gmt};
END;
julianDay()
BEGIN
LOCAL yy, mm, dd,dde, hh, mn, ss;
LOCAL aa,bb, jd, jde;
yy:= Y; mm:= m;
hh:= datetimelist(4);
mn:= datetimelist(5);
ss:= datetimelist(6);
dd:= D + hh/24 + mn/(60*24) + ss/(60*60*24);
dde:= D + hh/24 + mn/(60*24) + (ss+deltaT)/(60*60*24);
IF (m=1 OR m=2) THEN yy:=Y-1; mm:=m+12;END;
IF yy*100+mm+dd>=158225 THEN
aa:=IP(ABS(yy/100));
bb:=2-aa+IP(aa/4);
ELSE bb:=0; END;
jd:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dd+bb-1524.5;
jde:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dde+bb-1524.5;
RETURN {jd, jde};
END;
tempo()
BEGIN
LOCAL n, jd, TSid;
n:= calcN();
jd:= julianDay();
TSid:= calcTS;
PRINT;
PRINT("Date " + m + " " + D + " " + Y);
PRINT("Julian Day " + jd(1));
PRINT("Julian Day Effem. " + jd(2));
PRINT("N (from 1901jan1) " + n);
PRINT("Tempo Siderale 0h " + trunc_angle(TSid(1)));
PRINT("Tempo Siderale locale " + trunc_angle(TSid(2)));
PRINT("Tempo Siderale apparente " + trunc_angle(TSid(3)));
PRINT("Tempo Siderale gradi " + trunc_angle(TSid(4)));
PRINT("Time " + lstd + " GMT " + gmt);
RETURN({ROUND(jd(1),5), ROUND(jd(2),5), ROUND(n, 5), trunc_angle(TSid(1)), trunc_angle(TSid(2)),
trunc_angle(TSid(3)), trunc_angle(TSid(4)), 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, h0, H0, temp;
h0:= -0°50′; // 'standard' altitude of Sun
calcN();
julianDay();
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
// temp:= sin(h0)-sin(lat)*sin(δ)/cos(lat)*cos(δ);
// IF temp<0 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END;
// H0:= acos(temp);
// MSGBOX("H0: " + H0);
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)
//temp:= ((α*15 + long - ts(3)*15) / 360);
// IF temp<0 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END;
// culm:= (temp - Hd/360) *24;
// sorg:= (culm - H0/360) *24;
// tram:= (culm + H0/360) *24;
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, jd;
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(); //
jd:= julianDay();
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("at Long " + long + "Lat " + lat);
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("TS apparente " + trunc_angle(ts(3)));
PRINT("Equazione del tempo (min) " + trunc_angle(E));
PRINT("Julian day " + jd(1));
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));
PRINT("ε epsilon " + ε);
// RETURN {L, ω, M, v, l, r, x, y};
RETURN [[L, ω], [l, M], [α,δ], [az, h], [x,y], [sorg, tram], [culm, r], [H,ts(1)], [lt, ε]];
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;RE: Astronomy... - Marcel - 06-22-2015 06:26 PM
Hi Salvomic
We can calculate Delta T with some analytical way but this is not always precise. In fact we know the value of Delta T after the year. For instance, we will know the value for the year 2015, in 2016! This is a correction or a difference between two times scales.
In my app's, many data are in variable associate with the App. You can see these var if you click on the "var" tab in the conn kit. At this point this is more complex than your technique because many calculations most been done with these data in a FOR statement. If you check for the book of Jean Meeus on the internet, you will learn how to use that data.
Marcel
RE: Astronomy... - salvomic - 06-22-2015 06:34 PM
(06-22-2015 06:26 PM)Marcel Wrote: Hi Salvomic
We can calculate Delta T with some analytical way but this is not always precise. In fact we know the value of Delta T after the year. For instance, we will know the value for the year 2015, in 2016! This is a correction or a difference between two times scales.
In my app's, many data are in variable associate with the App. You can see these var if you click on the "var" tab in the conn kit. At this point this is more complex than your technique because many calculations most been done with these data in a FOR statement. If you check for the book of Jean Meeus on the internet, you will learn how to use that data.
Marcel
hi Marcel,
thank you for info.
I found the book of Meeus, yes, but just few minutes ago...
But I should rewrite all the program to conform 100% to it, and in this case there is already your app...
So, I'm trying to improve something adding formulas by Meeus to those by Bouiges. It isn't an easy job, hi...
Now, to improve Sun precision I should add at least 15 parameters, and solve the Keplero equation:
u - e sin(u) = M
in the var u, given e (eccentricity, by formula, ok) and M (anomaly), stopping it when the last iteration is no more great than previous for 1/10000, then calc v with u (v=acos(cos(u)-e/1-e*cos(u))) and so on...
Any hint to solve that equation?
I tried solve() but without luck...
thank you
Salvo
RE: Astronomy... - Marcel - 06-22-2015 06:48 PM
Hi,
One of the chapter in the book describe a way to solve the Kepler equation (Chapter 29). You will find in this chapter a program to do it.
Marcel
RE: Astronomy... - salvomic - 06-22-2015 07:29 PM
(06-22-2015 06:48 PM)Marcel Wrote: Hi,
One of the chapter in the book describe a way to solve the Kepler equation (Chapter 29). You will find in this chapter a program to do it.
Marcel
I'll read it!
Salvo
RE: Astronomy... - salvomic - 06-23-2015 03:16 PM
New version, greater precision (especially for the Moon) and rationalization...
Code:
data();
calcN();
calcTS();
julianDay();
julianDayAtZero();
tempo();
precession();
parametri();
sunCalc();
planetCalc();
planets();
sun();
moon();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendent();
zodiaco();
trunc_angle();
simp_angle();
atan2();
thisday:= 2000.0101;
lstd:=0; gmt:= 12;
long:= 0; lat:= 0;
m:= 1; dst:= 0; utc:=0;
deltaT:= 68; // 2015
datetimelist:= MAKELIST(X,X,1,6);
sunlist:= MAKELIST(X,X,1,25);
planet_const:= MAKELIST(X,X,1,10);
ε:= →HMS(23.4392911); // intial value, approx 23.45 (23°26'21")
EXPORT Effemeridi()
// from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015
// Credits: Serge Bouiges, Calcule astronomique pour amateurs
// Jean Meeus, Astronomical Algorithms
// Thanks Didier Lachieze, Marcel Pelletier, Eried
BEGIN
LOCAL ch;
HAngle:=1;
data(); // input data and calc Sun
CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "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 tempo(); END;
IF ch==5 THEN parametri(); 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;
LOCAL t, jd;
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}},
{deltaT, [0], {20,25,5}} },
"Data: Use Shift+9 for H°M′S″",
{"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",
"Long (-E):","Lat (+N)", "delta T"},
{"Enter month", "Enter day", "Enter year (1901-2099)",
"Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",
"Longitude", "Latitude", "Delta T (sec.)"},
{1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});
// 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);
datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };
// Greenwich Mean Time
gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)
thisday:= Y+m/100+D/10000;
jd:= julianDay();
t:=(jd(1)-2451545)/36525;
U:= t/100;
// ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
ε:= 23°26′21.448″ - 1°18′00.93″*U-1.55*U^2+1999.25*U^3-51.38*U^4-249.67*U^5;
ε:= ε-39.05*U^6+7.12*U^7+27.87*U^8+5.79*U^9+2.45*U^10; // Laskar, Meeus
ε:= →HMS(ε);
sunCalc(); // calculate Sun data
END;
planets()
BEGIN
LOCAL ch, nameplanet;
INPUT({{ch, nameplanet:={"Mercurius","Venus","Mars", "Jupiter",
"Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},
"Choose planet",
"planet: ", "Choose the planet to calculate an press OK", 0,5 );
planetCalc(nameplanet(ch));
END;
calcN()
BEGIN
N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT
RETURN N;
END;
calcTS()
BEGIN
LOCAL T, TS, TSd, N0, TSL, TSapp;
LOCAL θ0, θ00, jd, jd0;
calcN();
jd:= julianDay();
jd0:= julianDayAtZero();
T:= (jd(1)-2451545.0)/36525;
θ0:= 280.46061837+360.98564736629*(jd(1)-2451545.0)+0.000387933*T^2-T^3/38710000; // DEG at 0h
θ0:= θ0 MOD 360;
θ00:= 280.46061837+360.98564736629*(jd0(1)-2451545.0)+0.000387933*T^2-T^3/38710000; // DEG at 0h
θ00:= θ00 MOD 360;
N0:= DDAYS(1901.0101, thisday)+1+0; // N at 0 GMT (for TS)
TS:= ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds
TS:= →HMS(TS); // in HMS
TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST
TSL:= →HMS(TSL); // in HMS
TSapp:= (TS + gmt) MOD 24; // apparent ST (ST at Greenwitch at local our)
TSapp:= →HMS(TSapp); // in HMS
TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees
RETURN {TS, TSL, TSapp, TSd, θ0, θ00, T};
END;
julianDay()
BEGIN
LOCAL yy, mm, dd,dde, hh, mn, ss;
LOCAL aa,bb, jd, jde;
yy:= Y; mm:= m;
hh:= datetimelist(4);
mn:= datetimelist(5);
ss:= datetimelist(6);
dd:= D + hh/24 + mn/(60*24) + ss/(60*60*24);
dde:= D + hh/24 + mn/(60*24) + (ss+deltaT)/(60*60*24);
IF (m=1 OR m=2) THEN yy:=Y-1; mm:=m+12;END;
IF yy*100+mm+dd>=158225 THEN
aa:=IP(ABS(yy/100));
bb:=2-aa+IP(aa/4);
ELSE bb:=0; END;
jd:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dd+bb-1524.5;
jde:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dde+bb-1524.5;
RETURN {jd, jde};
END;
julianDayAtZero()
BEGIN
LOCAL yy, mm, dd,dde, hh, mn, ss;
LOCAL aa,bb, jd, jde;
yy:= Y; mm:= m;
hh:= 0;
mn:= 0;
ss:= 0;
dd:= D;
dde:= D + (ss+deltaT)/(60*60*24);
IF (m=1 OR m=2) THEN yy:=Y-1; mm:=m+12;END;
IF yy*100+mm+dd>=158225 THEN
aa:=IP(ABS(yy/100));
bb:=2-aa+IP(aa/4);
ELSE bb:=0; END;
jd:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dd+bb-1524.5;
jde:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dde+bb-1524.5;
RETURN {jd, jde};
END;
tempo()
BEGIN
LOCAL n, jd, TSid;
n:= calcN();
jd:= julianDay();
TSid:= calcTS();
PRINT;
PRINT("Date " + m + " " + D + " " + Y);
PRINT("Julian Day " + jd(1));
PRINT("Julian Day Effem. " + jd(2));
PRINT("N (from 1901jan1) " + n);
PRINT("T (from JD) " + TSid(7));
PRINT("Tempo Siderale 0h " + trunc_angle(TSid(1)));
PRINT("Tempo Siderale locale " + trunc_angle(TSid(2)));
PRINT("Tempo Siderale apparente " + trunc_angle(TSid(3)));
PRINT("Tempo Siderale (DEG) " + trunc_angle(TSid(4)));
PRINT("θ0 Mean ST 0h " + trunc_angle(TSid(5)/15));
PRINT("θ0 Mean ST 0h (DEG) " + trunc_angle(TSid(5)));
PRINT("θ00 MST 00h GMT " + trunc_angle(TSid(6)/15));
PRINT("θ00 MST 00h GMT (DEG) " + trunc_angle(TSid(6)));
PRINT("Time " + lstd + " GMT " + gmt);
RETURN({ROUND(jd(1),5), ROUND(jd(2),5), ROUND(n, 5), trunc_angle(TSid(1)), trunc_angle(TSid(2)),
trunc_angle(TSid(3)), trunc_angle(TSid(4)), trunc_angle(TSid(5)) });
END;
deltaEpsilonPsi()
BEGIN
// Nutazione Δψ (longitudine), Δε (obliquità)
LOCAL L1, deltaPsi, deltaEpsilon, Ω, jd;
jd:= julianDayAtZero();
T:= (jd(2)-2451545)/36525; // use JDE
Ω:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)
L:= (→HMS(280.4665)+36000.76983*T+→HMS(0.0003032*T^2)) MOD 360; // Long media Sole
L1:= →HMS(218.3164591) + →HMS(481267.88134236)*T-0.0013268*T^2; // Long media Luna
L1:= L1+(T^3)/538841-(T^4)/65194000;
L1:= L1 MOD 360;
deltaPsi:= -0°0′17.20″*sin(Ω)-0°0′1.32″*sin(2*L)-0°0′0.23″*sin(2*L1)+0°0′0.21″*sin(2*Ω);
deltaEpsilon:= 0°0′9.2″*cos(Ω)+0°0′0.57″*cos(2*L)+0°0′0.10″*cos(2*L1)-0°0′0.09″*cos(2*Ω);
RETURN({deltaPsi, deltaEpsilon});
END;
precession(alfa, delta)
BEGIN
// precessione dati asce retta e declinazione
LOCAL m,n,n1, deltaAlfa, deltaDelta, Tsec;
Tsec:= IP((Y-2000)/100)+1; // T century from 2000.0
m:= 0°0′3.07496″+0°0′0.00186″*Tsec;
n:= 0°0′1.33621″+0°0′0.00057″*Tsec; // sec
n1:= 0°0′20.0431″+0°0′0.0085″*Tsec; // ° ' "
deltaAlfa:= m+n*sin(alfa)*tan(delta); // Δα
deltaDelta:= n1*cos(alfa); // Δδ
RETURN({deltaAlfa, deltaDelta});
END;
parametri()
BEGIN
LOCAL dep, prec, α, δ;
dep:= deltaEpsilonPsi();
α:= sunlist(5); δ:= sunlist(6);
prec:= precession(α, δ);
// Nutazione Δψ (longitudine), Δε (obliquità), ε (inclinazione eclittica)...
PRINT;
PRINT("ε Inclinazione eclittica "+ ε);
PRINT("Δψ Nutazione longitudine "+ dep(1));
PRINT("Δε Nutazione longitudine "+ dep(2));
PRINT(" ");
PRINT("Δα Precessione AR "+ prec(1));
PRINT("Δδ Precessione dec "+ prec(2));
RETURN({ε, dep(1), dep(2), prec(1), prec(2)}); // ε, Δψ, Δε, Δα, Δδ
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;
simp_angle(ang)
BEGIN
LOCAL angle;
angle:=360*FP(ang/360);
IF angle<0
THEN
angle:=angle+360
END;
RETURN angle;
END;
kepler(ec, M)
// needs eccentricity and true anomaly M
BEGIN
LOCAL M1, j, u;
HAngle:=0; // RAD
M1:= HMS→(M)*PI/180; // anomalia RAD
F:= SIGN(M1); M1:= ABS(M1)/(2*PI);
M1:=(M1-IP(M1))*2*PI*F;
IF M1<0 THEN M1:= M1+2*PI; END;
F:= 1;
IF M1>PI THEN F:= -1; END;
IF M1>PI THEN M1:=2*PI-M; END;
u:= PI/2; D:=PI/4;
FOR j FROM 1 TO 33 DO // 33 iterations (Meeus, Roger Sinnot)
M2:= u-ec*sin(u); u:= u+D*SIGN(M1-M2); D:= D/2;
END; // for
u:= u * F; // anomalia eccentrica
HAngle:=1;
u:= u*180/PI;
RETURN u;
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, h0, H0, temp;
LOCAL ec, a, u, M1, M2, j, t;
LOCAL m0, m1, m2, jd, c, lapp, Ω;
calcN();
jd:= julianDay();
a:= 1.00000023; // semiasse maggiore
t:= N/36525;
T:= (jd(2)-2451545)/36525; // use JDE
// ec:= 0.01675062-0.0000418*t-0.000000137*t^2; // eccentricità (RAD)
ec:= 0.016708617-0.000042037*T-0.0000001236*T^2; // eccentricità Meeus
// ec:= ec*180/PI; // DEG
h0:= -0°50′; // 'standard' altitude of Sun
// L:= →HMS(278.9654+0.985647342*N + 0.0003*t^2) MOD 360; // long media
L:= (→HMS(280.4665)+36000.76983*T+→HMS(0.0003032*T^2)) MOD 360; // Long media Sole Meeus
ω:= →HMS(281.238 + 0.000047067*N + 0.00045*t^2) MOD 360; // long del perielio
// M:= (L - ω); // anomalia vera
//IF M<0 THEN M:= (360+M); END; // potrebbe essere negativa
Ω:= (→HMS(125.04)-→HMS(1934.136)*T) MOD 360; // Nutazione e aberrazione (nodo)
M:= →HMS(357.52910)+→HMS(35999.05030)*T+→HMS(0.0001559)*T^2-→HMS(0.00000048)*T^3; // anom vera Meeus
M:= M MOD 360;
u:=atan2(sin(M),(cos(M)-ec)); // only for small eccentricity
c:= (→HMS(1.914600)-→HMS(0.004817)*T-→HMS(0.000016)*T^2)*sin(M); // equation of center
c:= c+(→HMS(0.019993)-→HMS(0.000101)*T)*sin(2*M)+(→HMS(0.000290))*sin(3*M);
// v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio (Bouiges)
// v:= acos((cos(u)-ec)/(1-ec*cos(u))); // kepler
v:= M + c; // true anomaly, Meeus
// l:= (v + ω) MOD 360; // Longitudine Sole
l:= (L+c) MOD 360; // Ө Longitudine Sole, Meeus
lt:= l + 180; // Longitudine eliocentrica della Terra
lapp:= l - →HMS(0.00569)-→HMS(0.00478)*sin(Ω); // Longitudine apparente (ref true equinox)
lapp:= simp_angle(lapp);
// r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)
r:= (1.000001018*(1-ec^2))/(1+ec*cos(v)); // Distanza Sole-Terra (UA) Meeus
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))) );
// α:= (acos(cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta
// IF l>180 THEN α:= α +180; END;
α:= →HMS(atan2(cos(ε)*sin(l),cos(l))); // Ascensione retta, Meeus (b never > 1.2", l=Ө)
α:= →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
az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd)); // Azimuth
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, transito
sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)
// temp:= (sin(h0)-sin(lat)*sin(δ))/(cos(lat)*cos(δ));
//IF temp<-1 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END;
//H0:= acos(temp);
// m0:= ((α*15 + long - ts(6)) / 360);
//IF m0<0 THEN m0:= m0 +1; END; IF m0>1 THEN m0:=m0-1; END;
// culm:= m0*24;
//temp:= ts(3)*15+360.985647*m0;
// H1:= temp - long - α*15;
// IF temp>360 THEN temp:= temp MOD 360; END;
//culm:= (m0 - H/360) * 24 ;
// sorg:= (culm - H0/360) *24;
// tram:= (culm + H0/360) *24;
sunlist:= {L,ω,l, M, α,δ, az, h, x,y,
sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA,
SAD, E, lt, v, c, Ω, lapp, ec};
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, jd, v, c;
LOCAL Ω, lapp, ec;
L:= sunlist(1); ω:= sunlist(2);
l:= sunlist(3); M:= sunlist(4);
Ω:= sunlist(26); lapp:=sunlist(27);
α:= 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);
v:= sunlist(24); c:= sunlist(25);
ts:= calcTS(); //
jd:= julianDay();
DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);
ec:= sunlist(28);
segno:= zodiaco(l);
PRINT;
PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("at Long " + long + " Lat " + lat);
PRINT("");
PRINT("L longitudine media " + trunc_angle(L));
PRINT("ω long del perielio " + trunc_angle(ω));
PRINT("M anomalia " + trunc_angle(M));
PRINT("v anomalia vera " + v);
PRINT("c equazione del centro " + c);
PRINT(" ");
PRINT("Ө Longitudine Sole " + trunc_angle(l) + " " + segno);
PRINT("λ longitudine apparente " + trunc_angle(lapp));
PRINT("Ω aberr, nutaz nodo "+ Ω);
PRINT("r distanza (UA) " + r);
PRINT(" ");
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("TS apparente " + trunc_angle(ts(3)));
PRINT("θ0 Mean ST 0h " + trunc_angle(ts(5)/15));
PRINT("Equazione del tempo (min) " + trunc_angle(E));
PRINT("Julian day " + jd(1));
PRINT("Dif asc." + trunc_angle(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(" ");
PRINT("Long elioc. Terra "+ trunc_angle(lt));
PRINT("ε epsilon " + ε);
PRINT("ec eccentricità "+ ec);
// RETURN {L, ω, M, v, l, r, x, y};
RETURN [[L, ω], [l, lapp], [M,v], [Ω, c], [α,δ], [az, h], [x,y],
[sorg, tram], [culm, r], [H,ts(1)], [lt, 0], [ε, ec]];
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 d1, h1, m1, segno, jd, rr;
LOCAL lapp, dep, ill, ee, a1,a2,a3;
i:= 5.13333333; // inclinazione 5°8' sull'eclittica
calcN();
jd:= julianDay();
T:= (jd(2)-2451545)/36525; // use JDE
dep:= deltaEpsilonPsi();
// Sun
Ls:= sunlist(1); // Longitudine media
Ms:= sunlist(4); // Anomalia
// end Sun
//L:= (33.231 + 13.17639653*N) MOD 360; // long. media
L:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+(T^3)/538841-(T^4)/65194000);
// Long media Luna Meeus
// Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente
Ω:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)
Ω:= simp_angle(Ω);
// M:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera
M:= 134.9634114+477198.8676313*T+0.0089970*T^2+(T^3)/69699-(T^4)/14712000; // anomalia media
M:= simp_angle(M);
ω:= (L - M) MOD 360; // perigeo medio
// DD:= L - Ls; F:= L-Ω; // for the variations
DD:= 297.8502042+445267.1115168*T-0.0016300*T^2+(T^3)/545868-(T^4)/113065000; // mean elongation
DD:= simp_angle(DD);
F:= 93.2720993+483202.0175273*T-0.0034029*T^2-(T^3)/3526000+(T^4)/863310000; // dist from asc. node
F:= simp_angle(F);
ee:= 1-0.002516*T-0.0000074*T^2; // multiply those w/ Ms *ee
a1:= →HMS(119.75)+→HMS(131.849)*T; // action of Venus
a2:= →HMS(53.09)+→HMS(479264.290)*T; // action of Jupiter
a3:= →HMS(313.45)+→HMS(481266.484)*T;
LongM:= L + 6.288774 * sin(M); // calc Moon longitude
LongM:= LongM + 1.274027 * sin(2*DD-M); // evezione
LongM:= LongM + 0.658314 * sin(2*DD); // variazione
LongM:= LongM + 0.213618 * sin(2*M); // equazione del centro
LongM:= LongM - 0.185116 * sin(Ms) *ee; // equazione annuale
LongM:= LongM - 0.114332 * sin(2*F); // riduzione all'eclittica
LongM:= LongM + 0.058793 * sin(2*DD-2*M);
LongM:= LongM + 0.057066 * sin(2*DD-M-Ms) *ee;
LongM:= LongM + 0.053322 * sin(2*DD+M);
LongM:= LongM + 0.045758 * sin(2*DD-Ms) *ee;
LongM:= LongM - 0.040923 * sin(Ms-M) *ee;
LongM:= LongM - 0.034720 * sin(DD);
LongM:= LongM - 0.030383 * sin(M+Ms) *ee;
LongM:= LongM + 0.015327 * sin(2*DD-2*F);
LongM:= LongM - 0.012528 * sin(M+2*F);
LongM:= LongM + 0.010980 * sin(M-2*F);
LongM:= LongM + 0.010980 * sin(M-2*F);
LongM:= LongM + 0.010675 * sin(4*DD-M);
LongM:= LongM + 0.010034 * sin(3*M);
LongM:= LongM + 0.008548 * sin(4*DD-2*M);
LongM:= LongM - 0.007888 * sin(2*DD+Ms-M) *ee;
LongM:= LongM - 0.006766 * sin(2*DD+Ms) *ee;
LongM:= LongM - 0.005163 * sin(DD-M);
LongM:= LongM + 0.004987 * sin(DD+Ms) *ee;
LongM:= LongM + 0.004036 * sin(2*DD-Ms+M);
LongM:= LongM + 0.003994 * sin(2*DD+2*M);
LongM:= LongM + 0.003861 * sin(4*DD);
LongM:= LongM + 0.003665 * sin(2*DD-3*M);
LongM:= LongM - 0.002689 * sin(Ms-2*M) *ee;
LongM:= LongM - 0.002602 * sin(2*DD-M+2*F);
LongM:= LongM + 0.002390 * sin(2*DD-Ms-2*M) *ee;
LongM:= LongM - 0.002348 * sin(DD+M);
LongM:= LongM + 0.002236 * sin(2*DD-2*Ms) *ee^2;
LongM:= LongM - 0.002120 * sin(Ms+2*M) *ee;
LongM:= LongM - 0.002069 * sin(2*Ms) *ee^2;
LongM:= LongM + 0.002048 * sin(2*DD-2*Ms-M) *ee^2;
LongM:= LongM - 0.001773 * sin(2*DD+M-2*F);
LongM:= LongM - 0.001595 * sin(2*DD+2*F);
LongM:= LongM + 0.001215 * sin(4*DD-Ms-M) *ee;
LongM:= LongM - 0.001110 * sin(2*M+2*F);
LongM:= LongM - 0.000892 * sin(3*DD-M);
LongM:= LongM - 0.000810 * sin(2*DD+Ms+M) *ee;
LongM:= LongM + 0.000759 * sin(4*DD-Ms-2*M) *ee;
LongM:= LongM - 0.000713 * sin(2*Ms-M) * ee^2;
LongM:= LongM - 0.000700 * sin(2*DD+2*Ms-M) *ee^2;
LongM:= LongM + 0.000691 * sin(2*DD+Ms-2*M) *ee;
LongM:= LongM + 0.000596 * sin(2*DD-Ms-2*F) *ee;
LongM:= LongM + 0.000549 * sin(4*DD+M);
LongM:= LongM + 0.000537 * sin(4*M);
LongM:= LongM + 0.000520 * sin(4*DD-Ms) *ee;
LongM:= LongM - 0.000487 * sin(DD-2*M);
LongM:= LongM - 0.000399 * sin(2*DD+Ms-2*F) *ee;
LongM:= LongM - 0.000381 * sin(2*M-2*F);
LongM:= LongM + 0.000351 * sin(DD+Ms+M) *ee;
LongM:= LongM - 0.000340 * sin(3*DD-2*M);
LongM:= LongM + 0.000330 * sin(4*DD-3*M);
LongM:= LongM + 0.000327 * sin(2*DD-Ms+2*M) *ee;
LongM:= LongM - 0.000323 * sin(2*Ms+M) *ee^2;
LongM:= LongM + 0.000299 * sin(DD+Ms-M) *ee;
LongM:= LongM + 0.000294 * sin(2*DD+3*M);
LongM:= LongM + (3958*sin(a1) + 1962*sin(L-F)+318*sin(a2))/1000000;
LongM:= simp_angle(LongM);
lapp:= LongM + dep(1); // longitudine apparente (+ nutazione long, Δψ)
LatM := 5.128122 * sin(F); // termine principale
LatM := LatM + 0.280026 * sin(M+F);
LatM := LatM + 0.277693 * sin(M-F); // equazione del centro
LatM := LatM + 0.173237 * sin(2*DD-F); // grande ineguaglianza
LatM := LatM + 0.055413 * sin(2*DD-M+F); // evezione
LatM := LatM + 0.046271 * sin(2*DD-M-F); // evezione (2)
LatM := LatM + 0.032573 * sin(2*DD+F);
LatM := LatM + 0.017198 * sin(2*M+F);
LatM := LatM + 0.009266 * sin(2*D+M-F);
LatM := LatM + 0.008822 * sin(2*M-F);
LatM := LatM + 0.008216 * sin(2*DD-Ms-F) *ee;
LatM := LatM + 0.004324 * sin(2*DD-2*M-F);
LatM := LatM + 0.004200 * sin(2*DD+M-F);
LatM := LatM - 0.003359 * sin(2*DD+Ms-F) *ee;
LatM := LatM + 0.002463 * sin(2*DD-Ms-M+F) *ee;
LatM := LatM + 0.002211 * sin(2*DD-Ms+F) *ee;
LatM := LatM + 0.002065 * sin(2*DD-Ms-M-F) *ee;
LatM := LatM - 0.001870 * sin(Ms-M-F) *ee;
LatM := LatM + 0.001828 * sin(4*DD-M-F);
LatM := LatM - 0.001794 * sin(Ms-F) *ee;
LatM := LatM - 0.001749 * sin(3*F);
LatM := LatM - 0.001565 * sin(Ms-M+F) *ee;
LatM := LatM - 0.001491 * sin(DD+F);
LatM := LatM - 0.001475 * sin(Ms+M+F) *ee;
LatM := LatM - 0.001410 * sin(Ms+M-F) *ee;
LatM := LatM - 0.001344 * sin(Ms-F) *ee;
LatM := LatM - 0.001335 * sin(DD-F);
LatM := LatM + 0.001107 * sin(3*M+F);
LatM := LatM + 0.001021 * sin(4*DD-F);
LatM := LatM + 0.000833 * sin(4*DD-M+F);
LatM := LatM + 0.000777 * sin(M-3*F);
LatM := LatM + 0.000671 * sin(4*DD-2*M+F);
LatM := LatM + 0.000607 * sin(2*DD-3*F);
LatM := LatM + 0.000596 * sin(2*DD-F);
LatM := LatM + 0.000491 * sin(2*DD-M+F);
LatM := LatM - 0.000451 * sin(2*DD-M-F);
LatM := LatM + 0.000439 * sin(2*DD+F);
LatM := LatM + 0.000422 * sin(2*M-F);
LatM := LatM + 0.000421 * sin(2*DD+M-F);
LatM := LatM - 0.000366 * sin(2*M-F);
LatM := LatM - 0.000351 * sin(2*DD+Ms-F) *ee;
LatM := LatM + 0.000331 * sin(4*DD+F);
LatM := LatM + 0.000315 * sin(2*DD-Ms+M+F) *ee;
LatM := LatM + 0.000302 * sin(2*DD-2*Ms-F) *ee^2;
LatM := LatM - 0.000283 * sin(M+3*F);
LatM := LatM - 0.000229 * sin(2*DD+Ms+M-F) *ee;
LatM := LatM + 0.000223 * sin(DD+Ms-F) *ee;
LatM := LatM + 0.000223 * sin(DD+Ms+F) *ee;
LatM := LatM - 0.000220 * sin(Ms-2*M-F) *ee;
LatM := LatM - 0.000220 * sin(2*DD+Ms-M-F) *ee;
LatM := LatM - 0.000185 * sin(DD+M+F);
LatM := LatM + 0.000181 * sin(2*DD-Ms-2*M-F) *ee;
LatM := LatM - 0.000177 * sin(Ms+2*M+F);
LatM := LatM + 0.000176 * sin(4*DD-2*M-F);
LatM := LatM + 0.000166 * sin(4*DD-Ms-M+F) *ee;
LatM := LatM - 0.000164 * sin(DD+M-F);
LatM := LatM + 0.000132 * sin(4*DD+M-F);
LatM := LatM - 0.000119 * sin(DD-M-F);
LatM := LatM + 0.000115 * sin(4*DD-Ms-F) *ee;
LatM := LatM + 0.000107 * sin(2*DD-2*Ms+F) *ee^2;
LatM := LatM + (0-2235*sin(L)+382*sin(a3)+175*sin(a1-F)
+175*sin(a1+F)+127*sin(L-M)-115*sin(L+M)) / 1000000;
// correzioni per la distanza dalla Terra (Δ)
rr:=0;
rr:= rr -20.905355 * cos(M); // coseni
rr:= rr - 3.699111 * cos(2*DD-M); // evezione
rr:= rr - 2.955968 * cos(2*DD); // variazione
rr:= rr - 0.569925 * cos(2*M); // equazione del centro
rr:= rr + 0.048888 * cos(Ms) *ee; // equazione annuale
rr:= rr - 0.003149 * cos(2*F); // riduzione all'eclittica
rr:= rr + 0.246158 * cos(2*DD-2*M);
rr:= rr - 0.152138 * cos(2*DD-M-Ms) *ee;
rr:= rr - 0.170733 * cos(2*DD+M);
rr:= rr - 0.204586 * cos(2*DD-Ms) *ee;
rr:= rr - 0.129620 * cos(Ms-M);
rr:= rr + 0.108743 * cos(DD);
rr:= rr + 0.104755 * cos(M+Ms) *ee;
rr:= rr + 0.010321 * cos(2*DD-2*F);
rr:= rr - 0.000000 * cos(M+2*F); // null
rr:= rr + 0.079661 * cos(M-2*F);
rr:= rr - 0.034782 * cos(4*DD-M);
rr:= rr - 0.023210 * cos(3*M);
rr:= rr - 0.021636 * cos(4*DD-2*M);
rr:= rr + 0.024208 * cos(2*DD+Ms-M) *ee;
rr:= rr + 0.030824 * cos(2*DD+Ms) *ee;
rr:= rr - 0.008379 * cos(DD-M);
rr:= rr - 0.016675 * cos(DD+Ms) *ee;
rr:= rr - 0.012831 * cos(2*DD-Ms+M);
rr:= rr - 0.010445 * cos(2*DD+2*M);
rr:= rr - 0.011650 * cos(4*DD);
rr:= rr + 0.014403 * cos(2*DD-3*M);
rr:= rr - 0.007003 * cos(Ms-2*M) *ee;
rr:= rr - 0.000000 * cos(2*DD-M+2*F); // null
rr:= rr + 0.010056 * cos(2*DD-Ms-2*M) * ee;
rr:= rr + 0.006322 * cos(DD+M);
rr:= rr - 0.009884 * cos(2*DD-2*Ms) * ee;
rr:= rr + 0.005751 * cos(Ms+2*M) *ee;
rr:= rr - 0.000000 * cos(2*Ms) *ee^2; // null
rr:= rr - 0.004950 * cos(2*DD-2*Ms-M) *ee^2;
rr:= rr + 0.004130 * cos(2*DD+M-2*F);
rr:= rr - 0.000000 * cos(2*DD+2*F); // null
rr:= rr - 0.003958 * cos(4*DD-Ms-M) *ee;
rr:= rr - 0.000000 * cos(2*M+2*F); // null
rr:= rr + 0.003258 * cos(3*DD-M);
rr:= rr + 0.002616 * cos(2*DD+Ms+M) *ee;
rr:= rr - 0.001897 * cos(4*DD-Ms-2*M) *ee;
rr:= rr - 0.002117 * cos(2*Ms-M) * ee^2;
rr:= rr + 0.002354 * cos(2*DD+2*Ms-M) *ee^2;
rr:= rr + 0.000000 * cos(2*DD+Ms-2*M) *ee; // null
rr:= rr + 0.000000 * cos(2*DD-Ms-2*F) *ee; // null
rr:= rr - 0.001423 * cos(4*DD+M);
rr:= rr - 0.001117 * cos(4*M);
rr:= rr - 0.001571 * cos(4*DD-Ms) *ee;
rr:= rr - 0.001739 * cos(DD-2*M);
rr:= rr - 0.000000 * cos(2*DD+Ms-2*F) *ee; // null
rr:= rr - 0.004421 * cos(2*M-2*F);
rr:= rr + 0.000000 * cos(DD+Ms+M) *ee; // null
rr:= rr - 0.000000 * cos(3*DD-2*M); // null
rr:= rr + 0.000000 * cos(4*DD-3*M); // null
rr:= rr + 0.000000 * cos(2*DD-Ms+2*M) *ee; // null
rr:= rr + 0.001165 * cos(2*Ms+M) *ee^2;
rr:= rr + 0.000000 * cos(DD+Ms-M) *ee; // null
rr:= rr + 0.000000 * cos(2*DD+3*M); // null
rr:= rr + 0.008752 * cos(2*DD-M-2*F); // no sin
b:= LatM; // Latitudine
// 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)
// α:= (acos(cos(b)*cos(IFTE(LongM>180, LongM-180, LongM))/cos(δ))); // ascensione retta
//IF LongM>180 THEN α:= α +180; END;
α:= atan2(sin(LongM)*cos(ε)-tan(b)*sin(ε), cos(LongM));
α:= →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)
// Raggio quadratico medio della Terra: 6373,044737
d:= 385000.56 + rr*1000; // Δ distanza Terra-Luna in km (=1/sin(P)),rr in 10^-3
//P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse
P:= asin(6378.14/d); // Parallasse
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
d1:= IP(HMS→(ph/15));
h1:= IP(FP(HMS→(ph/15))*60);
IF (h1>24) THEN d1:= d1+1; h1:= h1 MOD 24; END;
m1:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);
eta:= →HMS(d1+h1/60+m1/3600); // Age of the Moon
segno:= zodiaco(LongM);
ill:= (1+cos(ph))/2; // illuminated fraction
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("λ Longitudine apparente "+ trunc_angle(lapp));
PRINT("Δψ Nutazione long "+ dep(1));
PRINT("Δε Nutazione obl. " + dep(2));
PRINT("ε inclinazione ecl. " + ε);
PRINT("ω Perigeo medio " + trunc_angle(ω));
PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
PRINT(" ");
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("TS (apparente) " + trunc_angle(ts(3)));
PRINT("θ0 Mean ST 0h " + trunc_angle(ts(5)/15));
PRINT("Julian day " + jd(1));
PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT("P Parallasse " + trunc_angle(P));
PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");
PRINT("Semidiametro "+ trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,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);
PRINT("Parte illuminata " + ill);
RETURN [[L, M],[LongM, LatM],[lapp,P], [ω,Ω],[α,δ], [az, h], [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 β
IF b> 360 THEN b:= b MOD 360; END;
// 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;RE: Astronomy... - salvomic - 06-25-2015 05:35 PM
Marcel (and Manuel and others),
please, help me to find this problem now: I'm using for Sun and Moon the Meeus formulae mostly... If I calc Sun I get longitude almost equal to your program, if I calculate Moon I get the same results only if I use your app "at 0h time" (like {2015, 6, 25.0}) and not for the "justnow" time...
I'm trying using formulas like in your program, but including parameters (for sin and so on) into program, not in a variable...
Every numbers seems to be ok, but...
At the real local hour my program returns for Moon now some degrees of difference in longitude (with the first version was almost ok) if your app is used with "justnow" parameter; the same or so, if your is calculated at 0h...
So, the program could be in Julian day (but it seems to be ok) or in a global variable (I'm trying to put if off the program, as far as I can)...
Salvo
Code:
data();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
tempo();
precession();
parametri();
transit();
sunCalc();
planetCalc();
planets();
sun();
sunAh();
moon();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendent();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();
thisday:= 2000.0101;
lstd:=0; gmt:= 12;
long:= 0; lat:= 0;
m:= 1; dst:= 0; utc:=0;
deltaT:= 68; // 2015
datetimelist:= MAKELIST(X,X,1,6);
dateshortlist:= MAKELIST(X,X,1,3);
dateshortlistgmt:= MAKELIST(X,X,1,3);
sunlist:= MAKELIST(X,X,1,25);
planet_const:= MAKELIST(X,X,1,10);
// ε:= epsilon(); // approx 23.45 (23°26'21")
EXPORT Effemeridi()
// from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015
// Credits: Serge Bouiges, Calcule astronomique pour amateurs
// Jean Meeus, Astronomical Algorithms
// Thanks Didier Lachieze, Marcel Pelletier, Eried
BEGIN
LOCAL ch, ε;
HAngle:=1;
data(); // input data and calc Sun
ε:=epsilon();
CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendent", "Quit");
CASE
IF ch==1 THEN sunCalc(); sun(); END;
IF ch==2 THEN moon(); END;
IF ch==3 THEN planets(); END;
IF ch==4 THEN tempo(); END;
IF ch==5 THEN parametri(); END;
IF ch==6 THEN ascendent(); END;
IF ch==7 THEN HAngle:=0; RETURN; END;
DEFAULT sunCalc(); sun();
END; // case
END;
data()
BEGIN
LOCAL dd, mm, yy, hh, mn, ss, loct;
LOCAL t, JD, JDE, ε, dep, eps;
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}},
{deltaT, [0], {20,25,5}} },
"Data: Use Shift+9 for H°M′S″",
{"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",
"Long (-E):","Lat (+N)", "delta T"},
{"Enter month", "Enter day", "Enter year (1901-2099)",
"Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",
"Longitude", "Latitude", "Delta T (sec.)"},
{1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});
// Ragusa (-14°43′41″, 36°55′15″) -
// Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)
gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)
thisday:= Y+m/100+D/10000;
datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };
dateshortlist:= {Y, m, D+HMS→(lstd)/24};
dateshortlistgmt:= {Y, m, D+HMS→(gmt)/24};
// Greenwich Mean Time
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
ε:=epsilon(); // (formula IAU J2000.0)
sunCalc();
END;
epsilon()
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
// ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
ε0:= 23°26′21.448″ - 1°18′00.93″*U-1.55*U^2+1999.25*U^3-51.38*U^4-249.67*U^5;
ε0:= ε0-39.05*U^6+7.12*U^7+27.87*U^8+5.79*U^9+2.45*U^10; // Laskar, Meeus
ε0:= →HMS(ε0);
dep:= deltaEpsilonPsi();
eps:= dep(2);
ε:= ε0+eps; // true obliquity ε
RETURN ε;
END;
deltaEpsilonPsi()
BEGIN
// Nutazione Δψ (longitudine), Δε (obliquità)
LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
LOCAL DD;
JD:= julianDayAtZero(Y,m,D);
JDE:= julianEphemerisDay(dateshortlist);
T:= (JDE-2451545)/36525; // use JDE
DD:=297.85036+445267.111480*T-0.0019142*T^2+(T^3)/189474; // elongazione
DD:=FP(DD/360)*360;
Ωmoon:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)
Lsun:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Sole
Lmoon:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);
// Long media Luna
deltaPsi:= (-0°0′17.20″)*sin(Ωmoon)-(0°0′1.32″)*sin(2*Lsun)+(-0°0′0.23″)*sin(2*Lmoon)+(0°0′0.21″)*sin(2*Ωmoon);
deltaEpsilon:= (0°0′9.2″)*cos(Ωmoon)+(0°0′0.57″)*cos(2*Lsun)+(0°0′0.10″)*cos(2*Lmoon)-(0°0′0.09″)*cos(2*Ωmoon);
RETURN({deltaPsi, deltaEpsilon});
END;
precession(alfa, delta)
BEGIN
// precessione dati asce retta e declinazione
LOCAL m,n,n1, deltaAlfa, deltaDelta, Tsec;
Tsec:= IP((Y-2000)/100)+1; // T century from 2000.0
m:= 0°0′3.07496″+0°0′0.00186″*Tsec;
n:= 0°0′1.33621″+0°0′0.00057″*Tsec; // sec
n1:= 0°0′20.0431″+0°0′0.0085″*Tsec; // ° ' "
deltaAlfa:= m+n*sin(alfa)*tan(delta); // Δα
deltaDelta:= n1*cos(alfa); // Δδ
RETURN({deltaAlfa, deltaDelta});
END;
planets()
BEGIN
LOCAL ch, nameplanet;
INPUT({{ch, nameplanet:={"Mercurius","Venus","Mars", "Jupiter",
"Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},
"Choose planet",
"planet: ", "Choose the planet to calculate an press OK", 0,5 );
planetCalc(nameplanet(ch));
END;
transformEclipticalToEquatorial(lambda, beta)
BEGIN
// trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
LOCAL ε,Epsilon,Lambda, Beta, alpha, delta;
Lambda:=HMS→(lambda);
Beta:=HMS→(beta);
Epsilon:= epsilon();
ε:= epsilon();
alpha:=atan2(sin(Lambda)*cos(Epsilon)-tan(Beta)*sin(Epsilon),cos(Lambda));
IF alpha>=360 THEN alpha:=alpha-360 END;
IF alpha<0 THEN alpha:=alpha+360 END;
delta:=asin(sin(Beta)*cos(Epsilon)+cos(Beta)*sin(Epsilon)*sin(Lambda));
alpha:=→HMS(alpha/15);
delta:=→HMS(delta);
RETURN ({alpha,delta});
END;
transformEquatorialToHorizontal(Lha, Phi, Delta)
BEGIN
LOCAL azimuth,altitude,lha,phi,delta;
lha:=HMS→(Lha);
phi:=HMS→(Phi); // latitude
delta:=HMS→(Delta);
azimuth:=atan2(sin(lha),cos(lha)*sin(phi)-tan(delta)*cos(phi));
altitude:=asin(sin(phi)*sin(delta)+cos(phi)*cos(delta)*cos(lha));
azimuth:=→HMS(azimuth);
altitude:=→HMS(altitude);
RETURN({azimuth,altitude});
END;
transformEclipticalToHorizontal(Lambda, Beta, Lha)
BEGIN
LOCAL ε;
// transform λ and β (long, lat) into azimuth and height
LOCAL lambda, beta, epsilon,phi, lha;
LOCAL alpha,delta,azimuth,altitude;
lambda:=HMS→(Lambda); // longitudine astro
beta:=HMS→(Beta); // latitudine astro
epsilon:= epsilon();
ε:= epsilon();
lha:=HMS→(Lha); // angolo orario
phi:=HMS→(lat); // lat geo
alpha:=atan2(SIN(lambda)*COS(epsilon)-TAN(beta)*SIN(epsilon),COS(lambda));
delta:=asin(sin(beta)*cos(epsilon)+cos(beta)*sin(epsilon)*sin(lambda));
azimuth:=atan2(SIN(lha),COS(lha)*SIN(phi)-TAN(delta)*COS(phi));
altitude:=ASIN(SIN(phi)*SIN(delta)+COS(phi)*COS(delta)*COS(lha));
azimuth:=→HMS(azimuth);
altitude:=→HMS(altitude);
RETURN ({azimuth,altitude});
END;
calcN()
BEGIN
N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT
RETURN N;
END;
calcTS()
BEGIN
LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
ε:= epsilon;
calcN();
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
JD0:= julianDayAtZero(Y,m,D);
dep:= deltaEpsilonPsi();
psi:=dep(1);
eps:=dep(2);
T:= (JD-2451545.0)/36525;
θ0:= meanSideralTime(dateshortlist); // MST DEG at 0h
θ0GMT:= meanSideralTime({Y,m,D}); // MST DEG at 0h at Greenwitch
// TS:= ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds
// TS:= →HMS(TS); // in HMS
TS:= θ0GMT;
// TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST
// TSL:= →HMS(TSL); // in HMS
TSapp:= apparentSideralTime(dateshortlist); // apparent ST (ST at Greenwitch at local our)
TSL:= apparentSideralTime(dateshortlist) - 1*dst;
TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees
RETURN {TS, TSL, TSapp, θ0, T};
END;
meanSideralTime(lista)
// lista like dateshortlist
BEGIN
LOCAL T,MST, JD, ε;
ε:= epsilon;
JD:= julianDay(lista);
T:=(julianDay(lista)-2451545.0)/36525;
MST:=280.46061837+360.98564736629*(JD-2451545.0)+0.000387933*T^2-(T^3)/38710000;
MST:=FP(MST/360)*360; // θ0
IF MST<0 THEN MST:=MST+360 END;
MST:=→HMS(MST/15);
RETURN MST;
END;
apparentSideralTime(lista)
// lista like dateshortlist
BEGIN
LOCAL MST,correction, deltaPsi, AST, dep, ε;
ε:= epsilon;
MST:= meanSideralTime(lista);
dep:= deltaEpsilonPsi();
deltaPsi:= dep(1);
correction:=(HMS→(deltaPsi)*3600*COS(HMS→(ε)))/15;
AST:=→HMS(HMS→(MST)+(correction/3600));
RETURN AST;
END;
julianDay(lista)
// lista like dateshortlist
BEGIN
LOCAL y,m,d,a,b, JD;
y:=lista(1);
m:=lista(2);
d:=lista(3);
IF (m=1 OR m=2) THEN
y:=y-1;
m:=m+12
END;
IF y*100+m+d>=158225 THEN
a:=IP(ABS(y/100));
b:=2-a+IP(a/4);
ELSE
b:=0
END;
JD:=IP(365.25*(y+4716))+IP(30.6001*(m+1))+d+b-1524.5;
RETURN JD;
END;
julianEphemerisDay(lista)
BEGIN
LOCAL JDE;
JDE:=julianDay(makeDateListShort(makeDateList(lista)+{0,0,0,0,0,deltaT}));
RETURN JDE;
END;
julianDayAtZero(Y,m,D)
BEGIN
// JD at 0h GMT
LOCAL y,mm,d,a,b, JD;
y:=Y;
mm:=m;
d:=dateshortlist(3);
IF (mm=1 OR mm=2) THEN
y:=y-1;
mm:=mm+12
END;
IF y*100+mm+d>=158225 THEN
a:=IP(ABS(y/100));
b:=2-a+IP(a/4);
ELSE
b:=0
END;
JD:=IP(365.25*(y+4716))+IP(30.6001*(mm+1))+d+b-1524.5;
RETURN JD;
END;
makeDateList(lista)
// lista {Y m D.d} (dateshortlist)
BEGIN
LOCAL y,o,m,d,h,s,t;
y:=lista(1);
o:=lista(2);
d:=IP(lista(3));
t:=FP(lista(3));
h:=IP(t*24);
t:=FP(t*24);
m:=IP(t*60);
t:=FP(t*60);
s:=IP(t*60);
lista:={y,o,d,h,m,s};
RETURN lista;
END;
makeDateListShort(lista)
BEGIN
// lista {Y m D h mn, s}
LOCAL y,m,d, listacorta;
y:=lista(1);
m:=lista(2);
d:=lista(3)+lista(4)/24+lista(5)/1440+lista(6)/86400;
listacorta:={y,m,d};
RETURN listacorta;
END;
tempo()
BEGIN
LOCAL n, JD, JDE, TSid, dep, psi, eps;
LOCAL ε;
ε:= epsilon();
n:= calcN();
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
TSid:= calcTS();
dep:= deltaEpsilonPsi();
psi:=dep(1);
eps:=dep(2);
PRINT;
PRINT("Date " + m + " " + D + " " + Y);
PRINT("Julian Day " + JD);
PRINT("Julian Day Effem. " + JDE);
PRINT("N (from 1901jan1) " + n);
PRINT("T (from JD) " + TSid(5));
PRINT(" ");
PRINT("Tempo Siderale 0h " + →HMS(TSid(1)));
PRINT("Tempo Siderale locale " + →HMS(TSid(2)));
PRINT("Tempo Siderale apparente " + →HMS(TSid(3)));
PRINT("θ0 Mean ST GMT " + →HMS(TSid(4)));
PRINT("Tempo Siderale (DEG) " + trunc_angle(TSid(1)*15));
PRINT("θ0 Mean ST GMT (DEG) " + trunc_angle(TSid(4)*15));
PRINT(" ");
PRINT("Δψ Nutazione (longit.) "+ psi);
PRINT("Δε Nutazione (obl.) "+ eps);
PRINT("ε Obliquità Ecl. "+ ε);
PRINT(" ");
PRINT("Time " + lstd + " GMT " + gmt);
RETURN({ROUND(JD,5), ROUND(JDE,5), trunc_angle(TSid(1)), trunc_angle(TSid(2)),
trunc_angle(TSid(3)), trunc_angle(TSid(4)), trunc_angle(TSid(5)) });
END;
transit(alfa, delta, h, h0)
BEGIN
// α, δ (AR, dec), height, h0 = standard height of body
LOCAL temp, H0, tsideral, Hor, m0, m1, m2;
LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
α:= alfa*15; // DEG
δ:= delta;
θ0:= apparentSideralTime(makeDateListShort({Y,m,D,0,0,0})) *15; // TS in DEG
ts:= calcTS();
// θ0:= ts(2)*15; // TS=θ00 GMT
temp:= (sin(h0)-sin(lat)*sin(δ))/(cos(lat)*cos(δ));
IF temp<-1 OR temp>1 THEN culm:=0; END; // above horizon?
H0:= acos(temp);
IF H0<-180 THEN H0:=H0+180; END; IF H0>180 THEN H0:= H0-180; END; // it should be from -180 to 180
Hor:=15*HMS→(ts(2)-α); // Hour angle, local ST vs apparent ST
Hor:=simp_angle(→HMS(Hor));
// Hor:= simp_angle(tsideral - long - α);
m0:= HMS→(α + long - θ0) / 360;
IF m0<0 THEN m0:= m0 +1; END; IF m0>1 THEN m0:=m0-1; END;
tsideral:= simp_angle(θ0 + 360.985647*m0); // TS transit DEG
deltaM:= -HMS→(H/360);
culm:= (m0 + deltaM)*24; // transit
m1:= m0-H0/360;
IF m1<0 THEN m1:= m1 +1; END; IF m1>1 THEN m1:=m1-1; END;
tsideral:= simp_angle(θ0 + 360.985647*m1); // TS rising
deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));
sorg:= ((m1 + deltaM)*24) MOD 24; // rising
m2:= m0+H0/360;
IF m2<0 THEN m2:= m2 +1; END; IF m2>1 THEN m2:=m2-1; END;
tsideral:= simp_angle(θ0 + 360.985647*m2); // TS setting
deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));
tram:= ((m2 + deltaM) * 24) MOD 24; // setting
RETURN ({culm, sorg, tram});
END;
parametri()
BEGIN
LOCAL dep, prec, α, δ, ε;
dep:= deltaEpsilonPsi();
α:= sunlist(5); δ:= sunlist(6);
prec:= precession(α, δ);
ε:= epsilon();
// Nutazione Δψ (longitudine), Δε (obliquità), ε (inclinazione eclittica)...
PRINT;
PRINT("ε Obliquità eclittica "+ ε);
PRINT("Δψ Nutazione longitudine "+ dep(1));
PRINT("Δε Nutazione longitudine "+ dep(2));
PRINT(" ");
PRINT("Δα Precessione AR "+ prec(1));
PRINT("Δδ Precessione dec "+ prec(2));
RETURN({ε, dep(1), dep(2), prec(1), prec(2)}); // ε, Δψ, Δε, Δα, Δδ
END;
ascendent()
BEGIN
LOCAL ASC, ts, segno, ε;
ε:= epsilon;
ts:= calcTS();
ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));
ASC:= simp_angle(ASC);
segno:= zodiaco(ASC);
RETURN {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), 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;
simp_angle(ang)
BEGIN
LOCAL angle;
angle:=360*FP(ang/360);
IF angle<0
THEN
angle:=angle+360
END;
RETURN angle;
END;
kepler(ec, M)
// needs eccentricity and true anomaly M
BEGIN
LOCAL M1, j, u;
HAngle:=0; // RAD
M1:= HMS→(M)*PI/180; // anomalia RAD
F:= SIGN(M1); M1:= ABS(M1)/(2*PI);
M1:=(M1-IP(M1))*2*PI*F;
IF M1<0 THEN M1:= M1+2*PI; END;
F:= 1;
IF M1>PI THEN F:= -1; END;
IF M1>PI THEN M1:=2*PI-M; END;
u:= PI/2; D:=PI/4;
FOR j FROM 1 TO 33 DO // 33 iterations (Meeus, Roger Sinnot)
M2:= u-ec*sin(u); u:= u+D*SIGN(M1-M2); D:= D/2;
END; // for
u:= u * F; // anomalia eccentrica
HAngle:=1;
u:= u*180/PI;
RETURN u;
END;
interpolation(lista,n)
BEGIN
LOCAL a,b,c,R,dif1,dif2;
dif1:=ΔLIST(lista);
dif2:=ΔLIST(dif1);
a:=dif1(1);
b:=dif1(2);
c:=dif2(1);
R:=lista(2)+(n/2)*(a+b+n*c);
RETURN R;
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, h0, H0, temp;
LOCAL ec, a, u, M1, M2, j, t, tau;
LOCAL JD, JDE, c, lapp, Ω, θ0, transito;
LOCAL dep, deltapsi, longitude;
ε:= epsilon;
N:= calcN();
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
a:= 1.00000023; // semiasse maggiore
dep:= deltaEpsilonPsi();
deltapsi:= dep(1);
t:= N/36525;
T:= (JD-2451545.0)/36525;
tau:= (JDE-2451545.0)/365250; // millenni - JDE
ec:= 0.016708617-0.000042037*T-0.0000001236*T^2; // eccentricità Meeus
// ec:= ec*180/PI; // DEG
// L:= →HMS(278.9654+0.985647342*N + 0.0003*t^2) MOD 360; // long media Bouiges
// L:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Meeus
L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
ω:= →HMS(281.238 + 0.000047067*N + 0.00045*t^2) MOD 360; // long del perielio
// M:= (L - ω); // anomalia vera
//IF M<0 THEN M:= (360+M); END; // potrebbe essere negativa
Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-0.00000048*T^3); // anom vera Meeus
u:=atan2(sin(M),(cos(M)-ec)); // only for small eccentricity
c:= (1.914600-0.004817*T-0.000016*T^2)*sin(M); // equation of center
c:= c+(0.019993-0.000101*T)*sin(2*M)+(0.000290)*sin(3*M);
// v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio (Bouiges)
// v:= acos((cos(u)-ec)/(1-ec*cos(u))); // kepler
v:= M + c; // true anomaly, Meeus
// l:= (v + ω) MOD 360; // Longitudine Sole
l:= simp_angle(L+c); // Ө Longitudine Sole, Meeus
lt:= l + 180; // Longitudine eliocentrica della Terra
lapp:= simp_angle(l - 0.00569-0.00478*sin(Ω)); // Longitudine apparente (ref true equinox)
// r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)
r:= (1.000001018*(1-ec^2))/(1+ec*cos(v)); // Distanza Sole-Terra (UA) Meeus
x:= r * cos(l); y:= r * sin(l);
δ:= →HMS((asin(sin(ε)*sin(l))) );
α:= →HMS(atan2(cos(ε)*sin(l),cos(l))); // Ascensione retta, Meeus (b never > 1.2", l=Ө)
α:= →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
// longitude:=-HMS→(long); // - if long -E
Hd:=15*HMS→(ts(2)-α); // local ST vs apparent ST
IF Hd<0 THEN Hd:=Hd+360 END;
IF Hd>=360 THEN Hd:=Hd-360 END;
Hd:=→HMS(Hd);
H:= Hd/15; // in HMS
temp:= sunAh(dateshortlist);
az:= temp(1);
h:= temp(2);
// E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)
E:= 4*( L - 0.0057183 - α*15 + deltapsi*cos(ε) ); // α in DEG
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, transito (Bouiges)
// sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
// tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)
// transit routine
h0:= -0°50′; // 'standard' altitude of Sun
transito:= transit(α, δ, h, h0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
sunlist:= {L,ω,l, M, α,δ, az, h, x,y,
sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA,
SAD, E, lt, v, c, Ω, lapp, ec};
RETURN sunlist;
END;
sun()
BEGIN
LOCAL ω, l, v, r, x, y, dz, ε;
LOCAL α, δ, dayY, ts, az, z, h;
LOCAL SAD, DA, sorg, tram, culm, Hd, lt;
LOCAL as, ac, zc, segno, JD, JDE, v, c;
LOCAL Ω, lapp, ec;
ε:= epsilon;
L:= sunlist(1); ω:= sunlist(2);
l:= sunlist(3); M:= sunlist(4);
Ω:= sunlist(26); lapp:=sunlist(27);
α:= 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);
v:= sunlist(24); c:= sunlist(25);
ts:= calcTS(); //
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);
ec:= sunlist(28);
segno:= zodiaco(l);
PRINT;
PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("at Long " + long + " Lat " + lat);
PRINT("");
PRINT("L longitudine media " + trunc_angle(L));
PRINT("ω long del perielio " + trunc_angle(ω));
PRINT("M anomalia media " + trunc_angle(M));
PRINT("v anomalia vera " + trunc_angle(v));
PRINT("c equazione del centro " + c);
PRINT(" ");
PRINT("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("λ longitudine apparente " + trunc_angle(lapp));
PRINT("Ω aberr, nutaz nodo "+ Ω);
PRINT("r distanza (UA) " + r);
PRINT(" ");
PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
PRINT("α Asc.r. (DEG) " + trunc_angle(α*15));
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("TS apparente " + ts(3));
PRINT("θ0 Mean ST GMT " + ts(4));
PRINT("Equazione del tempo (min) " + trunc_angle(E));
PRINT("Julian day " + JD);
PRINT("Dif asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT("");
PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
PRINT(" ");
PRINT("Long elioc. Terra "+ trunc_angle(lt));
PRINT("ε obliquità ecl. " + ε);
PRINT("ec eccentricità "+ ec);
// RETURN {L, ω, M, v, l, r, x, y};
RETURN [[L, ω], [l, lapp], [M,v], [Ω, c], [α,δ], [az, h], [x,y],
[sorg, tram], [culm, r], [H,ts(1)], [lt, 0], [ε, ec]];
END;
sunAh(lista)
// lista like dateshortlist
BEGIN
LOCAL lha,alpha,delta,phi, ah, α;
LOCAL temp,longitude, Azimuth, Altitude, ts;
alpha:=15*HMS→(sunlist(5));
α:= sunlist(5); // α in HMS
delta:=sunlist(6);
phi:=lat;
ts:= calcTS();
// longitude:= HMS→(long); // - if long -E
// lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
lha:= 15*HMS→(ts(2)-α); // local ST vs apparent ST;
IF lha<0 THEN lha:=lha+360 END;
IF lha>=360 THEN lha:=lha-360 END;
lha:=→HMS(lha);
temp:= transformEquatorialToHorizontal(lha, phi, delta);
Azimuth:= HMS→(temp(1))+180;
IF Azimuth>360 THEN Azimuth:=Azimuth-360 END;
Azimuth:= →HMS(Azimuth);
Altitude:= temp(2);
ah:= {Azimuth,Altitude};
RETURN ah;
END;
moon()
BEGIN
LOCAL ω,Ω, L1, M1, l, lsum, bsum, rsum;
LOCAL v, r, x, y, dep, ε;
LOCAL b, P, d, s, temp;
LOCAL α,δ, ts, Hd, zc, dz, h, az;
LOCAL i, as, ac, DA, SAD, DD;
LOCAL culm, sorg, tram, h0, transito;
LOCAL xs, ys, ph, fase, eta;
LOCAL d1, h1, m1, segno, JD, JDE, rr;
LOCAL lapp, ill, ee, a1,a2,a3, longitude;
i:= 5.13333333; // inclinazione 5°8' sull'eclittica
ε:= epsilon();
N:= calcN();
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:= (JDE-2451545.0)/36525; // use JDE
dep:= deltaEpsilonPsi();
// Sun
L:= simp_angle(280.4665+36000.76983*T+0.0003032*T^2); // Long media Sole Meeus
M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+T^3/24490000); // mean anomaly Sun
// end Sun
// L1:= (33.231 + 13.17639653*N) MOD 360; // long. media
L1:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);
// Long media Luna Meeus
// Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente
Ω:= simp_angle(125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000); // Nodo Luna (Moon)
// M1:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera
M1:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media
ω:= (L1 - M1) MOD 360; // perigeo medio
// DD:= L1 - L; F:= L1-Ω; // for the variations
DD:= simp_angle(297.8502042+445267.1115168*T-0.0016300*T^2+T^3/545868-T^4/113065000); // mean elongation
F:= simp_angle(93.2720993+483202.0175273*T-0.0034029*T^2-T^3/3526000+T^4/863310000); // dist from asc. node
// ee:= 1-0.002516*T-0.0000074*T^2; // multiply those w/ M *ee
a1:= simp_angle(119.75+131.849*T); // action of Venus
a2:= simp_angle(53.09+479264.290*T); // action of Jupiter
a3:= simp_angle(313.45+481266.484*T);
lsum:=0; // da sommare a L1
lsum:= lsum + 6288774 * sin(M1); // calc Moon longitude
lsum:= lsum + 1274027 * sin(2*DD-M1); // evezione
lsum:= lsum + 658314 * sin(2*DD); // variazione
lsum:= lsum + 213618 * sin(2*M1); // equazione del centro
lsum:= lsum - 185116 * sin(M) *ee; // equazione annuale
lsum:= lsum - 114332 * sin(2*F); // riduzione all'eclittica
lsum:= lsum + 58793 * sin(2*DD-2*M1);
lsum:= lsum + 57066 * sin(2*DD-M1-M) *ee;
lsum:= lsum + 53322 * sin(2*DD+M1);
lsum:= lsum + 45758 * sin(2*DD-M) *ee;
lsum:= lsum - 40923 * sin(M-M1) *ee;
lsum:= lsum - 34720 * sin(DD);
lsum:= lsum - 30383 * sin(M1+M) *ee;
lsum:= lsum + 15327 * sin(2*DD-2*F);
lsum:= lsum - 12528 * sin(M1+2*F);
lsum:= lsum + 10980 * sin(M1-2*F);
lsum:= lsum + 10675 * sin(4*DD-M1);
lsum:= lsum + 10034 * sin(3*M1);
lsum:= lsum + 8548 * sin(4*DD-2*M1);
lsum:= lsum - 7888 * sin(2*DD+M-M1) *ee;
lsum:= lsum - 6766 * sin(2*DD+M) *ee;
lsum:= lsum - 5163 * sin(DD-M1);
lsum:= lsum + 4987 * sin(DD+M) *ee;
lsum:= lsum + 4036 * sin(2*DD-M+M1) *ee;
lsum:= lsum + 3994 * sin(2*DD+2*M1);
lsum:= lsum + 3861 * sin(4*DD);
lsum:= lsum + 3665 * sin(2*DD-3*M1);
lsum:= lsum - 2689 * sin(M-2*M1) *ee;
lsum:= lsum - 2602 * sin(2*DD-M1+2*F);
lsum:= lsum + 2390 * sin(2*DD-M-2*M1) *ee;
lsum:= lsum - 2348 * sin(DD+M1);
lsum:= lsum + 2236 * sin(2*DD-2*M) *ee^2;
lsum:= lsum - 2120 * sin(M+2*M1) *ee;
lsum:= lsum - 2069 * sin(2*M) *ee^2;
lsum:= lsum + 2048 * sin(2*DD-2*M-M1) *ee^2;
lsum:= lsum - 1773 * sin(2*DD+M1-2*F);
lsum:= lsum - 1595 * sin(2*DD+2*F);
lsum:= lsum + 1215 * sin(4*DD-M-M1) *ee;
lsum:= lsum - 1110 * sin(2*M1+2*F);
lsum:= lsum - 892 * sin(3*DD-M1);
lsum:= lsum - 810 * sin(2*DD+M+M1) *ee;
lsum:= lsum + 759 * sin(4*DD-M-2*M1) *ee;
lsum:= lsum - 713 * sin(2*M-M1) * ee^2;
lsum:= lsum - 700 * sin(2*DD+2*M-M1) *ee^2;
lsum:= lsum + 691 * sin(2*DD+M-2*M1) *ee;
lsum:= lsum + 596 * sin(2*DD-M-2*F) *ee;
lsum:= lsum + 549 * sin(4*DD+M1);
lsum:= lsum + 537 * sin(4*M1);
lsum:= lsum + 520 * sin(4*DD-M) *ee;
lsum:= lsum - 487 * sin(DD-2*M1);
lsum:= lsum - 399 * sin(2*DD+M-2*F) *ee;
lsum:= lsum - 381 * sin(2*M1-2*F);
lsum:= lsum + 351 * sin(DD+M+M1) *ee;
lsum:= lsum - 340 * sin(3*DD-2*M1);
lsum:= lsum + 330 * sin(4*DD-3*M1);
lsum:= lsum + 327 * sin(2*DD-M+2*M1) *ee;
lsum:= lsum - 323 * sin(2*M+M1) *ee^2;
lsum:= lsum + 299 * sin(DD+M-M1) *ee;
lsum:= lsum + 294 * sin(2*DD+3*M1);
lsum:= lsum + 3958*sin(a1)+1962*sin(L1-F)+318*sin(a2);
l:= (simp_angle(L1 + lsum/1000000)); // λ divided by 1 milion
lapp:= →HMS(l + dep(1)); // longitudine apparente (+ nutazione long, Δψ)
// fattori per latitudine
bsum:= 5128122 * sin(F); // termine principale latitudine
bsum:= bsum + 280026 * sin(M1+F);
bsum:= bsum + 277693 * sin(M1-F); // equazione del centro
bsum:= bsum + 173237 * sin(2*DD-F); // grande ineguaglianza
bsum:= bsum + 55413 * sin(2*DD-M1+F); // evezione
bsum:= bsum + 46271 * sin(2*DD-M1-F); // evezione (2)
bsum:= bsum + 32573 * sin(2*DD+F);
bsum:= bsum + 17198 * sin(2*M1+F);
bsum:= bsum + 9266 * sin(2*D+M1-F);
bsum:= bsum + 8822 * sin(2*M1-F);
bsum:= bsum + 8216 * sin(2*DD-M-F) *ee;
bsum:= bsum + 4324 * sin(2*DD-2*M1-F);
bsum:= bsum + 4200 * sin(2*DD+M1+F);
bsum:= bsum - 3359 * sin(2*DD+M-F) *ee;
bsum:= bsum + 2463 * sin(2*DD-M-M1+F) *ee;
bsum:= bsum + 2211 * sin(2*DD-M+F) *ee;
bsum:= bsum + 2065 * sin(2*DD-M-M1-F) *ee;
bsum:= bsum - 1870 * sin(M-M1-F) *ee;
bsum:= bsum + 1828 * sin(4*DD-M1-F);
bsum:= bsum - 1794 * sin(M+F) *ee;
bsum:= bsum - 1749 * sin(3*F);
bsum:= bsum - 1565 * sin(M-M1+F) *ee;
bsum:= bsum - 1491 * sin(DD+F);
bsum:= bsum - 1475 * sin(M+M1+F) *ee;
bsum:= bsum - 1410 * sin(M+M1-F) *ee;
bsum:= bsum - 1344 * sin(M-F) *ee;
bsum:= bsum - 1335 * sin(DD-F);
bsum:= bsum + 1107 * sin(3*M1+F);
bsum:= bsum + 1021 * sin(4*DD-F);
bsum:= bsum + 833 * sin(4*DD-M1+F);
bsum:= bsum + 777 * sin(M1-3*F);
bsum:= bsum + 671 * sin(4*DD-2*M1+F);
bsum:= bsum + 607 * sin(2*DD-3*F);
bsum:= bsum + 596 * sin(2*DD+2*M1-F);
bsum:= bsum + 491 * sin(2*DD-M+M1-F) *ee;
bsum:= bsum - 451 * sin(2*DD-2*M1+F);
bsum:= bsum + 439 * sin(3*M1-F);
bsum:= bsum + 422 * sin(2*DD+2*M1+F);
bsum:= bsum + 421 * sin(2*DD-3*M1-F);
bsum:= bsum - 366 * sin(2*DD+M-M1+F) *ee;
bsum:= bsum - 351 * sin(2*DD+M+F) *ee;
bsum:= bsum + 331 * sin(4*DD+F);
bsum:= bsum + 315 * sin(2*DD-M+M1+F) *ee;
bsum:= bsum + 302 * sin(2*DD-2*M-F) *ee^2;
bsum:= bsum - 283 * sin(M1+3*F);
bsum:= bsum - 229 * sin(2*DD+M+M1-F) *ee;
bsum:= bsum + 223 * sin(DD+M-F) *ee;
bsum:= bsum + 223 * sin(DD+M+F) *ee;
bsum:= bsum - 220 * sin(M-2*M1-F) *ee;
bsum:= bsum - 220 * sin(2*DD+M-M1-F) *ee;
bsum:= bsum - 185 * sin(DD+M1+F);
bsum:= bsum + 181 * sin(2*DD-M-2*M1-F) *ee;
bsum:= bsum - 177 * sin(M+2*M1+F) *ee;
bsum:= bsum + 176 * sin(4*DD-2*M1-F);
bsum:= bsum + 166 * sin(4*DD-M-M1-F) *ee;
bsum:= bsum - 164 * sin(DD+M1-F);
bsum:= bsum + 132 * sin(4*DD+M1-F);
bsum:= bsum - 119 * sin(DD-M1-F);
bsum:= bsum + 115 * sin(4*DD-M-F) *ee;
bsum:= bsum + 107 * sin(2*DD-2*M+F) *ee^2;
bsum:= bsum-2235*sin(L1)+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(L1-M1)-115*sin(L1+M1);
b := →HMS(bsum/1000000); // β div by 1 milion
// correzioni per la distanza dalla Terra (Δ)
rsum:=0;
rsum:= rsum -20905355 * cos(M1); // coseni
rsum:= rsum - 3699111 * cos(2*DD-M1); // evezione
rsum:= rsum - 2955968 * cos(2*DD); // variazione
rsum:= rsum - 569925 * cos(2*M1); // equazione del centro
rsum:= rsum + 48888 * cos(M) *ee; // equazione annuale
rsum:= rsum - 3149 * cos(2*F); // riduzione all'eclittica
rsum:= rsum + 246158 * cos(2*DD-2*M1);
rsum:= rsum - 152138 * cos(2*DD-M1-M) *ee;
rsum:= rsum - 170733 * cos(2*DD+M1);
rsum:= rsum - 204586 * cos(2*DD-M) *ee;
rsum:= rsum - 129620 * cos(M-M1) *ee;
rsum:= rsum + 108743 * cos(DD);
rsum:= rsum + 104755 * cos(M1+M) *ee;
rsum:= rsum + 10321 * cos(2*DD-2*F);
rsum:= rsum + 79661 * cos(M1-2*F);
rsum:= rsum - 34782 * cos(4*DD-M1);
rsum:= rsum - 23210 * cos(3*M1);
rsum:= rsum - 21636 * cos(4*DD-2*M1);
rsum:= rsum + 24208 * cos(2*DD+M-M1) *ee;
rsum:= rsum + 30824 * cos(2*DD+M) *ee;
rsum:= rsum - 8379 * cos(DD-M1);
rsum:= rsum - 16675 * cos(DD+M) *ee;
rsum:= rsum - 12831 * cos(2*DD-M+M1) *ee;
rsum:= rsum - 10445 * cos(2*DD+2*M1);
rsum:= rsum - 11650 * cos(4*DD);
rsum:= rsum + 14403 * cos(2*DD-3*M1);
rsum:= rsum - 7003 * cos(M-2*M1) *ee;
rsum:= rsum + 10056 * cos(2*DD-M-2*M1) * ee;
rsum:= rsum + 6322 * cos(DD+M1);
rsum:= rsum - 9884 * cos(2*DD-2*M) * ee^2;
rsum:= rsum + 5751 * cos(M+2*M1) *ee;
rsum:= rsum - 4950 * cos(2*DD-2*M-M1) *ee^2;
rsum:= rsum + 4130 * cos(2*DD+M1-2*F);
rsum:= rsum - 3958 * cos(4*DD-M-M1) *ee;
rsum:= rsum + 3258 * cos(3*DD-M1);
rsum:= rsum + 2616 * cos(2*DD+M+M1) *ee;
rsum:= rsum - 1897 * cos(4*DD-M-2*M1) *ee;
rsum:= rsum - 2117 * cos(2*M-M1) * ee^2;
rsum:= rsum + 2354 * cos(2*DD+2*M-M1) *ee^2;
rsum:= rsum - 1423 * cos(4*DD+M1);
rsum:= rsum - 1117 * cos(4*M1);
rsum:= rsum - 1571 * cos(4*DD-M) *ee;
rsum:= rsum - 1739 * cos(DD-2*M1);
rsum:= rsum - 4421 * cos(2*M1-2*F);
rsum:= rsum + 1165 * cos(2*M+M1) *ee^2;
rsum:= rsum + 8752 * cos(2*DD-M1-2*F);
rsum:= rsum/1000; // div by 1000
temp:= transformEclipticalToEquatorial(l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
ts:= calcTS(); // calcola tempo siderale
//H:= (ts(2) - α) MOD 24; // Angolo orario
// // H:= asin(-tan(lat)*tan(δ));
//Hd:= →HMS(15*H); // Angolo orario in gradi
longitude:=-HMS→(long); // - if long -E
Hd:=15*HMS→(apparentSideralTime(dateshortlist))-longitude-α;
IF Hd<0 THEN Hd:=Hd+360 END;
IF Hd>=360 THEN Hd:=Hd-360 END;
Hd:=→HMS(Hd);
H:= Hd/15;
temp:= transformEquatorialToHorizontal(Hd, lat, δ);
//temp:= azimuthFromLBPhi(l, b, Hd);
az:= temp(1);
h:= temp(2);
// Raggio quadratico medio della Terra: 6373,044737
d:= (385000.56 + rsum); // Δ distanza Terra-Luna in km (=1/sin(P)),rr in 10^-3
//P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse
P:= asin(6378.14/d); // Parallasse
s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)
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 (Bouiges)
// sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
// tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)
// transit routine
h0:= 0.7275*P -0°34′; // 'standard' altitude of Moon (depend on parallax)
// h0:= 0°7′30″; // medium value: 0.125°
transito:= transit(α, δ, h, h0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
ph:= abs(l - L); // Phase Moon-Sun
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
d1:= IP(HMS→(ph/15));
h1:= IP(FP(HMS→(ph/15))*60);
IF (h1>24) THEN d1:= d1+1; h1:= h1 MOD 24; END;
m1:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);
eta:= →HMS(d1+h1/60+m1/3600); // Age of the Moon
segno:= zodiaco(l);
ill:= (1+cos(ph))/2; // illuminated fraction
PRINT;
PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("L Long media " + trunc_angle(L1));
PRINT("M anomalia vera " + trunc_angle(M1));
PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("β Latitudine " + trunc_angle(b));
PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
PRINT("Δψ Nutazione long "+ dep(1));
PRINT("Δε Nutazione obl. " + dep(2));
PRINT("ε inclinazione ecl. " + ε);
PRINT("ω Perigeo medio " + trunc_angle(ω));
PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
PRINT(" ");
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("TS (apparente) " + trunc_angle(ts(3)));
PRINT("θ0 Mean ST GMT " + trunc_angle(ts(5)));
PRINT("Julian day " + JD);
PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT("P Parallasse " + trunc_angle(P));
PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");
PRINT("Semidiametro "+ trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,3) + " km");
PRINT(" ");
PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
PRINT(" ");
PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");
PRINT("Età " + eta);
PRINT("Parte illuminata " + ill);
RETURN [[L1, M1],[l, b],[lapp,P], [ω,Ω],[α,δ], [az, h],
[ROUND(HMS→(d),6),s], [sorg, tram], [culm, ts(3)],[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 β
IF b> 360 THEN b:= b MOD 360; END;
// 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) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
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("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
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;RE: Astronomy... - salvomic - 06-27-2015 05:51 PM
new version, better handling of Sun and Moon RA, dec and Az,h...
Now should work almost well also the routine for transit, rising and setting of Sun, Moon, planets.
When the program will be complete, I'll put the actual version in the Software Library of the Forum.
Salvo
Code:
data();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
tempo();
precession();
parametri();
transit();
sunCalc();
planetCalc();
planets();
sun();
sunAh();
moon();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendent();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();
thisday:= 2000.0101;
lstd:=0; gmt:= 12;
long:= 0; lat:= 0;
m:= 1; dst:= 0; utc:=0;
deltaT:= 68; // 2015
datetimelist:= MAKELIST(X,X,1,6);
dateshortlist:= MAKELIST(X,X,1,3);
dateshortlistgmt:= MAKELIST(X,X,1,3);
sunlist:= MAKELIST(X,X,1,25);
planet_const:= MAKELIST(X,X,1,10);
// ε:= epsilon(); // approx 23.45 (23°26'21")
EXPORT Effemeridi()
// from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015
// Credits: Serge Bouiges, Calcule astronomique pour amateurs
// Jean Meeus, Astronomical Algorithms
// Thanks Didier Lachieze, Marcel Pelletier, Eried
BEGIN
LOCAL ch, ε;
HAngle:=1;
data(); // input data and calc Sun
ε:=epsilon();
CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendent", "Quit");
CASE
IF ch==1 THEN sunCalc(); sun(); END;
IF ch==2 THEN moon(); END;
IF ch==3 THEN planets(); END;
IF ch==4 THEN tempo(); END;
IF ch==5 THEN parametri(); END;
IF ch==6 THEN ascendent(); END;
IF ch==7 THEN HAngle:=0; RETURN; END;
DEFAULT sunCalc(); sun();
END; // case
END;
data()
BEGIN
LOCAL dd, mm, yy, hh, mn, ss, loct;
LOCAL t, JD, JDE, ε, dep, eps;
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}},
{deltaT, [0], {20,25,5}} },
"Data: Use Shift+9 for H°M′S″",
{"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",
"Long (-E):","Lat (+N)", "delta T"},
{"Enter month", "Enter day", "Enter year (1901-2099)",
"Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",
"Longitude", "Latitude", "Delta T (sec.)"},
{1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});
// Ragusa (-14°43′41″, 36°55′15″) -
// Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)
gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)
thisday:= Y+m/100+D/10000;
datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };
dateshortlist:= {Y, m, D+HMS→(lstd)/24};
dateshortlistgmt:= {Y, m, D+HMS→(gmt)/24};
// Greenwich Mean Time
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
ε:=epsilon(); // (formula IAU J2000.0)
sunCalc();
END;
epsilon()
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
// ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
ε0:= 23°26′21.448″ - 1°18′00.93″*U-1.55*U^2+1999.25*U^3-51.38*U^4-249.67*U^5;
ε0:= ε0-39.05*U^6+7.12*U^7+27.87*U^8+5.79*U^9+2.45*U^10; // Laskar, Meeus
ε0:= →HMS(ε0);
dep:= deltaEpsilonPsi();
eps:= dep(2);
ε:= ε0+eps; // true obliquity ε
RETURN ε;
END;
deltaEpsilonPsi()
BEGIN
// Nutazione Δψ (longitudine), Δε (obliquità)
LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
LOCAL DD;
JD:= julianDayAtZero(Y,m,D);
JDE:= julianEphemerisDay(dateshortlist);
T:= (JDE-2451545)/36525; // use JDE
DD:=297.85036+445267.111480*T-0.0019142*T^2+(T^3)/189474; // elongazione
DD:=FP(DD/360)*360;
Ωmoon:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)
Lsun:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Sole
Lmoon:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);
// Long media Luna
deltaPsi:= (-0°0′17.20″)*sin(Ωmoon)-(0°0′1.32″)*sin(2*Lsun)+(-0°0′0.23″)*sin(2*Lmoon)+(0°0′0.21″)*sin(2*Ωmoon);
deltaEpsilon:= (0°0′9.2″)*cos(Ωmoon)+(0°0′0.57″)*cos(2*Lsun)+(0°0′0.10″)*cos(2*Lmoon)-(0°0′0.09″)*cos(2*Ωmoon);
RETURN({deltaPsi, deltaEpsilon});
END;
precession(alfa, delta)
BEGIN
// precessione dati asce retta e declinazione
LOCAL m,n,n1, deltaAlfa, deltaDelta, Tsec;
Tsec:= IP((Y-2000)/100)+1; // T century from 2000.0
m:= 0°0′3.07496″+0°0′0.00186″*Tsec;
n:= 0°0′1.33621″+0°0′0.00057″*Tsec; // sec
n1:= 0°0′20.0431″+0°0′0.0085″*Tsec; // ° ' "
deltaAlfa:= m+n*sin(alfa)*tan(delta); // Δα
deltaDelta:= n1*cos(alfa); // Δδ
RETURN({deltaAlfa, deltaDelta});
END;
planets()
BEGIN
LOCAL ch, nameplanet;
INPUT({{ch, nameplanet:={"Mercurius","Venus","Mars", "Jupiter",
"Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},
"Choose planet",
"planet: ", "Choose the planet to calculate an press OK", 0,5 );
planetCalc(nameplanet(ch));
END;
transformEclipticalToEquatorial(lambda, beta)
BEGIN
// trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
LOCAL ε,Epsilon,Lambda, Beta, alpha, delta;
Lambda:=HMS→(lambda);
Beta:=HMS→(beta);
Epsilon:= epsilon();
ε:= epsilon();
alpha:=atan2(sin(Lambda)*cos(Epsilon)-tan(Beta)*sin(Epsilon),cos(Lambda));
IF alpha>=360 THEN alpha:=alpha-360 END;
IF alpha<0 THEN alpha:=alpha+360 END;
delta:=asin(sin(Beta)*cos(Epsilon)+cos(Beta)*sin(Epsilon)*sin(Lambda));
alpha:=→HMS(alpha/15);
delta:=→HMS(delta);
RETURN ({alpha,delta});
END;
transformEquatorialToHorizontal(Lha, Phi, Delta)
BEGIN
LOCAL azimuth,altitude,lha,phi,delta;
lha:=HMS→(Lha);
phi:=HMS→(Phi); // latitude
delta:=HMS→(Delta);
azimuth:=atan2(sin(lha),cos(lha)*sin(phi)-tan(delta)*cos(phi));
altitude:=asin(sin(phi)*sin(delta)+cos(phi)*cos(delta)*cos(lha));
azimuth:=→HMS(azimuth);
altitude:=→HMS(altitude);
RETURN({azimuth,altitude});
END;
transformEclipticalToHorizontal(Lambda, Beta, Lha)
BEGIN
LOCAL ε;
// transform λ and β (long, lat) into azimuth and height
LOCAL lambda, beta, epsilon,phi, lha;
LOCAL alpha,delta,azimuth,altitude;
lambda:=HMS→(Lambda); // longitudine astro
beta:=HMS→(Beta); // latitudine astro
epsilon:= epsilon();
ε:= epsilon();
lha:=HMS→(Lha); // angolo orario
phi:=HMS→(lat); // lat geo
alpha:=atan2(SIN(lambda)*COS(epsilon)-TAN(beta)*SIN(epsilon),COS(lambda));
delta:=asin(sin(beta)*cos(epsilon)+cos(beta)*sin(epsilon)*sin(lambda));
azimuth:=atan2(SIN(lha),COS(lha)*SIN(phi)-TAN(delta)*COS(phi));
altitude:=ASIN(SIN(phi)*SIN(delta)+COS(phi)*COS(delta)*COS(lha));
azimuth:=→HMS(azimuth);
altitude:=→HMS(altitude);
RETURN ({azimuth,altitude});
END;
calcN()
BEGIN
N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT
RETURN N;
END;
calcTS()
BEGIN
LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
ε:= epsilon;
calcN();
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
JD0:= julianDayAtZero(Y,m,D);
dep:= deltaEpsilonPsi();
psi:=dep(1);
eps:=dep(2);
T:= (JD-2451545.0)/36525;
θ0:= meanSideralTime(dateshortlist); // MST DEG at 0h
θ0GMT:= meanSideralTime({Y,m,D}); // MST DEG at 0h at Greenwitch
// TS:= ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds
// TS:= →HMS(TS); // in HMS
TS:= θ0GMT;
// TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST
// TSL:= →HMS(TSL); // in HMS
TSapp:= apparentSideralTime(dateshortlist); // apparent ST (ST at Greenwitch at local our)
TSL:= apparentSideralTime(dateshortlist) - 1*dst;
TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees
RETURN {TS, TSL, TSapp, θ0, T};
END;
meanSideralTime(lista)
// lista like dateshortlist
BEGIN
LOCAL T,MST, JD, ε;
ε:= epsilon;
JD:= julianDay(lista);
T:=(julianDay(lista)-2451545.0)/36525;
MST:=280.46061837+360.98564736629*(JD-2451545.0)+0.000387933*T^2-(T^3)/38710000;
MST:=FP(MST/360)*360; // θ0
IF MST<0 THEN MST:=MST+360 END;
MST:=→HMS(MST/15);
RETURN MST;
END;
apparentSideralTime(lista)
// lista like dateshortlist
BEGIN
LOCAL MST,correction, deltaPsi, AST, dep, ε;
ε:= epsilon;
MST:= meanSideralTime(lista);
dep:= deltaEpsilonPsi();
deltaPsi:= dep(1);
correction:=(HMS→(deltaPsi)*3600*COS(HMS→(ε)))/15;
AST:=→HMS(HMS→(MST)+(correction/3600));
RETURN AST;
END;
julianDay(lista)
// lista like dateshortlist
BEGIN
LOCAL y,m,d,a,b, JD;
y:=lista(1);
m:=lista(2);
d:=lista(3);
IF (m=1 OR m=2) THEN
y:=y-1;
m:=m+12
END;
IF y*100+m+d>=158225 THEN
a:=IP(ABS(y/100));
b:=2-a+IP(a/4);
ELSE
b:=0
END;
JD:=IP(365.25*(y+4716))+IP(30.6001*(m+1))+d+b-1524.5;
RETURN JD;
END;
julianEphemerisDay(lista)
BEGIN
LOCAL JDE;
JDE:=julianDay(makeDateListShort(makeDateList(lista)+{0,0,0,0,0,deltaT}));
RETURN JDE;
END;
julianDayAtZero(Y,m,D)
BEGIN
// JD at 0h GMT
LOCAL y,mm,d,a,b, JD;
y:=Y;
mm:=m;
d:=dateshortlist(3);
IF (mm=1 OR mm=2) THEN
y:=y-1;
mm:=mm+12
END;
IF y*100+mm+d>=158225 THEN
a:=IP(ABS(y/100));
b:=2-a+IP(a/4);
ELSE
b:=0
END;
JD:=IP(365.25*(y+4716))+IP(30.6001*(mm+1))+d+b-1524.5;
RETURN JD;
END;
makeDateList(lista)
// lista {Y m D.d} (dateshortlist)
BEGIN
LOCAL y,o,m,d,h,s,t;
y:=lista(1);
o:=lista(2);
d:=IP(lista(3));
t:=FP(lista(3));
h:=IP(t*24);
t:=FP(t*24);
m:=IP(t*60);
t:=FP(t*60);
s:=IP(t*60);
lista:={y,o,d,h,m,s};
RETURN lista;
END;
makeDateListShort(lista)
BEGIN
// lista {Y m D h mn, s}
LOCAL y,m,d, listacorta;
y:=lista(1);
m:=lista(2);
d:=lista(3)+lista(4)/24+lista(5)/1440+lista(6)/86400;
listacorta:={y,m,d};
RETURN listacorta;
END;
tempo()
BEGIN
LOCAL n, JD, JDE, TSid, dep, psi, eps;
LOCAL ε;
ε:= epsilon();
n:= calcN();
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
TSid:= calcTS();
dep:= deltaEpsilonPsi();
psi:=dep(1);
eps:=dep(2);
PRINT;
PRINT("Date " + m + " " + D + " " + Y);
PRINT("Julian Day " + JD);
PRINT("Julian Day Effem. " + JDE);
PRINT("N (from 1901jan1) " + n);
PRINT("T (from JD) " + TSid(5));
PRINT(" ");
PRINT("Tempo Siderale 0h " + →HMS(TSid(1)));
PRINT("Tempo Siderale locale " + →HMS(TSid(2)));
PRINT("Tempo Siderale apparente " + →HMS(TSid(3)));
PRINT("θ0 Mean ST GMT " + →HMS(TSid(4)));
PRINT("Tempo Siderale (DEG) " + trunc_angle(TSid(1)*15));
PRINT("θ0 Mean ST GMT (DEG) " + trunc_angle(TSid(4)*15));
PRINT(" ");
PRINT("Δψ Nutazione (longit.) "+ psi);
PRINT("Δε Nutazione (obl.) "+ eps);
PRINT("ε Obliquità Ecl. "+ ε);
PRINT(" ");
PRINT("Time " + lstd + " GMT " + gmt);
RETURN({ROUND(JD,5), ROUND(JDE,5), trunc_angle(TSid(1)), trunc_angle(TSid(2)),
trunc_angle(TSid(3)), trunc_angle(TSid(4)), trunc_angle(TSid(5)) });
END;
transit(alfa, delta, h0)
BEGIN
// α, δ (AR, dec), height, h0 = standard height of body
LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;
LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
α:= alfa; // AR
δ:= delta;
temp:= calcTS();
// θ0:= apparentSideralTime(makeDateListShort({Y,m,D,0,0,0})) *15; // TS in DEG
θ0 := temp(1) * 15;
temp:= (sin(h0)-sin(lat)*sin(δ))/(cos(lat)*cos(δ));
IF temp<-1 OR temp>1 THEN culm:=0; END; // above horizon?
H0:= acos(temp);
//IF H0<-180 THEN H0:=H0+180; END; IF H0>180 THEN H0:= H0-180; END; // it should be from -180 to 180
m0:= HMS→(α*15 + long - θ0) / 360;
IF m0<0 THEN m0:= m0 +1; END; IF m0>1 THEN m0:=m0-1; END;
tsideral:= simp_angle(θ0 + 360.985647*m0); // TS transit DEG
Hor:= HMS→(tsideral - long -α*15); // Hour angle, local ST vs apparent ST
Hor:=simp_angle(→HMS(Hor));
deltaM:= -HMS→(Hor/360);
culm:= (m0 + deltaM)*24; // transit
m1:= m0-H0/360;
IF m1<0 THEN m1:= m1 +1; END; IF m1>1 THEN m1:=m1-1; END;
tsideral:= simp_angle(θ0 + 360.985647*m1); // TS rising
Hor:= HMS→(tsideral - long - α*15);
Hor:=simp_angle(→HMS(Hor));
h:=asin(sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hor));
deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));
sorg:= ((m1 + deltaM)*24); // rising
m2:= m0+H0/360;
IF m2<0 THEN m2:= m2 +1; END; IF m2>1 THEN m2:=m2-1; END;
tsideral:= simp_angle(θ0 + 360.985647*m2); // TS setting
Hor:= HMS→(tsideral - long - α*15);
Hor:=simp_angle(→HMS(Hor));
h:=asin(sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hor));
deltaM:= (HMS→(h)-HMS→(h0))/(360*cos(δ)*cos(lat)*sin(Hor));
tram:= ((m2 + deltaM) * 24); // setting
RETURN ({culm, sorg, tram});
END;
parametri()
BEGIN
LOCAL dep, prec, α, δ, ε;
dep:= deltaEpsilonPsi();
α:= sunlist(5); δ:= sunlist(6);
prec:= precession(α, δ);
ε:= epsilon();
// Nutazione Δψ (longitudine), Δε (obliquità), ε (inclinazione eclittica)...
PRINT;
PRINT("ε Obliquità eclittica "+ ε);
PRINT("Δψ Nutazione longitudine "+ dep(1));
PRINT("Δε Nutazione longitudine "+ dep(2));
PRINT(" ");
PRINT("Δα Precessione AR "+ prec(1));
PRINT("Δδ Precessione dec "+ prec(2));
RETURN({ε, dep(1), dep(2), prec(1), prec(2)}); // ε, Δψ, Δε, Δα, Δδ
END;
ascendent()
BEGIN
LOCAL ASC, ts, segno, ε;
ε:= epsilon;
ts:= calcTS();
ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));
ASC:= simp_angle(ASC);
segno:= zodiaco(ASC);
RETURN {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), 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;
simp_angle(ang)
BEGIN
LOCAL angle;
angle:=360*FP(ang/360);
IF angle<0
THEN
angle:=angle+360
END;
RETURN angle;
END;
kepler(ec, M)
// needs eccentricity and true anomaly M
BEGIN
LOCAL M1, j, u;
HAngle:=0; // RAD
M1:= HMS→(M)*PI/180; // anomalia RAD
F:= SIGN(M1); M1:= ABS(M1)/(2*PI);
M1:=(M1-IP(M1))*2*PI*F;
IF M1<0 THEN M1:= M1+2*PI; END;
F:= 1;
IF M1>PI THEN F:= -1; END;
IF M1>PI THEN M1:=2*PI-M; END;
u:= PI/2; D:=PI/4;
FOR j FROM 1 TO 33 DO // 33 iterations (Meeus, Roger Sinnot)
M2:= u-ec*sin(u); u:= u+D*SIGN(M1-M2); D:= D/2;
END; // for
u:= u * F; // anomalia eccentrica
HAngle:=1;
u:= u*180/PI;
RETURN u;
END;
interpolation(lista,n)
BEGIN
LOCAL a,b,c,R,dif1,dif2;
dif1:=ΔLIST(lista);
dif2:=ΔLIST(dif1);
a:=dif1(1);
b:=dif1(2);
c:=dif2(1);
R:=lista(2)+(n/2)*(a+b+n*c);
RETURN R;
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, h0, H0, temp;
LOCAL ec, a, u, M1, M2, j, t, tau;
LOCAL JD, JDE, c, lapp, Ω, θ0, transito;
LOCAL dep, deltapsi, longitude;
ε:= epsilon;
N:= calcN();
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
a:= 1.00000023; // semiasse maggiore
dep:= deltaEpsilonPsi();
deltapsi:= dep(1);
t:= N/36525;
T:= (JD-2451545.0)/36525;
tau:= (JDE-2451545.0)/365250; // millenni - JDE
ec:= 0.016708617-0.000042037*T-0.0000001236*T^2; // eccentricità Meeus
// ec:= ec*180/PI; // DEG
// L:= →HMS(278.9654+0.985647342*N + 0.0003*t^2) MOD 360; // long media Bouiges
// L:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Meeus
L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
ω:= →HMS(281.238 + 0.000047067*N + 0.00045*t^2) MOD 360; // long del perielio
// M:= (L - ω); // anomalia vera
//IF M<0 THEN M:= (360+M); END; // potrebbe essere negativa
Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-0.00000048*T^3); // anom vera Meeus
u:=atan2(sin(M),(cos(M)-ec)); // only for small eccentricity
c:= (1.914600-0.004817*T-0.000016*T^2)*sin(M); // equation of center
c:= c+(0.019993-0.000101*T)*sin(2*M)+(0.000290)*sin(3*M);
// v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio (Bouiges)
// v:= acos((cos(u)-ec)/(1-ec*cos(u))); // kepler
v:= M + c; // true anomaly, Meeus
// l:= (v + ω) MOD 360; // Longitudine Sole
l:= simp_angle(L+c); // Ө Longitudine Sole, Meeus
lt:= l + 180; // Longitudine eliocentrica della Terra
lapp:= simp_angle(l - 0.00569-0.00478*sin(Ω)); // Longitudine apparente (ref true equinox)
// r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)
r:= (1.000001018*(1-ec^2))/(1+ec*cos(v)); // Distanza Sole-Terra (UA) Meeus
x:= r * cos(l); y:= r * sin(l);
δ:= →HMS((asin(sin(ε)*sin(l))) );
α:= →HMS(atan2(cos(ε)*sin(l),cos(l))); // Ascensione retta, Meeus (b never > 1.2", l=Ө)
α:= →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
Hd:=15*HMS→(ts(2)-α); // local ST vs apparent ST
Hd:=simp_angle(→HMS(Hd));
H:= Hd/15; // in HMS
temp:= sunAh(dateshortlist);
az:= temp(1); // Azimuth
h:= temp(2); // height
// E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)
E:= 4*( L - 0.0057183 - α*15 + deltapsi*cos(ε) ); // α in DEG
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, transito (Bouiges)
// sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
// tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)
// transit routine
h0:= -0°50′; // 'standard' altitude of Sun
transito:= transit(α, δ, h0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
sunlist:= {L,ω,l, M, α,δ, az, h, x,y,
sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA,
SAD, E, lt, v, c, Ω, lapp, ec};
RETURN sunlist;
END;
sun()
BEGIN
LOCAL ω, l, v, r, x, y, dz, ε;
LOCAL α, δ, dayY, ts, az, z, h;
LOCAL SAD, DA, sorg, tram, culm, Hd, lt;
LOCAL as, ac, zc, segno, JD, JDE, v, c;
LOCAL Ω, lapp, ec;
ε:= epsilon;
L:= sunlist(1); ω:= sunlist(2);
l:= sunlist(3); M:= sunlist(4);
Ω:= sunlist(26); lapp:=sunlist(27);
α:= 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);
v:= sunlist(24); c:= sunlist(25);
ts:= calcTS(); //
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);
ec:= sunlist(28);
segno:= zodiaco(l);
PRINT;
PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("at Long " + long + " Lat " + lat);
PRINT("");
PRINT("L longitudine media " + trunc_angle(L));
PRINT("ω long del perielio " + trunc_angle(ω));
PRINT("M anomalia media " + trunc_angle(M));
PRINT("v anomalia vera " + trunc_angle(v));
PRINT("c equazione del centro " + c);
PRINT(" ");
PRINT("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("λ longitudine apparente " + trunc_angle(lapp));
PRINT("Ω aberr, nutaz nodo "+ Ω);
PRINT("r distanza (UA) " + r);
PRINT(" ");
PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
PRINT("α Asc.r. (DEG) " + trunc_angle(α*15));
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("TS apparente " + ts(3));
PRINT("θ0 Mean ST GMT " + ts(4));
PRINT("Equazione del tempo (min) " + trunc_angle(E));
PRINT("Julian day " + JD);
PRINT("Dif asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT("");
PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
PRINT(" ");
PRINT("Long elioc. Terra "+ trunc_angle(lt));
PRINT("ε obliquità ecl. " + ε);
PRINT("ec eccentricità "+ ec);
// RETURN {L, ω, M, v, l, r, x, y};
RETURN [[L, ω], [l, lapp], [M,v], [Ω, c], [α,δ], [az, h], [x,y],
[sorg, tram], [culm, r], [H,ts(1)], [lt, 0], [ε, ec]];
END;
sunAh(lista)
// lista like dateshortlist
BEGIN
LOCAL lha,alpha,delta,phi, ah, α;
LOCAL temp,longitude, Azimuth, Altitude, ts;
alpha:=15*HMS→(sunlist(5));
α:= sunlist(5); // α in HMS
delta:=sunlist(6);
phi:=lat;
ts:= calcTS();
// longitude:= HMS→(long); // - if long -E
// lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
lha:= 15*HMS→(ts(2)-α); // local ST vs apparent ST;
// IF lha<0 THEN lha:=lha+360 END;
// IF lha>=360 THEN lha:=lha-360 END;
lha:=simp_angle(→HMS(lha));
temp:= transformEquatorialToHorizontal(lha, phi, delta);
Azimuth:= HMS→(temp(1))+180;
IF Azimuth>360 THEN Azimuth:=Azimuth-360 END;
Azimuth:= →HMS(Azimuth);
Altitude:= temp(2);
ah:= {Azimuth,Altitude};
RETURN ah;
END;
moon()
BEGIN
LOCAL ω,Ω, L1, M1, l, lsum, bsum, rsum;
LOCAL v, r, x, y, dep, ε;
LOCAL b, P, d, s, temp;
LOCAL α,δ, ts, Hd, zc, dz, h, az;
LOCAL i, as, ac, DA, SAD, DD;
LOCAL culm, sorg, tram, h0, transito;
LOCAL xs, ys, ph, fase, eta;
LOCAL d1, h1, m1, segno, JD, JDE, rr;
LOCAL lapp, ill, ee, a1,a2,a3, longitude;
i:= 5.13333333; // inclinazione 5°8' sull'eclittica
ε:= epsilon();
N:= calcN();
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:= (JDE-2451545.0)/36525; // use JDE
dep:= deltaEpsilonPsi();
// Sun
L:= simp_angle(280.4665+36000.76983*T+0.0003032*T^2); // Long media Sole Meeus
M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+T^3/24490000); // mean anomaly Sun
// end Sun
// L1:= (33.231 + 13.17639653*N) MOD 360; // long. media
L1:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);
// Long media Luna Meeus
// Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente
Ω:= simp_angle(125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000); // Nodo Luna (Moon)
// M1:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera
M1:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media
ω:= (L1 - M1) MOD 360; // perigeo medio
// DD:= L1 - L; F:= L1-Ω; // for the variations
DD:= simp_angle(297.8502042+445267.1115168*T-0.0016300*T^2+T^3/545868-T^4/113065000); // mean elongation
F:= simp_angle(93.2720993+483202.0175273*T-0.0034029*T^2-T^3/3526000+T^4/863310000); // dist from asc. node
// ee:= 1-0.002516*T-0.0000074*T^2; // multiply those w/ M *ee
a1:= simp_angle(119.75+131.849*T); // action of Venus
a2:= simp_angle(53.09+479264.290*T); // action of Jupiter
a3:= simp_angle(313.45+481266.484*T);
lsum:=0; // da sommare a L1
lsum:= lsum + 6288774 * sin(M1); // calc Moon longitude
lsum:= lsum + 1274027 * sin(2*DD-M1); // evezione
lsum:= lsum + 658314 * sin(2*DD); // variazione
lsum:= lsum + 213618 * sin(2*M1); // equazione del centro
lsum:= lsum - 185116 * sin(M) *ee; // equazione annuale
lsum:= lsum - 114332 * sin(2*F); // riduzione all'eclittica
lsum:= lsum + 58793 * sin(2*DD-2*M1);
lsum:= lsum + 57066 * sin(2*DD-M1-M) *ee;
lsum:= lsum + 53322 * sin(2*DD+M1);
lsum:= lsum + 45758 * sin(2*DD-M) *ee;
lsum:= lsum - 40923 * sin(M-M1) *ee;
lsum:= lsum - 34720 * sin(DD);
lsum:= lsum - 30383 * sin(M1+M) *ee;
lsum:= lsum + 15327 * sin(2*DD-2*F);
lsum:= lsum - 12528 * sin(M1+2*F);
lsum:= lsum + 10980 * sin(M1-2*F);
lsum:= lsum + 10675 * sin(4*DD-M1);
lsum:= lsum + 10034 * sin(3*M1);
lsum:= lsum + 8548 * sin(4*DD-2*M1);
lsum:= lsum - 7888 * sin(2*DD+M-M1) *ee;
lsum:= lsum - 6766 * sin(2*DD+M) *ee;
lsum:= lsum - 5163 * sin(DD-M1);
lsum:= lsum + 4987 * sin(DD+M) *ee;
lsum:= lsum + 4036 * sin(2*DD-M+M1) *ee;
lsum:= lsum + 3994 * sin(2*DD+2*M1);
lsum:= lsum + 3861 * sin(4*DD);
lsum:= lsum + 3665 * sin(2*DD-3*M1);
lsum:= lsum - 2689 * sin(M-2*M1) *ee;
lsum:= lsum - 2602 * sin(2*DD-M1+2*F);
lsum:= lsum + 2390 * sin(2*DD-M-2*M1) *ee;
lsum:= lsum - 2348 * sin(DD+M1);
lsum:= lsum + 2236 * sin(2*DD-2*M) *ee^2;
lsum:= lsum - 2120 * sin(M+2*M1) *ee;
lsum:= lsum - 2069 * sin(2*M) *ee^2;
lsum:= lsum + 2048 * sin(2*DD-2*M-M1) *ee^2;
lsum:= lsum - 1773 * sin(2*DD+M1-2*F);
lsum:= lsum - 1595 * sin(2*DD+2*F);
lsum:= lsum + 1215 * sin(4*DD-M-M1) *ee;
lsum:= lsum - 1110 * sin(2*M1+2*F);
lsum:= lsum - 892 * sin(3*DD-M1);
lsum:= lsum - 810 * sin(2*DD+M+M1) *ee;
lsum:= lsum + 759 * sin(4*DD-M-2*M1) *ee;
lsum:= lsum - 713 * sin(2*M-M1) * ee^2;
lsum:= lsum - 700 * sin(2*DD+2*M-M1) *ee^2;
lsum:= lsum + 691 * sin(2*DD+M-2*M1) *ee;
lsum:= lsum + 596 * sin(2*DD-M-2*F) *ee;
lsum:= lsum + 549 * sin(4*DD+M1);
lsum:= lsum + 537 * sin(4*M1);
lsum:= lsum + 520 * sin(4*DD-M) *ee;
lsum:= lsum - 487 * sin(DD-2*M1);
lsum:= lsum - 399 * sin(2*DD+M-2*F) *ee;
lsum:= lsum - 381 * sin(2*M1-2*F);
lsum:= lsum + 351 * sin(DD+M+M1) *ee;
lsum:= lsum - 340 * sin(3*DD-2*M1);
lsum:= lsum + 330 * sin(4*DD-3*M1);
lsum:= lsum + 327 * sin(2*DD-M+2*M1) *ee;
lsum:= lsum - 323 * sin(2*M+M1) *ee^2;
lsum:= lsum + 299 * sin(DD+M-M1) *ee;
lsum:= lsum + 294 * sin(2*DD+3*M1);
lsum:= lsum + 3958*sin(a1)+1962*sin(L1-F)+318*sin(a2);
l:= (simp_angle(L1 + lsum/1000000)); // λ divided by 1 milion
lapp:= →HMS(l + dep(1)); // longitudine apparente (+ nutazione long, Δψ)
// fattori per latitudine
bsum:= 5128122 * sin(F); // termine principale latitudine
bsum:= bsum + 280026 * sin(M1+F);
bsum:= bsum + 277693 * sin(M1-F); // equazione del centro
bsum:= bsum + 173237 * sin(2*DD-F); // grande ineguaglianza
bsum:= bsum + 55413 * sin(2*DD-M1+F); // evezione
bsum:= bsum + 46271 * sin(2*DD-M1-F); // evezione (2)
bsum:= bsum + 32573 * sin(2*DD+F);
bsum:= bsum + 17198 * sin(2*M1+F);
bsum:= bsum + 9266 * sin(2*D+M1-F);
bsum:= bsum + 8822 * sin(2*M1-F);
bsum:= bsum + 8216 * sin(2*DD-M-F) *ee;
bsum:= bsum + 4324 * sin(2*DD-2*M1-F);
bsum:= bsum + 4200 * sin(2*DD+M1+F);
bsum:= bsum - 3359 * sin(2*DD+M-F) *ee;
bsum:= bsum + 2463 * sin(2*DD-M-M1+F) *ee;
bsum:= bsum + 2211 * sin(2*DD-M+F) *ee;
bsum:= bsum + 2065 * sin(2*DD-M-M1-F) *ee;
bsum:= bsum - 1870 * sin(M-M1-F) *ee;
bsum:= bsum + 1828 * sin(4*DD-M1-F);
bsum:= bsum - 1794 * sin(M+F) *ee;
bsum:= bsum - 1749 * sin(3*F);
bsum:= bsum - 1565 * sin(M-M1+F) *ee;
bsum:= bsum - 1491 * sin(DD+F);
bsum:= bsum - 1475 * sin(M+M1+F) *ee;
bsum:= bsum - 1410 * sin(M+M1-F) *ee;
bsum:= bsum - 1344 * sin(M-F) *ee;
bsum:= bsum - 1335 * sin(DD-F);
bsum:= bsum + 1107 * sin(3*M1+F);
bsum:= bsum + 1021 * sin(4*DD-F);
bsum:= bsum + 833 * sin(4*DD-M1+F);
bsum:= bsum + 777 * sin(M1-3*F);
bsum:= bsum + 671 * sin(4*DD-2*M1+F);
bsum:= bsum + 607 * sin(2*DD-3*F);
bsum:= bsum + 596 * sin(2*DD+2*M1-F);
bsum:= bsum + 491 * sin(2*DD-M+M1-F) *ee;
bsum:= bsum - 451 * sin(2*DD-2*M1+F);
bsum:= bsum + 439 * sin(3*M1-F);
bsum:= bsum + 422 * sin(2*DD+2*M1+F);
bsum:= bsum + 421 * sin(2*DD-3*M1-F);
bsum:= bsum - 366 * sin(2*DD+M-M1+F) *ee;
bsum:= bsum - 351 * sin(2*DD+M+F) *ee;
bsum:= bsum + 331 * sin(4*DD+F);
bsum:= bsum + 315 * sin(2*DD-M+M1+F) *ee;
bsum:= bsum + 302 * sin(2*DD-2*M-F) *ee^2;
bsum:= bsum - 283 * sin(M1+3*F);
bsum:= bsum - 229 * sin(2*DD+M+M1-F) *ee;
bsum:= bsum + 223 * sin(DD+M-F) *ee;
bsum:= bsum + 223 * sin(DD+M+F) *ee;
bsum:= bsum - 220 * sin(M-2*M1-F) *ee;
bsum:= bsum - 220 * sin(2*DD+M-M1-F) *ee;
bsum:= bsum - 185 * sin(DD+M1+F);
bsum:= bsum + 181 * sin(2*DD-M-2*M1-F) *ee;
bsum:= bsum - 177 * sin(M+2*M1+F) *ee;
bsum:= bsum + 176 * sin(4*DD-2*M1-F);
bsum:= bsum + 166 * sin(4*DD-M-M1-F) *ee;
bsum:= bsum - 164 * sin(DD+M1-F);
bsum:= bsum + 132 * sin(4*DD+M1-F);
bsum:= bsum - 119 * sin(DD-M1-F);
bsum:= bsum + 115 * sin(4*DD-M-F) *ee;
bsum:= bsum + 107 * sin(2*DD-2*M+F) *ee^2;
bsum:= bsum-2235*sin(L1)+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(L1-M1)-115*sin(L1+M1);
b := →HMS(bsum/1000000); // β div by 1 milion
// correzioni per la distanza dalla Terra (Δ)
rsum:=0;
rsum:= rsum -20905355 * cos(M1); // coseni
rsum:= rsum - 3699111 * cos(2*DD-M1); // evezione
rsum:= rsum - 2955968 * cos(2*DD); // variazione
rsum:= rsum - 569925 * cos(2*M1); // equazione del centro
rsum:= rsum + 48888 * cos(M) *ee; // equazione annuale
rsum:= rsum - 3149 * cos(2*F); // riduzione all'eclittica
rsum:= rsum + 246158 * cos(2*DD-2*M1);
rsum:= rsum - 152138 * cos(2*DD-M1-M) *ee;
rsum:= rsum - 170733 * cos(2*DD+M1);
rsum:= rsum - 204586 * cos(2*DD-M) *ee;
rsum:= rsum - 129620 * cos(M-M1) *ee;
rsum:= rsum + 108743 * cos(DD);
rsum:= rsum + 104755 * cos(M1+M) *ee;
rsum:= rsum + 10321 * cos(2*DD-2*F);
rsum:= rsum + 79661 * cos(M1-2*F);
rsum:= rsum - 34782 * cos(4*DD-M1);
rsum:= rsum - 23210 * cos(3*M1);
rsum:= rsum - 21636 * cos(4*DD-2*M1);
rsum:= rsum + 24208 * cos(2*DD+M-M1) *ee;
rsum:= rsum + 30824 * cos(2*DD+M) *ee;
rsum:= rsum - 8379 * cos(DD-M1);
rsum:= rsum - 16675 * cos(DD+M) *ee;
rsum:= rsum - 12831 * cos(2*DD-M+M1) *ee;
rsum:= rsum - 10445 * cos(2*DD+2*M1);
rsum:= rsum - 11650 * cos(4*DD);
rsum:= rsum + 14403 * cos(2*DD-3*M1);
rsum:= rsum - 7003 * cos(M-2*M1) *ee;
rsum:= rsum + 10056 * cos(2*DD-M-2*M1) * ee;
rsum:= rsum + 6322 * cos(DD+M1);
rsum:= rsum - 9884 * cos(2*DD-2*M) * ee^2;
rsum:= rsum + 5751 * cos(M+2*M1) *ee;
rsum:= rsum - 4950 * cos(2*DD-2*M-M1) *ee^2;
rsum:= rsum + 4130 * cos(2*DD+M1-2*F);
rsum:= rsum - 3958 * cos(4*DD-M-M1) *ee;
rsum:= rsum + 3258 * cos(3*DD-M1);
rsum:= rsum + 2616 * cos(2*DD+M+M1) *ee;
rsum:= rsum - 1897 * cos(4*DD-M-2*M1) *ee;
rsum:= rsum - 2117 * cos(2*M-M1) * ee^2;
rsum:= rsum + 2354 * cos(2*DD+2*M-M1) *ee^2;
rsum:= rsum - 1423 * cos(4*DD+M1);
rsum:= rsum - 1117 * cos(4*M1);
rsum:= rsum - 1571 * cos(4*DD-M) *ee;
rsum:= rsum - 1739 * cos(DD-2*M1);
rsum:= rsum - 4421 * cos(2*M1-2*F);
rsum:= rsum + 1165 * cos(2*M+M1) *ee^2;
rsum:= rsum + 8752 * cos(2*DD-M1-2*F);
rsum:= rsum/1000; // div by 1000
temp:= transformEclipticalToEquatorial(l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
ts:= calcTS(); // calcola tempo siderale
Hd:=15*HMS→(ts(2)-α); // local ST vs apparent ST
Hd:=simp_angle(→HMS(Hd));
H:= Hd/15; // in HMS
temp:= transformEquatorialToHorizontal(Hd, lat, δ);
az:= HMS→(temp(1))+180; // azimuth
IF az>360 THEN az:=az-360 END;
az:= →HMS(az);
h:= temp(2); // height
// Raggio quadratico medio della Terra: 6373,044737
d:= (385000.56 + rsum); // Δ distanza Terra-Luna in km (=1/sin(P)),rr in 10^-3
//P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse
P:= asin(6378.14/d); // Parallasse
s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)
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 (Bouiges)
// sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
// tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)
// transit routine
h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)
// h0:= 0°7′30″; // medium value: 0.125°
transito:= transit(α, δ, h0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
ph:= abs(l - L); // Phase Moon-Sun
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
d1:= IP(HMS→(ph/15));
h1:= IP(FP(HMS→(ph/15))*60);
IF (h1>24) THEN d1:= d1+1; h1:= h1 MOD 24; END;
m1:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);
eta:= →HMS(d1+h1/60+m1/3600); // Age of the Moon
segno:= zodiaco(l);
ill:= 1-(1+cos(ph))/2; // illuminated fraction
PRINT;
PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("L Long media " + trunc_angle(L1));
PRINT("M anomalia vera " + trunc_angle(M1));
PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("β Latitudine " + trunc_angle(b));
PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
PRINT("Δψ Nutazione long "+ dep(1));
PRINT("Δε Nutazione obl. " + dep(2));
PRINT("ε inclinazione ecl. " + ε);
PRINT("ω Perigeo medio " + trunc_angle(ω));
PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
PRINT(" ");
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("TS apparente " + ts(3));
PRINT("θ0 Mean ST GMT " + ts(4));
PRINT("Julian day " + JD);
PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT("P Parallasse " + trunc_angle(P));
PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");
PRINT("Semidiametro "+ trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,3) + " km");
PRINT(" ");
PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
PRINT(" ");
PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");
PRINT("Età " + eta);
PRINT("Parte illuminata " + ill);
RETURN [[L1, M1],[l, b],[lapp,P], [ω,Ω],[α,δ], [az, h],
[ROUND(HMS→(d),6),s], [sorg, tram], [culm, ts(3)],[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, ε;
LOCAL transito, h0, JD, temp;
JD:= julianDay(dateshortlist);
// 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 β
IF b> 360 THEN b:= b MOD 360; END;
// 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;
IF l>180 THEN α:= α +180; END;
// δ:= →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
//α:= →HMS(α/15) MOD 24; // α in hms
temp:= transformEclipticalToEquatorial(l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
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)
// transit routine
h0:= -0°34′; // 'standard' altitude of a planet
transito:= transit(α, δ, h0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
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) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
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("TS apparente " + ts(3));
PRINT("θ0 Mean ST GMT " + ts(4));
PRINT("Julian day " + JD);
PRINT("Diff. ascensionale " + trunc_angle(DA));
PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));
PRINT(" ");
PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");
PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
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;