Astronomy: Effemeridi (Ephemeris)
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;
***
∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
|
Search
Member List
Calendar
Help
