Astronomy...

[salvomic]

New version, greater precision (especially for the Moon) and rationalization...

Code:


data();
calcN();
calcTS();
julianDay();
julianDayAtZero();
tempo();
precession();
parametri();
sunCalc();
planetCalc();
planets();
sun();
moon();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendent();
zodiaco();
trunc_angle();
simp_angle();
atan2();
thisday:= 2000.0101;
lstd:=0; gmt:= 12;
long:= 0; lat:=  0;
m:= 1; dst:= 0; utc:=0;
deltaT:= 68; // 2015
datetimelist:= MAKELIST(X,X,1,6);
sunlist:= MAKELIST(X,X,1,25);
planet_const:= MAKELIST(X,X,1,10);
ε:= →HMS(23.4392911);  // intial value, approx 23.45 (23°26'21")

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

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

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

  jd:= julianDay();
  t:=(jd(1)-2451545)/36525;
  U:= t/100; 
  // ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
  ε:= 23°26′21.448″ - 1°18′00.93″*U-1.55*U^2+1999.25*U^3-51.38*U^4-249.67*U^5;
  ε:= ε-39.05*U^6+7.12*U^7+27.87*U^8+5.79*U^9+2.45*U^10; // Laskar, Meeus
  ε:= →HMS(ε);
  sunCalc(); // calculate Sun data
END;

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

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

calcTS()
BEGIN
  LOCAL T, TS, TSd, N0, TSL, TSapp;
  LOCAL θ0, θ00, jd, jd0;
  calcN();
  jd:= julianDay();
  jd0:= julianDayAtZero();
  T:= (jd(1)-2451545.0)/36525;
  θ0:= 280.46061837+360.98564736629*(jd(1)-2451545.0)+0.000387933*T^2-T^3/38710000; // DEG at 0h
  θ0:= θ0 MOD 360;
  θ00:= 280.46061837+360.98564736629*(jd0(1)-2451545.0)+0.000387933*T^2-T^3/38710000; // DEG at 0h
  θ00:= θ00 MOD 360; 
  N0:= DDAYS(1901.0101, thisday)+1+0; // N at 0 GMT (for TS)
  TS:=  ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds
  TS:= →HMS(TS); // in HMS 
  TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST
  TSL:= →HMS(TSL); // in HMS
  TSapp:= (TS + gmt) MOD 24; // apparent ST (ST at Greenwitch at local our)
  TSapp:= →HMS(TSapp); // in HMS

  TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees
  RETURN {TS, TSL, TSapp, TSd, θ0, θ00, T};
END;

julianDay()
BEGIN
  LOCAL yy, mm, dd,dde, hh, mn, ss;
  LOCAL aa,bb, jd, jde;
  yy:= Y; mm:= m;
  hh:= datetimelist(4);
  mn:= datetimelist(5);
  ss:= datetimelist(6);
  dd:= D + hh/24 + mn/(60*24) + ss/(60*60*24);
  dde:= D + hh/24 + mn/(60*24) + (ss+deltaT)/(60*60*24);
  IF (m=1 OR m=2) THEN yy:=Y-1; mm:=m+12;END;
  IF yy*100+mm+dd>=158225 THEN
  aa:=IP(ABS(yy/100));
  bb:=2-aa+IP(aa/4);
  ELSE bb:=0; END;
  jd:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dd+bb-1524.5;
  jde:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dde+bb-1524.5;
  RETURN {jd, jde};
END;

julianDayAtZero()
BEGIN
  LOCAL yy, mm, dd,dde, hh, mn, ss;
  LOCAL aa,bb, jd, jde;
  yy:= Y; mm:= m;
  hh:= 0;
  mn:= 0;
  ss:= 0;
  dd:= D;
  dde:= D + (ss+deltaT)/(60*60*24);
  IF (m=1 OR m=2) THEN yy:=Y-1; mm:=m+12;END;
  IF yy*100+mm+dd>=158225 THEN
  aa:=IP(ABS(yy/100));
  bb:=2-aa+IP(aa/4);
  ELSE bb:=0; END;
  jd:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dd+bb-1524.5;
  jde:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dde+bb-1524.5;
  RETURN {jd, jde};
