Astronomy...

[salvomic]

new version, better handling of Sun and Moon RA, dec and Az,h...
Now should work almost well also the routine for transit, rising and setting of Sun, Moon, planets.

When the program will be complete, I'll put the actual version in the Software Library of the Forum.

Salvo

Code:


data();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
tempo();
precession();
parametri();
transit();
sunCalc();
planetCalc();
planets();
sun();
sunAh();
moon();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendent();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();
thisday:= 2000.0101;
lstd:=0; gmt:= 12;
long:= 0; lat:=  0;
m:= 1; dst:= 0; utc:=0;
deltaT:= 68; // 2015
datetimelist:= MAKELIST(X,X,1,6);
dateshortlist:= MAKELIST(X,X,1,3);
dateshortlistgmt:= MAKELIST(X,X,1,3);
sunlist:= MAKELIST(X,X,1,25);
planet_const:= MAKELIST(X,X,1,10);
// ε:= epsilon();  // approx 23.45 (23°26'21")

EXPORT Effemeridi()
// from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015
// Credits: Serge Bouiges, Calcule astronomique pour amateurs
// Jean Meeus, Astronomical Algorithms
// Thanks Didier Lachieze, Marcel Pelletier, Eried
BEGIN
  LOCAL ch, ε;
  HAngle:=1;
  data(); // input data and calc Sun
  ε:=epsilon();
  CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendent", "Quit");
  CASE
  IF ch==1 THEN sunCalc(); sun(); END;
  IF ch==2 THEN moon(); END;
  IF ch==3 THEN planets(); END;
  IF ch==4 THEN tempo(); END;
  IF ch==5 THEN parametri(); END;
  IF ch==6 THEN ascendent(); END;
  IF ch==7 THEN HAngle:=0; RETURN; END;
  DEFAULT sunCalc(); sun();
  END; // case
END;

data()
BEGIN
  LOCAL dd, mm, yy, hh, mn, ss, loct;
  LOCAL t, JD, JDE, ε, dep, eps;
  yy:= IP(Date);
  mm:= IP(FP(Date)*100);
  dd:= FP(Date*100)*100;
  hh:= IP(HMS→(Time));
  mn:= IP(FP(HMS→(Time))*60);
  ss:= IP(FP(FP(FP(HMS→(Time))*60))*60);
  loct:= →HMS(hh+mn/60+ss/3600);
  INPUT({ {m,[0],{15,15,1}}, {D,[0],{50,15,1}}, {Y,[0],{80,15,1}},
  {lstd,[0],{65,30,2}},
  {dst,0,{40,15,3}}, {utc,[0],{80,15,3}}, 
  {long,[0],{20,25,4}}, 
  {lat,[0],{70,25,4}},
  {deltaT, [0], {20,25,5}} },
  "Data: Use Shift+9 for H°M′S″",
  {"Month:","Day :","Year :","Local Time (24):", "DST", "UTC", 
  "Long (-E):","Lat (+N)", "delta T"},
  {"Enter month", "Enter day", "Enter year (1901-2099)", 
  "Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",
  "Longitude", "Latitude", "Delta T (sec.)"}, 
  {1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});
  // Ragusa (-14°43′41″, 36°55′15″) - 
  // Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)

  gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)
  thisday:= Y+m/100+D/10000;
  datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };
  dateshortlist:= {Y, m, D+HMS→(lstd)/24};
  dateshortlistgmt:= {Y, m, D+HMS→(gmt)/24};
  // Greenwich Mean Time

  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:=(JD-2451545)/36525;
  U:= T/100; 
  ε:=epsilon(); // (formula IAU J2000.0)
  sunCalc();
END;

epsilon()
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:=(JD-2451545)/36525;
  U:= T/100;
  // ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
  ε0:= 23°26′21.448″ - 1°18′00.93″*U-1.55*U^2+1999.25*U^3-51.38*U^4-249.67*U^5;
  ε0:= ε0-39.05*U^6+7.12*U^7+27.87*U^8+5.79*U^9+2.45*U^10; // Laskar, Meeus
  ε0:= →HMS(ε0);
  dep:= deltaEpsilonPsi();
  eps:= dep(2);
  ε:= ε0+eps; // true obliquity ε
  RETURN ε;
END;

deltaEpsilonPsi()
BEGIN
  // Nutazione Δψ (longitudine), Δε (obliquità)
  LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
  LOCAL DD;
  JD:= julianDayAtZero(Y,m,D);
  JDE:= julianEphemerisDay(dateshortlist);
  T:= (JDE-2451545)/36525; // use JDE
  DD:=297.85036+445267.111480*T-0.0019142*T^2+(T^3)/189474;  // elongazione
  DD:=FP(DD/360)*360;
  Ωmoon:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)
  Lsun:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Sole
  Lmoon:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000); 
  //  Long media Luna
  deltaPsi:= (-0°0′17.20″)*sin(Ωmoon)-(0°0′1.32″)*sin(2*Lsun)+(-0°0′0.23″)*sin(2*Lmoon)+(0°0′0.21″)*sin(2*Ωmoon);
  deltaEpsilon:= (0°0′9.2″)*cos(Ωmoon)+(0°0′0.57″)*cos(2*Lsun)+(0°0′0.10″)*cos(2*Lmoon)-(0°0′0.09″)*cos(2*Ωmoon);
  RETURN({deltaPsi, deltaEpsilon});   
END;

precession(alfa, delta)
BEGIN
  // precessione dati asce retta e declinazione
  LOCAL m,n,n1, deltaAlfa, deltaDelta, Tsec;
  Tsec:= IP((Y-2000)/100)+1; // T century from 2000.0
  m:= 0°0′3.07496″+0°0′0.00186″*Tsec;
  n:= 0°0′1.33621″+0°0′0.00057″*Tsec; // sec
  n1:= 0°0′20.0431″+0°0′0.0085″*Tsec; // ° ' "
  deltaAlfa:= m+n*sin(alfa)*tan(delta); // Δα
  deltaDelta:= n1*cos(alfa); // Δδ
  RETURN({deltaAlfa, deltaDelta});
END;

