Astronomy: Effemeridi (Ephemeris) - salvomic - 06-29-2015 08:24 PM
hi everybody,
I'm developing a program about Astronomy.
I developed it on a TI 74 basic many years ago, now I've revisited (and enhanced) it for the Prime; it's name is Effemeridi (Astrolabio).
It can calculates ephemeris for Sun, Moon and planets, transit, sidereal time, nutation and obliquity, julian day, ascendant and many others data.
I'll put here new version as they will be completed.
I advice to use the program together with the splendid Astro Lab 4 from Marcel
The program present data in Terminal (multi pages: for now use enter to go on, and *never* touch screen or the program go on without presenting the other pages) and export variables with data.
For now some data names are still in Italian, I'll translate those later.
Enjoy!
Salvo Micciché
EDIT: new version, reviewed, more precise transit calculation.
EDIT 2017 November 8: last version (firmware from 10077 to the actual beta 12951) in this post of the thread.
After the public release of the new FW I hope to confirm that version or put a new one.
***
The first code, however, was this one:
Code:
data();
changeData();
whatToDo();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
tempo();
precession();
parametri();
transit();
sunCalc();
sun();
sunAh();
sunAlphaDelta();
moon();
moonCalc();
moonAh();
moonAlphaDelta();
planets();
planetCalc();
planetDisplay();
planetAlphaDelta();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendant();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();
EXPORT Tempi;
EXPORT EpsilonPsi;
EXPORT Ascendant;
EXPORT Sun_Effemerids;
EXPORT Sun_Transit;
EXPORT Sun_Others;
EXPORT Moon_Effemerids;
EXPORT Moon_Transit;
EXPORT Moon_Others;
EXPORT Planet_Effemerids;
EXPORT Planet_Transit;
EXPORT Planet_Others;
thisday:= 1964.0529;
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,40);
moonlist:= MAKELIST(X,X,1,40);
planetlist:= MAKELIST(X,X,1,40);
planet_const:= MAKELIST(X,X,1,40);
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, DrD
BEGIN
HAngle:=1;
data(); // input data and calc Sun
END;
data()
BEGIN
// initial input
LOCAL dd, mm, yy, hh, mn, ss, loct;
LOCAL ε, JD, JDE;
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};
ε:=epsilon(dateshortlist);
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
sunCalc(dateshortlist); // set initial data for Sun
moonCalc(dateshortlist);
whatToDo(dateshortlist);
END;
changeData()
BEGIN
// input to change date, lat, long...
LOCAL ε, JD, JDE;
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},{m,D,Y,lstd, dst, utc, long, lat, deltaT});
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};
ε:=epsilon(dateshortlist);
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
whatToDo(dateshortlist);
END;
whatToDo(lista)
BEGIN
LOCAL ch;
CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendant",
"Quit");
CASE
IF ch==1 THEN sunCalc(lista); sun(lista); END;
IF ch==2 THEN moonCalc(lista); moon(lista); END;
IF ch==3 THEN planets(lista); END;
IF ch==4 THEN tempo(lista); END;
IF ch==5 THEN parametri(lista); END;
IF ch==6 THEN ascendant(lista); END;
IF ch==7 THEN HAngle:=0; RETURN; END;
DEFAULT
END; // case
END;
planets(lista)
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), lista);
planetDisplay(nameplanet(ch),lista);
END;
epsilon(lista)
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
T:=(JD-2451545)/36525;
U:= T/100;
// ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
// 4680.93″ = 1°18′00.93″
ε0:= 23°26′21.448″ - (HMS→(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(lista);
eps:= dep(2);
ε:= ε0+eps; // true obliquity ε
RETURN ε;
END;
deltaEpsilonPsi(lista)
BEGIN
// Nutazione Δψ (longitudine), Δε (obliquità)
LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
LOCAL DD;
// JD:= julianDayAtZero(Y,m,D);
JDE:= julianEphemerisDay(lista);
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;
transformEclipticalToEquatorial(lista, lambda, beta)
BEGIN
// trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
LOCAL ε,Lambda, Beta, alpha, delta;
Lambda:=HMS→(lambda);
Beta:=HMS→(beta);
ε:= epsilon(lista);
alpha:=atan2(sin(Lambda)*cos(ε)-tan(Beta)*sin(ε),cos(Lambda));
IF alpha>=360 THEN alpha:=alpha-360 END;
IF alpha<0 THEN alpha:=alpha+360 END;
delta:=asin(sin(Beta)*cos(ε)+cos(Beta)*sin(ε)*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(lista, Lambda, Beta, Lha)
BEGIN
// transform λ and β (long, lat) into azimuth and height
LOCAL lambda, beta,phi, lha, ε;
LOCAL alpha,delta,azimuth,altitude;
lambda:=HMS→(Lambda); // longitudine astro
beta:=HMS→(Beta); // latitudine astro
ε:= epsilon(lista);
lha:=HMS→(Lha); // angolo orario
phi:=HMS→(lat); // lat geo
alpha:=atan2(SIN(lambda)*COS(ε)-TAN(beta)*SIN(ε),COS(lambda));
delta:=asin(sin(beta)*cos(ε)+cos(beta)*sin(ε)*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(lista)
BEGIN
LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
JD0:= julianDayAtZero(Y,m,D);
dep:= deltaEpsilonPsi(lista);
psi:=dep(1);
eps:=dep(2);
T:= (JD-2451545.0)/36525;
θ0:= meanSideralTime(lista); // 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(lista); // apparent ST (ST at Greenwitch at local our)
TSL:= apparentSideralTime(lista) - 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(lista);
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(lista);
MST:= meanSideralTime(lista);
dep:= deltaEpsilonPsi(lista);
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(lista)
BEGIN
LOCAL JD, JDE, TSid, dep, psi, eps;
LOCAL ε, listatempi, ch;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
TSid:= calcTS(lista);
dep:= deltaEpsilonPsi(lista);
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(" ");
WAIT();
PRINT();
PRINT("Δψ Nutazione (longit.) "+ psi);
PRINT("Δε Nutazione (obl.) "+ eps);
PRINT("ε Obliquità Ecl. "+ ε);
PRINT(" ");
PRINT("Time " + lstd + " GMT " + gmt);
listatempi:= {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)) };
Tempi:= listatempi;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
transit(lista, listaARh0)
BEGIN
// lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
// α, δ (AR, dec), height, h0 = standard height of body
LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;
LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
LOCAL alfa1, alfa3, delta1, delta3, h0, alpha, del, n;
α:= listaARh0(2); // AR for day
δ:= listaARh0(5); // decl for day
alfa1:= listaARh0(1); // AR day-1
alfa3:= listaARh0(3); // AR day+1
delta1:= listaARh0(4); // decl day-1
delta3:= listaARh0(6); // decl day+1
h0:= listaARh0(7); // standard height
temp:= calcTS(lista);
θ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
n:= m0 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*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
n:= m1 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*15);
Hor:=simp_angle(→HMS(Hor));
h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
deltaM:= HMS→((h-h0)/(360*cos(del)*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
n:= m2 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*15);
Hor:=simp_angle(→HMS(Hor));
h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
deltaM:= HMS→((h-h0)/(360*cos(del)*cos(lat)*sin(Hor)));
tram:= ((m2 + deltaM) * 24); // setting
RETURN ({culm, sorg, tram});
END;
parametri(lista)
BEGIN
LOCAL dep, prec, α, δ, ε, listaparametri, ch, temp;
dep:= deltaEpsilonPsi(lista);
temp:= sunAlphaDelta(lista);
α:= sunlist(5); δ:= sunlist(6);
prec:= precession(α, δ);
ε:= epsilon(lista);
// 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));
listaparametri:= {ε, dep(1), dep(2), prec(1), prec(2)}; // ε, Δψ, Δε, Δα, Δδ
EpsilonPsi:= listaparametri;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
ascendant(lista)
BEGIN
LOCAL ASC, ts, segno, ε, listaascendente, ch, temp;
LOCAL tsb;
ε:= epsilon(lista);
ts:= calcTS(lista);
tsb:= ts(1)+lstd-dst; // ST of the birth
// temp:= sin(ts(2)*15)*cos(ε)+tan(lat)*sin(ε);
temp:= sin(tsb*15)*cos(ε)+tan(lat)*sin(ε);
ASC:= atan2(-cos(tsb*15), temp);
// IF temp<0 THEN ASC:= ASC+180 ELSE ASC:= ASC+360; END;
IF ASC<180 THEN ASC:=ASC+180 ELSE ASC:= ASC-180; END;
ASC:= simp_angle(ASC);
segno:= zodiaco(ASC);
listaascendente:= {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno};
Ascendant:= listaascendente;
RECT_P();
TEXTOUT_P("Ascendant " + trunc_angle(→HMS(ASC)), 25,50);
TEXTOUT_P(" " + trunc_angle(→HMS(ASC) MOD 30), 25,70);
TEXTOUT_P("in " + segno, 100,70);
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
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(lista)
// type dateshortlist
BEGIN
LOCAL ω, l, v, r, x, y, dz, ε;
LOCAL α, δ, dayY, ts, az, z, h;
LOCAL SAD, DA, Hd, lt;
LOCAL as, ac, zc, h0, H0, temp;
LOCAL ec, a, u, M1, M2, j, t, tau;
LOCAL JD, JDE, c, lapp, Ω, θ0;
LOCAL dep, deltapsi, longitude;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
a:= 1.00000023; // semiasse maggiore
dep:= deltaEpsilonPsi(lista);
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.003032028*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
//M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-(T^3)/24490000); // anom vera Meeus
M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // Meeus p 308
Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
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(lista); // 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:= transformEquatorialToHorizontal(Hd, lat, δ);
az:= HMS→(temp(1))+180; // Azimuth
IF az>360 THEN az:= az-360 END;
az:= →HMS(az);
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
sunlist:= {L,ω, l, lapp, α,δ, az, h, M, Ω, v, c,x,y,r, H, zc, dz, DA, SAD, E, lt, ec};
RETURN sunlist;
END;
sun(lista)
BEGIN
LOCAL ω, l, v, r, x, y, dz, ε;
LOCAL α, δ, dayY, ts, az, z, h;
LOCAL SAD, DA, sorg, tram, culm, Hd, h0, lt;
LOCAL as, ac, zc, segno, JD, JDE, v, c;
LOCAL Ω, lapp, ec, listasun, listasuntran, listasunothers;
LOCAL ch, transito, listaARh0, temp, Hor;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3;
ε:= epsilon(lista);
L:= sunlist(1); ω:= sunlist(2);
l:= sunlist(3); lapp:=sunlist(4);
α:= sunlist(5); δ:= sunlist(6);
az:= sunlist(7); h:= sunlist(8);
M:= sunlist(9); Ω:= sunlist(10);
v:= sunlist(11); c:= sunlist(12);
x:= sunlist(13); y:= sunlist(14);
r:= sunlist(15); Hor:= sunlist(16);
Hd:= 15*sunlist(16); zc:= sunlist(17);
dz:= sunlist(18);
DA:= sunlist(19); SAD:= sunlist(20);
E:= sunlist(21); lt:= sunlist(22);
ec:= sunlist(23);
ts:= calcTS(lista);
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
segno:= zodiaco(l);
// 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
// to pass lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
h0:= -0°50′; // 'standard' altitude of Sun
temp:= sunAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= sunAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= sunAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
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("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("λ longitudine apparente " + trunc_angle(lapp));
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(" ");
WAIT();
PRINT();
PRINT("Ω aberr, nutaz nodo "+ Ω);
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("Ang. orario "+ trunc_angle(Hor) + " " + 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(" ");
WAIT();
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);
l:= →HMS(l); lapp:= →HMS(lapp);
culm:= →HMS(culm); sorg:= →HMS(sorg); tram:= →HMS(tram);
listasun:= {l, α, δ, az, h, r};
listasuntran:= {culm, sorg, tram};
listasunothers:= {L, lapp, ω, Ω, M,v, c, x, y, Hor, lt, ε, ec};
Sun_Effemerids:= listasun;
Sun_Transit:= listasuntran;
Sun_Others:= listasunothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
sunAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= sunCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;
sunAh(lista)
// lista like dateshortlist
BEGIN
LOCAL lha,alpha,delta,phi, ah, alfa;
LOCAL temp,longitude, Azimuth, Altitude, ts;
temp:= sunAlphaDelta(lista);
alpha:=15*HMS→(temp(1));
alfa:= →HMS(temp(1)); // α in HMS
delta:= temp(2);
phi:=lat;
ts:= calcTS(lista);
longitude:= - HMS→(long); // - if long -E
lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
// lha:= 15*HMS→(ts(2)-alfa); // 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;
moonCalc(lista)
BEGIN
LOCAL ω,Ω, L1, M1, l, lsum, bsum, rsum;
LOCAL v, r, x, y, dep, ε;
LOCAL b, P, d, s, temp, tau;
LOCAL α,δ, ts, Hd, zc, dz, h, az;
LOCAL i, as, ac, DA, SAD, DD;
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(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
T:= (JDE-2451545.0)/36525; // use JDE
tau:= (JDE-2451545.0)/365250; // millenni - JDE
dep:= deltaEpsilonPsi(lista);
// Sun
L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // mean anomaly Sun, Meeus
// 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 + 280602 * 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*DD+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(lista, l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
ts:= calcTS(lista); // calcola tempo siderale
Hd:=15*HMS→(ts(2)-α); // local ST vs apparent ST
Hd:=simp_angle(→HMS(Hd));
H:= Hd/15; // in HMS
DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
SAD:= (90 + DA); // Semiarco Diurno
// SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
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)
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 Calante"; END;
IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
IF ph>100 AND ph<170 THEN fase:="Gibbosa Crescente"; 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
moonlist:= {L1, l, b, lapp, α, δ, az, h, M1, ω, Ω, d, DD, P, s, H, DA, SAD,ph, fase, eta, ill};
RETURN moonlist;
END;
moon(lista)
BEGIN
LOCAL ε, ts, JD, JDE, segno, L1, b;
LOCAL l, Ω, ω, M1, lapp, α, δ;
LOCAL az, h, sorg, tram, culm, Hd, DA, SAD;
LOCAL ph, fase, eta, DD, d, s, ill;
LOCAL dep, temp, listamoon, listamoontran, listamoonothers, ch;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
LOCAL listaARh0, transito, Hor;
ε:= epsilon(lista);
ts:= calcTS(lista); //
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
dep:= deltaEpsilonPsi(lista);
L1:= moonlist(1); l:= moonlist(2);
b:= moonlist(3);
α:= moonlist(5); δ:= moonlist(6);
az:= moonlist(7); h:= moonlist(8);
ω:= moonlist(10); Ω:= moonlist(11);
lapp:= moonlist(4); M1:= moonlist(9);
d:= moonlist(12); DD:= moonlist(13);
P:= moonlist(14); s:= moonlist(15);
Hor:= moonlist(16); Hd:= 15*moonlist(16);
DA:= moonlist(17); SAD:= moonlist(18);
ph:= moonlist(19); fase:= moonlist(20);
eta:= moonlist(21); ill:= moonlist(22);
segno:= zodiaco(l);
// transit routine
h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)
// h0:= 0°7′30″; // medium value: 0.125°
temp:= moonAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= moonAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= moonAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
PRINT;
PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("L Long media " + trunc_angle(L1));
PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("β Latitudine " + trunc_angle(b));
PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
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("P Parallasse " + trunc_angle(P));
PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");
PRINT(" ");
WAIT();
PRINT();
PRINT("ω Perigeo medio " + trunc_angle(ω));
PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
PRINT("M anomalia vera " + trunc_angle(M1));
PRINT("D Elongazione " + trunc_angle(DD));
PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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("Julian day of Ephemeris " + JDE);
PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT(" ");
WAIT();
PRINT();
PRINT("Δψ Nutazione long "+ dep(1));
PRINT("Δε Nutazione obl. " + dep(2));
PRINT("ε inclinazione ecl. " + ε);
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);
PRINT("Età " + eta);
PRINT("Parte illuminata " + ROUND(ill,5) + " " + ROUND(100*ill,2)+"%");
listamoon:= {→HMS(l), →HMS(b), ROUND(HMS→(d),6), →HMS(α), →HMS(δ), →HMS(az), →HMS(h)};
listamoonothers:= {L1, lapp, →HMS(α*15), M1,P, ω,Ω,s, Hor,ph, eta};
listamoontran:= →HMS({culm, sorg, tram});
Moon_Effemerids:= listamoon;
Moon_Transit:= listamoontran;
Moon_Others:= listamoonothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
moonAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= moonCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;
moonAh(lista)
// lista like dateshortlist
BEGIN
LOCAL lha,alpha,delta,phi, ah, alfa;
LOCAL temp,longitude, Azimuth, Altitude, ts;
temp:= sunAlphaDelta(lista);
alpha:=15*HMS→(temp(1));
alfa:= →HMS(temp(1)); // α in HMS
delta:= temp(2);
phi:=lat;
ts:= calcTS(lista);
longitude:= - HMS→(long); // - if long -E
lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
// lha:= 15*HMS→(ts(2)-alfa); // 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;
planetCalc(nameplanet, lista)
BEGIN
LOCAL ω,Ω, l, v, r, x, y, z;
LOCAL Ls, Ms, b, xs, ys, b_ec;
LOCAL x1, y1, z1, l_ec, 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, lapp;
LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
JD:= julianDay(lista);
ε:= epsilon(lista);
// calc of Sun (as basis)
sunCalc(lista);
xs:= sunlist(13); ys:= sunlist(14);
Ls:= sunlist(1); // Longitudine media
Ms:= sunlist(9); // 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(lista, l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
ts:= calcTS(lista); // 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
planetlist:= {L, l, b, lapp, α,δ, az, h, M, Ω, ω, v, x,y, z, x1, y1, z1,
r, R, H, zc, dz, DA, SAD, ec, in, a, θ, l_ec, b_ec, Ls, Ms, xs, ys, nameplanet};
RETURN planetlist;
END;
planetDisplay(nameplanet, lista)
BEGIN
LOCAL ω,Ω, l, v, r;
LOCAL x, y, z, x1, y1, z1;
LOCAL Ls, Ms, b, xs, ys, l_ec, b_ec;
LOCAL ec, in, a, δ, α, ts, az, h;
LOCAL DA, SAD, tram, sorg, culm, dayY;
LOCAL as, ac, zc, dz, θ, mov, segno, ε;
LOCAL transito, JD, JDE, temp, Hor, Hd;
LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
LOCAL listaARh0, lapp;
ts:= calcTS(lista);
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
L:= planetlist(1); l:= planetlist(2);
b:= planetlist(3); lapp:= planetlist(4);
α:= planetlist(5); δ:= planetlist(6);
az:= planetlist(7); h:= planetlist(8);
M:= planetlist(9); Ω:= planetlist(10);
ω:= planetlist(11); v:= planetlist(12);
x := planetlist(13); y:= planetlist(14); z:= planetlist(15);
x1:= planetlist(16); y1:= planetlist(17); z1:= planetlist(18);
r:= planetlist(19); R:= planetlist(20);
Hor:= planetlist(21); Hd:= 15*planetlist(21);
zc:= planetlist(22); dz:= planetlist(23);
DA:= planetlist(24); SAD:= planetlist(25);
ec:= planetlist(26); in:= planetlist(27);
a:= planetlist(28); θ:= planetlist(29);
l_ec:= planetlist(30); b_ec:= planetlist(31);
Ls:= planetlist(32); Ms:= planetlist(33); // Sun values
xs:= planetlist(34); ys:= planetlist(35);
// nameplanet:= planetlist(36);
segno:= zodiaco(l);
// transit routine
h0:= -0°34′; // 'standard' altitude of a planet
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
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("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("ω 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));
WAIT();
PRINT();
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(" ");
PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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(" ");
WAIT();
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);
listaplanet:= {l,b, α,δ,az,h, Hor, Hd, mov};
listaplanettran:= {culm, sorg, tram};
listaplanetothers:= {L,ω,Ω, M,l_ec,b_ec, v, x,y,z, x1,y1,z1, r, R, SAD, DA};
Planet_Effemerids:= listaplanet;
Planet_Transit:= listaplanettran;
Planet_Others:= listaplanetothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
planetAlphaDelta(nameplanet, lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= planetCalc(nameplanet, lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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: Effemeridi (Ephemeris) - salvomic - 07-01-2015 03:45 PM
New version, reviewed.
Corrected some bugs.
More user-friendly, I hope.
Transit routine returns more precise results, now.
Code:
data();
changeData();
whatToDo();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
tempo();
precession();
parametri();
transit();
sunCalc();
sun();
sunAh();
sunAlphaDelta();
moon();
moonCalc();
moonAh();
moonAlphaDelta();
planets();
planetCalc();
planetDisplay();
planetAlphaDelta();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendant();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();
EXPORT Tempi;
EXPORT EpsilonPsi;
EXPORT Ascendant;
EXPORT Sun_Effemerids;
EXPORT Sun_Transit;
EXPORT Sun_Others;
EXPORT Moon_Effemerids;
EXPORT Moon_Transit;
EXPORT Moon_Others;
EXPORT Planet_Effemerids;
EXPORT Planet_Transit;
EXPORT Planet_Others;
thisday:= 1964.0529;
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,40);
moonlist:= MAKELIST(X,X,1,40);
planetlist:= MAKELIST(X,X,1,40);
planet_const:= MAKELIST(X,X,1,40);
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, DrD
BEGIN
HAngle:=1;
data(); // input data and calc Sun
END;
data()
BEGIN
// initial input
LOCAL dd, mm, yy, hh, mn, ss, loct;
LOCAL ε, JD, JDE;
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};
ε:=epsilon(dateshortlist);
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
sunCalc(dateshortlist); // set initial data for Sun
moonCalc(dateshortlist);
whatToDo(dateshortlist);
END;
changeData()
BEGIN
// input to change date, lat, long...
LOCAL ε, JD, JDE;
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},{m,D,Y,lstd, dst, utc, long, lat, deltaT});
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};
ε:=epsilon(dateshortlist);
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
whatToDo(dateshortlist);
END;
whatToDo(lista)
BEGIN
LOCAL ch;
CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendant",
"Quit");
CASE
IF ch==1 THEN sunCalc(lista); sun(lista); END;
IF ch==2 THEN moonCalc(lista); moon(lista); END;
IF ch==3 THEN planets(lista); END;
IF ch==4 THEN tempo(lista); END;
IF ch==5 THEN parametri(lista); END;
IF ch==6 THEN ascendant(lista); END;
IF ch==7 THEN HAngle:=0; RETURN; END;
DEFAULT
END; // case
END;
planets(lista)
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), lista);
planetDisplay(nameplanet(ch),lista);
END;
epsilon(lista)
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
T:=(JD-2451545)/36525;
U:= T/100;
// ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
// 4680.93″ = 1°18′00.93″
ε0:= 23°26′21.448″ - (HMS→(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(lista);
eps:= dep(2);
ε:= ε0+eps; // true obliquity ε
RETURN ε;
END;
deltaEpsilonPsi(lista)
BEGIN
// Nutazione Δψ (longitudine), Δε (obliquità)
LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
LOCAL DD;
// JD:= julianDayAtZero(Y,m,D);
JDE:= julianEphemerisDay(lista);
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;
transformEclipticalToEquatorial(lista, lambda, beta)
BEGIN
// trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
LOCAL ε,Lambda, Beta, alpha, delta;
Lambda:=HMS→(lambda);
Beta:=HMS→(beta);
ε:= epsilon(lista);
alpha:=atan2(sin(Lambda)*cos(ε)-tan(Beta)*sin(ε),cos(Lambda));
IF alpha>=360 THEN alpha:=alpha-360 END;
IF alpha<0 THEN alpha:=alpha+360 END;
delta:=asin(sin(Beta)*cos(ε)+cos(Beta)*sin(ε)*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(lista, Lambda, Beta, Lha)
BEGIN
// transform λ and β (long, lat) into azimuth and height
LOCAL lambda, beta,phi, lha, ε;
LOCAL alpha,delta,azimuth,altitude;
lambda:=HMS→(Lambda); // longitudine astro
beta:=HMS→(Beta); // latitudine astro
ε:= epsilon(lista);
lha:=HMS→(Lha); // angolo orario
phi:=HMS→(lat); // lat geo
alpha:=atan2(SIN(lambda)*COS(ε)-TAN(beta)*SIN(ε),COS(lambda));
delta:=asin(sin(beta)*cos(ε)+cos(beta)*sin(ε)*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(lista)
BEGIN
LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
JD0:= julianDayAtZero(Y,m,D);
dep:= deltaEpsilonPsi(lista);
psi:=dep(1);
eps:=dep(2);
T:= (JD-2451545.0)/36525;
θ0:= meanSideralTime(lista); // 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(lista); // apparent ST (ST at Greenwitch at local our)
TSL:= apparentSideralTime(lista) - 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(lista);
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(lista);
MST:= meanSideralTime(lista);
dep:= deltaEpsilonPsi(lista);
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(lista)
BEGIN
LOCAL JD, JDE, TSid, dep, psi, eps;
LOCAL ε, listatempi, ch;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
TSid:= calcTS(lista);
dep:= deltaEpsilonPsi(lista);
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(" ");
WAIT();
PRINT();
PRINT("Δψ Nutazione (longit.) "+ psi);
PRINT("Δε Nutazione (obl.) "+ eps);
PRINT("ε Obliquità Ecl. "+ ε);
PRINT(" ");
PRINT("Time " + lstd + " GMT " + gmt);
listatempi:= {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)) };
Tempi:= listatempi;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
transit(lista, listaARh0)
BEGIN
// lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
// α, δ (AR, dec), height, h0 = standard height of body
LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;
LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
LOCAL alfa1, alfa3, delta1, delta3, h0, alpha, del, n;
α:= listaARh0(2); // AR for day
δ:= listaARh0(5); // decl for day
alfa1:= listaARh0(1); // AR day-1
alfa3:= listaARh0(3); // AR day+1
delta1:= listaARh0(4); // decl day-1
delta3:= listaARh0(6); // decl day+1
h0:= listaARh0(7); // standard height
temp:= calcTS(lista);
θ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
n:= m0 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*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
n:= m1 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*15);
Hor:=simp_angle(→HMS(Hor));
h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
deltaM:= HMS→((h-h0)/(360*cos(del)*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
n:= m2 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*15);
Hor:=simp_angle(→HMS(Hor));
h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
deltaM:= HMS→((h-h0)/(360*cos(del)*cos(lat)*sin(Hor)));
tram:= ((m2 + deltaM) * 24); // setting
RETURN ({culm, sorg, tram});
END;
parametri(lista)
BEGIN
LOCAL dep, prec, α, δ, ε, listaparametri, ch, temp;
dep:= deltaEpsilonPsi(lista);
temp:= sunAlphaDelta(lista);
α:= sunlist(5); δ:= sunlist(6);
prec:= precession(α, δ);
ε:= epsilon(lista);
// 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));
listaparametri:= {ε, dep(1), dep(2), prec(1), prec(2)}; // ε, Δψ, Δε, Δα, Δδ
EpsilonPsi:= listaparametri;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
ascendant(lista)
BEGIN
LOCAL ASC, ts, segno, ε, listaascendente, ch, temp;
LOCAL tsb;
ε:= epsilon(lista);
ts:= calcTS(lista);
tsb:= ts(1)+lstd-dst; // ST of the birth
// temp:= sin(ts(2)*15)*cos(ε)+tan(lat)*sin(ε);
temp:= sin(tsb*15)*cos(ε)+tan(lat)*sin(ε);
ASC:= atan2(-cos(tsb*15), temp);
// IF temp<0 THEN ASC:= ASC+180 ELSE ASC:= ASC+360; END;
IF ASC<180 THEN ASC:=ASC+180 ELSE ASC:= ASC-180; END;
ASC:= simp_angle(ASC);
segno:= zodiaco(ASC);
listaascendente:= {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno};
Ascendant:= listaascendente;
RECT_P();
TEXTOUT_P("Ascendant " + trunc_angle(→HMS(ASC)), 25,50);
TEXTOUT_P(" " + trunc_angle(→HMS(ASC) MOD 30), 25,70);
TEXTOUT_P("in " + segno, 100,70);
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
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(lista)
// type dateshortlist
BEGIN
LOCAL ω, l, v, r, x, y, dz, ε;
LOCAL α, δ, dayY, ts, az, z, h;
LOCAL SAD, DA, Hd, lt;
LOCAL as, ac, zc, h0, H0, temp;
LOCAL ec, a, u, M1, M2, j, t, tau;
LOCAL JD, JDE, c, lapp, Ω, θ0;
LOCAL dep, deltapsi, longitude;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
a:= 1.00000023; // semiasse maggiore
dep:= deltaEpsilonPsi(lista);
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.003032028*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
//M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-(T^3)/24490000); // anom vera Meeus
M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // Meeus p 308
Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
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(lista); // 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:= transformEquatorialToHorizontal(Hd, lat, δ);
az:= HMS→(temp(1))+180; // Azimuth
IF az>360 THEN az:= az-360 END;
az:= →HMS(az);
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
sunlist:= {L,ω, l, lapp, α,δ, az, h, M, Ω, v, c,x,y,r, H, zc, dz, DA, SAD, E, lt, ec};
RETURN sunlist;
END;
sun(lista)
BEGIN
LOCAL ω, l, v, r, x, y, dz, ε;
LOCAL α, δ, dayY, ts, az, z, h;
LOCAL SAD, DA, sorg, tram, culm, Hd, h0, lt;
LOCAL as, ac, zc, segno, JD, JDE, v, c;
LOCAL Ω, lapp, ec, listasun, listasuntran, listasunothers;
LOCAL ch, transito, listaARh0, temp, Hor, dep;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3;
ε:= epsilon(lista);
L:= sunlist(1); ω:= sunlist(2);
l:= sunlist(3); lapp:=sunlist(4);
α:= sunlist(5); δ:= sunlist(6);
az:= sunlist(7); h:= sunlist(8);
M:= sunlist(9); Ω:= sunlist(10);
v:= sunlist(11); c:= sunlist(12);
x:= sunlist(13); y:= sunlist(14);
r:= sunlist(15); Hor:= sunlist(16);
Hd:= 15*sunlist(16); zc:= sunlist(17);
dz:= sunlist(18);
DA:= sunlist(19); SAD:= sunlist(20);
E:= sunlist(21); lt:= sunlist(22);
ec:= sunlist(23);
ts:= calcTS(lista);
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
dep:= deltaEpsilonPsi(lista);
segno:= zodiaco(l);
// 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
// to pass lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
h0:= -0°50′; // 'standard' altitude of Sun
temp:= sunAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= sunAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= sunAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
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("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("λ longitudine apparente " + trunc_angle(lapp));
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(" ");
WAIT();
PRINT();
PRINT("Ω aberr, nutaz nodo "+ Ω);
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("Ang. orario "+ trunc_angle(Hor) + " " + 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(" ");
WAIT();
PRINT();
PRINT("Transita a: " + trunc_angle(culm MOD 24) + " 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("ε obliquità ecl. " + ε);
PRINT("ec eccentricità "+ ec);
PRINT("Δψ Nutazione long "+ dep(1));
PRINT("Δε Nutazione obl. " + dep(2));
PRINT(" ");
PRINT("Long elioc. Terra "+ trunc_angle(lt));
l:= →HMS(l); lapp:= →HMS(lapp);
culm:= →HMS(culm); sorg:= →HMS(sorg); tram:= →HMS(tram);
listasun:= {l, α, δ, az, h, r};
listasuntran:= {culm, sorg, tram};
listasunothers:= {L, lapp, ω, Ω, M,v, c, x, y, Hor, lt, ε, ec};
Sun_Effemerids:= listasun;
Sun_Transit:= listasuntran;
Sun_Others:= listasunothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
sunAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= sunCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;
sunAh(lista)
// lista like dateshortlist
BEGIN
LOCAL lha,alpha,delta,phi, ah, alfa;
LOCAL temp,longitude, Azimuth, Altitude, ts;
temp:= sunAlphaDelta(lista);
alpha:=15*HMS→(temp(1));
alfa:= →HMS(temp(1)); // α in HMS
delta:= temp(2);
phi:=lat;
ts:= calcTS(lista);
longitude:= HMS→(long); // - if long -E
lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
// lha:= 15*HMS→(ts(2)-alfa); // 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;
moonCalc(lista)
BEGIN
LOCAL ω,Ω, Lprime, Mprime, l, lsum, bsum, rsum;
LOCAL v, r, x, y, dep, ε;
LOCAL b, P, d, s, temp, tau;
LOCAL α,δ, ts, Hd, zc, dz, h, az;
LOCAL i, as, ac, DA, SAD, DD;
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(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
T:= (JDE-2451545.0)/36525; // use JDE
tau:= (JDE-2451545.0)/365250; // millenni - JDE
dep:= deltaEpsilonPsi(lista);
// Sun
L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // mean anomaly Sun, Meeus
// end Sun
// L1:= (33.231 + 13.17639653*N) MOD 360; // long. media
Lprime:= 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
Mprime:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media
ω:= (Lprime - Mprime) MOD 360; // perigeo medio
// DD:= Lprime - L; F:= Lprime-Ω; // 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(Mprime); // calc Moon longitude
lsum:= lsum + 1274027 * sin(2*DD-Mprime); // evezione
lsum:= lsum + 658314 * sin(2*DD); // variazione
lsum:= lsum + 213618 * sin(2*Mprime); // 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*Mprime);
lsum:= lsum + 57066 * sin(2*DD-Mprime-M) *ee;
lsum:= lsum + 53322 * sin(2*DD+Mprime);
lsum:= lsum + 45758 * sin(2*DD-M) *ee;
lsum:= lsum - 40923 * sin(M-Mprime) *ee;
lsum:= lsum - 34720 * sin(DD);
lsum:= lsum - 30383 * sin(Mprime+M) *ee;
lsum:= lsum + 15327 * sin(2*DD-2*F);
lsum:= lsum - 12528 * sin(Mprime+2*F);
lsum:= lsum + 10980 * sin(Mprime-2*F);
lsum:= lsum + 10675 * sin(4*DD-Mprime);
lsum:= lsum + 10034 * sin(3*Mprime);
lsum:= lsum + 8548 * sin(4*DD-2*Mprime);
lsum:= lsum - 7888 * sin(2*DD+M-Mprime) *ee;
lsum:= lsum - 6766 * sin(2*DD+M) *ee;
lsum:= lsum - 5163 * sin(DD-Mprime);
lsum:= lsum + 4987 * sin(DD+M) *ee;
lsum:= lsum + 4036 * sin(2*DD-M+Mprime) *ee;
lsum:= lsum + 3994 * sin(2*DD+2*Mprime);
lsum:= lsum + 3861 * sin(4*DD);
lsum:= lsum + 3665 * sin(2*DD-3*Mprime);
lsum:= lsum - 2689 * sin(M-2*Mprime) *ee;
lsum:= lsum - 2602 * sin(2*DD-Mprime+2*F);
lsum:= lsum + 2390 * sin(2*DD-M-2*Mprime) *ee;
lsum:= lsum - 2348 * sin(DD+Mprime);
lsum:= lsum + 2236 * sin(2*DD-2*M) *ee^2;
lsum:= lsum - 2120 * sin(M+2*Mprime) *ee;
lsum:= lsum - 2069 * sin(2*M) *ee^2;
lsum:= lsum + 2048 * sin(2*DD-2*M-Mprime) *ee^2;
lsum:= lsum - 1773 * sin(2*DD+Mprime-2*F);
lsum:= lsum - 1595 * sin(2*DD+2*F);
lsum:= lsum + 1215 * sin(4*DD-M-Mprime) *ee;
lsum:= lsum - 1110 * sin(2*Mprime+2*F);
lsum:= lsum - 892 * sin(3*DD-Mprime);
lsum:= lsum - 810 * sin(2*DD+M+Mprime) *ee;
lsum:= lsum + 759 * sin(4*DD-M-2*Mprime) *ee;
lsum:= lsum - 713 * sin(2*M-Mprime) * ee^2;
lsum:= lsum - 700 * sin(2*DD+2*M-Mprime) *ee^2;
lsum:= lsum + 691 * sin(2*DD+M-2*Mprime) *ee;
lsum:= lsum + 596 * sin(2*DD-M-2*F) *ee;
lsum:= lsum + 549 * sin(4*DD+Mprime);
lsum:= lsum + 537 * sin(4*Mprime);
lsum:= lsum + 520 * sin(4*DD-M) *ee;
lsum:= lsum - 487 * sin(DD-2*Mprime);
lsum:= lsum - 399 * sin(2*DD+M-2*F) *ee;
lsum:= lsum - 381 * sin(2*Mprime-2*F);
lsum:= lsum + 351 * sin(DD+M+Mprime) *ee;
lsum:= lsum - 340 * sin(3*DD-2*Mprime);
lsum:= lsum + 330 * sin(4*DD-3*Mprime);
lsum:= lsum + 327 * sin(2*DD-M+2*Mprime) *ee;
lsum:= lsum - 323 * sin(2*M+Mprime) *ee^2;
lsum:= lsum + 299 * sin(DD+M-Mprime) *ee;
lsum:= lsum + 294 * sin(2*DD+3*Mprime);
lsum:= lsum + 3958*sin(a1)+1962*sin(Lprime-F)+318*sin(a2);
l:= (simp_angle(Lprime + 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 + 280602 * sin(Mprime+F);
bsum:= bsum + 277693 * sin(Mprime-F); // equazione del centro
bsum:= bsum + 173237 * sin(2*DD-F); // grande ineguaglianza
bsum:= bsum + 55413 * sin(2*DD-Mprime+F); // evezione
bsum:= bsum + 46271 * sin(2*DD-Mprime-F); // evezione (2)
bsum:= bsum + 32573 * sin(2*DD+F);
bsum:= bsum + 17198 * sin(2*Mprime+F);
bsum:= bsum + 9266 * sin(2*DD+Mprime-F);
bsum:= bsum + 8822 * sin(2*Mprime-F);
bsum:= bsum + 8216 * sin(2*DD-M-F) *ee;
bsum:= bsum + 4324 * sin(2*DD-2*Mprime-F);
bsum:= bsum + 4200 * sin(2*DD+Mprime+F);
bsum:= bsum - 3359 * sin(2*DD+M-F) *ee;
bsum:= bsum + 2463 * sin(2*DD-M-Mprime+F) *ee;
bsum:= bsum + 2211 * sin(2*DD-M+F) *ee;
bsum:= bsum + 2065 * sin(2*DD-M-Mprime-F) *ee;
bsum:= bsum - 1870 * sin(M-Mprime-F) *ee;
bsum:= bsum + 1828 * sin(4*DD-Mprime-F);
bsum:= bsum - 1794 * sin(M+F) *ee;
bsum:= bsum - 1749 * sin(3*F);
bsum:= bsum - 1565 * sin(M-Mprime+F) *ee;
bsum:= bsum - 1491 * sin(DD+F);
bsum:= bsum - 1475 * sin(M+Mprime+F) *ee;
bsum:= bsum - 1410 * sin(M+Mprime-F) *ee;
bsum:= bsum - 1344 * sin(M-F) *ee;
bsum:= bsum - 1335 * sin(DD-F);
bsum:= bsum + 1107 * sin(3*Mprime+F);
bsum:= bsum + 1021 * sin(4*DD-F);
bsum:= bsum + 833 * sin(4*DD-Mprime+F);
bsum:= bsum + 777 * sin(Mprime-3*F);
bsum:= bsum + 671 * sin(4*DD-2*Mprime+F);
bsum:= bsum + 607 * sin(2*DD-3*F);
bsum:= bsum + 596 * sin(2*DD+2*Mprime-F);
bsum:= bsum + 491 * sin(2*DD-M+Mprime-F) *ee;
bsum:= bsum - 451 * sin(2*DD-2*Mprime+F);
bsum:= bsum + 439 * sin(3*Mprime-F);
bsum:= bsum + 422 * sin(2*DD+2*Mprime+F);
bsum:= bsum + 421 * sin(2*DD-3*Mprime-F);
bsum:= bsum - 366 * sin(2*DD+M-Mprime+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+Mprime+F) *ee;
bsum:= bsum + 302 * sin(2*DD-2*M-F) *ee^2;
bsum:= bsum - 283 * sin(Mprime+3*F);
bsum:= bsum - 229 * sin(2*DD+M+Mprime-F) *ee;
bsum:= bsum + 223 * sin(DD+M-F) *ee;
bsum:= bsum + 223 * sin(DD+M+F) *ee;
bsum:= bsum - 220 * sin(M-2*Mprime-F) *ee;
bsum:= bsum - 220 * sin(2*DD+M-Mprime-F) *ee;
bsum:= bsum - 185 * sin(DD+Mprime+F);
bsum:= bsum + 181 * sin(2*DD-M-2*Mprime-F) *ee;
bsum:= bsum - 177 * sin(M+2*Mprime+F) *ee;
bsum:= bsum + 176 * sin(4*DD-2*Mprime-F);
bsum:= bsum + 166 * sin(4*DD-M-Mprime-F) *ee;
bsum:= bsum - 164 * sin(DD+Mprime-F);
bsum:= bsum + 132 * sin(4*DD+Mprime-F);
bsum:= bsum - 119 * sin(DD-Mprime-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(Lprime))+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(Lprime-Mprime)-115*sin(Lprime+Mprime);
b := →HMS(bsum/1000000); // β div by 1 milion
// correzioni per la distanza dalla Terra (Δ)
rsum:=0;
rsum:= rsum -20905355 * cos(Mprime); // coseni
rsum:= rsum - 3699111 * cos(2*DD-Mprime); // evezione
rsum:= rsum - 2955968 * cos(2*DD); // variazione
rsum:= rsum - 569925 * cos(2*Mprime); // 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*Mprime);
rsum:= rsum - 152138 * cos(2*DD-Mprime-M) *ee;
rsum:= rsum - 170733 * cos(2*DD+Mprime);
rsum:= rsum - 204586 * cos(2*DD-M) *ee;
rsum:= rsum - 129620 * cos(M-Mprime) *ee;
rsum:= rsum + 108743 * cos(DD);
rsum:= rsum + 104755 * cos(Mprime+M) *ee;
rsum:= rsum + 10321 * cos(2*DD-2*F);
rsum:= rsum + 79661 * cos(Mprime-2*F);
rsum:= rsum - 34782 * cos(4*DD-Mprime);
rsum:= rsum - 23210 * cos(3*Mprime);
rsum:= rsum - 21636 * cos(4*DD-2*Mprime);
rsum:= rsum + 24208 * cos(2*DD+M-Mprime) *ee;
rsum:= rsum + 30824 * cos(2*DD+M) *ee;
rsum:= rsum - 8379 * cos(DD-Mprime);
rsum:= rsum - 16675 * cos(DD+M) *ee;
rsum:= rsum - 12831 * cos(2*DD-M+Mprime) *ee;
rsum:= rsum - 10445 * cos(2*DD+2*Mprime);
rsum:= rsum - 11650 * cos(4*DD);
rsum:= rsum + 14403 * cos(2*DD-3*Mprime);
rsum:= rsum - 7003 * cos(M-2*Mprime) *ee;
rsum:= rsum + 10056 * cos(2*DD-M-2*Mprime) * ee;
rsum:= rsum + 6322 * cos(DD+Mprime);
rsum:= rsum - 9884 * cos(2*DD-2*M) * ee^2;
rsum:= rsum + 5751 * cos(M+2*Mprime) *ee;
rsum:= rsum - 4950 * cos(2*DD-2*M-Mprime) *ee^2;
rsum:= rsum + 4130 * cos(2*DD+Mprime-2*F);
rsum:= rsum - 3958 * cos(4*DD-M-Mprime) *ee;
rsum:= rsum + 3258 * cos(3*DD-Mprime);
rsum:= rsum + 2616 * cos(2*DD+M+Mprime) *ee;
rsum:= rsum - 1897 * cos(4*DD-M-2*Mprime) *ee;
rsum:= rsum - 2117 * cos(2*M-Mprime) * ee^2;
rsum:= rsum + 2354 * cos(2*DD+2*M-Mprime) *ee^2;
rsum:= rsum - 1423 * cos(4*DD+Mprime);
rsum:= rsum - 1117 * cos(4*Mprime);
rsum:= rsum - 1571 * cos(4*DD-M) *ee;
rsum:= rsum - 1739 * cos(DD-2*Mprime);
rsum:= rsum - 4421 * cos(2*Mprime-2*F);
rsum:= rsum + 1165 * cos(2*M+Mprime) *ee^2;
rsum:= rsum + 8752 * cos(2*DD-Mprime-2*F);
rsum:= rsum/1000; // div by 1000
temp:= transformEclipticalToEquatorial(lista, l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
ts:= calcTS(lista); // calcola tempo siderale
Hd:=15*HMS→(ts(2)-α); // local ST vs apparent ST
Hd:=simp_angle(→HMS(Hd));
H:= Hd/15; // in HMS
DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
SAD:= (90 + DA); // Semiarco Diurno
// SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
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)
// s:= 358473400/d; // semidiametro (in " arcsec)
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 Calante"; END;
IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
IF ph>100 AND ph<170 THEN fase:="Gibbosa Crescente"; 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
moonlist:= {Lprime, l, b, lapp, α, δ, az, h, Mprime, ω, Ω, d, DD, P, s, H, DA, SAD,ph, fase, eta, ill};
RETURN moonlist;
END;
moon(lista)
BEGIN
LOCAL ε, ts, JD, JDE, segno, Lprime, b;
LOCAL l, Ω, ω, Mprime, lapp, α, δ;
LOCAL az, h, sorg, tram, culm, Hd, DA, SAD;
LOCAL ph, fase, eta, DD, d, s, ill;
LOCAL dep, temp, listamoon, listamoontran, listamoonothers, ch;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
LOCAL listaARh0, transito, Hor;
ε:= epsilon(lista);
ts:= calcTS(lista); //
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
dep:= deltaEpsilonPsi(lista);
Lprime:= moonlist(1); l:= moonlist(2);
b:= moonlist(3);
α:= moonlist(5); δ:= moonlist(6);
az:= moonlist(7); h:= moonlist(8);
ω:= moonlist(10); Ω:= moonlist(11);
lapp:= moonlist(4); Mprime:= moonlist(9);
d:= moonlist(12); DD:= moonlist(13);
P:= moonlist(14); s:= moonlist(15);
Hor:= moonlist(16); Hd:= 15*moonlist(16);
DA:= moonlist(17); SAD:= moonlist(18);
ph:= moonlist(19); fase:= moonlist(20);
eta:= moonlist(21); ill:= moonlist(22);
segno:= zodiaco(l);
// transit routine
h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)
// h0:= 0°7′30″; // medium value: 0.125°
temp:= moonAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= moonAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= moonAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
PRINT;
PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("L Long media " + trunc_angle(Lprime));
PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("β Latitudine " + trunc_angle(b));
PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
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("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(" ");
WAIT();
PRINT();
PRINT("ω Perigeo medio " + trunc_angle(ω));
PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
PRINT("M anomalia vera " + trunc_angle(Mprime));
PRINT("D Elongazione " + trunc_angle(DD));
PRINT(" ");
PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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("Julian day of Ephemeris " + JDE);
PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT(" ");
WAIT();
PRINT();
PRINT("Transita a: " + trunc_angle(culm MOD 24 ) + " 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("ε obliquità ecl. " + ε);
PRINT("Δψ Nutazione long "+ dep(1));
PRINT("Δε Nutazione obl. " + dep(2));
PRINT(" ");
PRINT("Fase " + trunc_angle(ph) + " " + fase);
PRINT("Età " + eta);
PRINT("Parte illuminata " + ROUND(ill,5) + " " + ROUND(100*ill,2)+"%");
listamoon:= {→HMS(l), →HMS(b), ROUND(HMS→(d),6), →HMS(α), →HMS(δ), →HMS(az), →HMS(h)};
listamoonothers:= {→HMS(Lprime), →HMS(lapp), →HMS(Mprime),→HMS(ω),→HMS(Ω), →HMS(P),
→HMS(s), →HMS(Hor), →HMS(α*15), ph, eta};
listamoontran:= →HMS({culm, sorg, tram});
Moon_Effemerids:= listamoon;
Moon_Transit:= listamoontran;
Moon_Others:= listamoonothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
moonAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= moonCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;
moonAh(lista)
// lista like dateshortlist
BEGIN
LOCAL lha,alpha,delta,phi, ah, alfa;
LOCAL temp,longitude, Azimuth, Altitude, ts;
temp:= sunAlphaDelta(lista);
alpha:=15*HMS→(temp(1));
alfa:= →HMS(temp(1)); // α in HMS
delta:= temp(2);
phi:=lat;
ts:= calcTS(lista);
longitude:= - HMS→(long); // - if long -E
lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
// lha:= 15*HMS→(ts(2)-alfa); // 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;
planetCalc(nameplanet, lista)
BEGIN
LOCAL ω,Ω, l, v, r, x, y, z;
LOCAL Ls, Ms, b, xs, ys, b_ec;
LOCAL x1, y1, z1, l_ec, 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, lapp, dep;
LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
JD:= julianDay(lista);
ε:= epsilon(lista);
dep:= deltaEpsilonPsi(lista);
// calc of Sun (as basis)
sunCalc(lista);
xs:= sunlist(13); ys:= sunlist(14);
Ls:= sunlist(1); // Longitudine media
Ms:= sunlist(9); // 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
lapp:= →HMS(l + dep(1)); // longitudine apparente (+ nutazione long, Δψ)
temp:= transformEclipticalToEquatorial(lista, l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
ts:= calcTS(lista); // 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
planetlist:= {L, l, b, lapp, α,δ, az, h, M, Ω, ω, v, x,y, z, x1, y1, z1,
r, R, H, zc, dz, DA, SAD, ec, in, a, θ, l_ec, b_ec, Ls, Ms, xs, ys, nameplanet};
RETURN planetlist;
END;
planetDisplay(nameplanet, lista)
BEGIN
LOCAL ω,Ω, l, v, r;
LOCAL x, y, z, x1, y1, z1;
LOCAL Ls, Ms, b, xs, ys, l_ec, b_ec;
LOCAL ec, in, a, δ, α, ts, az, h;
LOCAL DA, SAD, tram, sorg, culm, dayY;
LOCAL as, ac, zc, dz, θ, mov, segno, ε;
LOCAL transito, JD, JDE, temp, Hor, Hd;
LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
LOCAL listaARh0, lapp;
ts:= calcTS(lista);
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
L:= planetlist(1); l:= planetlist(2);
b:= planetlist(3); lapp:= planetlist(4);
α:= planetlist(5); δ:= planetlist(6);
az:= planetlist(7); h:= planetlist(8);
M:= planetlist(9); Ω:= planetlist(10);
ω:= planetlist(11); v:= planetlist(12);
x := planetlist(13); y:= planetlist(14); z:= planetlist(15);
x1:= planetlist(16); y1:= planetlist(17); z1:= planetlist(18);
r:= planetlist(19); R:= planetlist(20);
Hor:= planetlist(21); Hd:= 15*planetlist(21);
zc:= planetlist(22); dz:= planetlist(23);
DA:= planetlist(24); SAD:= planetlist(25);
ec:= planetlist(26); in:= planetlist(27);
a:= planetlist(28); θ:= planetlist(29);
l_ec:= planetlist(30); b_ec:= planetlist(31);
Ls:= planetlist(32); Ms:= planetlist(33); // Sun values
xs:= planetlist(34); ys:= planetlist(35);
// nameplanet:= planetlist(36);
segno:= zodiaco(l);
// transit routine
h0:= -0°34′; // 'standard' altitude of a planet
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
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(" ");
PRINT("L Long media " + trunc_angle(L));
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) + " (UA) ");
PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
PRINT(" ");
PRINT("ω Perielio " + trunc_angle(ω));
PRINT("Ω Nodo " + trunc_angle(Ω));
PRINT("M Anomalia vera " + trunc_angle(M));
WAIT();
PRINT();
PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));
PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(simp_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(" ");
PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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);
WAIT();
PRINT();
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);
listaplanet:= {l,b, R, α, δ, az, h};
listaplanettran:= {culm, sorg, tram};
listaplanetothers:= {L, lapp, l_ec,b_ec, M, ω,Ω, v, x,y,z, x1,y1,z1, r, Hor, Hd, SAD, DA, mov};
Planet_Effemerids:= listaplanet;
Planet_Transit:= listaplanettran;
Planet_Others:= listaplanetothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
planetAlphaDelta(nameplanet, lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= planetCalc(nameplanet, lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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: Effemeridi (Ephemeris) - Thomas_Sch - 07-01-2015 08:47 PM
Hi Salvo,
Many thanks for sharing your programs!
They are very useful, and in addition they show nicely the possibilities of the Prime.
edit:
May I express a few wishes or suggestions?
- it would be helpful, if in successive calculations the time in the input fields is updated
- is it possible to build an input dialog for user specifiv local coordinates?
Alternatively, variables or constants for lat and long would be useful at the beginning of the program,
so the 'customization' would be easier.
Again many thanks!
Thomas
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-01-2015 09:12 PM
(07-01-2015 08:47 PM)Thomas_Sch Wrote: Hi Salvo,
Many thanks for sharing your programs!
They are very useful, and in addition they show nicely the possibilities of the Prime.
Thomas
thank you Thomas,
I'm happy to share, and yes, the Prime are an "open book" of possibilities, very amazing to program
Salvo
RE: Astronomy: Effemeridi (Ephemeris) - Thomas_Sch - 07-01-2015 09:14 PM
(07-01-2015 09:12 PM)salvomic Wrote: thank you Thomas,
I'm happy to share, and yes, the Prime are an "open book" of possibilities, very amazing to program
Salvo
Sorry, I did some editing to my posting later:
edit:
May I express a few wishes or suggestions?
- it would be helpful, if in successive calculations the time in the input fields is updated
- is it possible to build an input dialog for user specifiv local coordinates?
Alternatively, variables or constants for lat and long would be useful at the beginning of the program,
so the 'customization' would be easier.
Again many thanks!
Thomas
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-01-2015 09:24 PM
(07-01-2015 09:14 PM)Thomas_Sch Wrote: Sorry, I did some editing to my posting later:
edit:
May I express a few wishes or suggestions?
- it would be helpful, if in successive calculations the time in the input fields is updated
- is it possible to build an input dialog for user specifiv local coordinates?
Alternatively, variables or constants for lat and long would be useful at the beginning of the program,
so the 'customization' would be easier.
ok,thank you; I'll think about these tips.
For now the input works so: at the start of program it calculates the actual date and time and the user can change the default lat/long (for now mine, while I'm still debugging), the program calculates Sun and Moon (to set some variables), then asks to the user for another calculation, and then save the old values (lat, long, date, time) or the other changes the user does.
This is useful to make different calculations with the same data.
I need help to make an option to choose automatically update, if the user chooses it.
It would be also nice to have automatic calculation of DST (but I've no idea how todo this in the Prime)...
I'll think also a way to save long/lat as in the app of Marcel, also if this is for now a program, not an app, that uses AVars...
Salvo
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-01-2015 09:48 PM
...in the meantime this version introduce other values for planets: parallax, illuminated fraction, magnitudo, phase.
EDIT: Thomas, try this with an option to save coordinates at the starting, the program should save them in Long_Lat var in the program variables and load them when used again...
I need a better safeguard, however, then only control if the size is ...2 ;-)
Code:
data();
changeData();
whatToDo();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEquatorialToEcliptical();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
transformHorizontalToEquatorial();
tempo();
precession();
parametri();
transit();
sunCalc();
sun();
sunAh();
sunAlphaDelta();
moon();
moonCalc();
moonAh();
moonAlphaDelta();
planets();
planetCalc();
planetDisplay();
planetAlphaDelta();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendant();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();
EXPORT Tempi;
EXPORT EpsilonPsi;
EXPORT Ascendant;
EXPORT Sun_Effemerids;
EXPORT Sun_Transit;
EXPORT Sun_Others;
EXPORT Moon_Effemerids;
EXPORT Moon_Transit;
EXPORT Moon_Others;
EXPORT Planet_Effemerids;
EXPORT Planet_Transit;
EXPORT Planet_Others;
EXPORT Long_Lat;
thisday:= 1964.0529;
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,40);
moonlist:= MAKELIST(X,X,1,40);
planetlist:= MAKELIST(X,X,1,40);
planet_const:= MAKELIST(X,X,1,40);
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, DrD
BEGIN
HAngle:=1;
data(); // input data and calc Sun
END;
data()
BEGIN
// initial input
LOCAL dd, mm, yy, hh, mn, ss, loct, save;
LOCAL ε, JD, JDE;
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);
IF (SIZE(Long_Lat)==2) THEN long:= Long_Lat(1); lat:= Long_Lat(2);
ELSE long:= →HMS(-14°43′41″); lat:= →HMS(36°55′15″); END;
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}},
{save,0,{85,15,5}} },
"Data: Use Shift+9 for H°M′S″",
{"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",
"Long (-E):","Lat (+N)", "delta T", "Save coord?"},
{"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.)", "Save Lat/Long in a var?"},
{1,1,1901,0, 0,0, 0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, long, lat, 69.3, 1});
// Ragusa (-14°43′41″, 36°55′15″) -
// Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)
IF save==1 THEN Long_Lat:= {long, lat}; END;
gmt:= (lstd-utc-(1*dst)) MOD 24; // 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};
ε:=epsilon(dateshortlist);
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
// sunCalc(dateshortlist); // set initial data for Sun
// moonCalc(dateshortlist);
whatToDo(dateshortlist);
END;
changeData()
BEGIN
// input to change date, lat, long...
LOCAL ε, JD, JDE;
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},{m,D,Y,lstd, dst, utc, long, lat, deltaT});
gmt:= (lstd-utc-(1*dst)) MOD 24; // 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};
ε:=epsilon(dateshortlist);
JD:= julianDay(dateshortlist);
JDE:= julianEphemerisDay(dateshortlist);
T:=(JD-2451545)/36525;
U:= T/100;
whatToDo(dateshortlist);
END;
whatToDo(lista)
BEGIN
LOCAL ch;
CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendant",
"Quit");
CASE
IF ch==1 THEN sunCalc(lista); sun(lista); END;
IF ch==2 THEN moonCalc(lista); moon(lista); END;
IF ch==3 THEN planets(lista); END;
IF ch==4 THEN tempo(lista); END;
IF ch==5 THEN parametri(lista); END;
IF ch==6 THEN ascendant(lista); END;
IF ch==7 THEN HAngle:=0; RETURN; END;
DEFAULT
END; // case
END;
planets(lista)
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), lista);
planetDisplay(nameplanet(ch),lista);
END;
epsilon(lista)
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
T:=(JD-2451545)/36525;
U:= T/100;
// ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
// 4680.93″ = 1°18′00.93″
ε0:= 23°26′21.448″ - (HMS→(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(lista);
eps:= dep(2);
ε:= ε0+eps; // true obliquity ε
RETURN ε;
END;
deltaEpsilonPsi(lista)
BEGIN
// Nutazione Δψ (longitudine), Δε (obliquità)
LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
LOCAL DD;
// JD:= julianDayAtZero(Y,m,D);
JDE:= julianEphemerisDay(lista);
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;
transformEquatorialToEcliptical(lista, alfa, delta)
BEGIN
LOCAL Alpha,Delta,ε,Lambda,Beta, lambda_beta;
Alpha:= HMS→(alfa)*15;
Delta:= HMS→(delta);
ε:= epsilon(lista);
Lambda:= atan2(sin(Alpha)*cos(ε)+tan(Delta)*sin(ε),cos(Alpha));
Beta:= asin(sin(Delta)*cos(ε)-cos(Delta)*sin(ε)*sin(Alpha));
lambda_beta:={→HMS(Lambda),→HMS(Beta)};
RETURN (lambda_beta);
END;
transformEclipticalToEquatorial(lista, lambda, beta)
BEGIN
// trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
LOCAL ε,Lambda, Beta, alpha, delta;
Lambda:=HMS→(lambda);
Beta:=HMS→(beta);
ε:= epsilon(lista);
alpha:=atan2(sin(Lambda)*cos(ε)-tan(Beta)*sin(ε),cos(Lambda));
IF alpha>=360 THEN alpha:=alpha-360 END;
IF alpha<0 THEN alpha:=alpha+360 END;
delta:=asin(sin(Beta)*cos(ε)+cos(Beta)*sin(ε)*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(lista, Lambda, Beta, Lha)
BEGIN
// transform λ and β (long, lat) into azimuth and height
LOCAL lambda, beta,phi, lha, ε;
LOCAL alpha,delta,azimuth,altitude;
lambda:=HMS→(Lambda); // longitudine astro
beta:=HMS→(Beta); // latitudine astro
ε:= epsilon(lista);
lha:=HMS→(Lha); // angolo orario
phi:=HMS→(lat); // lat geo
alpha:=atan2(SIN(lambda)*COS(ε)-TAN(beta)*SIN(ε),COS(lambda));
delta:=asin(sin(beta)*cos(ε)+cos(beta)*sin(ε)*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;
transformHorizontalToEquatorial(az, h)
BEGIN
LOCAL azimuth,altitude,phi,lha,delta, ang_delta;
azimuth:=HMS→(az);
altitude:=HMS→(h);
phi:=HMS→(lat);
lha:= atan2(sin(azimuth),cos(azimuth)*sin(phi)+tan(altitude)*cos(phi));
delta:= asin(sin(phi)*sin(altitude)-cos(phi)*cos(altitude)*cos(azimuth));
lha:=→HMS(lha);
delta:=→HMS(delta);
ang_delta:={lha,delta};
RETURN ang_delta;
END;
calcN()
BEGIN
N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT
RETURN N;
END;
calcTS(lista)
BEGIN
LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
JD0:= julianDayAtZero(Y,m,D);
dep:= deltaEpsilonPsi(lista);
psi:=dep(1);
eps:=dep(2);
T:= (JD-2451545.0)/36525;
θ0:= meanSideralTime(lista); // 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(lista); // apparent ST (ST at Greenwitch at local our)
TSL:= apparentSideralTime(lista) - 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(lista);
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(lista);
MST:= meanSideralTime(lista);
dep:= deltaEpsilonPsi(lista);
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(lista)
BEGIN
LOCAL JD, JDE, TSid, dep, psi, eps;
LOCAL ε, listatempi, ch;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
TSid:= calcTS(lista);
dep:= deltaEpsilonPsi(lista);
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(" ");
WAIT();
PRINT();
PRINT("Δψ Nutazione (longit.) "+ psi);
PRINT("Δε Nutazione (obl.) "+ eps);
PRINT("ε Obliquità Ecl. "+ ε);
PRINT(" ");
PRINT("Time " + lstd + " GMT " + gmt);
listatempi:= {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)) };
Tempi:= listatempi;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
transit(lista, listaARh0)
BEGIN
// lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
// α, δ (AR, dec), height, h0 = standard height of body
LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;
LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
LOCAL alfa1, alfa3, delta1, delta3, h0, alpha, del, n;
α:= listaARh0(2); // AR for day
δ:= listaARh0(5); // decl for day
alfa1:= listaARh0(1); // AR day-1
alfa3:= listaARh0(3); // AR day+1
delta1:= listaARh0(4); // decl day-1
delta3:= listaARh0(6); // decl day+1
h0:= listaARh0(7); // standard height
temp:= calcTS(lista);
θ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
n:= m0 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*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
n:= m1 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*15);
Hor:=simp_angle(→HMS(Hor));
h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
deltaM:= HMS→((h-h0)/(360*cos(del)*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
n:= m2 + deltaT/86400;
alpha := interpolation({alfa1, α, alfa3},n);
del := interpolation({delta1, δ, delta3},n);
Hor:= HMS→(tsideral - long - alpha*15);
Hor:=simp_angle(→HMS(Hor));
h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
deltaM:= HMS→((h-h0)/(360*cos(del)*cos(lat)*sin(Hor)));
tram:= ((m2 + deltaM) * 24); // setting
RETURN ({culm, sorg, tram});
END;
parametri(lista)
BEGIN
LOCAL dep, prec, α, δ, ε, listaparametri, ch, temp;
dep:= deltaEpsilonPsi(lista);
temp:= sunAlphaDelta(lista);
α:= sunlist(5); δ:= sunlist(6);
prec:= precession(α, δ);
ε:= epsilon(lista);
// 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));
listaparametri:= {ε, dep(1), dep(2), prec(1), prec(2)}; // ε, Δψ, Δε, Δα, Δδ
EpsilonPsi:= listaparametri;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
ascendant(lista)
BEGIN
LOCAL ASC, ts, segno, ε, listaascendente, ch, temp;
LOCAL tsb;
ε:= epsilon(lista);
ts:= calcTS(lista);
tsb:= ts(1)+lstd-dst; // ST of the birth
// temp:= sin(ts(2)*15)*cos(ε)+tan(lat)*sin(ε);
temp:= sin(tsb*15)*cos(ε)+tan(lat)*sin(ε);
ASC:= atan2(-cos(tsb*15), temp);
// IF temp<0 THEN ASC:= ASC+180 ELSE ASC:= ASC+360; END;
IF ASC<180 THEN ASC:=ASC+180 ELSE ASC:= ASC-180; END;
ASC:= simp_angle(ASC);
segno:= zodiaco(ASC);
listaascendente:= {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno};
Ascendant:= listaascendente;
RECT_P();
TEXTOUT_P("Ascendant " + trunc_angle(→HMS(ASC)), 25,50);
TEXTOUT_P(" " + trunc_angle(→HMS(ASC) MOD 30), 25,70);
TEXTOUT_P("in " + segno, 100,70);
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
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(lista)
// type dateshortlist
BEGIN
LOCAL ω, l, v, r, x, y, dz, ε;
LOCAL α, δ, dayY, ts, az, z, h;
LOCAL SAD, DA, Hd, lt;
LOCAL as, ac, zc, h0, H0, temp;
LOCAL ec, a, u, M1, M2, j, t, tau;
LOCAL JD, JDE, c, lapp, Ω, θ0;
LOCAL dep, deltapsi, longitude, b;
LOCAL lambdaprime, delta_long, delta_beta, correction;
ε:= epsilon(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
a:= 1.00000023; // semiasse maggiore
dep:= deltaEpsilonPsi(lista);
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.003032028*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
//M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-(T^3)/24490000); // anom vera Meeus
M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // Meeus p 308
Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
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
// 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(lista); // calcola tempo siderale
temp:= transformEquatorialToEcliptical(lista, α, δ); // β latitudine
b:= →HMS(temp(2));
lambdaprime:= l-1.397*T-0.00031*T^2;
delta_long=-(0.09033/3600);
delta_beta:=0.03916/3600*(COS(lambdaprime)-SIN(lambdaprime));
b:= b+delta_beta;
l:= l+delta_long;
correction:=-20.4898/r;
//lapp:= simp_angle(l - 0.00569-0.00478*sin(Ω)); // λ Longitudine apparente (ref true equinox)
lapp:=l+deltapsi+correction/3600;
lt:= simp_angle(l + 180); // Longitudine eliocentrica della Terra
// 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:= transformEquatorialToHorizontal(Hd, lat, δ);
az:= HMS→(temp(1))+180; // Azimuth
IF az>360 THEN az:= az-360 END;
az:= →HMS(az);
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
sunlist:= {L,ω, l, lapp, α,δ, az, h, M, Ω, v, c,x,y,r, H, zc, dz, DA, SAD, E, lt, ec, b};
RETURN sunlist;
END;
sun(lista)
BEGIN
LOCAL ω, l, b, v, r, x, y, dz;
LOCAL ε;
LOCAL α, δ, dayY, ts, az, z, h;
LOCAL SAD, DA, sorg, tram, culm, Hd, h0, lt;
LOCAL as, ac, zc, segno, JD, JDE, v, c;
LOCAL Ω, lapp, ec, listasun, listasuntran, listasunothers;
LOCAL ch, transito, listaARh0, temp, Hor, dep;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3;
ε:= epsilon(lista);
L:= sunlist(1); ω:= sunlist(2);
l:= sunlist(3); lapp:= sunlist(4);
α:= sunlist(5); δ:= sunlist(6);
az:= sunlist(7); h:= sunlist(8);
M:= sunlist(9); Ω:= sunlist(10);
v:= sunlist(11); c:= sunlist(12);
x:= sunlist(13); y:= sunlist(14);
r:= sunlist(15); Hor:= sunlist(16);
Hd:= 15*sunlist(16);
zc:= sunlist(17); dz:= sunlist(18);
DA:= sunlist(19); SAD:= sunlist(20);
E:= sunlist(21); lt:= sunlist(22);
ec:= sunlist(23);
b:= sunlist(24);
ts:= calcTS(lista);
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
dep:= deltaEpsilonPsi(lista);
segno:= zodiaco(l);
// 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
// to pass lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
h0:= -0°50′; // 'standard' altitude of Sun
temp:= sunAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= sunAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= sunAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
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("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("λ longitudine apparente " + trunc_angle(lapp));
PRINT("β latitudine " + →HMS(b)); // b never > 1.2" (Meeus)
PRINT("Δ 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(" ");
WAIT();
PRINT();
PRINT("Ω aberr, nutaz nodo "+ Ω);
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("Ang. orario "+ trunc_angle(Hor) + " " + 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(" ");
WAIT();
PRINT();
PRINT("Transita a: " + trunc_angle(culm MOD 24) + " 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("ε obliquità ecl. " + ε);
PRINT("ec eccentricità "+ ec);
PRINT("Δψ Nutazione long "+ dep(1));
PRINT("Δε Nutazione obl. " + dep(2));
PRINT(" ");
PRINT("Long elioc. Terra "+ trunc_angle(lt));
l:= →HMS(l); lapp:= →HMS(lapp);
culm:= →HMS(culm); sorg:= →HMS(sorg); tram:= →HMS(tram);
listasun:= {l, lapp, b, α, δ, az, h, r};
listasuntran:= {culm, sorg, tram};
listasunothers:= {L, ω, Ω, M,v, c, x, y, Hor, lt, ε, ec};
Sun_Effemerids:= listasun;
Sun_Transit:= listasuntran;
Sun_Others:= listasunothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
sunAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= sunCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;
sunAh(lista)
// lista like dateshortlist
BEGIN
LOCAL lha,alpha,delta,phi, ah, alfa;
LOCAL temp,longitude, Azimuth, Altitude, ts;
temp:= sunAlphaDelta(lista);
alpha:=15*HMS→(temp(1));
alfa:= →HMS(temp(1)); // α in HMS
delta:= temp(2);
phi:=lat;
ts:= calcTS(lista);
longitude:= HMS→(long); // - if long -E
lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
// lha:= 15*HMS→(ts(2)-alfa); // 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;
moonCalc(lista)
BEGIN
LOCAL ω,Ω, Lprime, Mprime, l, lsum, bsum, rsum;
LOCAL v, r, x, y, dep, ε;
LOCAL b, P, d, s, temp, tau;
LOCAL α,δ, ts, Hd, zc, dz, h, az;
LOCAL i, as, ac, DA, SAD, DD;
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(lista);
N:= calcN();
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
T:= (JDE-2451545.0)/36525; // use JDE
tau:= (JDE-2451545.0)/365250; // millenni - JDE
dep:= deltaEpsilonPsi(lista);
// Sun
L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // mean anomaly Sun, Meeus
// end Sun
// L1:= (33.231 + 13.17639653*N) MOD 360; // long. media
Lprime:= 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
Mprime:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media
ω:= (Lprime - Mprime) MOD 360; // perigeo medio
// DD:= Lprime - L; F:= Lprime-Ω; // 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(Mprime); // calc Moon longitude
lsum:= lsum + 1274027 * sin(2*DD-Mprime); // evezione
lsum:= lsum + 658314 * sin(2*DD); // variazione
lsum:= lsum + 213618 * sin(2*Mprime); // 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*Mprime);
lsum:= lsum + 57066 * sin(2*DD-Mprime-M) *ee;
lsum:= lsum + 53322 * sin(2*DD+Mprime);
lsum:= lsum + 45758 * sin(2*DD-M) *ee;
lsum:= lsum - 40923 * sin(M-Mprime) *ee;
lsum:= lsum - 34720 * sin(DD);
lsum:= lsum - 30383 * sin(Mprime+M) *ee;
lsum:= lsum + 15327 * sin(2*DD-2*F);
lsum:= lsum - 12528 * sin(Mprime+2*F);
lsum:= lsum + 10980 * sin(Mprime-2*F);
lsum:= lsum + 10675 * sin(4*DD-Mprime);
lsum:= lsum + 10034 * sin(3*Mprime);
lsum:= lsum + 8548 * sin(4*DD-2*Mprime);
lsum:= lsum - 7888 * sin(2*DD+M-Mprime) *ee;
lsum:= lsum - 6766 * sin(2*DD+M) *ee;
lsum:= lsum - 5163 * sin(DD-Mprime);
lsum:= lsum + 4987 * sin(DD+M) *ee;
lsum:= lsum + 4036 * sin(2*DD-M+Mprime) *ee;
lsum:= lsum + 3994 * sin(2*DD+2*Mprime);
lsum:= lsum + 3861 * sin(4*DD);
lsum:= lsum + 3665 * sin(2*DD-3*Mprime);
lsum:= lsum - 2689 * sin(M-2*Mprime) *ee;
lsum:= lsum - 2602 * sin(2*DD-Mprime+2*F);
lsum:= lsum + 2390 * sin(2*DD-M-2*Mprime) *ee;
lsum:= lsum - 2348 * sin(DD+Mprime);
lsum:= lsum + 2236 * sin(2*DD-2*M) *ee^2;
lsum:= lsum - 2120 * sin(M+2*Mprime) *ee;
lsum:= lsum - 2069 * sin(2*M) *ee^2;
lsum:= lsum + 2048 * sin(2*DD-2*M-Mprime) *ee^2;
lsum:= lsum - 1773 * sin(2*DD+Mprime-2*F);
lsum:= lsum - 1595 * sin(2*DD+2*F);
lsum:= lsum + 1215 * sin(4*DD-M-Mprime) *ee;
lsum:= lsum - 1110 * sin(2*Mprime+2*F);
lsum:= lsum - 892 * sin(3*DD-Mprime);
lsum:= lsum - 810 * sin(2*DD+M+Mprime) *ee;
lsum:= lsum + 759 * sin(4*DD-M-2*Mprime) *ee;
lsum:= lsum - 713 * sin(2*M-Mprime) * ee^2;
lsum:= lsum - 700 * sin(2*DD+2*M-Mprime) *ee^2;
lsum:= lsum + 691 * sin(2*DD+M-2*Mprime) *ee;
lsum:= lsum + 596 * sin(2*DD-M-2*F) *ee;
lsum:= lsum + 549 * sin(4*DD+Mprime);
lsum:= lsum + 537 * sin(4*Mprime);
lsum:= lsum + 520 * sin(4*DD-M) *ee;
lsum:= lsum - 487 * sin(DD-2*Mprime);
lsum:= lsum - 399 * sin(2*DD+M-2*F) *ee;
lsum:= lsum - 381 * sin(2*Mprime-2*F);
lsum:= lsum + 351 * sin(DD+M+Mprime) *ee;
lsum:= lsum - 340 * sin(3*DD-2*Mprime);
lsum:= lsum + 330 * sin(4*DD-3*Mprime);
lsum:= lsum + 327 * sin(2*DD-M+2*Mprime) *ee;
lsum:= lsum - 323 * sin(2*M+Mprime) *ee^2;
lsum:= lsum + 299 * sin(DD+M-Mprime) *ee;
lsum:= lsum + 294 * sin(2*DD+3*Mprime);
lsum:= lsum + 3958*sin(a1)+1962*sin(Lprime-F)+318*sin(a2);
l:= (simp_angle(Lprime + 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 + 280602 * sin(Mprime+F);
bsum:= bsum + 277693 * sin(Mprime-F); // equazione del centro
bsum:= bsum + 173237 * sin(2*DD-F); // grande ineguaglianza
bsum:= bsum + 55413 * sin(2*DD-Mprime+F); // evezione
bsum:= bsum + 46271 * sin(2*DD-Mprime-F); // evezione (2)
bsum:= bsum + 32573 * sin(2*DD+F);
bsum:= bsum + 17198 * sin(2*Mprime+F);
bsum:= bsum + 9266 * sin(2*DD+Mprime-F);
bsum:= bsum + 8822 * sin(2*Mprime-F);
bsum:= bsum + 8216 * sin(2*DD-M-F) *ee;
bsum:= bsum + 4324 * sin(2*DD-2*Mprime-F);
bsum:= bsum + 4200 * sin(2*DD+Mprime+F);
bsum:= bsum - 3359 * sin(2*DD+M-F) *ee;
bsum:= bsum + 2463 * sin(2*DD-M-Mprime+F) *ee;
bsum:= bsum + 2211 * sin(2*DD-M+F) *ee;
bsum:= bsum + 2065 * sin(2*DD-M-Mprime-F) *ee;
bsum:= bsum - 1870 * sin(M-Mprime-F) *ee;
bsum:= bsum + 1828 * sin(4*DD-Mprime-F);
bsum:= bsum - 1794 * sin(M+F) *ee;
bsum:= bsum - 1749 * sin(3*F);
bsum:= bsum - 1565 * sin(M-Mprime+F) *ee;
bsum:= bsum - 1491 * sin(DD+F);
bsum:= bsum - 1475 * sin(M+Mprime+F) *ee;
bsum:= bsum - 1410 * sin(M+Mprime-F) *ee;
bsum:= bsum - 1344 * sin(M-F) *ee;
bsum:= bsum - 1335 * sin(DD-F);
bsum:= bsum + 1107 * sin(3*Mprime+F);
bsum:= bsum + 1021 * sin(4*DD-F);
bsum:= bsum + 833 * sin(4*DD-Mprime+F);
bsum:= bsum + 777 * sin(Mprime-3*F);
bsum:= bsum + 671 * sin(4*DD-2*Mprime+F);
bsum:= bsum + 607 * sin(2*DD-3*F);
bsum:= bsum + 596 * sin(2*DD+2*Mprime-F);
bsum:= bsum + 491 * sin(2*DD-M+Mprime-F) *ee;
bsum:= bsum - 451 * sin(2*DD-2*Mprime+F);
bsum:= bsum + 439 * sin(3*Mprime-F);
bsum:= bsum + 422 * sin(2*DD+2*Mprime+F);
bsum:= bsum + 421 * sin(2*DD-3*Mprime-F);
bsum:= bsum - 366 * sin(2*DD+M-Mprime+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+Mprime+F) *ee;
bsum:= bsum + 302 * sin(2*DD-2*M-F) *ee^2;
bsum:= bsum - 283 * sin(Mprime+3*F);
bsum:= bsum - 229 * sin(2*DD+M+Mprime-F) *ee;
bsum:= bsum + 223 * sin(DD+M-F) *ee;
bsum:= bsum + 223 * sin(DD+M+F) *ee;
bsum:= bsum - 220 * sin(M-2*Mprime-F) *ee;
bsum:= bsum - 220 * sin(2*DD+M-Mprime-F) *ee;
bsum:= bsum - 185 * sin(DD+Mprime+F);
bsum:= bsum + 181 * sin(2*DD-M-2*Mprime-F) *ee;
bsum:= bsum - 177 * sin(M+2*Mprime+F) *ee;
bsum:= bsum + 176 * sin(4*DD-2*Mprime-F);
bsum:= bsum + 166 * sin(4*DD-M-Mprime-F) *ee;
bsum:= bsum - 164 * sin(DD+Mprime-F);
bsum:= bsum + 132 * sin(4*DD+Mprime-F);
bsum:= bsum - 119 * sin(DD-Mprime-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(Lprime))+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(Lprime-Mprime)-115*sin(Lprime+Mprime);
b := →HMS(bsum/1000000); // β div by 1 milion
// correzioni per la distanza dalla Terra (Δ)
rsum:=0;
rsum:= rsum -20905355 * cos(Mprime); // coseni
rsum:= rsum - 3699111 * cos(2*DD-Mprime); // evezione
rsum:= rsum - 2955968 * cos(2*DD); // variazione
rsum:= rsum - 569925 * cos(2*Mprime); // 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*Mprime);
rsum:= rsum - 152138 * cos(2*DD-Mprime-M) *ee;
rsum:= rsum - 170733 * cos(2*DD+Mprime);
rsum:= rsum - 204586 * cos(2*DD-M) *ee;
rsum:= rsum - 129620 * cos(M-Mprime) *ee;
rsum:= rsum + 108743 * cos(DD);
rsum:= rsum + 104755 * cos(Mprime+M) *ee;
rsum:= rsum + 10321 * cos(2*DD-2*F);
rsum:= rsum + 79661 * cos(Mprime-2*F);
rsum:= rsum - 34782 * cos(4*DD-Mprime);
rsum:= rsum - 23210 * cos(3*Mprime);
rsum:= rsum - 21636 * cos(4*DD-2*Mprime);
rsum:= rsum + 24208 * cos(2*DD+M-Mprime) *ee;
rsum:= rsum + 30824 * cos(2*DD+M) *ee;
rsum:= rsum - 8379 * cos(DD-Mprime);
rsum:= rsum - 16675 * cos(DD+M) *ee;
rsum:= rsum - 12831 * cos(2*DD-M+Mprime) *ee;
rsum:= rsum - 10445 * cos(2*DD+2*Mprime);
rsum:= rsum - 11650 * cos(4*DD);
rsum:= rsum + 14403 * cos(2*DD-3*Mprime);
rsum:= rsum - 7003 * cos(M-2*Mprime) *ee;
rsum:= rsum + 10056 * cos(2*DD-M-2*Mprime) * ee;
rsum:= rsum + 6322 * cos(DD+Mprime);
rsum:= rsum - 9884 * cos(2*DD-2*M) * ee^2;
rsum:= rsum + 5751 * cos(M+2*Mprime) *ee;
rsum:= rsum - 4950 * cos(2*DD-2*M-Mprime) *ee^2;
rsum:= rsum + 4130 * cos(2*DD+Mprime-2*F);
rsum:= rsum - 3958 * cos(4*DD-M-Mprime) *ee;
rsum:= rsum + 3258 * cos(3*DD-Mprime);
rsum:= rsum + 2616 * cos(2*DD+M+Mprime) *ee;
rsum:= rsum - 1897 * cos(4*DD-M-2*Mprime) *ee;
rsum:= rsum - 2117 * cos(2*M-Mprime) * ee^2;
rsum:= rsum + 2354 * cos(2*DD+2*M-Mprime) *ee^2;
rsum:= rsum - 1423 * cos(4*DD+Mprime);
rsum:= rsum - 1117 * cos(4*Mprime);
rsum:= rsum - 1571 * cos(4*DD-M) *ee;
rsum:= rsum - 1739 * cos(DD-2*Mprime);
rsum:= rsum - 4421 * cos(2*Mprime-2*F);
rsum:= rsum + 1165 * cos(2*M+Mprime) *ee^2;
rsum:= rsum + 8752 * cos(2*DD-Mprime-2*F);
rsum:= rsum/1000; // div by 1000
temp:= transformEclipticalToEquatorial(lista, l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
ts:= calcTS(lista); // calcola tempo siderale
Hd:=15*HMS→(ts(2)-α); // local ST vs apparent ST
Hd:=simp_angle(→HMS(Hd));
H:= Hd/15; // in HMS
DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
SAD:= (90 + DA); // Semiarco Diurno
// SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
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)
// s:= 358473400/d; // semidiametro (in " arcsec)
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 Calante"; END;
IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
IF ph>100 AND ph<170 THEN fase:="Gibbosa Crescente"; 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
moonlist:= {Lprime, l, b, lapp, α, δ, az, h, Mprime, ω, Ω, d, DD, P, s, H, DA, SAD,ph, fase, eta, ill};
RETURN moonlist;
END;
moon(lista)
BEGIN
LOCAL ε, ts, JD, JDE, segno, Lprime, b;
LOCAL l, Ω, ω, Mprime, lapp, α, δ;
LOCAL az, h, sorg, tram, culm, Hd, DA, SAD;
LOCAL ph, fase, eta, DD, d, s, ill;
LOCAL dep, temp, listamoon, listamoontran, listamoonothers, ch;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
LOCAL listaARh0, transito, Hor;
ε:= epsilon(lista);
ts:= calcTS(lista); //
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
dep:= deltaEpsilonPsi(lista);
Lprime:= moonlist(1); l:= moonlist(2);
b:= moonlist(3);
α:= moonlist(5); δ:= moonlist(6);
az:= moonlist(7); h:= moonlist(8);
ω:= moonlist(10); Ω:= moonlist(11);
lapp:= moonlist(4); Mprime:= moonlist(9);
d:= moonlist(12); DD:= moonlist(13);
P:= moonlist(14); s:= moonlist(15);
Hor:= moonlist(16); Hd:= 15*moonlist(16);
DA:= moonlist(17); SAD:= moonlist(18);
ph:= moonlist(19); fase:= moonlist(20);
eta:= moonlist(21); ill:= moonlist(22);
segno:= zodiaco(l);
// transit routine
h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)
// h0:= 0°7′30″; // medium value: 0.125°
temp:= moonAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= moonAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= moonAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
culm:= transito(1);
sorg:= transito(2);
tram:= transito(3);
// end transit
PRINT;
PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
PRINT("L Long media " + trunc_angle(Lprime));
PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("β Latitudine " + trunc_angle(b));
PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
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("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");
PRINT("Semidiametro "+ trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,3) + " km");
PRINT(" ");
WAIT();
PRINT();
PRINT("ω Perigeo medio " + trunc_angle(ω));
PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
PRINT("M anomalia vera " + trunc_angle(Mprime));
PRINT("D Elongazione " + trunc_angle(DD));
PRINT(" ");
PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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("Julian day of Ephemeris " + JDE);
PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
PRINT(" ");
WAIT();
PRINT();
PRINT("Transita a: " + trunc_angle(culm MOD 24 ) + " 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("ε obliquità ecl. " + ε);
PRINT("Δψ Nutazione long "+ dep(1));
PRINT("Δε Nutazione obl. " + dep(2));
PRINT(" ");
PRINT("π Parallasse " + trunc_angle(P));
PRINT("Fase " + trunc_angle(ph) + " " + fase);
PRINT("Età " + eta);
PRINT("Parte illuminata " + ROUND(ill,5) + " " + ROUND(100*ill,2)+"%");
listamoon:= {→HMS(l), →HMS(b), ROUND(HMS→(d),6), →HMS(α), →HMS(δ), →HMS(az), →HMS(h)};
listamoonothers:= {→HMS(Lprime), →HMS(lapp), →HMS(Mprime),→HMS(ω),→HMS(Ω), →HMS(P),
→HMS(s), →HMS(Hor), →HMS(α*15), ph, eta};
listamoontran:= →HMS({culm, sorg, tram});
Moon_Effemerids:= listamoon;
Moon_Transit:= listamoontran;
Moon_Others:= listamoonothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
moonAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= moonCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;
moonAh(lista)
// lista like dateshortlist
BEGIN
LOCAL lha,alpha,delta,phi, ah, alfa;
LOCAL temp,longitude, Azimuth, Altitude, ts;
temp:= sunAlphaDelta(lista);
alpha:=15*HMS→(temp(1));
alfa:= →HMS(temp(1)); // α in HMS
delta:= temp(2);
phi:=lat;
ts:= calcTS(lista);
longitude:= - HMS→(long); // - if long -E
lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
// lha:= 15*HMS→(ts(2)-alfa); // 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;
planetCalc(nameplanet, lista)
BEGIN
LOCAL ω,Ω, l, v, r, x, y, z;
LOCAL Ls, Ms, b, xs, ys, b_ec;
LOCAL x1, y1, z1, l_ec, 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, lapp, dep;
LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
JD:= julianDay(lista);
ε:= epsilon(lista);
dep:= deltaEpsilonPsi(lista);
// calc of Sun (as basis)
sunCalc(lista);
xs:= sunlist(13); ys:= sunlist(14);
Ls:= sunlist(1); // Longitudine media
Ms:= sunlist(9); // 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 eclittic heliocentric
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 (dist Sun)
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); // distanza Terra
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
lapp:= →HMS(l + dep(1)); // longitudine apparente (+ nutazione long, Δψ)
temp:= transformEclipticalToEquatorial(lista, l,b); // transform into AR, decl
α:= temp(1);
δ:= temp(2);
ts:= calcTS(lista); // 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
planetlist:= {L, l, b, lapp, α,δ, az, h, M, Ω, ω, v, x,y, z, x1, y1, z1,
r, R, H, zc, dz, DA, SAD, ec, in, a, θ, l_ec, b_ec, Ls, Ms, xs, ys, nameplanet};
RETURN planetlist;
END;
planetDisplay(nameplanet, lista)
BEGIN
LOCAL ω,Ω, l, v, r;
LOCAL x, y, z, x1, y1, z1;
LOCAL Ls, Ms, b, xs, ys, l_ec, b_ec;
LOCAL ec, in, a, δ, α, ts, az, h;
LOCAL DA, SAD, tram, sorg, culm, dayY;
LOCAL as, ac, zc, dz, θ, mov, segno, ε;
LOCAL transito, JD, JDE, temp, Hor, Hd;
LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
LOCAL ee, magni, m0, aa, bb,cc, listamagn, ill;
LOCAL d0,dapp;
LOCAL listaARh0, lapp;
ts:= calcTS(lista);
JD:= julianDay(lista);
JDE:= julianEphemerisDay(lista);
L:= planetlist(1); l:= planetlist(2);
b:= planetlist(3); lapp:= planetlist(4);
α:= planetlist(5); δ:= planetlist(6);
az:= planetlist(7); h:= planetlist(8);
M:= planetlist(9); Ω:= planetlist(10);
ω:= planetlist(11); v:= planetlist(12);
x := planetlist(13); y:= planetlist(14); z:= planetlist(15);
x1:= planetlist(16); y1:= planetlist(17); z1:= planetlist(18);
r:= planetlist(19); R:= planetlist(20);
Hor:= planetlist(21); Hd:= 15*planetlist(21);
zc:= planetlist(22); dz:= planetlist(23);
DA:= planetlist(24); SAD:= planetlist(25);
ec:= planetlist(26); in:= planetlist(27);
a:= planetlist(28); θ:= planetlist(29);
l_ec:= planetlist(30); b_ec:= planetlist(31);
Ls:= planetlist(32); Ms:= planetlist(33); // Sun values
xs:= planetlist(34); ys:= planetlist(35);
// nameplanet:= planetlist(36);
segno:= zodiaco(l);
ee:= abs(l - l_ec); // E (phase) to calc phase and magnitude
ill:= 0.5*(1+cos(ee));
CASE
IF nameplanet=="Mercurius" THEN listamagn:= {6.7, -0.21,3.8,-3.4,2}; END;
IF nameplanet=="Venus" THEN listamagn:= {16.7, -0.21, 3.8, -3.4, 2};END;
IF nameplanet=="Mars" THEN listamagn:= {9.4, -1.36, 1.5,0, 0}; END;
IF nameplanet=="Jupiter" THEN listamagn:= {196.9, -9, 1.48, 0,0}; END;
IF nameplanet=="Saturn" THEN listamagn:= {165.5, -8.7, 1.7, 0,0}; END;
IF nameplanet== "Uranus" THEN listamagn:= {71.4, -7, 0,0,0}; END;
IF nameplanet=="Neptunus" THEN listamagn:= {68.3, -7, 0,0,0}; END;
IF nameplanet=="Pluto" THEN listamagn:= {8.3, -1.5, 0, 0, 0}; END;
DEFAULT
END; // case
d0:= listamagn(1);
m0:= listamagn(2);
aa:= listamagn(3); bb:= listamagn(4); cc:= listamagn(5);
dapp:= d0/R; // apparent diameter
magni:= m0 + 5*LOG(r*R)+aa*ee/100+bb*(ee/100)^2+cc*(ee/100)^3;
P:= →HMS(HMS→(0°0′8.794″)/R); // Parallasse
// transit routine
h0:= -0°34′; // 'standard' altitude of a planet
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D-1,0,0,0}));
alfa1:= temp(1); delta1:= temp(2);
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D,0,0,0}));
alfa2:= temp(1); delta2:= temp(2);
temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D+1,0,0,0}));
alfa3:= temp(1); delta3:= temp(2);
listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
transito:= transit(lista, listaARh0);
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(" ");
PRINT("L Long media " + trunc_angle(L));
PRINT("λ Longitudine " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
PRINT("β Latitudine " + trunc_angle(b));
PRINT("Δ Distanza " + ROUND(R,4) + " (UA) from Earth ");
PRINT("α Asc. retta "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
PRINT(" ");
PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));
PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(simp_angle(v)));
PRINT("r Raggio vettore " + ROUND(r,4) + " (UA) from Sun "); // distanza eliocentrica
WAIT();
PRINT();
PRINT("ω Perielio " + trunc_angle(ω));
PRINT("Ω Nodo " + trunc_angle(Ω));
PRINT("M Anomalia vera " + trunc_angle(M));
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(" ");
PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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);
WAIT();
PRINT();
PRINT("Diff. ascensionale " + trunc_angle(DA));
PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));
PRINT(" ");
PRINT("π Parallasse " + ROUND(P, 5));
PRINT("Diametro apparente "+ ROUND(dapp, 5));
PRINT("Magnitudo " + ROUND(magni, 5));
PRINT("Frazione illuninata " + ROUND(ill*100, 2)+"%");
PRINT("Fase " + →HMS(ee));
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);
listaplanet:= {l,b, R, α, δ, az, h};
listaplanettran:= {culm, sorg, tram};
listaplanetothers:= {L, lapp, l_ec,b_ec, M, ω,Ω, v, x,y,z, x1,y1,z1, r, Hor, Hd, SAD, DA, mov, dapp, magni};
Planet_Effemerids:= listaplanet;
Planet_Transit:= listaplanettran;
Planet_Others:= listaplanetothers;
FREEZE();
WAIT();
ch:= MSGBOX("Another calculation?", 1);
IF (ch<>0) THEN changeData(); END;
END;
planetAlphaDelta(nameplanet, lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= planetCalc(nameplanet, lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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
N:= 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: Effemeridi (Ephemeris) - salvomic - 07-02-2015 11:10 PM
New version:
corrected some bugs, added Δα and Δδ correction for daytime parallax, exported a function to calc precession (given RA -DEG- and declination), "graphic" restyling of printed values...
Enjoy!
EDIT: corrected a but for age of Moon, added "epatta" for Moon (old calculation for Easter date...)
RE: Astronomy: Effemeridi (Ephemeris) - Thomas_Sch - 07-03-2015 06:27 AM
(07-01-2015 09:48 PM)salvomic Wrote: ...in the meantime this version introduce other values for planets: parallax, illuminated fraction, magnitudo, phase.
EDIT: Thomas, try this with an option to save coordinates at the starting, the program should save them in Long_Lat var in the program variables and load them when used again...
I need a better safeguard, however, then only control if the size is ...2 ;-)
Hello Salvo,
many thanks for your message & your efforts! I'll try it as soon as possible.
I've seen your other questions about checking the list, and the discussion about type(), very interesting! Maybe the special characters could be checked? Don't know... Range of values? (-90 >= latitude <= +90; -180 >= longitude <= +180)?
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-03-2015 06:53 AM
(07-03-2015 06:27 AM)Thomas_Sch Wrote: Hello Salvo,
many thanks for your message & your efforts! I'll try it as soon as possible.
I've seen your other questions about checking the list, and the discussion about type(), very interesting! Maybe the special characters could be checked? Don't know... Range of values? (-90 >= latitude <= +90; -180 >= longitude <= +180)?
I'm investigating, Thomas,
for now is enough to check for size(): if there are two (or more elements is a list...); why the user should change that variable? ;-)
Do you think is important save only coordinates or also time zone (and / or DST)?
Time and date I think no, better to leave it as is: the user can change them easily...
Salvo
RE: Astronomy: Effemeridi (Ephemeris) - Thomas_Sch - 07-03-2015 07:05 AM
(07-03-2015 06:53 AM)salvomic Wrote: I'm investigating, Thomas,
for now is enough to check for size(): if there are two (or more elements is a list...); why the user should change that variable? ;-)
Do you think is important save only coordinates or also time zone (and / or DST)?
Time and date I think no, better to leave it as is: the user can change them easily...
Salvo
You are right, time zone should also be saved, (also there is a rough correlation between longitude and time zone).
Maybe an option or checkbox ("Save?" yes / no ) could avoid uncertainties ..
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-03-2015 02:31 PM
(07-03-2015 07:05 AM)Thomas_Sch Wrote: You are right, time zone should also be saved, (also there is a rough correlation between longitude and time zone).
Maybe an option or checkbox ("Save?" yes / no ) could avoid uncertainties ..
ok, now the program saves long, lat, time zone, DST and deltaT (like the app by Marcel, in the same way, so we can use both program and app together).
Exported also some useful function (Ascendant, making lists, Julian Day, epsilon, precession, epatta,Sun and Moon alpha/delta, Sun and Moon Ah...)
Salvo
File's date: July 5, 2015
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-06-2015 03:01 PM
New version:
rewritten some parts: now calculations for Sun, Moon, Planets are a bit slow, but more precise (rewritten the routine for nutation, deltapsi, deltaepsilon, more precise but require a loop of 63 linear combination of parameters, used also in others part of program, for precision)...
More precise transit routine for Sun and Planets and finally also the Moon (thanks Marcel!).
Reviewed Ascendant calculation (more precise) with indication of Descendant...
More precise latitude for Sun.
Corrected a but in Moon phase.
Salvo
Date last version: 2015, July, 8
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-09-2015 12:30 PM
New version:
further improved Sun precision (about 100 parameters, precision about 0.01").
Exported functions Sun_Coordinates, Moon_Coordinates, Planets_Coordinates.
They gives however only few important coordinates (lambda, beta, alfa, delta, R, parallax); to get more use Effemeridi that gives (for Sun, Moon, Planets) other data: anomaly, node, perigee/apogee, sidereal time, azimuth, height, cartesian coordinates, age of Moon..., and transit (exported also as independent function).
Date: 2015, July, 9
RE: Astronomy: Effemeridi (Ephemeris) - Marcel - 07-09-2015 01:50 PM
Hi Salvomic!
All that DATA (anomaly, node,...) is a great addition to my program.
Thank you.
Marcel
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-09-2015 04:17 PM
(07-09-2015 01:50 PM)Marcel Wrote: Hi Salvomic!
All that DATA (anomaly, node,...) is a great addition to my program.
Thank you.
Marcel
thank you Marcel,
I think that the best experience for our users is to use both app and program: they give the same data for the principal aspects (Sun, Moon, Planets, time, most important parameters), then different approach (graphic and detailed). They are becoming "complementary"...
And we must always thank J. Meeus for the great input and his precise calculations!
Eventually, the Prime is a great machine to do Astronomy tools, and our programs demonstrate that.
Salvo
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-10-2015 08:15 PM
hi,
in this new version I added a function dedicated to the Earth (ellipsoid).
Earth(datelist);
With it the user can input a latitude (or use that already saved with the program) -and (or not) heigh above the sea level- and calculate some geodetic parameters: geocentric latitude, flattening, equatorial and meridian radius, eccentricity of the meridian, rho•sin(phi) and rho•cos(phi), difference between geographic and geocentric latitude, radius of parallel of latitude, rotational angular velocity, radius of curvature of the meridian...
Giving lon/lat for two points the program returns the angular and linear distance (with an error or ±50 meters!).
Then the program output actual heliocentric longitude and latitude for the Earth.
Examples:
• to input φ, h (height above sea level): ρ, ρ·sin(φ'), ρ·sin(φ') :: useful data to get diurnal parallax, eclipses, occultations)
• to input φ: Rp, Rm, λ°, φ° (radius of parallel of latitude, radius of curvature, degree of longitude or latitude length
• to input φ1, φ2: s, d (distance from two points given their coordinates: both linear and angular)
Enjoy!
Salvo M.
Date: 2015, July, 13
Version: 2.1
RE: Astronomy: Effemeridi (Ephemeris) - Marcel - 07-12-2015 12:38 PM
Hi Salvomic!
Great to see that this chapter on the Earth globe is in your program!
The book is covered from the start to the end..
If the probe New Horizons interest you, I sugest to you the download of the program "NASA's Eyes". This is a beta application from the JPL Lab. If you go with the probe New Horizons, you will have a "Live" option and a "Preview" option. This "Preview" option will show you what the probe will do tomorrow and the day after. Is you look at the right of the screen, you will see what is camera will do... This is really a great animation... I think Astro Lab 8 or may be 9 will show this view!
Marcel
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-12-2015 08:57 PM
(07-12-2015 12:38 PM)Marcel Wrote: Hi Salvomic!
Great to see that this chapter on the Earth globe is in your program!
The book is covered from the start to the end..
yes, Marcel, the book is covered all! and it deserve: it's a great book!
Quote:If the probe New Horizons interest you, I sugest to you the download of the program "NASA's Eyes". This is a beta application from the JPL Lab. If you go with the probe New Horizons, you will have a "Live" option and a "Preview" option. This "Preview" option will show you what the probe will do tomorrow and the day after. Is you look at the right of the screen, you will see what is camera will do... This is really a great animation...
yes, I'm interested enough, that's a historic event with Pluto, we'll know many new things about this fascinating planet.
I'll download the program tomorrow.
Quote:I think Astro Lab 8 or may be 9 will show this view!
Marcel
sure! :-)
Thanks always for your great app!
Salvo
RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-14-2015 10:17 AM
New version:
implemented a new method to calculate planets coordinates (for more precise result -max ±30', mostly ±2'- use the correspondent functions in Astro Lab 4 that uses another precise method to calculates values);
exported a function to calculate angular separation for two bodies, giving right ascension and declination for both (d = angular separation from the Earth, p = relative separation);
exported also the function to interpolate from three values, giving the "n" factor, useful to get also angular separation (and many other purpose) from a table of RA and dec;
exported Anomaly (mean, eccentric, true) and Moon elongation;
translated other parts;
corrected printing of perihelion and added "argument of perihelion";
corrected (and exported) the tool to get "eccentric anomaly" via Kepler equation...
I must also write a guide: soon.
Enjoy!
Salvo
Version: 2.3
Date: 2015, July, 18
|