END;


tempo()
BEGIN
  LOCAL n, jd, TSid;
  n:= calcN();
  jd:= julianDay();
  TSid:= calcTS();
  PRINT;
  PRINT("Date " + m + " " + D + " " + Y);
  PRINT("Julian Day " + jd(1));
  PRINT("Julian Day Effem. " + jd(2));
  PRINT("N (from 1901jan1) " + n);
  PRINT("T (from JD) " + TSid(7));
  PRINT("Tempo Siderale 0h " + trunc_angle(TSid(1)));
  PRINT("Tempo Siderale locale " + trunc_angle(TSid(2)));
  PRINT("Tempo Siderale apparente " + trunc_angle(TSid(3)));
  PRINT("Tempo Siderale (DEG) " + trunc_angle(TSid(4)));
  PRINT("θ0 Mean ST 0h " + trunc_angle(TSid(5)/15));
  PRINT("θ0 Mean ST 0h (DEG) " + trunc_angle(TSid(5)));
  PRINT("θ00 MST 00h GMT " + trunc_angle(TSid(6)/15));
  PRINT("θ00 MST 00h GMT (DEG) " + trunc_angle(TSid(6)));
  PRINT("Time " + lstd + " GMT " + gmt);
  RETURN({ROUND(jd(1),5), ROUND(jd(2),5), ROUND(n, 5), trunc_angle(TSid(1)), trunc_angle(TSid(2)), 
    trunc_angle(TSid(3)), trunc_angle(TSid(4)), trunc_angle(TSid(5)) });
END;

deltaEpsilonPsi()
BEGIN
  // Nutazione Δψ (longitudine), Δε (obliquità)
  LOCAL L1, deltaPsi, deltaEpsilon, Ω, jd;
  jd:= julianDayAtZero();
  T:= (jd(2)-2451545)/36525; // use JDE
  Ω:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)
  L:= (→HMS(280.4665)+36000.76983*T+→HMS(0.0003032*T^2)) MOD 360; // Long media Sole
  L1:= →HMS(218.3164591) + →HMS(481267.88134236)*T-0.0013268*T^2; //  Long media Luna
  L1:= L1+(T^3)/538841-(T^4)/65194000;
  L1:= L1 MOD 360;
  deltaPsi:= -0°0′17.20″*sin(Ω)-0°0′1.32″*sin(2*L)-0°0′0.23″*sin(2*L1)+0°0′0.21″*sin(2*Ω);
  deltaEpsilon:= 0°0′9.2″*cos(Ω)+0°0′0.57″*cos(2*L)+0°0′0.10″*cos(2*L1)-0°0′0.09″*cos(2*Ω);
  RETURN({deltaPsi, deltaEpsilon});   
END;

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

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

ascendent()
BEGIN
  LOCAL ASC, ts, segno;
  ts:= calcTS();
  ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));
  IF ASC<0 THEN ASC:= ASC+360; END; IF ASC>360 THEN ASC:=ASC-360; END;
  segno:= zodiaco(ASC);
RETURN {→HMS(ASC), segno};
END;

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

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

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

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

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