planets()
BEGIN
  LOCAL ch, nameplanet;
  INPUT({{ch, nameplanet:={"Mercurius","Venus","Mars", "Jupiter", 
  "Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},
  "Choose planet", 
  "planet: ", "Choose the planet to calculate an press OK", 0,5 );
  planetCalc(nameplanet(ch));
END;

transformEclipticalToEquatorial(lambda, beta)
BEGIN
  // trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
  LOCAL ε,Epsilon,Lambda, Beta, alpha, delta;
  Lambda:=HMS→(lambda);
  Beta:=HMS→(beta);
  Epsilon:= epsilon();
  ε:= epsilon();
  alpha:=atan2(sin(Lambda)*cos(Epsilon)-tan(Beta)*sin(Epsilon),cos(Lambda));
  IF alpha>=360 THEN alpha:=alpha-360 END;
  IF alpha<0 THEN alpha:=alpha+360 END;
  delta:=asin(sin(Beta)*cos(Epsilon)+cos(Beta)*sin(Epsilon)*sin(Lambda));
  alpha:=→HMS(alpha/15);
  delta:=→HMS(delta);
  RETURN ({alpha,delta});
END;

transformEquatorialToHorizontal(Lha, Phi, Delta)
BEGIN
  LOCAL azimuth,altitude,lha,phi,delta;
  lha:=HMS→(Lha);
  phi:=HMS→(Phi);   // latitude
  delta:=HMS→(Delta);
  azimuth:=atan2(sin(lha),cos(lha)*sin(phi)-tan(delta)*cos(phi));
  altitude:=asin(sin(phi)*sin(delta)+cos(phi)*cos(delta)*cos(lha));
  azimuth:=→HMS(azimuth);
  altitude:=→HMS(altitude);
  RETURN({azimuth,altitude});
END;

transformEclipticalToHorizontal(Lambda, Beta, Lha)
BEGIN
  LOCAL ε;
  // transform λ and β (long, lat) into azimuth and height
  LOCAL lambda, beta, epsilon,phi, lha;
  LOCAL alpha,delta,azimuth,altitude;
  lambda:=HMS→(Lambda); // longitudine astro
  beta:=HMS→(Beta); // latitudine astro
  epsilon:= epsilon();
  ε:= epsilon();
  lha:=HMS→(Lha);  // angolo orario
  phi:=HMS→(lat); // lat geo
  alpha:=atan2(SIN(lambda)*COS(epsilon)-TAN(beta)*SIN(epsilon),COS(lambda));
  delta:=asin(sin(beta)*cos(epsilon)+cos(beta)*sin(epsilon)*sin(lambda));
  azimuth:=atan2(SIN(lha),COS(lha)*SIN(phi)-TAN(delta)*COS(phi));
  altitude:=ASIN(SIN(phi)*SIN(delta)+COS(phi)*COS(delta)*COS(lha));
  azimuth:=→HMS(azimuth);
  altitude:=→HMS(altitude);
  RETURN ({azimuth,altitude});
END;

calcN()
BEGIN
  N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT
  RETURN N;
END;

calcTS()
BEGIN
  LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
  LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
  ε:= epsilon;
  calcN();
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  JD0:= julianDayAtZero(Y,m,D);
  dep:= deltaEpsilonPsi();
  psi:=dep(1);
  eps:=dep(2);
  T:= (JD-2451545.0)/36525;
  θ0:= meanSideralTime(dateshortlist); // MST DEG at 0h
  θ0GMT:= meanSideralTime({Y,m,D}); // MST DEG at 0h at Greenwitch
  // TS:=  ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds
  // TS:= →HMS(TS); // in HMS
  TS:= θ0GMT;
  // TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST
  // TSL:= →HMS(TSL); // in HMS
  TSapp:= apparentSideralTime(dateshortlist); // apparent ST (ST at Greenwitch at local our)
  TSL:= apparentSideralTime(dateshortlist) - 1*dst;
  TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees

  RETURN {TS, TSL, TSapp, θ0, T};
END;

meanSideralTime(lista)
// lista like dateshortlist
BEGIN
  LOCAL T,MST, JD, ε;
  ε:= epsilon;
  JD:= julianDay(lista);
  T:=(julianDay(lista)-2451545.0)/36525;
  MST:=280.46061837+360.98564736629*(JD-2451545.0)+0.000387933*T^2-(T^3)/38710000;
  MST:=FP(MST/360)*360; // θ0
  IF MST<0 THEN MST:=MST+360 END;
  MST:=→HMS(MST/15);
  RETURN MST;
END;

apparentSideralTime(lista)
// lista like dateshortlist
BEGIN
  LOCAL MST,correction, deltaPsi, AST, dep, ε;
  ε:= epsilon;
  MST:= meanSideralTime(lista);
  dep:= deltaEpsilonPsi();
  deltaPsi:= dep(1);
  correction:=(HMS→(deltaPsi)*3600*COS(HMS→(ε)))/15;
  AST:=→HMS(HMS→(MST)+(correction/3600));
  RETURN AST;
END;

julianDay(lista)
// lista like dateshortlist
BEGIN
  LOCAL y,m,d,a,b, JD;
  y:=lista(1);
  m:=lista(2);
  d:=lista(3);
  IF (m=1 OR m=2) THEN
  y:=y-1;
  m:=m+12
  END;
  IF y*100+m+d>=158225 THEN
  a:=IP(ABS(y/100));
  b:=2-a+IP(a/4);
  ELSE
  b:=0
  END;
  JD:=IP(365.25*(y+4716))+IP(30.6001*(m+1))+d+b-1524.5;
  RETURN JD;
END;

julianEphemerisDay(lista)
BEGIN
  LOCAL JDE;
  JDE:=julianDay(makeDateListShort(makeDateList(lista)+{0,0,0,0,0,deltaT}));
  RETURN JDE;
END;

julianDayAtZero(Y,m,D)
BEGIN
  // JD at 0h GMT
  LOCAL y,mm,d,a,b, JD;
  y:=Y;
  mm:=m;
  d:=dateshortlist(3);
  IF (mm=1 OR mm=2) THEN
  y:=y-1;
  mm:=mm+12
  END;
  IF y*100+mm+d>=158225 THEN
  a:=IP(ABS(y/100));
  b:=2-a+IP(a/4);
  ELSE
  b:=0
  END;
  JD:=IP(365.25*(y+4716))+IP(30.6001*(mm+1))+d+b-1524.5;
  RETURN JD;
END;

makeDateList(lista)
// lista {Y m D.d} (dateshortlist)
BEGIN
LOCAL y,o,m,d,h,s,t;
y:=lista(1);
o:=lista(2);
d:=IP(lista(3));
t:=FP(lista(3));
h:=IP(t*24);
t:=FP(t*24);
m:=IP(t*60);
t:=FP(t*60);
s:=IP(t*60);
lista:={y,o,d,h,m,s};
RETURN lista;
END;

makeDateListShort(lista)
BEGIN
// lista {Y m D h mn, s}
LOCAL y,m,d, listacorta;
y:=lista(1);
m:=lista(2);
d:=lista(3)+lista(4)/24+lista(5)/1440+lista(6)/86400;
listacorta:={y,m,d};
RETURN listacorta;
END;

tempo()
BEGIN
  LOCAL n, JD, JDE, TSid, dep, psi, eps;
  LOCAL ε;
  ε:= epsilon();
  n:= calcN();
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  TSid:= calcTS();
  dep:= deltaEpsilonPsi();
  psi:=dep(1);
  eps:=dep(2);
  PRINT;
  PRINT("Date " + m + " " + D + " " + Y);
  PRINT("Julian Day " + JD);
  PRINT("Julian Day Effem. " + JDE);
  PRINT("N (from 1901jan1) " + n);
  PRINT("T (from JD) " + TSid(5));
  PRINT(" ");
  PRINT("Tempo Siderale 0h " + →HMS(TSid(1)));
  PRINT("Tempo Siderale locale " + →HMS(TSid(2)));
  PRINT("Tempo Siderale apparente " + →HMS(TSid(3)));
  PRINT("θ0 Mean ST GMT " + →HMS(TSid(4)));
  PRINT("Tempo Siderale (DEG) " + trunc_angle(TSid(1)*15));
  PRINT("θ0 Mean ST GMT (DEG) " + trunc_angle(TSid(4)*15));
  PRINT(" ");
  PRINT("Δψ Nutazione (longit.) "+ psi);
  PRINT("Δε Nutazione (obl.) "+ eps);
  PRINT("ε Obliquità Ecl. "+ ε);
  PRINT(" ");
  PRINT("Time " + lstd + " GMT " + gmt);
  RETURN({ROUND(JD,5), ROUND(JDE,5), trunc_angle(TSid(1)), trunc_angle(TSid(2)), 
    trunc_angle(TSid(3)), trunc_angle(TSid(4)), trunc_angle(TSid(5)) });
END;

transit(alfa, delta, h0)
BEGIN
// α, δ (AR, dec), height, h0 = standard height of body
  LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;
  LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
  α:= alfa; // AR
  δ:= delta;
  temp:= calcTS();
  // θ0:= apparentSideralTime(makeDateListShort({Y,m,D,0,0,0})) *15; // TS in DEG
  θ0 := temp(1) * 15;
  temp:= (sin(h0)-sin(lat)*sin(δ))/(cos(lat)*cos(δ));
  IF temp<-1 OR temp>1 THEN culm:=0; END; // above horizon?
  H0:= acos(temp);
  //IF H0<-180 THEN H0:=H0+180; END; IF H0>180 THEN H0:= H0-180; END; // it should be from -180 to 180
 
  m0:= HMS→(α*15 + long - θ0) / 360;
  IF m0<0 THEN m0:= m0 +1; END; IF m0>1 THEN m0:=m0-1; END;
  tsideral:= simp_angle(θ0 + 360.985647*m0); // TS transit DEG
  Hor:= HMS→(tsideral - long -α*15);  // Hour angle, local ST vs apparent ST
  Hor:=simp_angle(→HMS(Hor));
  deltaM:= -HMS→(Hor/360);
  culm:= (m0 + deltaM)*24;  // transit

  m1:= m0-H0/360;
  IF m1<0 THEN m1:= m1 +1; END; IF m1>1 THEN m1:=m1-1; END;
  tsideral:= simp_angle(θ0 + 360.985647*m1); // TS rising
  Hor:= HMS→(tsideral - long - α*15);
  Hor:=simp_angle(→HMS(Hor));
  h:=asin(sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hor));
  deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));
  sorg:= ((m1 + deltaM)*24); // rising

  m2:= m0+H0/360;
  IF m2<0 THEN m2:= m2 +1; END; IF m2>1 THEN m2:=m2-1; END;
  tsideral:= simp_angle(θ0 + 360.985647*m2); // TS setting
  Hor:= HMS→(tsideral - long - α*15); 
  Hor:=simp_angle(→HMS(Hor));
  h:=asin(sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hor));
  deltaM:= (HMS→(h)-HMS→(h0))/(360*cos(δ)*cos(lat)*sin(Hor));
  tram:= ((m2 + deltaM) * 24); // setting

