Astronomy...

[salvomic]

some improvements (calculation of times, sidereal times, julian day, epsilon, data for Pluto...)

Code:


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

EXPORT Effemeridi()
// from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015
// Credits: from a book, Calcule astronomique pour amateurs, by Serge Bouiges
// Thanks Didier Lachieze, Marcel, 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 ascendent(); END;
  IF ch==6 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;
ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3;
ε:= →HMS(ε);
  sunCalc(); // calculate Sun data
END;

planets()
BEGIN
  LOCAL ch, nameplanet;
  INPUT({{ch,{"Mercurius","Venus","Mars", "Jupiter", "Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},
  "Choose planet", 
  "planet: ", "Choose the planet to calculate an press OK", 0,5 );
  CASE
  IF ch==1 THEN nameplanet:="Mercurius"; planetCalc(nameplanet); END;
  IF ch==2 THEN nameplanet:="Venus"; planetCalc(nameplanet); END;
  IF ch==3 THEN nameplanet:="Mars"; planetCalc(nameplanet); END;
  IF ch==4 THEN nameplanet:="Jupiter"; planetCalc(nameplanet); END;
  IF ch==5 THEN nameplanet:="Saturn"; planetCalc(nameplanet); END;
  IF ch==6 THEN nameplanet:="Uranus"; planetCalc(nameplanet); END;
  IF ch==7 THEN nameplanet:="Neptunus"; planetCalc(nameplanet); END;
  IF ch==8 THEN nameplanet:="Pluto"; planetCalc(nameplanet); END;
  DEFAULT
  END; // case
END;

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

calcTS()
BEGIN
  LOCAL TS, TSd, N0, TSL, TSapp;
  calcN();
  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 {trunc_angle(TS), trunc_angle(TSL), trunc_angle(TSapp), trunc_angle(TSd), gmt};
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;


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("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 gradi " + trunc_angle(TSid(4)));
  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)), gmt });
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;

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;
  h0:= -0°50′; // 'standard' altitude of Sun
  calcN();
  julianDay();
  L:= →HMS((278.965+0.985647342*N)) MOD 360; // long media
  ω:= →HMS((281.235 + 0.0000469*N)) MOD 360; // long del perielio
  M:= (L - ω); // anomalia vera
  v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio
  l:= v + ω; // Longitudine Sole
  lt:= l + 180; // Longitudine eliocentrica della Terra
  r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)
  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))) );
  // α:= →HMS((acos(cos(l)/cos(δ))/15) MOD 24); // ascensione retta
  α:= (acos(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
// temp:= sin(h0)-sin(lat)*sin(δ)/cos(lat)*cos(δ);
// IF temp<0 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END;
// H0:= acos(temp);
// MSGBOX("H0: " + H0);
  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
  // 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;
  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
  sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)
  tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)

//temp:= ((α*15 + long - ts(3)*15) / 360);
// IF temp<0 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END;
// culm:= (temp - Hd/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};
  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;
  L:= sunlist(1); ω:= sunlist(2); 
  l:= sunlist(3); M:= sunlist(4); 
  α:= 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); 
  ts:= calcTS(); //
  jd:= julianDay();
  DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);
  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 vera " + trunc_angle(M));
  PRINT("λ longitudine Sole " + trunc_angle(l) + " " + segno);
  PRINT("r Distanza (UA) " + r);
  PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
  PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
  PRINT(" ");
  PRINT("Ang. orario "+ 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("Equazione del tempo (min) " + trunc_angle(E));
  PRINT("Julian day " + jd(1));
  PRINT("Dif asc." + 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("Long elioc. Terra "+ trunc_angle(lt));
  PRINT("ε epsilon " + ε);
  // RETURN {L, ω, M, v, l, r, x, y};
  RETURN [[L, ω], [l, M], [α,δ], [az, h], [x,y], [sorg, tram], [culm, r], [H,ts(1)], [lt, ε]];
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 dd, hh, mm, segno;
  i:= 5.13333333; // inclinazione 5°8' sull'eclittica
  calcN();
  // Sun
  Ls:= sunlist(1); // Longitudine media
  Ms:= sunlist(4); // Anomalia
  // end Sun

  L:= (33.231 + 13.17639653*N) MOD 360; // long. media
  Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente
  M:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera
  ω:= (L - M) MOD 360; // perigeo medio
  DD:= L - Ls; F:= L-Ω; // for the variations

  LongM:= L + 6.28875 * sin(M); // calculus (h/ and above) of Moon longitude
  LongM:= LongM + 0.2136 * sin(2*M); // equazione del centro
  LongM:= LongM + 0.6583 * sin(2*DD); // variazione
  LongM:= LongM - 0.1856 * sin(Ms); // equazione annuale
  LongM:= LongM + 1.2740 * sin(2*DD-M); // evezione
  LongM:= LongM - 0.1143 * sin(2*F); // riduzione all'eclittica
  LongM:= LongM + 0.0588 * sin(2*DD-2*M);
  LongM:= LongM + 0.0572 * sin(2*DD-M-Ms);
  LongM:= LongM + 0.0533 * sin(2*DD+M);
  LongM:= LongM + 0.0459 * sin(2*DD-Ms);
  LongM:= LongM + 0.0410 * sin(M-Ms);
  LongM:= LongM - 0.0305 * sin(M+Ms);
  LongM:= LongM - 0.0348 * sin(DD);
  LongM:= LongM MOD 360;
  LatM := 5.1280 * sin(F); // termine principale
  LatM := LatM + 0.2806 * sin(M+F);
  LatM := LatM + 0.2777 * sin(M-F); // equazione del centro
  LatM := LatM + 0.1732 * sin(2*DD-F); // grande ineguaglianza
  LatM := LatM + 0.0554 * sin(2*DD-M+F); // evezione
  LatM := LatM + 0.0462 * sin(2*DD-M-F); // evezione (2)
  LatM := LatM + 0.0326 * sin(2*DD+F);
  LatM := LatM MOD 360;

b:= LatM;
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)
  // α:= →HMS((acos(cos(b)*cos(LongM)/cos(δ))/15) MOD 24); // ascensione retta
  α:= (acos(cos(b)*cos(IFTE(LongM>180, LongM-180, LongM))/cos(δ))); // ascensione retta
  IF LongM>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
  // 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)

  P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse
  // Raggio quadratico medio della Terra: 6373,044737
  d:= 1/sin(P); // distanza Terra-Luna in raggi terrestri
  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
  dd:= IP(HMS→(ph/15));
  hh:= IP(FP(HMS→(ph/15))*60);
  IF (hh>24) THEN dd:= dd+1; hh:= hh MOD 24; END;
  mm:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);
  eta:= →HMS(dd+hh/60+mm/3600); // Age of the Moon
  segno:= zodiaco(LongM);

  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("ω Perigeo medio " + trunc_angle(ω));
  PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
  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("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
  PRINT("P Parallasse " + trunc_angle(P));
  PRINT("d Distanza Terra - Luna "+ ROUND(d*6373.044737,3) + " km");
  PRINT("Semidiametro "+  trunc_angle(→HMS(s)) + " :: " + ROUND(s*6373.044737,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);
  RETURN [[L, M],[LongM, LatM],[ω,Ω],[α,δ], [az, h], [P,0],[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  β
  // 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;