sunCalc()
BEGIN
  LOCAL ω, l, v, r, x, y, dz;
  LOCAL α, δ, dayY, H, ts, az, z, h;
  LOCAL SAD, DA, sorg, tram, culm, Hd, lt;
  LOCAL as, ac, zc, h0, H0, temp;
  LOCAL ec, a, u, M1, M2, j, t;
  LOCAL m0, m1, m2, jd, c, lapp, Ω;
  calcN();
  jd:= julianDay();
  a:= 1.00000023; // semiasse maggiore
  t:= N/36525;
  T:= (jd(2)-2451545)/36525; // use JDE
  // ec:= 0.01675062-0.0000418*t-0.000000137*t^2; // eccentricità (RAD)
  ec:= 0.016708617-0.000042037*T-0.0000001236*T^2; // eccentricità Meeus
  // ec:= ec*180/PI; // DEG
  h0:= -0°50′; // 'standard' altitude of Sun
  // L:= →HMS(278.9654+0.985647342*N + 0.0003*t^2) MOD 360; // long media
  L:= (→HMS(280.4665)+36000.76983*T+→HMS(0.0003032*T^2)) MOD 360; // Long media Sole Meeus
  ω:= →HMS(281.238 + 0.000047067*N + 0.00045*t^2) MOD 360; // long del perielio
  // M:= (L - ω); // anomalia vera
  //IF M<0 THEN M:= (360+M); END; // potrebbe essere negativa
  Ω:= (→HMS(125.04)-→HMS(1934.136)*T) MOD 360; // Nutazione e aberrazione (nodo)
  M:= →HMS(357.52910)+→HMS(35999.05030)*T+→HMS(0.0001559)*T^2-→HMS(0.00000048)*T^3; // anom vera Meeus
  M:= M MOD 360;
  u:=atan2(sin(M),(cos(M)-ec)); // only for small eccentricity
  c:= (→HMS(1.914600)-→HMS(0.004817)*T-→HMS(0.000016)*T^2)*sin(M);  // equation of center
  c:= c+(→HMS(0.019993)-→HMS(0.000101)*T)*sin(2*M)+(→HMS(0.000290))*sin(3*M);
  // v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio (Bouiges)
  // v:= acos((cos(u)-ec)/(1-ec*cos(u))); // kepler
  v:= M + c; // true anomaly, Meeus
  // l:= (v + ω) MOD 360; // Longitudine Sole
  l:= (L+c) MOD 360; // Ө Longitudine Sole, Meeus
  lt:= l + 180; // Longitudine eliocentrica della Terra
  lapp:= l - →HMS(0.00569)-→HMS(0.00478)*sin(Ω); // Longitudine apparente (ref true equinox)
  lapp:= simp_angle(lapp);
  // r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)
  r:= (1.000001018*(1-ec^2))/(1+ec*cos(v)); // Distanza Sole-Terra (UA) Meeus
  x:= r * cos(l); y:= r * sin(l);
  dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number
  // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper
  δ:= →HMS((asin(sin(ε)*sin(l))) );
  // α:= (acos(cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta
  // IF l>180 THEN α:= α +180; END;
  α:= →HMS(atan2(cos(ε)*sin(l),cos(l))); // Ascensione retta, Meeus (b never > 1.2", l=Ө)
  α:= →HMS(α/15) MOD 24; // α in hms
  ts:= calcTS();  // calcola tempo siderale
  H:= (ts(2) - α) MOD 24; // Angolo orario
  Hd:= →HMS(15*H); // Angolo orario in gradi
  zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt
  dz:= →HMS(acos(zc)) ; // distanza zenitale
  // alt:= →HMS(90 - z); // Altezza
  h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza
  az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd)); // Azimuth
  IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;
  IF (az<0) THEN az:= 360+az; END;
  E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)
  DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
  SAD:= (90 + DA); // Semiarco Diurno
  //   SAD:= →HMS(acos(-tan(lat)*tan(δ))); // Semiarco Diurno
 culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione, transito
 sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
 tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)

// temp:= (sin(h0)-sin(lat)*sin(δ))/(cos(lat)*cos(δ));
//IF temp<-1 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END;
//H0:= acos(temp);
// m0:= ((α*15 + long - ts(6)) / 360);
//IF m0<0 THEN m0:= m0 +1; END; IF m0>1 THEN m0:=m0-1; END;
// culm:= m0*24;


//temp:= ts(3)*15+360.985647*m0;
// H1:= temp - long - α*15;
//        IF temp>360 THEN temp:= temp MOD 360; END;
//culm:= (m0 - H/360) * 24 ;

// sorg:= (culm - H0/360) *24;
// tram:= (culm + H0/360) *24;

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


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

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