RETURN ({culm, sorg, tram});
END;

parametri()
BEGIN
  LOCAL dep, prec, α, δ, ε;
  dep:= deltaEpsilonPsi();
  α:= sunlist(5); δ:= sunlist(6);
  prec:= precession(α, δ);
  ε:= epsilon();
  // Nutazione Δψ (longitudine), Δε (obliquità), ε (inclinazione eclittica)...
  PRINT;
  PRINT("ε Obliquità eclittica "+ ε);
  PRINT("Δψ Nutazione longitudine "+ dep(1));
  PRINT("Δε Nutazione longitudine "+ dep(2));
  PRINT(" ");
  PRINT("Δα Precessione AR "+ prec(1));
  PRINT("Δδ Precessione dec "+ prec(2));
  RETURN({ε, dep(1), dep(2), prec(1), prec(2)}); // ε, Δψ, Δε, Δα, Δδ
END;

ascendent()
BEGIN
  LOCAL ASC, ts, segno, ε;
  ε:= epsilon;
  ts:= calcTS();
  ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));
  ASC:= simp_angle(ASC);
  segno:= zodiaco(ASC);
RETURN {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno};
END;

zodiaco(astro)
BEGIN
  LOCAL segno;
  CASE  // Segni e case
  IF astro>=0   AND astro <30  THEN segno:="Ariete"; END;
  IF astro>=30  AND astro <60  THEN segno:="Toro"; END;
  IF astro>=60  AND astro <90  THEN segno:="Gemelli"; END;
  IF astro>=90  AND astro <120 THEN segno:="Cancro"; END;
  IF astro>=120 AND astro <150 THEN segno:="Leone"; END;
  IF astro>=150 AND astro <180 THEN segno:="Vergine"; END;
  IF astro>=180 AND astro <210 THEN segno:="Bilancia"; END;
  IF astro>=210 AND astro <240 THEN segno:="Scorpione"; END;
  IF astro>=240 AND astro <270 THEN segno:="Sagittario"; END;
  IF astro>=270 AND astro <300 THEN segno:="Capricorno"; END;
  IF astro>=300 AND astro <330 THEN segno:="Acquario"; END;
  IF astro>=330 AND astro <360 THEN segno:="Pesci"; END;
  END; // case
  RETURN segno;
END;

atan2(y,x)
BEGIN
// atan2 returns coordinates in correct angle for all four quadrants (DEG)  
  CASE
    IF x==0 AND y==0 then return "undefined"; end;         // x=0, y=0  
    IF x>0 then return atan(y/x); end;                     // x>0
    IF x<0 AND y>=0 then return atan(y/x)+180; end;         // x<0, y>=0
    IF x==0 AND y>0 then return 90; end;                 // x=0, y>0
    IF x<0 AND y<=0 then return atan(y/x)-180; end;         // x<0, y<0
    IF x==0 AND y<0 then return -90; end;                // x=0, y<0 
  END;
  RETURN; 
END;

trunc_angle(ang)
// trunc decimal from a DEG form angle
BEGIN
→HMS((IP(HMS→(ang)*3600)/3600));
END;

simp_angle(ang)
BEGIN
LOCAL angle;
angle:=360*FP(ang/360);
IF angle<0
THEN
angle:=angle+360
END;
RETURN angle;
END;

kepler(ec, M) 
// needs eccentricity and true anomaly M
BEGIN
  LOCAL M1, j, u;
  HAngle:=0; // RAD
  M1:= HMS→(M)*PI/180; // anomalia RAD
  F:= SIGN(M1); M1:= ABS(M1)/(2*PI);
  M1:=(M1-IP(M1))*2*PI*F;
  IF M1<0 THEN M1:= M1+2*PI; END;
  F:= 1; 
  IF M1>PI THEN F:= -1; END;
  IF M1>PI THEN M1:=2*PI-M; END;
  u:= PI/2; D:=PI/4;
  FOR j FROM 1 TO 33 DO // 33 iterations (Meeus, Roger Sinnot)
  M2:= u-ec*sin(u); u:= u+D*SIGN(M1-M2); D:= D/2;
  END; // for
  u:= u * F; // anomalia eccentrica 
  HAngle:=1;
  u:= u*180/PI;
  RETURN u;
END;

interpolation(lista,n)
BEGIN
LOCAL a,b,c,R,dif1,dif2;
dif1:=ΔLIST(lista);
dif2:=ΔLIST(dif1);
a:=dif1(1);
b:=dif1(2);
c:=dif2(1);
R:=lista(2)+(n/2)*(a+b+n*c);
RETURN R;
END;

