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;

[Marcel]

Hi Salvomic
We don't use the same technique for calculating the planetary position. I use the VSOP technique (Variation séculaire des orbites planétaires). In this technique, we use polynomials and series of DATA to obtain the heliocentric position of a planet...
I complete today two new functions for the location of Pluto. This is the year of Pluto because in 25 days, the probe New Horizon will arrive at Pluto. For this reason, I will post my astronomy app's just before that moment.
Good lock in your project, this is fun and we can use both programs.

Marcel

[salvomic]

(06-20-2015 07:18 PM)Marcel Wrote:  Hi Salvomic
We don't use the same technique for calculating the planetary position. I use the VSOP technique (Variation séculaire des orbites planétaires). In this technique, we use polynomials and series of DATA to obtain the heliocentric position of a planet...
I complete today two new functions for the location of Pluto. This is the year of Pluto because in 25 days, the probe New Horizon will arrive at Pluto. For this reason, I will post my astronomy app's just before that moment.
Good lock in your project, this is fun and we can use both programs.

Marcel

hi Marcel,
yes, I'm using a different approach, following Bouiges: using the constants for planets, Sun and Moon...
About Pluto, my approach is very few, for now, as I'm using the old data (they were written in 1973...) and then Pluto was not completely known about Longitude, for example, so my program for now use some values as constant, like perigeo or longitude and anomaly...
You are right: we could now give more attention to Pluto.

However we could perhaps to search if some functions are common to the two programs and to find a common point of view...

Beside Pluto, in your program, I'm looking forward to see also the Moon functions, phases and so on...

Salvo

[Marcel]

Hi,
I send you my app's with Pluto, so you will be able to compare with your program. If you have Mathematica or WolframAlpha, check for precision. It is very good.
Marcel

[salvomic]

(06-20-2015 10:29 PM)Marcel Wrote:  Hi,
I send you my app's with Pluto, so you will be able to compare with your program. If you have Mathematica or WolframAlpha, check for precision. It is very good.
Marcel

thank Marcel,
I'm trying it!

for now with my program I get the same declination, but not the same RA, but my program has a problem to get right longitude of Pluto for now, your does it better...

I've both Mathematica and Wolfram, I'll check there also.

Salvo

EDIT: I reduced difference for Pluto between the to programs about (1' - 3'), for now: there was an error in a formula of mine...