moon()
BEGIN
  LOCAL ω,Ω, l, v, r, x, y, i;
  LOCAL Ls, Ms, LongM, LatM, P, d, s;
  LOCAL α,δ, ts, Hd, zc, dz, h, az;
  LOCAL as, ac, DA, SAD, DD, culm, sorg, tram;
  LOCAL b, dayY, xs, ys, ph, fase, eta;
  LOCAL d1, h1, m1, segno, jd, rr;
  LOCAL lapp, dep, ill, ee, a1,a2,a3;
  i:= 5.13333333; // inclinazione 5°8' sull'eclittica
  calcN();
  jd:= julianDay();
  T:= (jd(2)-2451545)/36525; // use JDE
  dep:= deltaEpsilonPsi();
  // Sun
  Ls:= sunlist(1); // Longitudine media
  Ms:= sunlist(4); // Anomalia
  // end Sun

  //L:= (33.231 + 13.17639653*N) MOD 360; // long. media
  L:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+(T^3)/538841-(T^4)/65194000);
  //  Long media Luna Meeus

  // Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente
  Ω:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)
  Ω:= simp_angle(Ω);
  // M:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera
  M:= 134.9634114+477198.8676313*T+0.0089970*T^2+(T^3)/69699-(T^4)/14712000; // anomalia media
  M:= simp_angle(M);
  ω:= (L - M) MOD 360; // perigeo medio
  // DD:= L - Ls; F:= L-Ω; // for the variations
  DD:= 297.8502042+445267.1115168*T-0.0016300*T^2+(T^3)/545868-(T^4)/113065000; // mean elongation
  DD:= simp_angle(DD);
  F:= 93.2720993+483202.0175273*T-0.0034029*T^2-(T^3)/3526000+(T^4)/863310000; // dist from asc. node
  F:= simp_angle(F);
  ee:= 1-0.002516*T-0.0000074*T^2;  // multiply those w/ Ms *ee
  a1:= →HMS(119.75)+→HMS(131.849)*T; // action of Venus
  a2:= →HMS(53.09)+→HMS(479264.290)*T; // action of Jupiter
  a3:= →HMS(313.45)+→HMS(481266.484)*T; 

  LongM:= L     + 6.288774 * sin(M); // calc Moon longitude
  LongM:= LongM + 1.274027 * sin(2*DD-M); // evezione
  LongM:= LongM + 0.658314 * sin(2*DD); // variazione
  LongM:= LongM + 0.213618 * sin(2*M); // equazione del centro
  LongM:= LongM - 0.185116 * sin(Ms) *ee; // equazione annuale
  LongM:= LongM - 0.114332 * sin(2*F); // riduzione all'eclittica
  LongM:= LongM + 0.058793 * sin(2*DD-2*M);
  LongM:= LongM + 0.057066 * sin(2*DD-M-Ms) *ee;
  LongM:= LongM + 0.053322 * sin(2*DD+M);
  LongM:= LongM + 0.045758 * sin(2*DD-Ms) *ee;
  LongM:= LongM - 0.040923 * sin(Ms-M) *ee;
  LongM:= LongM - 0.034720 * sin(DD);
  LongM:= LongM - 0.030383 * sin(M+Ms) *ee;
  LongM:= LongM + 0.015327 * sin(2*DD-2*F);
  LongM:= LongM - 0.012528 * sin(M+2*F);
  LongM:= LongM + 0.010980 * sin(M-2*F);
  LongM:= LongM + 0.010980 * sin(M-2*F);
  LongM:= LongM + 0.010675 * sin(4*DD-M);
  LongM:= LongM + 0.010034 * sin(3*M);
  LongM:= LongM + 0.008548 * sin(4*DD-2*M);
  LongM:= LongM - 0.007888 * sin(2*DD+Ms-M) *ee;
  LongM:= LongM - 0.006766 * sin(2*DD+Ms) *ee;
  LongM:= LongM - 0.005163 * sin(DD-M);
  LongM:= LongM + 0.004987 * sin(DD+Ms) *ee;
  LongM:= LongM + 0.004036 * sin(2*DD-Ms+M);
  LongM:= LongM + 0.003994 * sin(2*DD+2*M);
  LongM:= LongM + 0.003861 * sin(4*DD);
  LongM:= LongM + 0.003665 * sin(2*DD-3*M);
  LongM:= LongM - 0.002689 * sin(Ms-2*M) *ee;
  LongM:= LongM - 0.002602 * sin(2*DD-M+2*F);
  LongM:= LongM + 0.002390 * sin(2*DD-Ms-2*M) *ee;
  LongM:= LongM - 0.002348 * sin(DD+M);
  LongM:= LongM + 0.002236 * sin(2*DD-2*Ms) *ee^2;
  LongM:= LongM - 0.002120 * sin(Ms+2*M) *ee;
  LongM:= LongM - 0.002069 * sin(2*Ms) *ee^2;
  LongM:= LongM + 0.002048 * sin(2*DD-2*Ms-M) *ee^2;
  LongM:= LongM - 0.001773 * sin(2*DD+M-2*F);
  LongM:= LongM - 0.001595 * sin(2*DD+2*F);
  LongM:= LongM + 0.001215 * sin(4*DD-Ms-M) *ee;
  LongM:= LongM - 0.001110 * sin(2*M+2*F);
  LongM:= LongM - 0.000892 * sin(3*DD-M);
  LongM:= LongM - 0.000810 * sin(2*DD+Ms+M) *ee;
  LongM:= LongM + 0.000759 * sin(4*DD-Ms-2*M) *ee;
  LongM:= LongM - 0.000713 * sin(2*Ms-M) * ee^2;
  LongM:= LongM - 0.000700 * sin(2*DD+2*Ms-M) *ee^2;
  LongM:= LongM + 0.000691 * sin(2*DD+Ms-2*M) *ee;
  LongM:= LongM + 0.000596 * sin(2*DD-Ms-2*F) *ee;
  LongM:= LongM + 0.000549 * sin(4*DD+M);
  LongM:= LongM + 0.000537 * sin(4*M);
  LongM:= LongM + 0.000520 * sin(4*DD-Ms) *ee;
  LongM:= LongM - 0.000487 * sin(DD-2*M);
  LongM:= LongM - 0.000399 * sin(2*DD+Ms-2*F) *ee;
  LongM:= LongM - 0.000381 * sin(2*M-2*F);
  LongM:= LongM + 0.000351 * sin(DD+Ms+M) *ee;
  LongM:= LongM - 0.000340 * sin(3*DD-2*M);
  LongM:= LongM + 0.000330 * sin(4*DD-3*M);
  LongM:= LongM + 0.000327 * sin(2*DD-Ms+2*M) *ee;
  LongM:= LongM - 0.000323 * sin(2*Ms+M) *ee^2;
  LongM:= LongM + 0.000299 * sin(DD+Ms-M) *ee;
  LongM:= LongM + 0.000294 * sin(2*DD+3*M);

  LongM:= LongM + (3958*sin(a1) + 1962*sin(L-F)+318*sin(a2))/1000000;

  LongM:= simp_angle(LongM);

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

  LatM :=        5.128122 * sin(F); // termine principale
  LatM := LatM + 0.280026 * sin(M+F);
  LatM := LatM + 0.277693 * sin(M-F); // equazione del centro
  LatM := LatM + 0.173237 * sin(2*DD-F); // grande ineguaglianza
  LatM := LatM + 0.055413 * sin(2*DD-M+F); // evezione
  LatM := LatM + 0.046271 * sin(2*DD-M-F); // evezione (2)
  LatM := LatM + 0.032573 * sin(2*DD+F);
  LatM := LatM + 0.017198 * sin(2*M+F);
  LatM := LatM + 0.009266 * sin(2*D+M-F);
  LatM := LatM + 0.008822 * sin(2*M-F);
  LatM := LatM + 0.008216 * sin(2*DD-Ms-F) *ee;
  LatM := LatM + 0.004324 * sin(2*DD-2*M-F);
  LatM := LatM + 0.004200 * sin(2*DD+M-F);
  LatM := LatM - 0.003359 * sin(2*DD+Ms-F) *ee;
  LatM := LatM + 0.002463 * sin(2*DD-Ms-M+F) *ee;
  LatM := LatM + 0.002211 * sin(2*DD-Ms+F) *ee;
  LatM := LatM + 0.002065 * sin(2*DD-Ms-M-F) *ee;
  LatM := LatM - 0.001870 * sin(Ms-M-F) *ee;
  LatM := LatM + 0.001828 * sin(4*DD-M-F);
  LatM := LatM - 0.001794 * sin(Ms-F) *ee;
  LatM := LatM - 0.001749 * sin(3*F);
  LatM := LatM - 0.001565 * sin(Ms-M+F) *ee;
  LatM := LatM - 0.001491 * sin(DD+F);
  LatM := LatM - 0.001475 * sin(Ms+M+F) *ee;
  LatM := LatM - 0.001410 * sin(Ms+M-F) *ee;
  LatM := LatM - 0.001344 * sin(Ms-F) *ee;
  LatM := LatM - 0.001335 * sin(DD-F);
  LatM := LatM + 0.001107 * sin(3*M+F);
  LatM := LatM + 0.001021 * sin(4*DD-F);
  LatM := LatM + 0.000833 * sin(4*DD-M+F);
  LatM := LatM + 0.000777 * sin(M-3*F);
  LatM := LatM + 0.000671 * sin(4*DD-2*M+F);
  LatM := LatM + 0.000607 * sin(2*DD-3*F);
  LatM := LatM + 0.000596 * sin(2*DD-F);
  LatM := LatM + 0.000491 * sin(2*DD-M+F);
  LatM := LatM - 0.000451 * sin(2*DD-M-F);
  LatM := LatM + 0.000439 * sin(2*DD+F);
  LatM := LatM + 0.000422 * sin(2*M-F);
  LatM := LatM + 0.000421 * sin(2*DD+M-F);
  LatM := LatM - 0.000366 * sin(2*M-F);
  LatM := LatM - 0.000351 * sin(2*DD+Ms-F) *ee;
  LatM := LatM + 0.000331 * sin(4*DD+F);
  LatM := LatM + 0.000315 * sin(2*DD-Ms+M+F) *ee;
  LatM := LatM + 0.000302 * sin(2*DD-2*Ms-F) *ee^2;
  LatM := LatM - 0.000283 * sin(M+3*F);
  LatM := LatM - 0.000229 * sin(2*DD+Ms+M-F) *ee;
  LatM := LatM + 0.000223 * sin(DD+Ms-F) *ee;
  LatM := LatM + 0.000223 * sin(DD+Ms+F) *ee;
  LatM := LatM - 0.000220 * sin(Ms-2*M-F) *ee;
  LatM := LatM - 0.000220 * sin(2*DD+Ms-M-F) *ee;
  LatM := LatM - 0.000185 * sin(DD+M+F);
  LatM := LatM + 0.000181 * sin(2*DD-Ms-2*M-F) *ee;
  LatM := LatM - 0.000177 * sin(Ms+2*M+F);
  LatM := LatM + 0.000176 * sin(4*DD-2*M-F);
  LatM := LatM + 0.000166 * sin(4*DD-Ms-M+F) *ee;
  LatM := LatM - 0.000164 * sin(DD+M-F);
  LatM := LatM + 0.000132 * sin(4*DD+M-F);
  LatM := LatM - 0.000119 * sin(DD-M-F);
  LatM := LatM + 0.000115 * sin(4*DD-Ms-F) *ee;
  LatM := LatM + 0.000107 * sin(2*DD-2*Ms+F) *ee^2;

  LatM := LatM + (0-2235*sin(L)+382*sin(a3)+175*sin(a1-F)
    +175*sin(a1+F)+127*sin(L-M)-115*sin(L+M)) / 1000000;

  // correzioni per la distanza dalla Terra (Δ)
  rr:=0;
  rr:= rr -20.905355 * cos(M); // coseni
  rr:= rr - 3.699111 * cos(2*DD-M); // evezione
  rr:= rr - 2.955968 * cos(2*DD); // variazione
  rr:= rr - 0.569925 * cos(2*M); // equazione del centro
  rr:= rr + 0.048888 * cos(Ms) *ee; // equazione annuale
  rr:= rr - 0.003149 * cos(2*F); // riduzione all'eclittica
  rr:= rr + 0.246158 * cos(2*DD-2*M);
  rr:= rr - 0.152138 * cos(2*DD-M-Ms) *ee;
  rr:= rr - 0.170733 * cos(2*DD+M);
  rr:= rr - 0.204586 * cos(2*DD-Ms) *ee;
  rr:= rr - 0.129620 * cos(Ms-M);
  rr:= rr + 0.108743 * cos(DD);
  rr:= rr + 0.104755 * cos(M+Ms) *ee;
  rr:= rr + 0.010321 * cos(2*DD-2*F);
  rr:= rr - 0.000000 * cos(M+2*F); // null
  rr:= rr + 0.079661 * cos(M-2*F);
  rr:= rr - 0.034782 * cos(4*DD-M);
  rr:= rr - 0.023210 * cos(3*M);
  rr:= rr - 0.021636 * cos(4*DD-2*M);
  rr:= rr + 0.024208 * cos(2*DD+Ms-M) *ee;
  rr:= rr + 0.030824 * cos(2*DD+Ms) *ee;
  rr:= rr - 0.008379 * cos(DD-M);
  rr:= rr - 0.016675 * cos(DD+Ms) *ee;
  rr:= rr - 0.012831 * cos(2*DD-Ms+M);
  rr:= rr - 0.010445 * cos(2*DD+2*M);
  rr:= rr - 0.011650 * cos(4*DD);
  rr:= rr + 0.014403 * cos(2*DD-3*M);
  rr:= rr - 0.007003 * cos(Ms-2*M) *ee;
  rr:= rr - 0.000000 * cos(2*DD-M+2*F); // null
  rr:= rr + 0.010056 * cos(2*DD-Ms-2*M) * ee;
  rr:= rr + 0.006322 * cos(DD+M);
  rr:= rr - 0.009884 * cos(2*DD-2*Ms) * ee;
  rr:= rr + 0.005751 * cos(Ms+2*M) *ee;
  rr:= rr - 0.000000 * cos(2*Ms) *ee^2; // null
  rr:= rr - 0.004950 * cos(2*DD-2*Ms-M) *ee^2;
  rr:= rr + 0.004130 * cos(2*DD+M-2*F);
  rr:= rr - 0.000000 * cos(2*DD+2*F); // null
  rr:= rr - 0.003958 * cos(4*DD-Ms-M) *ee;
  rr:= rr - 0.000000 * cos(2*M+2*F); // null
  rr:= rr + 0.003258 * cos(3*DD-M);
  rr:= rr + 0.002616 * cos(2*DD+Ms+M) *ee;
  rr:= rr - 0.001897 * cos(4*DD-Ms-2*M) *ee;
  rr:= rr - 0.002117 * cos(2*Ms-M) * ee^2;
  rr:= rr + 0.002354 * cos(2*DD+2*Ms-M) *ee^2;
  rr:= rr + 0.000000 * cos(2*DD+Ms-2*M) *ee; // null
  rr:= rr + 0.000000 * cos(2*DD-Ms-2*F) *ee; // null
  rr:= rr - 0.001423 * cos(4*DD+M);
  rr:= rr - 0.001117 * cos(4*M);
  rr:= rr - 0.001571 * cos(4*DD-Ms) *ee;
  rr:= rr - 0.001739 * cos(DD-2*M);
  rr:= rr - 0.000000 * cos(2*DD+Ms-2*F) *ee; // null
  rr:= rr - 0.004421 * cos(2*M-2*F);
  rr:= rr + 0.000000 * cos(DD+Ms+M) *ee; // null
  rr:= rr - 0.000000 * cos(3*DD-2*M); // null
  rr:= rr + 0.000000 * cos(4*DD-3*M); // null
  rr:= rr + 0.000000 * cos(2*DD-Ms+2*M) *ee; // null
  rr:= rr + 0.001165 * cos(2*Ms+M) *ee^2;
  rr:= rr + 0.000000 * cos(DD+Ms-M) *ee; // null
  rr:= rr + 0.000000 * cos(2*DD+3*M); // null
  rr:= rr + 0.008752 * cos(2*DD-M-2*F); // no sin

  b:= LatM; // Latitudine
  // dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number
  δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(LongM))) ); // declinazione (coord locali)
  //  α:= (acos(cos(b)*cos(IFTE(LongM>180, LongM-180, LongM))/cos(δ))); // ascensione retta
  //IF LongM>180 THEN α:= α +180; END;
  α:= atan2(sin(LongM)*cos(ε)-tan(b)*sin(ε), cos(LongM));
  α:= →HMS(α/15) MOD 24; // α in hms
  ts:= calcTS();  // calcola tempo siderale
  H:= (ts(2) - α) MOD 24; // Angolo orario
  Hd:= →HMS(15*H); // Angolo orario in gradi
  zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt
  dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza
  h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)
  as:= cos(δ)*sin(Hd)/sin(dz); // calc az
  ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);
  // az:= →HMS(ATAN(as/ac) ); // Azimuth
  // IF (ac<0) THEN az:= ((360 + az) MOD 360); END;
  az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));
  IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;
  DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
  SAD:= (90 + DA); // Semiarco Diurno
  // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
  culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione
  sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
  tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)

  // Raggio quadratico medio della Terra: 6373,044737
  d:= 385000.56 + rr*1000; // Δ distanza Terra-Luna in km (=1/sin(P)),rr in 10^-3
  //P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse
  P:= asin(6378.14/d); // Parallasse
  s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)
  ph:= abs(LongM - Ls); // Phase
  CASE
  // IF ph == 0 THEN fase:="New Moon"; END;
  IF ph>280 AND ph<350 THEN fase:="Decrescente"; END;
  IF ph>=260 AND ph<=280 THEN fase:="Last Quart"; END;
  IF ph>190 AND ph<260 THEN fase:="Gibbosa Cal."; END;
  IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
  IF ph>100 AND ph<170 THEN fase:="Gibbosa Cresc."; END;
  IF ph>=80 AND ph<=100 THEN fase:="1st Quart"; END;
  IF ph>10 AND ph<80 THEN fase:="Crescente"; END;
  DEFAULT fase:="New Moon";
  END; // case
  // età della Luna
  d1:= IP(HMS→(ph/15));
  h1:= IP(FP(HMS→(ph/15))*60);
  IF (h1>24) THEN d1:= d1+1; h1:= h1 MOD 24; END;
  m1:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);
  eta:= →HMS(d1+h1/60+m1/3600); // Age of the Moon
  segno:= zodiaco(LongM);
  ill:= (1+cos(ph))/2; // illuminated fraction

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

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

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

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

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

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

  dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number
  // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper
  δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(l))) ); // declinazione (coord locali)
  α:= (acos(cos(b)*cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta
  IF l>180 THEN α:= α +180; END;
  α:= →HMS(α/15) MOD 24; // α in hms
  ts:= calcTS();  // calcola tempo siderale
  H:= (ts(2) - α) MOD 24; // Angolo orario
  Hd:= →HMS(15*H); // Angolo orario in gradi
  zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt
  dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza
  h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)
  as:= cos(δ)*sin(Hd)/sin(dz); // calc az
  ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);
  // az:= →HMS(ATAN(as/ac) ); // Azimuth
  az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));
  IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;
  DA:= →HMS(asin(tan(lat)*tan(δ))); // Differenza ascensionale
  SAD:= →HMS(90 + DA); // Semiarco Diurno
  // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
  culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) MOD 24 ; // Culminazione
  sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
  tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)
  IF abs(l-220.4)>θ THEN mov:="diretto"; ELSE mov:="retrogrado"; END;
  segno:= zodiaco(l);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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