sunCalc()
BEGIN
  LOCAL ω, l, v, r, x, y, dz, ε;
  LOCAL α, δ, dayY, H, ts, az, z, h;
  LOCAL SAD, DA, sorg, tram, culm, Hd, lt;
  LOCAL as, ac, zc, h0, H0, temp;
  LOCAL ec, a, u, M1, M2, j, t, tau;
  LOCAL JD, JDE, c, lapp, Ω, θ0, transito;
  LOCAL dep, deltapsi, longitude;
  ε:= epsilon;
  N:= calcN();
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  a:= 1.00000023; // semiasse maggiore
  dep:= deltaEpsilonPsi();
  deltapsi:= dep(1);
  t:= N/36525;
  T:=   (JD-2451545.0)/36525; 
  tau:= (JDE-2451545.0)/365250; // millenni  - JDE
  ec:= 0.016708617-0.000042037*T-0.0000001236*T^2; // eccentricità Meeus
  // ec:= ec*180/PI; // DEG
  // L:= →HMS(278.9654+0.985647342*N + 0.0003*t^2) MOD 360; // long media Bouiges
  // L:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Meeus
  L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
  ω:= →HMS(281.238 + 0.000047067*N + 0.00045*t^2) MOD 360; // long del perielio
  // M:= (L - ω); // anomalia vera
  //IF M<0 THEN M:= (360+M); END; // potrebbe essere negativa
  Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
  M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-0.00000048*T^3); // anom vera Meeus
  u:=atan2(sin(M),(cos(M)-ec)); // only for small eccentricity
  c:= (1.914600-0.004817*T-0.000016*T^2)*sin(M);  // equation of center
  c:= c+(0.019993-0.000101*T)*sin(2*M)+(0.000290)*sin(3*M);
  // v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio (Bouiges)
  // v:= acos((cos(u)-ec)/(1-ec*cos(u))); // kepler
  v:= M + c; // true anomaly, Meeus
  // l:= (v + ω) MOD 360; // Longitudine Sole
  l:= simp_angle(L+c); // Ө Longitudine Sole, Meeus
  lt:= l + 180; // Longitudine eliocentrica della Terra
  lapp:= simp_angle(l - 0.00569-0.00478*sin(Ω)); // Longitudine apparente (ref true equinox)
  // r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)
  r:= (1.000001018*(1-ec^2))/(1+ec*cos(v)); // Distanza Sole-Terra (UA) Meeus
  x:= r * cos(l); y:= r * sin(l);

  δ:= →HMS((asin(sin(ε)*sin(l))) );
  α:= →HMS(atan2(cos(ε)*sin(l),cos(l))); // Ascensione retta, Meeus (b never > 1.2", l=Ө)
  α:= →HMS(α/15) MOD 24; // α in hms
  ts:= calcTS();  // calcola tempo siderale
 // H:= (ts(2) - α) MOD 24; // Angolo orario
 // Hd:= →HMS(15*H); // Angolo orario in gradi

  Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST
  Hd:=simp_angle(→HMS(Hd));
  H:= Hd/15; // in HMS
  
  temp:= sunAh(dateshortlist);
  az:= temp(1);  // Azimuth
  h:=  temp(2);  // height

  // E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)
  E:= 4*( L - 0.0057183 - α*15 + deltapsi*cos(ε) );  // α in DEG
  DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
  SAD:= (90 + DA); // Semiarco Diurno
  //   SAD:= →HMS(acos(-tan(lat)*tan(δ))); // Semiarco Diurno
 
   // culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione, transito (Bouiges)
   // sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
   // tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)

  // transit routine
  h0:= -0°50′; // 'standard' altitude of Sun
  transito:= transit(α, δ, h0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  sunlist:= {L,ω,l, M, α,δ, az, h, x,y,
             sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA,
             SAD, E, lt, v, c, Ω, lapp, ec};
  RETURN sunlist;
END;

sun()
BEGIN
  LOCAL ω, l, v, r, x, y, dz, ε;
  LOCAL α, δ, dayY, ts, az, z, h;
  LOCAL SAD, DA, sorg, tram, culm, Hd, lt;
  LOCAL as, ac, zc, segno, JD, JDE, v, c;
  LOCAL Ω, lapp, ec;
  ε:= epsilon;
  L:= sunlist(1); ω:= sunlist(2); 
  l:= sunlist(3); M:= sunlist(4);
  Ω:= sunlist(26); lapp:=sunlist(27); 
  α:= sunlist(5); δ:= sunlist(6); 
  az:= sunlist(7); h:= sunlist(8); 
  x:= sunlist(9); y:= sunlist(10);
  sorg:= sunlist(11); tram:= sunlist(12); culm:= sunlist(13);
  r:= sunlist(14); H:= sunlist(15); Hd:= sunlist(16);
  zc:= sunlist(17); dz:= sunlist(18);
  v:= sunlist(24); c:= sunlist(25); 
  ts:= calcTS(); //
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);
  ec:= sunlist(28);
  segno:= zodiaco(l);

  PRINT;
  PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);
  PRINT("at Long " + long + " Lat " + lat);
  PRINT("");
  PRINT("L longitudine media " + trunc_angle(L));
  PRINT("ω long del perielio " + trunc_angle(ω));
  PRINT("M anomalia media " + trunc_angle(M));
  PRINT("v anomalia vera " + trunc_angle(v));
  PRINT("c equazione del centro " + c);
  PRINT(" ");
  PRINT("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("λ longitudine apparente " + trunc_angle(lapp));
  PRINT("Ω aberr, nutaz nodo "+ Ω);
  PRINT("r distanza (UA) " + r);
  PRINT(" ");
  PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
  PRINT("α Asc.r. (DEG) " + trunc_angle(α*15));
  PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
  PRINT(" ");
  PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));
  PRINT("TS (UTC) " + trunc_angle(ts(1)));
  PRINT("TS locale " + trunc_angle(ts(2)));
  PRINT("TS apparente " + ts(3));
  PRINT("θ0 Mean ST GMT " + ts(4));
  PRINT("Equazione del tempo (min) " + trunc_angle(E));
  PRINT("Julian day " + JD);
  PRINT("Dif asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
  PRINT("");
  PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
  PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)"); 
  PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
  PRINT(" ");
  PRINT("Long elioc. Terra "+ trunc_angle(lt));
  PRINT("ε obliquità ecl. " + ε);
  PRINT("ec eccentricità "+ ec);
  // RETURN {L, ω, M, v, l, r, x, y};
  RETURN [[L, ω], [l, lapp], [M,v], [Ω, c], [α,δ], [az, h], [x,y], 
  [sorg, tram], [culm, r], [H,ts(1)], [lt, 0], [ε, ec]];
END;

sunAh(lista)
// lista like dateshortlist
BEGIN
  LOCAL lha,alpha,delta,phi, ah, α;
  LOCAL temp,longitude, Azimuth, Altitude, ts;
  alpha:=15*HMS→(sunlist(5));
  α:= sunlist(5); // α in HMS
  delta:=sunlist(6);
  phi:=lat;
  ts:= calcTS();
  // longitude:= HMS→(long); // - if long -E
  // lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
  lha:= 15*HMS→(ts(2)-α);  // local ST vs apparent ST;
  // IF lha<0 THEN lha:=lha+360 END;
  // IF lha>=360 THEN lha:=lha-360 END;
  lha:=simp_angle(→HMS(lha));

  temp:= transformEquatorialToHorizontal(lha, phi, delta);
  Azimuth:= HMS→(temp(1))+180;
  IF Azimuth>360 THEN Azimuth:=Azimuth-360 END;
  Azimuth:= →HMS(Azimuth);
  Altitude:= temp(2);
  ah:= {Azimuth,Altitude};
  RETURN ah;
END;

