Astronomy...

[salvomic]

hi everybody,
waiting for the evolution of the very interesting app by Marcel (I'm already using!), I'm developing also this program to get astronomical results for Sun, Moon and Planets...
Effemeridi (Astrolabio) is a program that I'm adapted from an old version I made many years ago for a TI pocket calculator.
The program has still some incorrect output (very few), so I'm glad to have your advice and hints to improve it, thank you.
The program calculates Sidereal Time, Longitude, Latitude (various coordinates), Rect Ascension and Declination, Azimuth and Height, phases of the Moon, Retrogradation of planets, Ascendent, after a book of S. Bouiges...
I'm working to add eclipses, and improving program also in "astrological" sense (I was asked for this) and Satellites of Jupiter and Saturn.
Also data for Pluto are still not so precise as those in the app of Marcel...

Some routines and data are still in Italian, sorry.
I need your help, please, test, and thank you in advance!

Salvo

The code (I'll post here the new versions for now):

Code:


data();
calcN();
calcTS();
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;
sunlist:= MAKELIST(X,X,1,25);
planet_const:= MAKELIST(X,X,1,10);
ε:= →HMS(23.4392911);  // approx 23.45

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 N", "Calc TS", "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 calcN(); END;
  IF ch==5 THEN calcTS(); 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;
  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}} },
  "Data: Use Shift+9 for H°M′S″",
  {"Month:","Day :","Year :","Local Time (24):", "DST", "GMT", 
  "Long (-E):","Lat (+N)"},
  {"Enter month", "Enter day", "Enter year (1901-2099)", 
  "Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC lag",
  "Longitude", "Latitude"}, 
  {1,1,1901,0, 0,0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″});
  // 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);
  lstd:= lstd; 
  // Greenwich Mean Time
  gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)
  thisday:= Y+m/100+D/10000;
  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;
  calcN();
  N0:= DDAYS(1901.0101, thisday)+1+0; // N at 0 GMT (for TS)
  TS:= ((23750.3 + 236.555362*(N0+gmt/24)) / 3600) MOD 24; // Sideral Time (GMT) in seconds
  TS:= →HMS(TS); // in HMS
  TSL:= (((23750.3 + 236.555362*N+lstd/24) / 3600) MOD 24); // local ST
  TSL:= (TSL + gmt - 4*(HMS→(long))/60) MOD 24;
  TSL:= →HMS(TSL); // in HMS
  TSd:= →HMS((98.965 + 0.985647342*N0) MOD 360); // ST at 0 GMT in degrees
  RETURN {TS, TSL, TSd, 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;
  calcN();
  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
  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)

  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;
  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(); // no 19
  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("");
  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("Equazione del tempo (min) " + trunc_angle(E));
  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));
  // RETURN {L, ω, M, v, l, r, x, y};
  RETURN [[L, ω], [l, M], [α,δ], [az, h], [x,y], [sorg, tram], [culm, r], [H,ts(1)]];
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;