[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;

[Marcel]

Hi Salvomic
We can calculate Delta T with some analytical way but this is not always precise. In fact we know the value of Delta T after the year. For instance, we will know the value for the year 2015, in 2016! This is a correction or a difference between two times scales.

In my app's, many data are in variable associate with the App. You can see these var if you click on the "var" tab in the conn kit. At this point this is more complex than your technique because many calculations most been done with these data in a FOR statement. If you check for the book of Jean Meeus on the internet, you will learn how to use that data.
Marcel

[salvomic]

(06-22-2015 06:26 PM)Marcel Wrote:  Hi Salvomic
We can calculate Delta T with some analytical way but this is not always precise. In fact we know the value of Delta T after the year. For instance, we will know the value for the year 2015, in 2016! This is a correction or a difference between two times scales.

In my app's, many data are in variable associate with the App. You can see these var if you click on the "var" tab in the conn kit. At this point this is more complex than your technique because many calculations most been done with these data in a FOR statement. If you check for the book of Jean Meeus on the internet, you will learn how to use that data.
Marcel

hi Marcel,
thank you for info.
I found the book of Meeus, yes, but just few minutes ago...
But I should rewrite all the program to conform 100% to it, and in this case there is already your app...
So, I'm trying to improve something adding formulas by Meeus to those by Bouiges. It isn't an easy job, hi...

Now, to improve Sun precision I should add at least 15 parameters, and solve the Keplero equation:
u - e sin(u) = M
in the var u, given e (eccentricity, by formula, ok) and M (anomaly), stopping it when the last iteration is no more great than previous for 1/10000, then calc v with u (v=acos(cos(u)-e/1-e*cos(u))) and so on...

Any hint to solve that equation?
I tried solve() but without luck...

thank you
Salvo

[Marcel]

Hi,
One of the chapter in the book describe a way to solve the Kepler equation (Chapter 29). You will find in this chapter a program to do it.
Marcel

[salvomic]

(06-22-2015 06:48 PM)Marcel Wrote:  Hi,
One of the chapter in the book describe a way to solve the Kepler equation (Chapter 29). You will find in this chapter a program to do it.
Marcel

I'll read it!

Salvo

[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;

[salvomic]

Marcel (and Manuel and others),
please, help me to find this problem now: I'm using for Sun and Moon the Meeus formulae mostly... If I calc Sun I get longitude almost equal to your program, if I calculate Moon I get the same results only if I use your app "at 0h time" (like {2015, 6, 25.0}) and not for the "justnow" time...
I'm trying using formulas like in your program, but including parameters (for sin and so on) into program, not in a variable...
Every numbers seems to be ok, but...
At the real local hour my program returns for Moon now some degrees of difference in longitude (with the first version was almost ok) if your app is used with "justnow" parameter; the same or so, if your is calculated at 0h...
So, the program could be in Julian day (but it seems to be ok) or in a global variable (I'm trying to put if off the program, as far as I can)...


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, h, h0)
BEGIN
// α, δ (AR, dec), height, h0 = standard height of body
  LOCAL temp, H0, tsideral, Hor, m0, m1, m2;
  LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
  α:= alfa*15; // DEG
  δ:= delta;
  θ0:= apparentSideralTime(makeDateListShort({Y,m,D,0,0,0})) *15; // TS in DEG
  ts:= calcTS();
  // θ0:= ts(2)*15; // TS=θ00 GMT
  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
 
  Hor:=15*HMS→(ts(2)-α);  // Hour angle, local ST vs apparent ST
  Hor:=simp_angle(→HMS(Hor));
  // Hor:= simp_angle(tsideral - long - α);
 
  m0:= HMS→(α + 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
  deltaM:= -HMS→(H/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
  deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));
  sorg:= ((m1 + deltaM)*24) MOD 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
  deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));
  tram:= ((m2 + deltaM) * 24) MOD 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

  // longitude:=-HMS→(long); // - if long -E
  Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST
  IF Hd<0 THEN Hd:=Hd+360 END;
  IF Hd>=360 THEN Hd:=Hd-360 END;
  Hd:=→HMS(Hd);
  H:= Hd/15; // in HMS
  
  temp:= sunAh(dateshortlist);
  az:= temp(1);
  h:=  temp(2);

  // 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(α, δ, h, 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:=→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
  //H:= (ts(2) - α) MOD 24; // Angolo orario
  //    // H:= asin(-tan(lat)*tan(δ));
  //Hd:= →HMS(15*H); // Angolo orario in gradi

  longitude:=-HMS→(long); // - if long -E
  Hd:=15*HMS→(apparentSideralTime(dateshortlist))-longitude-α;
  IF Hd<0 THEN Hd:=Hd+360 END;
  IF Hd>=360 THEN Hd:=Hd-360 END;
  Hd:=→HMS(Hd);
  H:= Hd/15;

  temp:= transformEquatorialToHorizontal(Hd, lat, δ);  
  //temp:= azimuthFromLBPhi(l, b, Hd);
  az:= temp(1);
  h:= temp(2);

  // 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 -0°34′; // 'standard' altitude of Moon (depend on parallax)
  // h0:= 0°7′30″; // medium value: 0.125°
  transito:= transit(α, δ, h, 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+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) " + trunc_angle(ts(3)));
  PRINT("θ0 Mean ST GMT " + trunc_angle(ts(5)));
  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, ε;

  // 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) + " (" +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("Diff. ascensionale " + trunc_angle(DA));
  PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));
  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;

[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;