moon()
BEGIN
  LOCAL ω,Ω, L1, M1, l, lsum, bsum, rsum; 
  LOCAL v, r, x, y, dep, ε;
  LOCAL b, P, d, s, temp;
  LOCAL α,δ, ts, Hd, zc, dz, h, az;
  LOCAL i, as, ac, DA, SAD, DD;
  LOCAL culm, sorg, tram, h0, transito;
  LOCAL xs, ys, ph, fase, eta;
  LOCAL d1, h1, m1, segno, JD, JDE, rr;
  LOCAL lapp, ill, ee, a1,a2,a3, longitude;
  i:= 5.13333333; // inclinazione 5°8' sull'eclittica
  ε:= epsilon();
  N:= calcN();
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:= (JDE-2451545.0)/36525; // use JDE
  dep:= deltaEpsilonPsi();
  // Sun
  L:= simp_angle(280.4665+36000.76983*T+0.0003032*T^2); // Long media Sole Meeus
  M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+T^3/24490000); // mean anomaly Sun   
  // end Sun
  // L1:= (33.231 + 13.17639653*N) MOD 360; // long. media
  L1:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);
  //  Long media Luna Meeus
  // Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente
  Ω:= simp_angle(125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000); // Nodo Luna (Moon)
  // M1:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera
  M1:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media
  ω:= (L1 - M1) MOD 360; // perigeo medio
  // DD:= L1 - L; F:= L1-Ω; // for the variations
  DD:= simp_angle(297.8502042+445267.1115168*T-0.0016300*T^2+T^3/545868-T^4/113065000); // mean elongation
  F:= simp_angle(93.2720993+483202.0175273*T-0.0034029*T^2-T^3/3526000+T^4/863310000); // dist from asc. node
  // ee:= 1-0.002516*T-0.0000074*T^2;  // multiply those w/ M *ee
  a1:= simp_angle(119.75+131.849*T); // action of Venus
  a2:= simp_angle(53.09+479264.290*T); // action of Jupiter
  a3:= simp_angle(313.45+481266.484*T); 

  lsum:=0;  // da sommare a L1
  lsum:= lsum + 6288774 * sin(M1); // calc Moon longitude
  lsum:= lsum + 1274027 * sin(2*DD-M1); // evezione
  lsum:= lsum +  658314 * sin(2*DD); // variazione
  lsum:= lsum +  213618 * sin(2*M1); // equazione del centro
  lsum:= lsum -  185116 * sin(M) *ee; // equazione annuale
  lsum:= lsum -  114332 * sin(2*F); // riduzione all'eclittica
  lsum:= lsum +   58793 * sin(2*DD-2*M1);
  lsum:= lsum +   57066 * sin(2*DD-M1-M) *ee;
  lsum:= lsum +   53322 * sin(2*DD+M1);
  lsum:= lsum +   45758 * sin(2*DD-M) *ee;
  lsum:= lsum -   40923 * sin(M-M1) *ee;
  lsum:= lsum -   34720 * sin(DD);
  lsum:= lsum -   30383 * sin(M1+M) *ee;
  lsum:= lsum +   15327 * sin(2*DD-2*F);
  lsum:= lsum -   12528 * sin(M1+2*F);
  lsum:= lsum +   10980 * sin(M1-2*F);
  lsum:= lsum +   10675 * sin(4*DD-M1);
  lsum:= lsum +   10034 * sin(3*M1);
  lsum:= lsum +    8548 * sin(4*DD-2*M1);
  lsum:= lsum -    7888 * sin(2*DD+M-M1) *ee;
  lsum:= lsum -    6766 * sin(2*DD+M) *ee;
  lsum:= lsum -    5163 * sin(DD-M1);
  lsum:= lsum +    4987 * sin(DD+M) *ee;
  lsum:= lsum +    4036 * sin(2*DD-M+M1) *ee;
  lsum:= lsum +    3994 * sin(2*DD+2*M1);
  lsum:= lsum +    3861 * sin(4*DD);
  lsum:= lsum +    3665 * sin(2*DD-3*M1);
  lsum:= lsum -    2689 * sin(M-2*M1) *ee;
  lsum:= lsum -    2602 * sin(2*DD-M1+2*F);
  lsum:= lsum +    2390 * sin(2*DD-M-2*M1) *ee;
  lsum:= lsum -    2348 * sin(DD+M1);
  lsum:= lsum +    2236 * sin(2*DD-2*M) *ee^2;
  lsum:= lsum -    2120 * sin(M+2*M1) *ee;
  lsum:= lsum -    2069 * sin(2*M) *ee^2;
  lsum:= lsum +    2048 * sin(2*DD-2*M-M1) *ee^2;
  lsum:= lsum -    1773 * sin(2*DD+M1-2*F);
  lsum:= lsum -    1595 * sin(2*DD+2*F);
  lsum:= lsum +    1215 * sin(4*DD-M-M1) *ee;
  lsum:= lsum -    1110 * sin(2*M1+2*F);
  lsum:= lsum -     892 * sin(3*DD-M1);
  lsum:= lsum -     810 * sin(2*DD+M+M1) *ee;
  lsum:= lsum +     759 * sin(4*DD-M-2*M1) *ee;
  lsum:= lsum -     713 * sin(2*M-M1) * ee^2;
  lsum:= lsum -     700 * sin(2*DD+2*M-M1) *ee^2;
  lsum:= lsum +     691 * sin(2*DD+M-2*M1) *ee;
  lsum:= lsum +     596 * sin(2*DD-M-2*F) *ee;
  lsum:= lsum +     549 * sin(4*DD+M1);
  lsum:= lsum +     537 * sin(4*M1);
  lsum:= lsum +     520 * sin(4*DD-M) *ee;
  lsum:= lsum -     487 * sin(DD-2*M1);
  lsum:= lsum -     399 * sin(2*DD+M-2*F) *ee;
  lsum:= lsum -     381 * sin(2*M1-2*F);
  lsum:= lsum +     351 * sin(DD+M+M1) *ee;
  lsum:= lsum -     340 * sin(3*DD-2*M1);
  lsum:= lsum +     330 * sin(4*DD-3*M1);
  lsum:= lsum +     327 * sin(2*DD-M+2*M1) *ee;
  lsum:= lsum -     323 * sin(2*M+M1) *ee^2;
  lsum:= lsum +     299 * sin(DD+M-M1) *ee;
  lsum:= lsum +     294 * sin(2*DD+3*M1);

  lsum:= lsum + 3958*sin(a1)+1962*sin(L1-F)+318*sin(a2);

  l:= (simp_angle(L1 + lsum/1000000)); // λ divided by 1 milion

  lapp:= →HMS(l + dep(1)); // longitudine apparente (+ nutazione long, Δψ)

  // fattori per latitudine
  bsum:=        5128122 * sin(F); // termine principale latitudine
  bsum:= bsum +  280026 * sin(M1+F);
  bsum:= bsum +  277693 * sin(M1-F); // equazione del centro
  bsum:= bsum +  173237 * sin(2*DD-F); // grande ineguaglianza
  bsum:= bsum +   55413 * sin(2*DD-M1+F); // evezione
  bsum:= bsum +   46271 * sin(2*DD-M1-F); // evezione (2)
  bsum:= bsum +   32573 * sin(2*DD+F);
  bsum:= bsum +   17198 * sin(2*M1+F);
  bsum:= bsum +    9266 * sin(2*D+M1-F);
  bsum:= bsum +    8822 * sin(2*M1-F);
  bsum:= bsum +    8216 * sin(2*DD-M-F) *ee;
  bsum:= bsum +    4324 * sin(2*DD-2*M1-F);
  bsum:= bsum +    4200 * sin(2*DD+M1+F);
  bsum:= bsum -    3359 * sin(2*DD+M-F) *ee;
  bsum:= bsum +    2463 * sin(2*DD-M-M1+F) *ee;
  bsum:= bsum +    2211 * sin(2*DD-M+F) *ee;
  bsum:= bsum +    2065 * sin(2*DD-M-M1-F) *ee;
  bsum:= bsum -    1870 * sin(M-M1-F) *ee;
  bsum:= bsum +    1828 * sin(4*DD-M1-F);
  bsum:= bsum -    1794 * sin(M+F) *ee;
  bsum:= bsum -    1749 * sin(3*F);
  bsum:= bsum -    1565 * sin(M-M1+F) *ee;
  bsum:= bsum -    1491 * sin(DD+F);
  bsum:= bsum -    1475 * sin(M+M1+F) *ee;
  bsum:= bsum -    1410 * sin(M+M1-F) *ee;
  bsum:= bsum -    1344 * sin(M-F) *ee;
  bsum:= bsum -    1335 * sin(DD-F);
  bsum:= bsum +    1107 * sin(3*M1+F);
  bsum:= bsum +    1021 * sin(4*DD-F);
  bsum:= bsum +     833 * sin(4*DD-M1+F);
  bsum:= bsum +     777 * sin(M1-3*F);
  bsum:= bsum +     671 * sin(4*DD-2*M1+F);
  bsum:= bsum +     607 * sin(2*DD-3*F);
  bsum:= bsum +     596 * sin(2*DD+2*M1-F);
  bsum:= bsum +     491 * sin(2*DD-M+M1-F) *ee;
  bsum:= bsum -     451 * sin(2*DD-2*M1+F);
  bsum:= bsum +     439 * sin(3*M1-F);
  bsum:= bsum +     422 * sin(2*DD+2*M1+F);
  bsum:= bsum +     421 * sin(2*DD-3*M1-F);
  bsum:= bsum -     366 * sin(2*DD+M-M1+F) *ee;
  bsum:= bsum -     351 * sin(2*DD+M+F) *ee;
  bsum:= bsum +     331 * sin(4*DD+F);
  bsum:= bsum +     315 * sin(2*DD-M+M1+F) *ee;
  bsum:= bsum +     302 * sin(2*DD-2*M-F) *ee^2;
  bsum:= bsum -     283 * sin(M1+3*F);
  bsum:= bsum -     229 * sin(2*DD+M+M1-F) *ee;
  bsum:= bsum +     223 * sin(DD+M-F) *ee;
  bsum:= bsum +     223 * sin(DD+M+F) *ee;
  bsum:= bsum -     220 * sin(M-2*M1-F) *ee;
  bsum:= bsum -     220 * sin(2*DD+M-M1-F) *ee;
  bsum:= bsum -     185 * sin(DD+M1+F);
  bsum:= bsum +     181 * sin(2*DD-M-2*M1-F) *ee;
  bsum:= bsum -     177 * sin(M+2*M1+F) *ee;
  bsum:= bsum +     176 * sin(4*DD-2*M1-F);
  bsum:= bsum +     166 * sin(4*DD-M-M1-F) *ee;
  bsum:= bsum -     164 * sin(DD+M1-F);
  bsum:= bsum +     132 * sin(4*DD+M1-F);
  bsum:= bsum -     119 * sin(DD-M1-F);
  bsum:= bsum +     115 * sin(4*DD-M-F) *ee;
  bsum:= bsum +     107 * sin(2*DD-2*M+F) *ee^2;

  bsum:= bsum-2235*sin(L1)+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(L1-M1)-115*sin(L1+M1);
  b := →HMS(bsum/1000000); // β div by 1 milion

  // correzioni per la distanza dalla Terra (Δ)
  rsum:=0;
  rsum:= rsum -20905355 * cos(M1); // coseni
  rsum:= rsum - 3699111 * cos(2*DD-M1); // evezione
  rsum:= rsum - 2955968 * cos(2*DD); // variazione
  rsum:= rsum -  569925 * cos(2*M1); // equazione del centro
  rsum:= rsum +   48888 * cos(M) *ee; // equazione annuale
  rsum:= rsum -    3149 * cos(2*F); // riduzione all'eclittica
  rsum:= rsum +  246158 * cos(2*DD-2*M1);
  rsum:= rsum -  152138 * cos(2*DD-M1-M) *ee;
  rsum:= rsum -  170733 * cos(2*DD+M1);
  rsum:= rsum -  204586 * cos(2*DD-M) *ee;
  rsum:= rsum -  129620 * cos(M-M1) *ee;
  rsum:= rsum +  108743 * cos(DD);
  rsum:= rsum +  104755 * cos(M1+M) *ee;
  rsum:= rsum +   10321 * cos(2*DD-2*F);
  rsum:= rsum +   79661 * cos(M1-2*F);
  rsum:= rsum -   34782 * cos(4*DD-M1);
  rsum:= rsum -   23210 * cos(3*M1);
  rsum:= rsum -   21636 * cos(4*DD-2*M1);
  rsum:= rsum +   24208 * cos(2*DD+M-M1) *ee;
  rsum:= rsum +   30824 * cos(2*DD+M) *ee;
  rsum:= rsum -    8379 * cos(DD-M1);
  rsum:= rsum -   16675 * cos(DD+M) *ee;
  rsum:= rsum -   12831 * cos(2*DD-M+M1) *ee;
  rsum:= rsum -   10445 * cos(2*DD+2*M1);
  rsum:= rsum -   11650 * cos(4*DD);
  rsum:= rsum +   14403 * cos(2*DD-3*M1);
  rsum:= rsum -    7003 * cos(M-2*M1) *ee;
  rsum:= rsum +   10056 * cos(2*DD-M-2*M1) * ee;
  rsum:= rsum +    6322 * cos(DD+M1);
  rsum:= rsum -    9884 * cos(2*DD-2*M) * ee^2;
  rsum:= rsum +    5751 * cos(M+2*M1) *ee;
  rsum:= rsum -    4950 * cos(2*DD-2*M-M1) *ee^2;
  rsum:= rsum +    4130 * cos(2*DD+M1-2*F);
  rsum:= rsum -    3958 * cos(4*DD-M-M1) *ee;
  rsum:= rsum +    3258 * cos(3*DD-M1);
  rsum:= rsum +    2616 * cos(2*DD+M+M1) *ee;
  rsum:= rsum -    1897 * cos(4*DD-M-2*M1) *ee;
  rsum:= rsum -    2117 * cos(2*M-M1) * ee^2;
  rsum:= rsum +    2354 * cos(2*DD+2*M-M1) *ee^2;
  rsum:= rsum -    1423 * cos(4*DD+M1);
  rsum:= rsum -    1117 * cos(4*M1);
  rsum:= rsum -    1571 * cos(4*DD-M) *ee;
  rsum:= rsum -    1739 * cos(DD-2*M1);
  rsum:= rsum -    4421 * cos(2*M1-2*F);
  rsum:= rsum +    1165 * cos(2*M+M1) *ee^2;
  rsum:= rsum +    8752 * cos(2*DD-M1-2*F);

  rsum:= rsum/1000;  // div by 1000
 
  temp:= transformEclipticalToEquatorial(l,b);  // transform into AR, decl
  α:= temp(1);
  δ:= temp(2);

  ts:= calcTS();  // calcola tempo siderale
  Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST
  Hd:=simp_angle(→HMS(Hd));
  H:= Hd/15; // in HMS
  temp:= transformEquatorialToHorizontal(Hd, lat, δ);  
  az:= HMS→(temp(1))+180; // azimuth
  IF az>360 THEN az:=az-360 END;
  az:= →HMS(az);
  h:= temp(2);  // height

  // Raggio quadratico medio della Terra: 6373,044737
  d:= (385000.56 + rsum); // Δ distanza Terra-Luna in km (=1/sin(P)),rr in 10^-3
  //P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse
  P:= asin(6378.14/d); // Parallasse
  s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)

  DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
  SAD:= (90 + DA); // Semiarco Diurno
  // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
    // culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione (Bouiges)
    // sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
    // tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)

  // transit routine
  h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)
  // h0:= 0°7′30″; // medium value: 0.125°
  transito:= transit(α, δ, h0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  ph:= abs(l - L); // Phase  Moon-Sun
  CASE
  // IF ph == 0 THEN fase:="New Moon"; END;
  IF ph>280 AND ph<350 THEN fase:="Decrescente"; END;
  IF ph>=260 AND ph<=280 THEN fase:="Last Quart"; END;
  IF ph>190 AND ph<260 THEN fase:="Gibbosa Cal."; END;
  IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
  IF ph>100 AND ph<170 THEN fase:="Gibbosa Cresc."; END;
  IF ph>=80 AND ph<=100 THEN fase:="1st Quart"; END;
  IF ph>10 AND ph<80 THEN fase:="Crescente"; END;
  DEFAULT fase:="New Moon";
  END; // case
  // età della Luna
  d1:= IP(HMS→(ph/15));
  h1:= IP(FP(HMS→(ph/15))*60);
  IF (h1>24) THEN d1:= d1+1; h1:= h1 MOD 24; END;
  m1:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);
  eta:= →HMS(d1+h1/60+m1/3600); // Age of the Moon
  segno:= zodiaco(l);
  ill:= 1-(1+cos(ph))/2; // illuminated fraction

  PRINT;
  PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
  PRINT("L Long media " + trunc_angle(L1));
  PRINT("M anomalia vera " + trunc_angle(M1));
  PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("β Latitudine " + trunc_angle(b));
  PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
  PRINT("Δψ Nutazione long "+ dep(1));
  PRINT("Δε Nutazione obl. " + dep(2));
  PRINT("ε inclinazione ecl. " + ε);
  PRINT("ω Perigeo medio " + trunc_angle(ω));
  PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
  PRINT(" ");
  PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
  PRINT("α Ascensione retta (gradi) " + trunc_angle((α*15) MOD 360));
  PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
  PRINT("");
  PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));
  PRINT("TS (UTC) " + trunc_angle(ts(1)));
  PRINT("TS locale " + trunc_angle(ts(2)));
  PRINT("TS apparente " + ts(3));
  PRINT("θ0 Mean ST GMT " + ts(4));
  PRINT("Julian day " + JD);
  PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
  PRINT("P Parallasse " + trunc_angle(P));
  PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");
  PRINT("Semidiametro "+  trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,3) + " km");
  PRINT(" ");
  PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
  PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)"); 
  PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
  PRINT(" ");
  PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");
  PRINT("Età " + eta);
  PRINT("Parte illuminata " + ill);
  RETURN [[L1, M1],[l, b],[lapp,P], [ω,Ω],[α,δ], [az, h], 
    [ROUND(HMS→(d),6),s], [sorg, tram], [culm, ts(3)],[ph, eta] ];
