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