END;

planetCalc(nameplanet)
BEGIN
  LOCAL ω,Ω, l, v, r, x, y, z;
  LOCAL Ls, Ms, b, xs, ys, b_ec;
  LOCAL x1, y1, z1, R, l_ec, H, Hd;
  LOCAL ec, in, a, δ, α, ts, az, h;
  LOCAL DA, SAD, tram, sorg, culm, dayY;
  LOCAL as, ac, zc, dz, θ, mov, segno, ε;
  LOCAL transito, h0, JD, temp;

  JD:= julianDay(dateshortlist);

  // calc of Sun (as basis)
  xs:= sunlist(9); ys:= sunlist(10);
  Ls:= sunlist(1); // Longitudine media
  Ms:= sunlist(4); // Anomalia
  // end Sun

 EVAL(EXPR(nameplanet+"()")); // call planet-name function

  L:= planet_const(1); // Longitudine media
  ω:= planet_const(2); // perielio
  Ω:= planet_const(3); // nodo
  M:= planet_const(4); // anomalia
  v:= planet_const(5); 
  ec:= planet_const(6); // eccentricità
  in:= planet_const(7); // inclinazione
  a:= planet_const(8);  // semiasse maggiore
  θ:= planet_const(9);  // θs to calc if retrogade

  b_ec:= asin(sin(in)*sin(v+ω-Ω)); // argomento di latitudine del perielio
  l_ec:= acos(cos(v+ω-Ω)/cos(b_ec));  // long from eclittic
  IF b_ec<0 THEN l_ec:= 360 - l_ec; END; // argomento secondario
  l_ec:= →HMS((l_ec + Ω)) MOD 360; // aggiunge Nodo e riporta a 2PI
  r:= a*(1-ec^2)/(1+ec*cos(v)); // 9.524899/(1+ 0.0559*cos(v)); raggio vettore
  x1:= r*cos(b_ec)*cos(l_ec); y1:= r*cos(b_ec)*sin(l_ec); z1:=r*sin(b_ec); // cartesian
  x:= x1+xs; y:= y1+ys; z:=z1; // cohordinates cartesian->spheric
  R:= sqrt(x^2+y^2+z^2); 
  b:=asin(z/R); // latitudine  β
  IF b> 360 THEN b:= b MOD 360; END;
  // l:=→HMS(atan(y/x)); // cart->spheric longtitudine λ
  l:= atan2(y,x);
  IF l<0 THEN l:= l+360; END; IF l>360 THEN l:=l-360; END;
  IF l>180 THEN α:= α +180; END;
  // δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(l))) ); // declinazione (coord locali)
  // α:= (acos(cos(b)*cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta
  //α:= →HMS(α/15) MOD 24; // α in hms

  temp:= transformEclipticalToEquatorial(l,b);  // transform into AR, decl
  α:= temp(1);
  δ:= temp(2);

  ts:= calcTS();  // calcola tempo siderale
  H:= (ts(2) - α) MOD 24; // Angolo orario
  Hd:= →HMS(15*H); // Angolo orario in gradi
  zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt
  dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza
  h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)
  as:= cos(δ)*sin(Hd)/sin(dz); // calc az
  ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);
  // az:= →HMS(ATAN(as/ac) ); // Azimuth
  az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));
  IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;
  DA:= →HMS(asin(tan(lat)*tan(δ))); // Differenza ascensionale
  SAD:= →HMS(90 + DA); // Semiarco Diurno
  // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
  // culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) MOD 24 ; // Culminazione
  // sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
  // tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)

  // transit routine
  h0:= -0°34′; // 'standard' altitude of a planet
  transito:= transit(α, δ, h0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  IF abs(l-220.4)>θ THEN mov:="diretto"; ELSE mov:="retrogrado"; END;
  segno:= zodiaco(l);

  PRINT;
  PRINT("Effemeridi di " + nameplanet + " a " + m+" "+D+" "+Y+" "+ lstd);
  PRINT("L Long media " + trunc_angle(L));
  PRINT("ω Perielio " + trunc_angle(ω));
  PRINT("Ω Nodo " + trunc_angle(Ω));
  PRINT("M Anomalia vera " + trunc_angle(M));
  PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));
  PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(v));
  PRINT("r Raggio vettore " + r);
  PRINT("Cartes. elioc.: x' " + ROUND(x1,3) + " y' " + ROUND(y1,3) + " z' " + ROUND(z1,3));
  PRINT("Cartes. geoc.: x " + ROUND(x,3) + " y " + ROUND(y,3) + " z " + ROUND(z,3));
  PRINT("Spher.: λ Longitudine " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("Spher.: β Latitudine " + trunc_angle(b));
  PRINT("Spher.: R Distanza " + ROUND(R,3));
  PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
  PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
  PRINT("");
  PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));
  PRINT("TS (UTC) " + trunc_angle(ts(1)));
  PRINT("TS locale " + trunc_angle(ts(2)));
  PRINT("TS apparente " + ts(3));
  PRINT("θ0 Mean ST GMT " + ts(4));
  PRINT("Julian day " + JD);
  PRINT("Diff. ascensionale " + trunc_angle(DA));
  PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));
  PRINT(" ");
  PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");
  PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)"); 
  PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");
  PRINT("Movimento " + mov);
  RETURN [[L,ω,Ω ],[M,l_ec,b_ec], [v,r,0], [x1,y1,z1],[x,y,z],[l,b,R], [α,δ,H], [az,h,ts(1)],[sorg, tram, culm]];
END;

Mercurius()
BEGIN
  // Costanti per Mercurio
  LOCAL ω,Ω, L, M, v;
  LOCAL ec, in, a, θ;

  ec:= 0.205615; // eccentricità
  in:= 7.003; // inclinazione
  a:= 0.387098; // semiasse maggiore
  θ:= 35.6;  // θs to calc if retrogade
  calcN();
  // costanti planetarie
  L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media
  ω:= →HMS((75.913 + 0.00004253*N) MOD 360); // perielio
  Ω:= →HMS((47.157 + 0.00003244*N) MOD 360); // nodo
  M:= (L - ω); // Anomalia vera
  v:= M + 23.4372*sin(M)+2.9810*sin(2*M)+0.5396*sin(3*M)+0.1098*sin(4*M); 
  planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
  RETURN planet_const;
END;

Venus()
BEGIN
  // Costanti per Venere
  LOCAL ω,Ω, L, M, v;
  LOCAL ec, in, a, θ;

  ec:= 0.006816; // eccentricità
  in:= 3.393636; // inclinazione
  a:= 0.72333; // semiasse maggiore
  θ:= 13.0;  // θs to calc if retrogade
  calcN();
  // costanti planetarie
  L:= →HMS((206.758 + 1.60216873*N) MOD 360); // Longitudine media
  ω:= →HMS((130.154 + 0.00003757*N) MOD 360); // perielio
  Ω:= →HMS((75.797 + 0.00002503*N) MOD 360); // nodo
  M:= (L - ω); // Anomalia vera
  v:= M + 0.7811*sin(M); 
  planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
  RETURN planet_const;
END;

Mars()
BEGIN
  // Costanti per Marte
  LOCAL ω,Ω, L, M, v;
  LOCAL ec, in, a, θ;

  ec:= 0.093309; // eccentricità
  in:= 1.850303; // inclinazione
  a:= 1.523678; // semiasse maggiore
  θ:= 16.8;  // θs to calc if retrogade
  calcN();
  // costanti planetarie
  L:= →HMS((124.765 + 0.52407108*N) MOD 360); // Longitudine media
  ω:= →HMS((334.251 + 0.00005038*N) MOD 360); // perielio
  Ω:= →HMS((48.794 + 0.00002127*N) MOD 360); // nodo
  M:= (L - ω); // Anomalia vera
  v:= M + 10.608*sin(M)+0.6216*sin(2*M)+0.0504*sin(3*M)+0.0047*sin(4*M); 
  planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
  RETURN planet_const;
END;

Jupiter()
BEGIN
  // Costanti per Giove
  LOCAL ω,Ω, L, M, v;
  LOCAL ec, in, a, θ;

  ec:= 0.048335; // eccentricità
  in:= 1.308736; // inclinazione
  a:= 5.202561; // semiasse maggiore
  θ:= 54.43;  // θs to calc if retrogade
  calcN();
  // costanti planetarie
  L:= →HMS((268.350 + 0.08312942*N) MOD 360); // Longitudine media
  ω:= →HMS((12.737 + 0.00004408*N) MOD 360); // perielio
  Ω:= →HMS((99.453 + 0.00002767*N) MOD 360); // nodo
  M:= (L - ω); // Anomalia vera
  v:= M + 5372*sin(M)+0.1662*sin(2*M)+0.007*sin(3*M); 
  planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
  RETURN planet_const;
END;

Saturn()
BEGIN
  // Costanti per Saturno
  LOCAL ω,Ω, L, M, v;
  LOCAL ec, in, a, θ;

  ec:= 0.055892; // eccentricità
  in:= 2.492519; // inclinazione
  a:= 9.554747; // semiasse maggiore
  θ:= 65.53;  // θs to calc if retrogade
  calcN();
  // costanti planetarie
  L:= →HMS((278.774 + 0.03349788*N) MOD 360); // Longitudine media
  ω:= →HMS((91.117 + 0.00005362*N) MOD 360); // perielio
  Ω:= →HMS((112.79 + 0.00002391*N) MOD 360); // nodo
  M:= (L - ω); // Anomalia vera
  v:= M + 6.4023*sin(M)+0.2235*sin(2*M)+0.011*sin(3*M); 
  planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
  RETURN planet_const;
END;

Uranus()
BEGIN
  // Costanti per Urano
  LOCAL ω,Ω, L, M, v;
  LOCAL ec, in, a, θ;

  ec:= 0.046344; // eccentricità
  in:= 0.772464; // inclinazione
  a:= 19.21814; // semiasse maggiore
  θ:= 73.92;  // θs to calc if retrogade
  calcN();
  // costanti planetarie
  L:= →HMS((248.487 + 0.01176902*N) MOD 360); // Longitudine media
  ω:= →HMS((171.563 + 0.00004064*N) MOD 360); // perielio
  Ω:= →HMS((73.482 + 0.00001365*N) MOD 360); // nodo
  M:= (L - ω); // Anomalia vera
  v:= M + 5.3092*sin(M)+0.1537*sin(2*M)+0.0062*sin(3*M); 
  planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
  RETURN planet_const;
END;

Neptunus()
BEGIN
  // Costanti per Nettuno
  LOCAL ω,Ω, L, M, v;
  LOCAL ec, in, a, θ;

  ec:= 0.008997; // eccentricità
  in:= 1.779242; // inclinazione
  a:= 30.10957; // semiasse maggiore
  θ:= 77.63;  // θs to calc if retrogade
  calcN();
  // costanti planetarie
  L:= →HMS((86.652 + 0.00602015*N) MOD 360); // Longitudine media
  ω:= →HMS((46.742 + 0.000039*N) MOD 360); // perielio
  Ω:= →HMS((130.693 + 0.00003009*N) MOD 360); // nodo
  M:= (L - ω); // Anomalia vera
  v:= M + 1.031*sin(M)+0.0058*sin(2*M); 
  planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
  RETURN planet_const;
END;

Pluto()
BEGIN
  // Costanti per Plutone
  LOCAL ω,Ω, L, M, v;
  LOCAL ec, in, a, θ;

  ec:= 0.250236; // eccentricità
  in:= 17.17047; // inclinazione
  a:= 39.438712; // semiasse maggiore
  θ:= 79.41;  // θs to calc if retrogade
  calcN();
  // costanti planetarie
  // Long, perielio, nodo, non precisi, considerati costanti
  // L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media
  M:= →HMS((230.67 + 0.00397943*N) MOD 360); // Anomalia vera (calcolata direttamente)
  Ω:= →HMS((109.06 + 0.00003823*N) MOD 360); // nodo
  ω:= →HMS((114.27 + Ω) MOD 360); // perielio
  L:= →HMS((M + ω) MOD 360); // L (non sicuro)
  
  v:= M + 28.45*sin(M)+4.38*sin(2*M)+0.97*sin(3*M)+0.24*sin(4*M); 
  planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};
  RETURN planet_const;
END;