HP Forums
Astronomy: Effemeridi (Ephemeris) - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP Prime Software Library (/forum-15.html)
+--- Thread: Astronomy: Effemeridi (Ephemeris) (/thread-4262.html)

Pages: 1 2


Astronomy: Effemeridi (Ephemeris) - salvomic - 06-29-2015 08:24 PM

hi everybody,
I'm developing a program about Astronomy.
I developed it on a TI 74 basic many years ago, now I've revisited (and enhanced) it for the Prime; it's name is Effemeridi (Astrolabio).
It can calculates ephemeris for Sun, Moon and planets, transit, sidereal time, nutation and obliquity, julian day, ascendant and many others data.
I'll put here new version as they will be completed.
I advice to use the program together with the splendid Astro Lab 4 from Marcel
The program present data in Terminal (multi pages: for now use enter to go on, and *never* touch screen or the program go on without presenting the other pages) and export variables with data.
For now some data names are still in Italian, I'll translate those later.

Enjoy!
Salvo Micciché

EDIT: new version, reviewed, more precise transit calculation.

EDIT 2017 November 8: last version (firmware from 10077 to the actual beta 12951) in this post of the thread.
After the public release of the new FW I hope to confirm that version or put a new one.

***
The first code, however, was this one:
Code:

data();
changeData();
whatToDo();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
tempo();
precession();
parametri();
transit();
sunCalc();
sun();
sunAh();
sunAlphaDelta();
moon();
moonCalc();
moonAh();
moonAlphaDelta();
planets();
planetCalc();
planetDisplay();
planetAlphaDelta();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendant();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();

EXPORT Tempi;
EXPORT EpsilonPsi;
EXPORT Ascendant;
EXPORT Sun_Effemerids;
EXPORT Sun_Transit;
EXPORT Sun_Others;
EXPORT Moon_Effemerids;
EXPORT Moon_Transit;
EXPORT Moon_Others;
EXPORT Planet_Effemerids;
EXPORT Planet_Transit;
EXPORT Planet_Others;

thisday:= 1964.0529;
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,40);
moonlist:= MAKELIST(X,X,1,40);
planetlist:= MAKELIST(X,X,1,40);
planet_const:= MAKELIST(X,X,1,40);

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, DrD
BEGIN
  HAngle:=1;
  data(); // input data and calc Sun
END;

data()
BEGIN
  // initial input
  LOCAL dd, mm, yy, hh, mn, ss, loct;
  LOCAL ε, JD, JDE;
  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};

  ε:=epsilon(dateshortlist);
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:=(JD-2451545)/36525;
  U:= T/100;
  sunCalc(dateshortlist);  // set initial data for Sun
  moonCalc(dateshortlist);
  whatToDo(dateshortlist);
END;


changeData()
BEGIN
  // input to change date, lat, long...
  LOCAL ε, JD, JDE;

  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},{m,D,Y,lstd, dst, utc, long, lat, deltaT});

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

  ε:=epsilon(dateshortlist);
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:=(JD-2451545)/36525;
  U:= T/100;

  whatToDo(dateshortlist);
END;

whatToDo(lista)
BEGIN
  LOCAL ch;
  CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendant",
    "Quit");
  CASE
  IF ch==1 THEN sunCalc(lista); sun(lista); END;
  IF ch==2 THEN moonCalc(lista); moon(lista); END;
  IF ch==3 THEN planets(lista); END;
  IF ch==4 THEN tempo(lista); END;
  IF ch==5 THEN parametri(lista); END;
  IF ch==6 THEN ascendant(lista); END;
  IF ch==7 THEN HAngle:=0; RETURN; END;
  DEFAULT 
  END; // case
END;


planets(lista)
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), lista);
  planetDisplay(nameplanet(ch),lista); 
END;

epsilon(lista)
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  T:=(JD-2451545)/36525;
  U:= T/100;
  // ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
  // 4680.93″ = 1°18′00.93″
  ε0:= 23°26′21.448″ - (HMS→(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(lista);
  eps:= dep(2);
  ε:= ε0+eps; // true obliquity ε
  RETURN ε;
END;

deltaEpsilonPsi(lista)
BEGIN
  // Nutazione Δψ (longitudine), Δε (obliquità)
  LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
  LOCAL DD;
  // JD:= julianDayAtZero(Y,m,D);
  JDE:= julianEphemerisDay(lista);
  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;

transformEclipticalToEquatorial(lista, lambda, beta)
BEGIN
  // trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
  LOCAL ε,Lambda, Beta, alpha, delta;
  Lambda:=HMS→(lambda);
  Beta:=HMS→(beta);
  ε:= epsilon(lista);
  alpha:=atan2(sin(Lambda)*cos(ε)-tan(Beta)*sin(ε),cos(Lambda));
  IF alpha>=360 THEN alpha:=alpha-360 END;
  IF alpha<0 THEN alpha:=alpha+360 END;
  delta:=asin(sin(Beta)*cos(ε)+cos(Beta)*sin(ε)*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(lista, Lambda, Beta, Lha)
BEGIN
  // transform λ and β (long, lat) into azimuth and height
  LOCAL lambda, beta,phi, lha, ε;
  LOCAL alpha,delta,azimuth,altitude;
  lambda:=HMS→(Lambda); // longitudine astro
  beta:=HMS→(Beta); // latitudine astro
  ε:= epsilon(lista);
  lha:=HMS→(Lha);  // angolo orario
  phi:=HMS→(lat); // lat geo
  alpha:=atan2(SIN(lambda)*COS(ε)-TAN(beta)*SIN(ε),COS(lambda));
  delta:=asin(sin(beta)*cos(ε)+cos(beta)*sin(ε)*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(lista)
BEGIN
  LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
  LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  JD0:= julianDayAtZero(Y,m,D);
  dep:= deltaEpsilonPsi(lista);
  psi:=dep(1);
  eps:=dep(2);
  T:= (JD-2451545.0)/36525;
  θ0:= meanSideralTime(lista); // 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(lista); // apparent ST (ST at Greenwitch at local our)
  TSL:= apparentSideralTime(lista) - 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(lista);
  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(lista);
  MST:= meanSideralTime(lista);
  dep:= deltaEpsilonPsi(lista);
  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(lista)
BEGIN
  LOCAL JD, JDE, TSid, dep, psi, eps;
  LOCAL ε, listatempi, ch;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  TSid:= calcTS(lista);
  dep:= deltaEpsilonPsi(lista);
  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(" ");
WAIT();
PRINT();
  PRINT("Δψ Nutazione (longit.) "+ psi);
  PRINT("Δε Nutazione (obl.) "+ eps);
  PRINT("ε Obliquità Ecl. "+ ε);
  PRINT(" ");
  PRINT("Time " + lstd + " GMT " + gmt);
  listatempi:= {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)) };
  Tempi:= listatempi;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

transit(lista, listaARh0)
BEGIN
  // lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
  // α, δ (AR, dec), height, h0 = standard height of body
  LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;
  LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
  LOCAL alfa1, alfa3, delta1, delta3, h0, alpha, del, n;
  α:= listaARh0(2); // AR for day
  δ:= listaARh0(5); // decl for day
  alfa1:= listaARh0(1); // AR day-1
  alfa3:= listaARh0(3); // AR day+1
  delta1:= listaARh0(4); // decl day-1
  delta3:= listaARh0(6); // decl day+1
  h0:= listaARh0(7); // standard height
  temp:= calcTS(lista);
  θ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
  n:= m0 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*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
  n:= m1 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*15);
  Hor:=simp_angle(→HMS(Hor));
  h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
  deltaM:= HMS→((h-h0)/(360*cos(del)*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
  n:= m2 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*15); 
  Hor:=simp_angle(→HMS(Hor));
  h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
  deltaM:= HMS→((h-h0)/(360*cos(del)*cos(lat)*sin(Hor)));
  tram:= ((m2 + deltaM) * 24); // setting

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

parametri(lista)
BEGIN
  LOCAL dep, prec, α, δ, ε, listaparametri, ch, temp;
  dep:= deltaEpsilonPsi(lista);
  temp:= sunAlphaDelta(lista);
  α:= sunlist(5); δ:= sunlist(6);
  prec:= precession(α, δ);
  ε:= epsilon(lista);
  // 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));
  listaparametri:= {ε, dep(1), dep(2), prec(1), prec(2)}; // ε, Δψ, Δε, Δα, Δδ
  EpsilonPsi:= listaparametri;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

ascendant(lista)
BEGIN
  LOCAL ASC, ts, segno, ε, listaascendente, ch, temp;
  LOCAL tsb;
  ε:= epsilon(lista);
  ts:= calcTS(lista);
  tsb:= ts(1)+lstd-dst;  // ST of the birth
  // temp:= sin(ts(2)*15)*cos(ε)+tan(lat)*sin(ε);
  temp:= sin(tsb*15)*cos(ε)+tan(lat)*sin(ε);
  ASC:= atan2(-cos(tsb*15), temp);
  // IF temp<0 THEN ASC:= ASC+180 ELSE ASC:= ASC+360; END;
  IF ASC<180 THEN ASC:=ASC+180 ELSE ASC:= ASC-180; END;
  ASC:= simp_angle(ASC);
  segno:= zodiaco(ASC);
  listaascendente:= {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno};
  Ascendant:= listaascendente;
    RECT_P();
   TEXTOUT_P("Ascendant " + trunc_angle(→HMS(ASC)), 25,50);
   TEXTOUT_P(" " + trunc_angle(→HMS(ASC) MOD 30), 25,70);
   TEXTOUT_P("in " + segno, 100,70);

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
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(lista)
// type dateshortlist
BEGIN
  LOCAL ω, l, v, r, x, y, dz, ε;
  LOCAL α, δ, dayY, ts, az, z, h;
  LOCAL SAD, DA, Hd, lt;
  LOCAL as, ac, zc, h0, H0, temp;
  LOCAL ec, a, u, M1, M2, j, t, tau;
  LOCAL JD, JDE, c, lapp, Ω, θ0;
  LOCAL dep, deltapsi, longitude;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  a:= 1.00000023; // semiasse maggiore
  dep:= deltaEpsilonPsi(lista);
  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.003032028*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
  //M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-(T^3)/24490000); // anom vera Meeus
  M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // Meeus p 308
  Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
  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(lista);  // 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:= transformEquatorialToHorizontal(Hd, lat, δ);
  az:= HMS→(temp(1))+180;  // Azimuth
  IF az>360 THEN az:= az-360 END;
  az:= →HMS(az);
  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

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

sun(lista)
BEGIN
  LOCAL ω, l, v, r, x, y, dz, ε;
  LOCAL α, δ, dayY, ts, az, z, h;
  LOCAL SAD, DA, sorg, tram, culm, Hd, h0, lt;
  LOCAL as, ac, zc, segno, JD, JDE, v, c;
  LOCAL Ω, lapp, ec, listasun, listasuntran, listasunothers; 
  LOCAL ch, transito, listaARh0, temp, Hor;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3;
  ε:= epsilon(lista);
  L:= sunlist(1); ω:= sunlist(2); 
  l:= sunlist(3); lapp:=sunlist(4);
  α:= sunlist(5); δ:= sunlist(6); 
  az:= sunlist(7); h:= sunlist(8);
  M:= sunlist(9); Ω:= sunlist(10);   
  v:= sunlist(11); c:= sunlist(12); 
  x:= sunlist(13); y:= sunlist(14);
  r:= sunlist(15); Hor:= sunlist(16); 
  Hd:= 15*sunlist(16); zc:= sunlist(17); 
  dz:= sunlist(18);
  DA:= sunlist(19); SAD:= sunlist(20); 
  E:= sunlist(21); lt:= sunlist(22);
  ec:= sunlist(23);
  ts:= calcTS(lista);
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  segno:= zodiaco(l);

   // 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
  // to pass lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
  h0:= -0°50′; // 'standard' altitude of Sun
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  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("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("λ longitudine apparente " + trunc_angle(lapp));
  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(" ");
WAIT();
PRINT();
  PRINT("Ω aberr, nutaz nodo "+ Ω);
  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("Ang. orario "+ trunc_angle(Hor) + " " + 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(" ");
WAIT();
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);
  l:= →HMS(l); lapp:= →HMS(lapp);
  culm:= →HMS(culm); sorg:= →HMS(sorg); tram:= →HMS(tram);
  listasun:= {l, α, δ, az, h, r};
  listasuntran:= {culm, sorg, tram};
  listasunothers:= {L, lapp, ω, Ω, M,v, c, x, y, Hor, lt, ε, ec};
  Sun_Effemerids:= listasun;
  Sun_Transit:= listasuntran;
  Sun_Others:= listasunothers;
  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

sunAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= sunCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;

sunAh(lista)
// lista like dateshortlist
BEGIN
  LOCAL lha,alpha,delta,phi, ah, alfa;
  LOCAL temp,longitude, Azimuth, Altitude, ts;
  temp:= sunAlphaDelta(lista);
  alpha:=15*HMS→(temp(1));
  alfa:= →HMS(temp(1)); // α in HMS
  delta:= temp(2);
  phi:=lat;
  ts:= calcTS(lista);
  longitude:= - HMS→(long); // - if long -E
  lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
  // lha:= 15*HMS→(ts(2)-alfa);  // 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;

moonCalc(lista)
BEGIN
  LOCAL ω,Ω, L1, M1, l, lsum, bsum, rsum; 
  LOCAL v, r, x, y, dep, ε;
  LOCAL b, P, d, s, temp, tau;
  LOCAL α,δ, ts, Hd, zc, dz, h, az;
  LOCAL i, as, ac, DA, SAD, DD;
  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(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  T:= (JDE-2451545.0)/36525; // use JDE
  tau:= (JDE-2451545.0)/365250; // millenni  - JDE
  dep:= deltaEpsilonPsi(lista);
  // Sun
  L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
  M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // mean anomaly Sun, Meeus  
  // 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 +  280602 * 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*DD+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(lista, l,b);  // transform into AR, decl
  α:= temp(1);
  δ:= temp(2);

  ts:= calcTS(lista);  // calcola tempo siderale
  Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST
  Hd:=simp_angle(→HMS(Hd));
  H:= Hd/15; // in HMS

  DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
  SAD:= (90 + DA); // Semiarco Diurno
  // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
  
  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)

  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 Calante"; END;
  IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
  IF ph>100 AND ph<170 THEN fase:="Gibbosa Crescente"; 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
  moonlist:= {L1, l, b, lapp, α, δ, az, h, M1, ω, Ω, d, DD, P, s, H, DA, SAD,ph, fase, eta, ill};
RETURN moonlist;
END;

moon(lista)
BEGIN
  LOCAL ε, ts, JD, JDE, segno, L1, b;
  LOCAL l, Ω, ω, M1, lapp, α, δ;
  LOCAL az, h, sorg, tram, culm, Hd, DA, SAD;
  LOCAL ph, fase, eta, DD, d, s, ill;
  LOCAL dep, temp, listamoon, listamoontran, listamoonothers, ch;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
  LOCAL listaARh0, transito, Hor;

  ε:= epsilon(lista);
  ts:= calcTS(lista); //
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  dep:= deltaEpsilonPsi(lista);
  L1:= moonlist(1); l:= moonlist(2);
  b:= moonlist(3);
  α:= moonlist(5); δ:= moonlist(6); 
  az:= moonlist(7); h:= moonlist(8); 
  ω:= moonlist(10); Ω:= moonlist(11); 
  lapp:= moonlist(4);  M1:= moonlist(9);
  d:= moonlist(12); DD:= moonlist(13);
  P:= moonlist(14); s:= moonlist(15);
  Hor:= moonlist(16); Hd:= 15*moonlist(16);
  DA:= moonlist(17); SAD:= moonlist(18);
  ph:= moonlist(19); fase:= moonlist(20);

  eta:= moonlist(21); ill:= moonlist(22);
  segno:= zodiaco(l);

  // transit routine
  h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)
  // h0:= 0°7′30″; // medium value: 0.125°
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  PRINT;
  PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
  PRINT("L Long media " + trunc_angle(L1));
  PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("β Latitudine " + trunc_angle(b));
  PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
  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("P Parallasse " + trunc_angle(P));
  PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");
  PRINT(" ");
WAIT();
PRINT();
  PRINT("ω Perigeo medio " + trunc_angle(ω));
  PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
  PRINT("M anomalia vera " + trunc_angle(M1));
  PRINT("D Elongazione " + trunc_angle(DD));
  PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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("Julian day of Ephemeris " + JDE);
  PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
  PRINT(" ");
WAIT();
PRINT();
  PRINT("Δψ Nutazione long "+ dep(1));
  PRINT("Δε Nutazione obl. " + dep(2));
  PRINT("ε inclinazione ecl. " + ε);
  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);
  PRINT("Età " + eta);
  PRINT("Parte illuminata " + ROUND(ill,5) + " " + ROUND(100*ill,2)+"%");
  listamoon:= {→HMS(l), →HMS(b), ROUND(HMS→(d),6), →HMS(α), →HMS(δ), →HMS(az), →HMS(h)};
  listamoonothers:= {L1, lapp, →HMS(α*15), M1,P, ω,Ω,s, Hor,ph, eta};
  listamoontran:= →HMS({culm, sorg, tram});
  Moon_Effemerids:= listamoon;
  Moon_Transit:= listamoontran;
  Moon_Others:= listamoonothers;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

moonAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= moonCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;

moonAh(lista)
// lista like dateshortlist
BEGIN
  LOCAL lha,alpha,delta,phi, ah, alfa;
  LOCAL temp,longitude, Azimuth, Altitude, ts;
  temp:= sunAlphaDelta(lista);
  alpha:=15*HMS→(temp(1));
  alfa:= →HMS(temp(1)); // α in HMS
  delta:= temp(2);
  phi:=lat;
  ts:= calcTS(lista);
  longitude:= - HMS→(long); // - if long -E
  lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
  // lha:= 15*HMS→(ts(2)-alfa);  // 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;

planetCalc(nameplanet, lista)
BEGIN
  LOCAL ω,Ω, l, v, r, x, y, z;
  LOCAL Ls, Ms, b, xs, ys, b_ec;
  LOCAL x1, y1, z1, l_ec, 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, lapp;
  LOCAL listaplanet, listaplanettran, listaplanetothers, ch;

  JD:= julianDay(lista);
  ε:= epsilon(lista);

  // calc of Sun (as basis)
  sunCalc(lista);
  xs:= sunlist(13); ys:= sunlist(14);
  Ls:= sunlist(1); // Longitudine media
  Ms:= sunlist(9); // 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(lista, l,b);  // transform into AR, decl
  α:= temp(1);
  δ:= temp(2);

  ts:= calcTS(lista);  // 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

  planetlist:= {L, l, b, lapp, α,δ, az, h, M, Ω, ω, v, x,y, z, x1, y1, z1, 
    r, R, H, zc, dz, DA, SAD, ec, in, a, θ, l_ec, b_ec, Ls, Ms, xs, ys, nameplanet};
  RETURN planetlist;
END;

planetDisplay(nameplanet, lista)
BEGIN
  LOCAL ω,Ω, l, v, r;
  LOCAL x, y, z, x1, y1, z1;
  LOCAL Ls, Ms, b, xs, ys, l_ec, b_ec;
  LOCAL ec, in, a, δ, α, ts, az, h;
  LOCAL DA, SAD, tram, sorg, culm, dayY;
  LOCAL as, ac, zc, dz, θ, mov, segno, ε;
  LOCAL transito, JD, JDE, temp, Hor, Hd;
  LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
  LOCAL listaARh0, lapp;

  ts:= calcTS(lista);
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);

  L:= planetlist(1); l:= planetlist(2);
  b:= planetlist(3); lapp:= planetlist(4);
  α:= planetlist(5); δ:= planetlist(6);
  az:= planetlist(7); h:= planetlist(8);
  M:= planetlist(9); Ω:= planetlist(10);
  ω:= planetlist(11); v:= planetlist(12);
  x := planetlist(13);  y:= planetlist(14);  z:= planetlist(15);
  x1:= planetlist(16); y1:= planetlist(17); z1:= planetlist(18);
  r:= planetlist(19);   R:= planetlist(20);
  Hor:= planetlist(21); Hd:= 15*planetlist(21);
  zc:= planetlist(22); dz:= planetlist(23);
  DA:= planetlist(24); SAD:= planetlist(25);
  ec:= planetlist(26); in:= planetlist(27);
  a:= planetlist(28); θ:= planetlist(29);
  l_ec:= planetlist(30); b_ec:= planetlist(31);
  Ls:= planetlist(32); Ms:= planetlist(33);  // Sun values
  xs:= planetlist(34); ys:= planetlist(35);
  // nameplanet:= planetlist(36);

  segno:= zodiaco(l);

  // transit routine
  h0:= -0°34′; // 'standard' altitude of a planet
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  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("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("ω 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));
WAIT();
PRINT();
  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(" ");
  PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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(" ");
WAIT();
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);
  listaplanet:= {l,b, α,δ,az,h, Hor, Hd, mov};
  listaplanettran:= {culm, sorg, tram};
  listaplanetothers:= {L,ω,Ω, M,l_ec,b_ec, v, x,y,z, x1,y1,z1, r, R, SAD, DA};
  Planet_Effemerids:= listaplanet;
  Planet_Transit:= listaplanettran;
  Planet_Others:= listaplanetothers;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;

END;

planetAlphaDelta(nameplanet, lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= planetCalc(nameplanet, lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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;
***


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-01-2015 03:45 PM

New version, reviewed.
Corrected some bugs.
More user-friendly, I hope.
Transit routine returns more precise results, now.

Code:

data();
changeData();
whatToDo();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
tempo();
precession();
parametri();
transit();
sunCalc();
sun();
sunAh();
sunAlphaDelta();
moon();
moonCalc();
moonAh();
moonAlphaDelta();
planets();
planetCalc();
planetDisplay();
planetAlphaDelta();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendant();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();

EXPORT Tempi;
EXPORT EpsilonPsi;
EXPORT Ascendant;
EXPORT Sun_Effemerids;
EXPORT Sun_Transit;
EXPORT Sun_Others;
EXPORT Moon_Effemerids;
EXPORT Moon_Transit;
EXPORT Moon_Others;
EXPORT Planet_Effemerids;
EXPORT Planet_Transit;
EXPORT Planet_Others;

thisday:= 1964.0529;
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,40);
moonlist:= MAKELIST(X,X,1,40);
planetlist:= MAKELIST(X,X,1,40);
planet_const:= MAKELIST(X,X,1,40);

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, DrD
BEGIN
  HAngle:=1;
  data(); // input data and calc Sun
END;

data()
BEGIN
  // initial input
  LOCAL dd, mm, yy, hh, mn, ss, loct;
  LOCAL ε, JD, JDE;
  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};

  ε:=epsilon(dateshortlist);
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:=(JD-2451545)/36525;
  U:= T/100;
  sunCalc(dateshortlist);  // set initial data for Sun
  moonCalc(dateshortlist);
  whatToDo(dateshortlist);
END;


changeData()
BEGIN
  // input to change date, lat, long...
  LOCAL ε, JD, JDE;

  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},{m,D,Y,lstd, dst, utc, long, lat, deltaT});

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

  ε:=epsilon(dateshortlist);
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:=(JD-2451545)/36525;
  U:= T/100;

  whatToDo(dateshortlist);
END;

whatToDo(lista)
BEGIN
  LOCAL ch;
  CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendant",
    "Quit");
  CASE
  IF ch==1 THEN sunCalc(lista); sun(lista); END;
  IF ch==2 THEN moonCalc(lista); moon(lista); END;
  IF ch==3 THEN planets(lista); END;
  IF ch==4 THEN tempo(lista); END;
  IF ch==5 THEN parametri(lista); END;
  IF ch==6 THEN ascendant(lista); END;
  IF ch==7 THEN HAngle:=0; RETURN; END;
  DEFAULT 
  END; // case
END;


planets(lista)
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), lista);
  planetDisplay(nameplanet(ch),lista); 
END;

epsilon(lista)
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  T:=(JD-2451545)/36525;
  U:= T/100;
  // ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
  // 4680.93″ = 1°18′00.93″
  ε0:= 23°26′21.448″ - (HMS→(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(lista);
  eps:= dep(2);
  ε:= ε0+eps; // true obliquity ε
  RETURN ε;
END;

deltaEpsilonPsi(lista)
BEGIN
  // Nutazione Δψ (longitudine), Δε (obliquità)
  LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
  LOCAL DD;
  // JD:= julianDayAtZero(Y,m,D);
  JDE:= julianEphemerisDay(lista);
  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;

transformEclipticalToEquatorial(lista, lambda, beta)
BEGIN
  // trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
  LOCAL ε,Lambda, Beta, alpha, delta;
  Lambda:=HMS→(lambda);
  Beta:=HMS→(beta);
  ε:= epsilon(lista);
  alpha:=atan2(sin(Lambda)*cos(ε)-tan(Beta)*sin(ε),cos(Lambda));
  IF alpha>=360 THEN alpha:=alpha-360 END;
  IF alpha<0 THEN alpha:=alpha+360 END;
  delta:=asin(sin(Beta)*cos(ε)+cos(Beta)*sin(ε)*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(lista, Lambda, Beta, Lha)
BEGIN
  // transform λ and β (long, lat) into azimuth and height
  LOCAL lambda, beta,phi, lha, ε;
  LOCAL alpha,delta,azimuth,altitude;
  lambda:=HMS→(Lambda); // longitudine astro
  beta:=HMS→(Beta); // latitudine astro
  ε:= epsilon(lista);
  lha:=HMS→(Lha);  // angolo orario
  phi:=HMS→(lat); // lat geo
  alpha:=atan2(SIN(lambda)*COS(ε)-TAN(beta)*SIN(ε),COS(lambda));
  delta:=asin(sin(beta)*cos(ε)+cos(beta)*sin(ε)*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(lista)
BEGIN
  LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
  LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  JD0:= julianDayAtZero(Y,m,D);
  dep:= deltaEpsilonPsi(lista);
  psi:=dep(1);
  eps:=dep(2);
  T:= (JD-2451545.0)/36525;
  θ0:= meanSideralTime(lista); // 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(lista); // apparent ST (ST at Greenwitch at local our)
  TSL:= apparentSideralTime(lista) - 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(lista);
  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(lista);
  MST:= meanSideralTime(lista);
  dep:= deltaEpsilonPsi(lista);
  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(lista)
BEGIN
  LOCAL JD, JDE, TSid, dep, psi, eps;
  LOCAL ε, listatempi, ch;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  TSid:= calcTS(lista);
  dep:= deltaEpsilonPsi(lista);
  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(" ");
WAIT();
PRINT();
  PRINT("Δψ Nutazione (longit.) "+ psi);
  PRINT("Δε Nutazione (obl.) "+ eps);
  PRINT("ε Obliquità Ecl. "+ ε);
  PRINT(" ");
  PRINT("Time " + lstd + " GMT " + gmt);
  listatempi:= {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)) };
  Tempi:= listatempi;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

transit(lista, listaARh0)
BEGIN
  // lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
  // α, δ (AR, dec), height, h0 = standard height of body
  LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;
  LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
  LOCAL alfa1, alfa3, delta1, delta3, h0, alpha, del, n;
  α:= listaARh0(2); // AR for day
  δ:= listaARh0(5); // decl for day
  alfa1:= listaARh0(1); // AR day-1
  alfa3:= listaARh0(3); // AR day+1
  delta1:= listaARh0(4); // decl day-1
  delta3:= listaARh0(6); // decl day+1
  h0:= listaARh0(7); // standard height
  temp:= calcTS(lista);
  θ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
  n:= m0 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*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
  n:= m1 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*15);
  Hor:=simp_angle(→HMS(Hor));
  h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
  deltaM:= HMS→((h-h0)/(360*cos(del)*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
  n:= m2 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*15); 
  Hor:=simp_angle(→HMS(Hor));
  h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
  deltaM:= HMS→((h-h0)/(360*cos(del)*cos(lat)*sin(Hor)));
  tram:= ((m2 + deltaM) * 24); // setting

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

parametri(lista)
BEGIN
  LOCAL dep, prec, α, δ, ε, listaparametri, ch, temp;
  dep:= deltaEpsilonPsi(lista);
  temp:= sunAlphaDelta(lista);
  α:= sunlist(5); δ:= sunlist(6);
  prec:= precession(α, δ);
  ε:= epsilon(lista);
  // 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));
  listaparametri:= {ε, dep(1), dep(2), prec(1), prec(2)}; // ε, Δψ, Δε, Δα, Δδ
  EpsilonPsi:= listaparametri;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

ascendant(lista)
BEGIN
  LOCAL ASC, ts, segno, ε, listaascendente, ch, temp;
  LOCAL tsb;
  ε:= epsilon(lista);
  ts:= calcTS(lista);
  tsb:= ts(1)+lstd-dst;  // ST of the birth
  // temp:= sin(ts(2)*15)*cos(ε)+tan(lat)*sin(ε);
  temp:= sin(tsb*15)*cos(ε)+tan(lat)*sin(ε);
  ASC:= atan2(-cos(tsb*15), temp);
  // IF temp<0 THEN ASC:= ASC+180 ELSE ASC:= ASC+360; END;
  IF ASC<180 THEN ASC:=ASC+180 ELSE ASC:= ASC-180; END;
  ASC:= simp_angle(ASC);
  segno:= zodiaco(ASC);
  listaascendente:= {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno};
  Ascendant:= listaascendente;
    RECT_P();
   TEXTOUT_P("Ascendant " + trunc_angle(→HMS(ASC)), 25,50);
   TEXTOUT_P(" " + trunc_angle(→HMS(ASC) MOD 30), 25,70);
   TEXTOUT_P("in " + segno, 100,70);

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
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(lista)
// type dateshortlist
BEGIN
  LOCAL ω, l, v, r, x, y, dz, ε;
  LOCAL α, δ, dayY, ts, az, z, h;
  LOCAL SAD, DA, Hd, lt;
  LOCAL as, ac, zc, h0, H0, temp;
  LOCAL ec, a, u, M1, M2, j, t, tau;
  LOCAL JD, JDE, c, lapp, Ω, θ0;
  LOCAL dep, deltapsi, longitude;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  a:= 1.00000023; // semiasse maggiore
  dep:= deltaEpsilonPsi(lista);
  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.003032028*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
  //M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-(T^3)/24490000); // anom vera Meeus
  M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // Meeus p 308
  Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
  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(lista);  // 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:= transformEquatorialToHorizontal(Hd, lat, δ);
  az:= HMS→(temp(1))+180;  // Azimuth
  IF az>360 THEN az:= az-360 END;
  az:= →HMS(az);
  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

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

sun(lista)
BEGIN
  LOCAL ω, l, v, r, x, y, dz, ε;
  LOCAL α, δ, dayY, ts, az, z, h;
  LOCAL SAD, DA, sorg, tram, culm, Hd, h0, lt;
  LOCAL as, ac, zc, segno, JD, JDE, v, c;
  LOCAL Ω, lapp, ec, listasun, listasuntran, listasunothers; 
  LOCAL ch, transito, listaARh0, temp, Hor, dep;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3;
  ε:= epsilon(lista);
  L:= sunlist(1); ω:= sunlist(2); 
  l:= sunlist(3); lapp:=sunlist(4);
  α:= sunlist(5); δ:= sunlist(6); 
  az:= sunlist(7); h:= sunlist(8);
  M:= sunlist(9); Ω:= sunlist(10);   
  v:= sunlist(11); c:= sunlist(12); 
  x:= sunlist(13); y:= sunlist(14);
  r:= sunlist(15); Hor:= sunlist(16); 
  Hd:= 15*sunlist(16); zc:= sunlist(17); 
  dz:= sunlist(18);
  DA:= sunlist(19); SAD:= sunlist(20); 
  E:= sunlist(21); lt:= sunlist(22);
  ec:= sunlist(23);
  ts:= calcTS(lista);
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  dep:= deltaEpsilonPsi(lista);
  segno:= zodiaco(l);

   // 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
  // to pass lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
  h0:= -0°50′; // 'standard' altitude of Sun
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  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("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("λ longitudine apparente " + trunc_angle(lapp));
  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(" ");
WAIT();
PRINT();
  PRINT("Ω aberr, nutaz nodo "+ Ω);
  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("Ang. orario "+ trunc_angle(Hor) + " " + 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(" ");
WAIT();
PRINT();
  PRINT("Transita a: " + trunc_angle(culm MOD 24) + " 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("ε obliquità ecl. " + ε);
  PRINT("ec eccentricità "+ ec);
  PRINT("Δψ Nutazione long "+ dep(1));
  PRINT("Δε Nutazione obl. " + dep(2));
  PRINT(" ");
  PRINT("Long elioc. Terra "+ trunc_angle(lt));

  l:= →HMS(l); lapp:= →HMS(lapp);
  culm:= →HMS(culm); sorg:= →HMS(sorg); tram:= →HMS(tram);
  listasun:= {l, α, δ, az, h, r};
  listasuntran:= {culm, sorg, tram};
  listasunothers:= {L, lapp, ω, Ω, M,v, c, x, y, Hor, lt, ε, ec};
  Sun_Effemerids:= listasun;
  Sun_Transit:= listasuntran;
  Sun_Others:= listasunothers;
  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

sunAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= sunCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;

sunAh(lista)
// lista like dateshortlist
BEGIN
  LOCAL lha,alpha,delta,phi, ah, alfa;
  LOCAL temp,longitude, Azimuth, Altitude, ts;
  temp:= sunAlphaDelta(lista);
  alpha:=15*HMS→(temp(1));
  alfa:= →HMS(temp(1)); // α in HMS
  delta:= temp(2);
  phi:=lat;
  ts:= calcTS(lista);
  longitude:= HMS→(long); // - if long -E
  lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
  // lha:= 15*HMS→(ts(2)-alfa);  // 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;

moonCalc(lista)
BEGIN
  LOCAL ω,Ω, Lprime, Mprime, l, lsum, bsum, rsum; 
  LOCAL v, r, x, y, dep, ε;
  LOCAL b, P, d, s, temp, tau;
  LOCAL α,δ, ts, Hd, zc, dz, h, az;
  LOCAL i, as, ac, DA, SAD, DD;
  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(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  T:= (JDE-2451545.0)/36525; // use JDE
  tau:= (JDE-2451545.0)/365250; // millenni  - JDE
  dep:= deltaEpsilonPsi(lista);
  // Sun
  L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
  M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // mean anomaly Sun, Meeus  
  // end Sun
  // L1:= (33.231 + 13.17639653*N) MOD 360; // long. media
  Lprime:= 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
  Mprime:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media
  ω:= (Lprime - Mprime) MOD 360; // perigeo medio
  // DD:= Lprime - L; F:= Lprime-Ω; // 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(Mprime); // calc Moon longitude
  lsum:= lsum + 1274027 * sin(2*DD-Mprime); // evezione
  lsum:= lsum +  658314 * sin(2*DD); // variazione
  lsum:= lsum +  213618 * sin(2*Mprime); // 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*Mprime);
  lsum:= lsum +   57066 * sin(2*DD-Mprime-M) *ee;
  lsum:= lsum +   53322 * sin(2*DD+Mprime);
  lsum:= lsum +   45758 * sin(2*DD-M) *ee;
  lsum:= lsum -   40923 * sin(M-Mprime) *ee;
  lsum:= lsum -   34720 * sin(DD);
  lsum:= lsum -   30383 * sin(Mprime+M) *ee;
  lsum:= lsum +   15327 * sin(2*DD-2*F);
  lsum:= lsum -   12528 * sin(Mprime+2*F);
  lsum:= lsum +   10980 * sin(Mprime-2*F);
  lsum:= lsum +   10675 * sin(4*DD-Mprime);
  lsum:= lsum +   10034 * sin(3*Mprime);
  lsum:= lsum +    8548 * sin(4*DD-2*Mprime);
  lsum:= lsum -    7888 * sin(2*DD+M-Mprime) *ee;
  lsum:= lsum -    6766 * sin(2*DD+M) *ee;
  lsum:= lsum -    5163 * sin(DD-Mprime);
  lsum:= lsum +    4987 * sin(DD+M) *ee;
  lsum:= lsum +    4036 * sin(2*DD-M+Mprime) *ee;
  lsum:= lsum +    3994 * sin(2*DD+2*Mprime);
  lsum:= lsum +    3861 * sin(4*DD);
  lsum:= lsum +    3665 * sin(2*DD-3*Mprime);
  lsum:= lsum -    2689 * sin(M-2*Mprime) *ee;
  lsum:= lsum -    2602 * sin(2*DD-Mprime+2*F);
  lsum:= lsum +    2390 * sin(2*DD-M-2*Mprime) *ee;
  lsum:= lsum -    2348 * sin(DD+Mprime);
  lsum:= lsum +    2236 * sin(2*DD-2*M) *ee^2;
  lsum:= lsum -    2120 * sin(M+2*Mprime) *ee;
  lsum:= lsum -    2069 * sin(2*M) *ee^2;
  lsum:= lsum +    2048 * sin(2*DD-2*M-Mprime) *ee^2;
  lsum:= lsum -    1773 * sin(2*DD+Mprime-2*F);
  lsum:= lsum -    1595 * sin(2*DD+2*F);
  lsum:= lsum +    1215 * sin(4*DD-M-Mprime) *ee;
  lsum:= lsum -    1110 * sin(2*Mprime+2*F);
  lsum:= lsum -     892 * sin(3*DD-Mprime);
  lsum:= lsum -     810 * sin(2*DD+M+Mprime) *ee;
  lsum:= lsum +     759 * sin(4*DD-M-2*Mprime) *ee;
  lsum:= lsum -     713 * sin(2*M-Mprime) * ee^2;
  lsum:= lsum -     700 * sin(2*DD+2*M-Mprime) *ee^2;
  lsum:= lsum +     691 * sin(2*DD+M-2*Mprime) *ee;
  lsum:= lsum +     596 * sin(2*DD-M-2*F) *ee;
  lsum:= lsum +     549 * sin(4*DD+Mprime);
  lsum:= lsum +     537 * sin(4*Mprime);
  lsum:= lsum +     520 * sin(4*DD-M) *ee;
  lsum:= lsum -     487 * sin(DD-2*Mprime);
  lsum:= lsum -     399 * sin(2*DD+M-2*F) *ee;
  lsum:= lsum -     381 * sin(2*Mprime-2*F);
  lsum:= lsum +     351 * sin(DD+M+Mprime) *ee;
  lsum:= lsum -     340 * sin(3*DD-2*Mprime);
  lsum:= lsum +     330 * sin(4*DD-3*Mprime);
  lsum:= lsum +     327 * sin(2*DD-M+2*Mprime) *ee;
  lsum:= lsum -     323 * sin(2*M+Mprime) *ee^2;
  lsum:= lsum +     299 * sin(DD+M-Mprime) *ee;
  lsum:= lsum +     294 * sin(2*DD+3*Mprime);

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

  l:= (simp_angle(Lprime + 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 +  280602 * sin(Mprime+F);
  bsum:= bsum +  277693 * sin(Mprime-F); // equazione del centro
  bsum:= bsum +  173237 * sin(2*DD-F); // grande ineguaglianza
  bsum:= bsum +   55413 * sin(2*DD-Mprime+F); // evezione
  bsum:= bsum +   46271 * sin(2*DD-Mprime-F); // evezione (2)
  bsum:= bsum +   32573 * sin(2*DD+F);
  bsum:= bsum +   17198 * sin(2*Mprime+F);
  bsum:= bsum +    9266 * sin(2*DD+Mprime-F);
  bsum:= bsum +    8822 * sin(2*Mprime-F);
  bsum:= bsum +    8216 * sin(2*DD-M-F) *ee;
  bsum:= bsum +    4324 * sin(2*DD-2*Mprime-F);
  bsum:= bsum +    4200 * sin(2*DD+Mprime+F);
  bsum:= bsum -    3359 * sin(2*DD+M-F) *ee;
  bsum:= bsum +    2463 * sin(2*DD-M-Mprime+F) *ee;
  bsum:= bsum +    2211 * sin(2*DD-M+F) *ee;
  bsum:= bsum +    2065 * sin(2*DD-M-Mprime-F) *ee;
  bsum:= bsum -    1870 * sin(M-Mprime-F) *ee;
  bsum:= bsum +    1828 * sin(4*DD-Mprime-F);
  bsum:= bsum -    1794 * sin(M+F) *ee;
  bsum:= bsum -    1749 * sin(3*F);
  bsum:= bsum -    1565 * sin(M-Mprime+F) *ee;
  bsum:= bsum -    1491 * sin(DD+F);
  bsum:= bsum -    1475 * sin(M+Mprime+F) *ee;
  bsum:= bsum -    1410 * sin(M+Mprime-F) *ee;
  bsum:= bsum -    1344 * sin(M-F) *ee;
  bsum:= bsum -    1335 * sin(DD-F);
  bsum:= bsum +    1107 * sin(3*Mprime+F);
  bsum:= bsum +    1021 * sin(4*DD-F);
  bsum:= bsum +     833 * sin(4*DD-Mprime+F);
  bsum:= bsum +     777 * sin(Mprime-3*F);
  bsum:= bsum +     671 * sin(4*DD-2*Mprime+F);
  bsum:= bsum +     607 * sin(2*DD-3*F);
  bsum:= bsum +     596 * sin(2*DD+2*Mprime-F);
  bsum:= bsum +     491 * sin(2*DD-M+Mprime-F) *ee;
  bsum:= bsum -     451 * sin(2*DD-2*Mprime+F);
  bsum:= bsum +     439 * sin(3*Mprime-F);
  bsum:= bsum +     422 * sin(2*DD+2*Mprime+F);
  bsum:= bsum +     421 * sin(2*DD-3*Mprime-F);
  bsum:= bsum -     366 * sin(2*DD+M-Mprime+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+Mprime+F) *ee;
  bsum:= bsum +     302 * sin(2*DD-2*M-F) *ee^2;
  bsum:= bsum -     283 * sin(Mprime+3*F);
  bsum:= bsum -     229 * sin(2*DD+M+Mprime-F) *ee;
  bsum:= bsum +     223 * sin(DD+M-F) *ee;
  bsum:= bsum +     223 * sin(DD+M+F) *ee;
  bsum:= bsum -     220 * sin(M-2*Mprime-F) *ee;
  bsum:= bsum -     220 * sin(2*DD+M-Mprime-F) *ee;
  bsum:= bsum -     185 * sin(DD+Mprime+F);
  bsum:= bsum +     181 * sin(2*DD-M-2*Mprime-F) *ee;
  bsum:= bsum -     177 * sin(M+2*Mprime+F) *ee;
  bsum:= bsum +     176 * sin(4*DD-2*Mprime-F);
  bsum:= bsum +     166 * sin(4*DD-M-Mprime-F) *ee;
  bsum:= bsum -     164 * sin(DD+Mprime-F);
  bsum:= bsum +     132 * sin(4*DD+Mprime-F);
  bsum:= bsum -     119 * sin(DD-Mprime-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(Lprime))+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(Lprime-Mprime)-115*sin(Lprime+Mprime);
  b := →HMS(bsum/1000000); // β div by 1 milion

  // correzioni per la distanza dalla Terra (Δ)
  rsum:=0;
  rsum:= rsum -20905355 * cos(Mprime); // coseni
  rsum:= rsum - 3699111 * cos(2*DD-Mprime); // evezione
  rsum:= rsum - 2955968 * cos(2*DD); // variazione
  rsum:= rsum -  569925 * cos(2*Mprime); // 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*Mprime);
  rsum:= rsum -  152138 * cos(2*DD-Mprime-M) *ee;
  rsum:= rsum -  170733 * cos(2*DD+Mprime);
  rsum:= rsum -  204586 * cos(2*DD-M) *ee;
  rsum:= rsum -  129620 * cos(M-Mprime) *ee;
  rsum:= rsum +  108743 * cos(DD);
  rsum:= rsum +  104755 * cos(Mprime+M) *ee;
  rsum:= rsum +   10321 * cos(2*DD-2*F);
  rsum:= rsum +   79661 * cos(Mprime-2*F);
  rsum:= rsum -   34782 * cos(4*DD-Mprime);
  rsum:= rsum -   23210 * cos(3*Mprime);
  rsum:= rsum -   21636 * cos(4*DD-2*Mprime);
  rsum:= rsum +   24208 * cos(2*DD+M-Mprime) *ee;
  rsum:= rsum +   30824 * cos(2*DD+M) *ee;
  rsum:= rsum -    8379 * cos(DD-Mprime);
  rsum:= rsum -   16675 * cos(DD+M) *ee;
  rsum:= rsum -   12831 * cos(2*DD-M+Mprime) *ee;
  rsum:= rsum -   10445 * cos(2*DD+2*Mprime);
  rsum:= rsum -   11650 * cos(4*DD);
  rsum:= rsum +   14403 * cos(2*DD-3*Mprime);
  rsum:= rsum -    7003 * cos(M-2*Mprime) *ee;
  rsum:= rsum +   10056 * cos(2*DD-M-2*Mprime) * ee;
  rsum:= rsum +    6322 * cos(DD+Mprime);
  rsum:= rsum -    9884 * cos(2*DD-2*M) * ee^2;
  rsum:= rsum +    5751 * cos(M+2*Mprime) *ee;
  rsum:= rsum -    4950 * cos(2*DD-2*M-Mprime) *ee^2;
  rsum:= rsum +    4130 * cos(2*DD+Mprime-2*F);
  rsum:= rsum -    3958 * cos(4*DD-M-Mprime) *ee;
  rsum:= rsum +    3258 * cos(3*DD-Mprime);
  rsum:= rsum +    2616 * cos(2*DD+M+Mprime) *ee;
  rsum:= rsum -    1897 * cos(4*DD-M-2*Mprime) *ee;
  rsum:= rsum -    2117 * cos(2*M-Mprime) * ee^2;
  rsum:= rsum +    2354 * cos(2*DD+2*M-Mprime) *ee^2;
  rsum:= rsum -    1423 * cos(4*DD+Mprime);
  rsum:= rsum -    1117 * cos(4*Mprime);
  rsum:= rsum -    1571 * cos(4*DD-M) *ee;
  rsum:= rsum -    1739 * cos(DD-2*Mprime);
  rsum:= rsum -    4421 * cos(2*Mprime-2*F);
  rsum:= rsum +    1165 * cos(2*M+Mprime) *ee^2;
  rsum:= rsum +    8752 * cos(2*DD-Mprime-2*F);

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

  ts:= calcTS(lista);  // calcola tempo siderale
  Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST
  Hd:=simp_angle(→HMS(Hd));
  H:= Hd/15; // in HMS

  DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
  SAD:= (90 + DA); // Semiarco Diurno
  // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
  
  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)
  // s:= 358473400/d;  // semidiametro (in " arcsec)

  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 Calante"; END;
  IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
  IF ph>100 AND ph<170 THEN fase:="Gibbosa Crescente"; 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
  moonlist:= {Lprime, l, b, lapp, α, δ, az, h, Mprime, ω, Ω, d, DD, P, s, H, DA, SAD,ph, fase, eta, ill};
RETURN moonlist;
END;

moon(lista)
BEGIN
  LOCAL ε, ts, JD, JDE, segno, Lprime, b;
  LOCAL l, Ω, ω, Mprime, lapp, α, δ;
  LOCAL az, h, sorg, tram, culm, Hd, DA, SAD;
  LOCAL ph, fase, eta, DD, d, s, ill;
  LOCAL dep, temp, listamoon, listamoontran, listamoonothers, ch;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
  LOCAL listaARh0, transito, Hor;

  ε:= epsilon(lista);
  ts:= calcTS(lista); //
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  dep:= deltaEpsilonPsi(lista);
  Lprime:= moonlist(1); l:= moonlist(2);
  b:= moonlist(3);
  α:= moonlist(5); δ:= moonlist(6); 
  az:= moonlist(7); h:= moonlist(8); 
  ω:= moonlist(10); Ω:= moonlist(11); 
  lapp:= moonlist(4);  Mprime:= moonlist(9);
  d:= moonlist(12); DD:= moonlist(13);
  P:= moonlist(14); s:= moonlist(15);
  Hor:= moonlist(16); Hd:= 15*moonlist(16);
  DA:= moonlist(17); SAD:= moonlist(18);
  ph:= moonlist(19); fase:= moonlist(20);

  eta:= moonlist(21); ill:= moonlist(22);
  segno:= zodiaco(l);

  // transit routine
  h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)
  // h0:= 0°7′30″; // medium value: 0.125°
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  PRINT;
  PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
  PRINT("L Long media " + trunc_angle(Lprime));
  PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("β Latitudine " + trunc_angle(b));
  PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
  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("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(" ");
WAIT();
PRINT();
  PRINT("ω Perigeo medio " + trunc_angle(ω));
  PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
  PRINT("M anomalia vera " + trunc_angle(Mprime));
  PRINT("D Elongazione " + trunc_angle(DD));
  PRINT(" ");
  PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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("Julian day of Ephemeris " + JDE);
  PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
  PRINT(" ");
WAIT();
PRINT();
  PRINT("Transita a: " + trunc_angle(culm MOD 24 ) + " 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("ε obliquità ecl. " + ε);
  PRINT("Δψ Nutazione long "+ dep(1));
  PRINT("Δε Nutazione obl. " + dep(2));
  PRINT(" ");
  PRINT("Fase " + trunc_angle(ph) + " " + fase);
  PRINT("Età " + eta);
  PRINT("Parte illuminata " + ROUND(ill,5) + " " + ROUND(100*ill,2)+"%");
  listamoon:= {→HMS(l), →HMS(b), ROUND(HMS→(d),6), →HMS(α), →HMS(δ), →HMS(az), →HMS(h)};
  listamoonothers:= {→HMS(Lprime), →HMS(lapp), →HMS(Mprime),→HMS(ω),→HMS(Ω), →HMS(P), 
→HMS(s), →HMS(Hor), →HMS(α*15), ph, eta};
  listamoontran:= →HMS({culm, sorg, tram});
  Moon_Effemerids:= listamoon;
  Moon_Transit:= listamoontran;
  Moon_Others:= listamoonothers;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

moonAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= moonCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;

moonAh(lista)
// lista like dateshortlist
BEGIN
  LOCAL lha,alpha,delta,phi, ah, alfa;
  LOCAL temp,longitude, Azimuth, Altitude, ts;
  temp:= sunAlphaDelta(lista);
  alpha:=15*HMS→(temp(1));
  alfa:= →HMS(temp(1)); // α in HMS
  delta:= temp(2);
  phi:=lat;
  ts:= calcTS(lista);
  longitude:= - HMS→(long); // - if long -E
  lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
  // lha:= 15*HMS→(ts(2)-alfa);  // 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;

planetCalc(nameplanet, lista)
BEGIN
  LOCAL ω,Ω, l, v, r, x, y, z;
  LOCAL Ls, Ms, b, xs, ys, b_ec;
  LOCAL x1, y1, z1, l_ec, 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, lapp, dep;
  LOCAL listaplanet, listaplanettran, listaplanetothers, ch;

  JD:= julianDay(lista);
  ε:= epsilon(lista);
  dep:= deltaEpsilonPsi(lista);

  // calc of Sun (as basis)
  sunCalc(lista);
  xs:= sunlist(13); ys:= sunlist(14);
  Ls:= sunlist(1); // Longitudine media
  Ms:= sunlist(9); // 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

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

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

  ts:= calcTS(lista);  // 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

  planetlist:= {L, l, b, lapp, α,δ, az, h, M, Ω, ω, v, x,y, z, x1, y1, z1, 
    r, R, H, zc, dz, DA, SAD, ec, in, a, θ, l_ec, b_ec, Ls, Ms, xs, ys, nameplanet};
  RETURN planetlist;
END;

planetDisplay(nameplanet, lista)
BEGIN
  LOCAL ω,Ω, l, v, r;
  LOCAL x, y, z, x1, y1, z1;
  LOCAL Ls, Ms, b, xs, ys, l_ec, b_ec;
  LOCAL ec, in, a, δ, α, ts, az, h;
  LOCAL DA, SAD, tram, sorg, culm, dayY;
  LOCAL as, ac, zc, dz, θ, mov, segno, ε;
  LOCAL transito, JD, JDE, temp, Hor, Hd;
  LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
  LOCAL listaARh0, lapp;

  ts:= calcTS(lista);
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);

  L:= planetlist(1); l:= planetlist(2);
  b:= planetlist(3); lapp:= planetlist(4);
  α:= planetlist(5); δ:= planetlist(6);
  az:= planetlist(7); h:= planetlist(8);
  M:= planetlist(9); Ω:= planetlist(10);
  ω:= planetlist(11); v:= planetlist(12);
  x := planetlist(13);  y:= planetlist(14);  z:= planetlist(15);
  x1:= planetlist(16); y1:= planetlist(17); z1:= planetlist(18);
  r:= planetlist(19);   R:= planetlist(20);
  Hor:= planetlist(21); Hd:= 15*planetlist(21);
  zc:= planetlist(22); dz:= planetlist(23);
  DA:= planetlist(24); SAD:= planetlist(25);
  ec:= planetlist(26); in:= planetlist(27);
  a:= planetlist(28); θ:= planetlist(29);
  l_ec:= planetlist(30); b_ec:= planetlist(31);
  Ls:= planetlist(32); Ms:= planetlist(33);  // Sun values
  xs:= planetlist(34); ys:= planetlist(35);
  // nameplanet:= planetlist(36);

  segno:= zodiaco(l);

  // transit routine
  h0:= -0°34′; // 'standard' altitude of a planet
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  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(" ");
  PRINT("L Long media " + trunc_angle(L));
  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) + " (UA) ");
  PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
  PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
  PRINT(" ");
  PRINT("ω Perielio " + trunc_angle(ω));
  PRINT("Ω Nodo " + trunc_angle(Ω));
  PRINT("M Anomalia vera " + trunc_angle(M));
WAIT();
PRINT();
  PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));
  PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(simp_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(" ");
  PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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);
WAIT();
PRINT();
  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);
  listaplanet:= {l,b, R, α, δ, az, h};
  listaplanettran:= {culm, sorg, tram};
  listaplanetothers:= {L, lapp, l_ec,b_ec, M, ω,Ω, v, x,y,z, x1,y1,z1, r, Hor, Hd, SAD, DA, mov};
  Planet_Effemerids:= listaplanet;
  Planet_Transit:= listaplanettran;
  Planet_Others:= listaplanetothers;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;

END;

planetAlphaDelta(nameplanet, lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= planetCalc(nameplanet, lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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;



RE: Astronomy: Effemeridi (Ephemeris) - Thomas_Sch - 07-01-2015 08:47 PM

Hi Salvo,

Many thanks for sharing your programs!
They are very useful, and in addition they show nicely the possibilities of the Prime.

edit:
May I express a few wishes or suggestions?
- it would be helpful, if in successive calculations the time in the input fields is updated
- is it possible to build an input dialog for user specifiv local coordinates?
Alternatively, variables or constants for lat and long would be useful at the beginning of the program,
so the 'customization' would be easier.

Again many thanks!
Thomas


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-01-2015 09:12 PM

(07-01-2015 08:47 PM)Thomas_Sch Wrote:  Hi Salvo,

Many thanks for sharing your programs!
They are very useful, and in addition they show nicely the possibilities of the Prime.

Thomas

thank you Thomas,
I'm happy to share, and yes, the Prime are an "open book" of possibilities, very amazing to program

Salvo


RE: Astronomy: Effemeridi (Ephemeris) - Thomas_Sch - 07-01-2015 09:14 PM

(07-01-2015 09:12 PM)salvomic Wrote:  thank you Thomas,
I'm happy to share, and yes, the Prime are an "open book" of possibilities, very amazing to program

Salvo

Sorry, I did some editing to my posting later:
edit:
May I express a few wishes or suggestions?
- it would be helpful, if in successive calculations the time in the input fields is updated
- is it possible to build an input dialog for user specifiv local coordinates?
Alternatively, variables or constants for lat and long would be useful at the beginning of the program,
so the 'customization' would be easier.

Again many thanks!
Thomas


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-01-2015 09:24 PM

(07-01-2015 09:14 PM)Thomas_Sch Wrote:  Sorry, I did some editing to my posting later:
edit:
May I express a few wishes or suggestions?
- it would be helpful, if in successive calculations the time in the input fields is updated
- is it possible to build an input dialog for user specifiv local coordinates?
Alternatively, variables or constants for lat and long would be useful at the beginning of the program,
so the 'customization' would be easier.

ok,thank you; I'll think about these tips.
For now the input works so: at the start of program it calculates the actual date and time and the user can change the default lat/long (for now mine, while I'm still debugging), the program calculates Sun and Moon (to set some variables), then asks to the user for another calculation, and then save the old values (lat, long, date, time) or the other changes the user does.
This is useful to make different calculations with the same data.
I need help to make an option to choose automatically update, if the user chooses it.

It would be also nice to have automatic calculation of DST (but I've no idea how todo this in the Prime)...
I'll think also a way to save long/lat as in the app of Marcel, also if this is for now a program, not an app, that uses AVars...

Salvo


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-01-2015 09:48 PM

...in the meantime this version introduce other values for planets: parallax, illuminated fraction, magnitudo, phase.

EDIT: Thomas, try this with an option to save coordinates at the starting, the program should save them in Long_Lat var in the program variables and load them when used again...
I need a better safeguard, however, then only control if the size is ...2 ;-)

Code:

data();
changeData();
whatToDo();
calcN();
calcTS();
meanSideralTime();
apparentSideralTime();
julianDay();
julianEphemerisDay();
julianDayAtZero();
makeDateList();
makeDateListShort();
deltaEpsilonPsi();
epsilon();
transformEquatorialToEcliptical();
transformEclipticalToEquatorial();
transformEquatorialToHorizontal();
transformEclipticalToHorizontal();
transformHorizontalToEquatorial();
tempo();
precession();
parametri();
transit();
sunCalc();
sun();
sunAh();
sunAlphaDelta();
moon();
moonCalc();
moonAh();
moonAlphaDelta();
planets();
planetCalc();
planetDisplay();
planetAlphaDelta();
Mercurius();
Venus();
Mars();
Jupiter();
Saturn();
Uranus();
Neptunus();
Pluto();
ascendant();
zodiaco();
trunc_angle();
simp_angle();
atan2();
interpolation();

EXPORT Tempi;
EXPORT EpsilonPsi;
EXPORT Ascendant;
EXPORT Sun_Effemerids;
EXPORT Sun_Transit;
EXPORT Sun_Others;
EXPORT Moon_Effemerids;
EXPORT Moon_Transit;
EXPORT Moon_Others;
EXPORT Planet_Effemerids;
EXPORT Planet_Transit;
EXPORT Planet_Others;
EXPORT Long_Lat;

thisday:= 1964.0529;
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,40);
moonlist:= MAKELIST(X,X,1,40);
planetlist:= MAKELIST(X,X,1,40);
planet_const:= MAKELIST(X,X,1,40);

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, DrD
BEGIN
  HAngle:=1;
  data(); // input data and calc Sun
END;

data()
BEGIN
  // initial input
  LOCAL dd, mm, yy, hh, mn, ss, loct, save;
  LOCAL ε, JD, JDE;
  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);
  IF (SIZE(Long_Lat)==2) THEN long:= Long_Lat(1); lat:= Long_Lat(2);
  ELSE long:= →HMS(-14°43′41″); lat:= →HMS(36°55′15″); END;
  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}},
  {save,0,{85,15,5}} },
  "Data: Use Shift+9 for H°M′S″",
  {"Month:","Day :","Year :","Local Time (24):", "DST", "UTC", 
  "Long (-E):","Lat (+N)", "delta T", "Save coord?"},
  {"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.)", "Save Lat/Long in a var?"}, 
  {1,1,1901,0, 0,0, 0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, long, lat, 69.3, 1});
  // Ragusa (-14°43′41″, 36°55′15″) - 
  // Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)

  IF save==1 THEN Long_Lat:= {long, lat}; END;

  gmt:= (lstd-utc-(1*dst)) MOD 24; // 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};

  ε:=epsilon(dateshortlist);
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:=(JD-2451545)/36525;
  U:= T/100;
  // sunCalc(dateshortlist);  // set initial data for Sun
  // moonCalc(dateshortlist);
  whatToDo(dateshortlist);
END;


changeData()
BEGIN
  // input to change date, lat, long...
  LOCAL ε, JD, JDE;

  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},{m,D,Y,lstd, dst, utc, long, lat, deltaT});

  gmt:= (lstd-utc-(1*dst)) MOD 24; // 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};

  ε:=epsilon(dateshortlist);
  JD:= julianDay(dateshortlist);
  JDE:= julianEphemerisDay(dateshortlist);
  T:=(JD-2451545)/36525;
  U:= T/100;

  whatToDo(dateshortlist);
END;

whatToDo(lista)
BEGIN
  LOCAL ch;
  CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendant",
    "Quit");
  CASE
  IF ch==1 THEN sunCalc(lista); sun(lista); END;
  IF ch==2 THEN moonCalc(lista); moon(lista); END;
  IF ch==3 THEN planets(lista); END;
  IF ch==4 THEN tempo(lista); END;
  IF ch==5 THEN parametri(lista); END;
  IF ch==6 THEN ascendant(lista); END;
  IF ch==7 THEN HAngle:=0; RETURN; END;
  DEFAULT 
  END; // case
END;


planets(lista)
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), lista);
  planetDisplay(nameplanet(ch),lista); 
END;

epsilon(lista)
BEGIN
LOCAL ε, ε0, dep, eps, JD, JDE, U;
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  T:=(JD-2451545)/36525;
  U:= T/100;
  // ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)
  // 4680.93″ = 1°18′00.93″
  ε0:= 23°26′21.448″ - (HMS→(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(lista);
  eps:= dep(2);
  ε:= ε0+eps; // true obliquity ε
  RETURN ε;
END;

deltaEpsilonPsi(lista)
BEGIN
  // Nutazione Δψ (longitudine), Δε (obliquità)
  LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;
  LOCAL DD;
  // JD:= julianDayAtZero(Y,m,D);
  JDE:= julianEphemerisDay(lista);
  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;

transformEquatorialToEcliptical(lista, alfa, delta)
BEGIN
  LOCAL Alpha,Delta,ε,Lambda,Beta, lambda_beta;
  Alpha:= HMS→(alfa)*15;
  Delta:= HMS→(delta);
  ε:= epsilon(lista);
  Lambda:= atan2(sin(Alpha)*cos(ε)+tan(Delta)*sin(ε),cos(Alpha));
  Beta:= asin(sin(Delta)*cos(ε)-cos(Delta)*sin(ε)*sin(Alpha));
  lambda_beta:={→HMS(Lambda),→HMS(Beta)};
  RETURN (lambda_beta);
END;

transformEclipticalToEquatorial(lista, lambda, beta)
BEGIN
  // trasnform λ and β (long, lat) into α and δ (Asc retta, decl)
  LOCAL ε,Lambda, Beta, alpha, delta;
  Lambda:=HMS→(lambda);
  Beta:=HMS→(beta);
  ε:= epsilon(lista);
  alpha:=atan2(sin(Lambda)*cos(ε)-tan(Beta)*sin(ε),cos(Lambda));
  IF alpha>=360 THEN alpha:=alpha-360 END;
  IF alpha<0 THEN alpha:=alpha+360 END;
  delta:=asin(sin(Beta)*cos(ε)+cos(Beta)*sin(ε)*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(lista, Lambda, Beta, Lha)
BEGIN
  // transform λ and β (long, lat) into azimuth and height
  LOCAL lambda, beta,phi, lha, ε;
  LOCAL alpha,delta,azimuth,altitude;
  lambda:=HMS→(Lambda); // longitudine astro
  beta:=HMS→(Beta); // latitudine astro
  ε:= epsilon(lista);
  lha:=HMS→(Lha);  // angolo orario
  phi:=HMS→(lat); // lat geo
  alpha:=atan2(SIN(lambda)*COS(ε)-TAN(beta)*SIN(ε),COS(lambda));
  delta:=asin(sin(beta)*cos(ε)+cos(beta)*sin(ε)*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;

transformHorizontalToEquatorial(az, h)
BEGIN
  LOCAL azimuth,altitude,phi,lha,delta, ang_delta;
  azimuth:=HMS→(az);
  altitude:=HMS→(h);
  phi:=HMS→(lat);
  lha:= atan2(sin(azimuth),cos(azimuth)*sin(phi)+tan(altitude)*cos(phi));
  delta:= asin(sin(phi)*sin(altitude)-cos(phi)*cos(altitude)*cos(azimuth));
  lha:=→HMS(lha);
  delta:=→HMS(delta);
  ang_delta:={lha,delta};
  RETURN ang_delta;
END;

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

calcTS(lista)
BEGIN
  LOCAL T, TS, TSd, N0, TSL, TSapp, ε;
  LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  JD0:= julianDayAtZero(Y,m,D);
  dep:= deltaEpsilonPsi(lista);
  psi:=dep(1);
  eps:=dep(2);
  T:= (JD-2451545.0)/36525;
  θ0:= meanSideralTime(lista); // 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(lista); // apparent ST (ST at Greenwitch at local our)
  TSL:= apparentSideralTime(lista) - 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(lista);
  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(lista);
  MST:= meanSideralTime(lista);
  dep:= deltaEpsilonPsi(lista);
  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(lista)
BEGIN
  LOCAL JD, JDE, TSid, dep, psi, eps;
  LOCAL ε, listatempi, ch;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  TSid:= calcTS(lista);
  dep:= deltaEpsilonPsi(lista);
  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(" ");
WAIT();
PRINT();
  PRINT("Δψ Nutazione (longit.) "+ psi);
  PRINT("Δε Nutazione (obl.) "+ eps);
  PRINT("ε Obliquità Ecl. "+ ε);
  PRINT(" ");
  PRINT("Time " + lstd + " GMT " + gmt);
  listatempi:= {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)) };
  Tempi:= listatempi;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

transit(lista, listaARh0)
BEGIN
  // lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
  // α, δ (AR, dec), height, h0 = standard height of body
  LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;
  LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;
  LOCAL alfa1, alfa3, delta1, delta3, h0, alpha, del, n;
  α:= listaARh0(2); // AR for day
  δ:= listaARh0(5); // decl for day
  alfa1:= listaARh0(1); // AR day-1
  alfa3:= listaARh0(3); // AR day+1
  delta1:= listaARh0(4); // decl day-1
  delta3:= listaARh0(6); // decl day+1
  h0:= listaARh0(7); // standard height
  temp:= calcTS(lista);
  θ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
  n:= m0 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*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
  n:= m1 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*15);
  Hor:=simp_angle(→HMS(Hor));
  h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
  deltaM:= HMS→((h-h0)/(360*cos(del)*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
  n:= m2 + deltaT/86400;
  alpha := interpolation({alfa1, α, alfa3},n);
  del   := interpolation({delta1, δ, delta3},n);
  Hor:= HMS→(tsideral - long - alpha*15); 
  Hor:=simp_angle(→HMS(Hor));
  h:=asin(sin(lat)*sin(del)+cos(lat)*cos(del)*cos(Hor));
  deltaM:= HMS→((h-h0)/(360*cos(del)*cos(lat)*sin(Hor)));
  tram:= ((m2 + deltaM) * 24); // setting

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

parametri(lista)
BEGIN
  LOCAL dep, prec, α, δ, ε, listaparametri, ch, temp;
  dep:= deltaEpsilonPsi(lista);
  temp:= sunAlphaDelta(lista);
  α:= sunlist(5); δ:= sunlist(6);
  prec:= precession(α, δ);
  ε:= epsilon(lista);
  // 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));
  listaparametri:= {ε, dep(1), dep(2), prec(1), prec(2)}; // ε, Δψ, Δε, Δα, Δδ
  EpsilonPsi:= listaparametri;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

ascendant(lista)
BEGIN
  LOCAL ASC, ts, segno, ε, listaascendente, ch, temp;
  LOCAL tsb;
  ε:= epsilon(lista);
  ts:= calcTS(lista);
  tsb:= ts(1)+lstd-dst;  // ST of the birth
  // temp:= sin(ts(2)*15)*cos(ε)+tan(lat)*sin(ε);
  temp:= sin(tsb*15)*cos(ε)+tan(lat)*sin(ε);
  ASC:= atan2(-cos(tsb*15), temp);
  // IF temp<0 THEN ASC:= ASC+180 ELSE ASC:= ASC+360; END;
  IF ASC<180 THEN ASC:=ASC+180 ELSE ASC:= ASC-180; END;
  ASC:= simp_angle(ASC);
  segno:= zodiaco(ASC);
  listaascendente:= {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno};
  Ascendant:= listaascendente;
    RECT_P();
   TEXTOUT_P("Ascendant " + trunc_angle(→HMS(ASC)), 25,50);
   TEXTOUT_P(" " + trunc_angle(→HMS(ASC) MOD 30), 25,70);
   TEXTOUT_P("in " + segno, 100,70);

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
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(lista)
// type dateshortlist
BEGIN
  LOCAL ω, l, v, r, x, y, dz, ε;
  LOCAL α, δ, dayY, ts, az, z, h;
  LOCAL SAD, DA, Hd, lt;
  LOCAL as, ac, zc, h0, H0, temp;
  LOCAL ec, a, u, M1, M2, j, t, tau;
  LOCAL JD, JDE, c, lapp, Ω, θ0;
  LOCAL dep, deltapsi, longitude, b;
  LOCAL lambdaprime, delta_long, delta_beta, correction;
  ε:= epsilon(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  a:= 1.00000023; // semiasse maggiore
  dep:= deltaEpsilonPsi(lista);
  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.003032028*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
  //M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-(T^3)/24490000); // anom vera Meeus
  M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // Meeus p 308
  Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)
  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
  // 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(lista);  // calcola tempo siderale

  temp:= transformEquatorialToEcliptical(lista, α, δ); // β latitudine
  b:= →HMS(temp(2));

  lambdaprime:= l-1.397*T-0.00031*T^2;
  delta_long=-(0.09033/3600);
  delta_beta:=0.03916/3600*(COS(lambdaprime)-SIN(lambdaprime));
  b:= b+delta_beta;
  l:= l+delta_long;

  correction:=-20.4898/r;
  //lapp:= simp_angle(l - 0.00569-0.00478*sin(Ω)); // λ Longitudine apparente (ref true equinox)
  lapp:=l+deltapsi+correction/3600;

  lt:= simp_angle(l + 180); // Longitudine eliocentrica della Terra

  // 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:= transformEquatorialToHorizontal(Hd, lat, δ);
  az:= HMS→(temp(1))+180;  // Azimuth
  IF az>360 THEN az:= az-360 END;
  az:= →HMS(az);
  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

  sunlist:= {L,ω, l, lapp, α,δ, az, h,  M, Ω, v, c,x,y,r, H, zc, dz, DA, SAD, E, lt, ec, b};
  RETURN sunlist;
END;

sun(lista)
BEGIN
  LOCAL ω, l, b, v, r, x, y, dz; 
  LOCAL ε;
  LOCAL α, δ, dayY, ts, az, z, h;
  LOCAL SAD, DA, sorg, tram, culm, Hd, h0, lt;
  LOCAL as, ac, zc, segno, JD, JDE, v, c;
  LOCAL Ω, lapp, ec, listasun, listasuntran, listasunothers; 
  LOCAL ch, transito, listaARh0, temp, Hor, dep;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3;
  ε:=  epsilon(lista);
  L:=  sunlist(1);    ω:= sunlist(2); 
  l:=  sunlist(3); lapp:= sunlist(4);
  α:=  sunlist(5);    δ:= sunlist(6); 
  az:= sunlist(7);    h:= sunlist(8);
  M:=  sunlist(9);    Ω:= sunlist(10);   
  v:=  sunlist(11);   c:= sunlist(12); 
  x:=  sunlist(13);   y:= sunlist(14);
  r:=  sunlist(15); Hor:= sunlist(16); 
  Hd:= 15*sunlist(16); 
  zc:= sunlist(17);  dz:= sunlist(18);
  DA:= sunlist(19); SAD:= sunlist(20); 
  E:=  sunlist(21);  lt:= sunlist(22);
  ec:= sunlist(23);
  b:=  sunlist(24);

  ts:= calcTS(lista);
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  dep:= deltaEpsilonPsi(lista);
  segno:= zodiaco(l);

   // 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
  // to pass lista like dateshortlis, listaARh0 {α1,α2,α3,δ1,δ2,δ3,h0}
  h0:= -0°50′; // 'standard' altitude of Sun
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= sunAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  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("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("λ longitudine apparente " + trunc_angle(lapp));
  PRINT("β latitudine " + →HMS(b));  // b never > 1.2" (Meeus) 
  PRINT("Δ 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(" ");
WAIT();
PRINT();
  PRINT("Ω aberr, nutaz nodo "+ Ω);
  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("Ang. orario "+ trunc_angle(Hor) + " " + 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(" ");
WAIT();
PRINT();
  PRINT("Transita a: " + trunc_angle(culm MOD 24) + " 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("ε obliquità ecl. " + ε);
  PRINT("ec eccentricità "+ ec);
  PRINT("Δψ Nutazione long "+ dep(1));
  PRINT("Δε Nutazione obl. " + dep(2));
  PRINT(" ");
  PRINT("Long elioc. Terra "+ trunc_angle(lt));

  l:= →HMS(l); lapp:= →HMS(lapp);
  culm:= →HMS(culm); sorg:= →HMS(sorg); tram:= →HMS(tram);
  listasun:= {l, lapp, b, α, δ, az, h, r};
  listasuntran:= {culm, sorg, tram};
  listasunothers:= {L, ω, Ω, M,v, c, x, y, Hor, lt, ε, ec};
  Sun_Effemerids:= listasun;
  Sun_Transit:= listasuntran;
  Sun_Others:= listasunothers;
  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

sunAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= sunCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;

sunAh(lista)
// lista like dateshortlist
BEGIN
  LOCAL lha,alpha,delta,phi, ah, alfa;
  LOCAL temp,longitude, Azimuth, Altitude, ts;
  temp:= sunAlphaDelta(lista);
  alpha:=15*HMS→(temp(1));
  alfa:= →HMS(temp(1)); // α in HMS
  delta:= temp(2);
  phi:=lat;
  ts:= calcTS(lista);
  longitude:= HMS→(long); // - if long -E
  lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
  // lha:= 15*HMS→(ts(2)-alfa);  // 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;

moonCalc(lista)
BEGIN
  LOCAL ω,Ω, Lprime, Mprime, l, lsum, bsum, rsum; 
  LOCAL v, r, x, y, dep, ε;
  LOCAL b, P, d, s, temp, tau;
  LOCAL α,δ, ts, Hd, zc, dz, h, az;
  LOCAL i, as, ac, DA, SAD, DD;
  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(lista);
  N:= calcN();
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  T:= (JDE-2451545.0)/36525; // use JDE
  tau:= (JDE-2451545.0)/365250; // millenni  - JDE
  dep:= deltaEpsilonPsi(lista);
  // Sun
  L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);
  M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+(T^3)/24490000); // mean anomaly Sun, Meeus  
  // end Sun
  // L1:= (33.231 + 13.17639653*N) MOD 360; // long. media
  Lprime:= 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
  Mprime:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media
  ω:= (Lprime - Mprime) MOD 360; // perigeo medio
  // DD:= Lprime - L; F:= Lprime-Ω; // 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(Mprime); // calc Moon longitude
  lsum:= lsum + 1274027 * sin(2*DD-Mprime); // evezione
  lsum:= lsum +  658314 * sin(2*DD); // variazione
  lsum:= lsum +  213618 * sin(2*Mprime); // 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*Mprime);
  lsum:= lsum +   57066 * sin(2*DD-Mprime-M) *ee;
  lsum:= lsum +   53322 * sin(2*DD+Mprime);
  lsum:= lsum +   45758 * sin(2*DD-M) *ee;
  lsum:= lsum -   40923 * sin(M-Mprime) *ee;
  lsum:= lsum -   34720 * sin(DD);
  lsum:= lsum -   30383 * sin(Mprime+M) *ee;
  lsum:= lsum +   15327 * sin(2*DD-2*F);
  lsum:= lsum -   12528 * sin(Mprime+2*F);
  lsum:= lsum +   10980 * sin(Mprime-2*F);
  lsum:= lsum +   10675 * sin(4*DD-Mprime);
  lsum:= lsum +   10034 * sin(3*Mprime);
  lsum:= lsum +    8548 * sin(4*DD-2*Mprime);
  lsum:= lsum -    7888 * sin(2*DD+M-Mprime) *ee;
  lsum:= lsum -    6766 * sin(2*DD+M) *ee;
  lsum:= lsum -    5163 * sin(DD-Mprime);
  lsum:= lsum +    4987 * sin(DD+M) *ee;
  lsum:= lsum +    4036 * sin(2*DD-M+Mprime) *ee;
  lsum:= lsum +    3994 * sin(2*DD+2*Mprime);
  lsum:= lsum +    3861 * sin(4*DD);
  lsum:= lsum +    3665 * sin(2*DD-3*Mprime);
  lsum:= lsum -    2689 * sin(M-2*Mprime) *ee;
  lsum:= lsum -    2602 * sin(2*DD-Mprime+2*F);
  lsum:= lsum +    2390 * sin(2*DD-M-2*Mprime) *ee;
  lsum:= lsum -    2348 * sin(DD+Mprime);
  lsum:= lsum +    2236 * sin(2*DD-2*M) *ee^2;
  lsum:= lsum -    2120 * sin(M+2*Mprime) *ee;
  lsum:= lsum -    2069 * sin(2*M) *ee^2;
  lsum:= lsum +    2048 * sin(2*DD-2*M-Mprime) *ee^2;
  lsum:= lsum -    1773 * sin(2*DD+Mprime-2*F);
  lsum:= lsum -    1595 * sin(2*DD+2*F);
  lsum:= lsum +    1215 * sin(4*DD-M-Mprime) *ee;
  lsum:= lsum -    1110 * sin(2*Mprime+2*F);
  lsum:= lsum -     892 * sin(3*DD-Mprime);
  lsum:= lsum -     810 * sin(2*DD+M+Mprime) *ee;
  lsum:= lsum +     759 * sin(4*DD-M-2*Mprime) *ee;
  lsum:= lsum -     713 * sin(2*M-Mprime) * ee^2;
  lsum:= lsum -     700 * sin(2*DD+2*M-Mprime) *ee^2;
  lsum:= lsum +     691 * sin(2*DD+M-2*Mprime) *ee;
  lsum:= lsum +     596 * sin(2*DD-M-2*F) *ee;
  lsum:= lsum +     549 * sin(4*DD+Mprime);
  lsum:= lsum +     537 * sin(4*Mprime);
  lsum:= lsum +     520 * sin(4*DD-M) *ee;
  lsum:= lsum -     487 * sin(DD-2*Mprime);
  lsum:= lsum -     399 * sin(2*DD+M-2*F) *ee;
  lsum:= lsum -     381 * sin(2*Mprime-2*F);
  lsum:= lsum +     351 * sin(DD+M+Mprime) *ee;
  lsum:= lsum -     340 * sin(3*DD-2*Mprime);
  lsum:= lsum +     330 * sin(4*DD-3*Mprime);
  lsum:= lsum +     327 * sin(2*DD-M+2*Mprime) *ee;
  lsum:= lsum -     323 * sin(2*M+Mprime) *ee^2;
  lsum:= lsum +     299 * sin(DD+M-Mprime) *ee;
  lsum:= lsum +     294 * sin(2*DD+3*Mprime);

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

  l:= (simp_angle(Lprime + 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 +  280602 * sin(Mprime+F);
  bsum:= bsum +  277693 * sin(Mprime-F); // equazione del centro
  bsum:= bsum +  173237 * sin(2*DD-F); // grande ineguaglianza
  bsum:= bsum +   55413 * sin(2*DD-Mprime+F); // evezione
  bsum:= bsum +   46271 * sin(2*DD-Mprime-F); // evezione (2)
  bsum:= bsum +   32573 * sin(2*DD+F);
  bsum:= bsum +   17198 * sin(2*Mprime+F);
  bsum:= bsum +    9266 * sin(2*DD+Mprime-F);
  bsum:= bsum +    8822 * sin(2*Mprime-F);
  bsum:= bsum +    8216 * sin(2*DD-M-F) *ee;
  bsum:= bsum +    4324 * sin(2*DD-2*Mprime-F);
  bsum:= bsum +    4200 * sin(2*DD+Mprime+F);
  bsum:= bsum -    3359 * sin(2*DD+M-F) *ee;
  bsum:= bsum +    2463 * sin(2*DD-M-Mprime+F) *ee;
  bsum:= bsum +    2211 * sin(2*DD-M+F) *ee;
  bsum:= bsum +    2065 * sin(2*DD-M-Mprime-F) *ee;
  bsum:= bsum -    1870 * sin(M-Mprime-F) *ee;
  bsum:= bsum +    1828 * sin(4*DD-Mprime-F);
  bsum:= bsum -    1794 * sin(M+F) *ee;
  bsum:= bsum -    1749 * sin(3*F);
  bsum:= bsum -    1565 * sin(M-Mprime+F) *ee;
  bsum:= bsum -    1491 * sin(DD+F);
  bsum:= bsum -    1475 * sin(M+Mprime+F) *ee;
  bsum:= bsum -    1410 * sin(M+Mprime-F) *ee;
  bsum:= bsum -    1344 * sin(M-F) *ee;
  bsum:= bsum -    1335 * sin(DD-F);
  bsum:= bsum +    1107 * sin(3*Mprime+F);
  bsum:= bsum +    1021 * sin(4*DD-F);
  bsum:= bsum +     833 * sin(4*DD-Mprime+F);
  bsum:= bsum +     777 * sin(Mprime-3*F);
  bsum:= bsum +     671 * sin(4*DD-2*Mprime+F);
  bsum:= bsum +     607 * sin(2*DD-3*F);
  bsum:= bsum +     596 * sin(2*DD+2*Mprime-F);
  bsum:= bsum +     491 * sin(2*DD-M+Mprime-F) *ee;
  bsum:= bsum -     451 * sin(2*DD-2*Mprime+F);
  bsum:= bsum +     439 * sin(3*Mprime-F);
  bsum:= bsum +     422 * sin(2*DD+2*Mprime+F);
  bsum:= bsum +     421 * sin(2*DD-3*Mprime-F);
  bsum:= bsum -     366 * sin(2*DD+M-Mprime+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+Mprime+F) *ee;
  bsum:= bsum +     302 * sin(2*DD-2*M-F) *ee^2;
  bsum:= bsum -     283 * sin(Mprime+3*F);
  bsum:= bsum -     229 * sin(2*DD+M+Mprime-F) *ee;
  bsum:= bsum +     223 * sin(DD+M-F) *ee;
  bsum:= bsum +     223 * sin(DD+M+F) *ee;
  bsum:= bsum -     220 * sin(M-2*Mprime-F) *ee;
  bsum:= bsum -     220 * sin(2*DD+M-Mprime-F) *ee;
  bsum:= bsum -     185 * sin(DD+Mprime+F);
  bsum:= bsum +     181 * sin(2*DD-M-2*Mprime-F) *ee;
  bsum:= bsum -     177 * sin(M+2*Mprime+F) *ee;
  bsum:= bsum +     176 * sin(4*DD-2*Mprime-F);
  bsum:= bsum +     166 * sin(4*DD-M-Mprime-F) *ee;
  bsum:= bsum -     164 * sin(DD+Mprime-F);
  bsum:= bsum +     132 * sin(4*DD+Mprime-F);
  bsum:= bsum -     119 * sin(DD-Mprime-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(Lprime))+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(Lprime-Mprime)-115*sin(Lprime+Mprime);
  b := →HMS(bsum/1000000); // β div by 1 milion

  // correzioni per la distanza dalla Terra (Δ)
  rsum:=0;
  rsum:= rsum -20905355 * cos(Mprime); // coseni
  rsum:= rsum - 3699111 * cos(2*DD-Mprime); // evezione
  rsum:= rsum - 2955968 * cos(2*DD); // variazione
  rsum:= rsum -  569925 * cos(2*Mprime); // 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*Mprime);
  rsum:= rsum -  152138 * cos(2*DD-Mprime-M) *ee;
  rsum:= rsum -  170733 * cos(2*DD+Mprime);
  rsum:= rsum -  204586 * cos(2*DD-M) *ee;
  rsum:= rsum -  129620 * cos(M-Mprime) *ee;
  rsum:= rsum +  108743 * cos(DD);
  rsum:= rsum +  104755 * cos(Mprime+M) *ee;
  rsum:= rsum +   10321 * cos(2*DD-2*F);
  rsum:= rsum +   79661 * cos(Mprime-2*F);
  rsum:= rsum -   34782 * cos(4*DD-Mprime);
  rsum:= rsum -   23210 * cos(3*Mprime);
  rsum:= rsum -   21636 * cos(4*DD-2*Mprime);
  rsum:= rsum +   24208 * cos(2*DD+M-Mprime) *ee;
  rsum:= rsum +   30824 * cos(2*DD+M) *ee;
  rsum:= rsum -    8379 * cos(DD-Mprime);
  rsum:= rsum -   16675 * cos(DD+M) *ee;
  rsum:= rsum -   12831 * cos(2*DD-M+Mprime) *ee;
  rsum:= rsum -   10445 * cos(2*DD+2*Mprime);
  rsum:= rsum -   11650 * cos(4*DD);
  rsum:= rsum +   14403 * cos(2*DD-3*Mprime);
  rsum:= rsum -    7003 * cos(M-2*Mprime) *ee;
  rsum:= rsum +   10056 * cos(2*DD-M-2*Mprime) * ee;
  rsum:= rsum +    6322 * cos(DD+Mprime);
  rsum:= rsum -    9884 * cos(2*DD-2*M) * ee^2;
  rsum:= rsum +    5751 * cos(M+2*Mprime) *ee;
  rsum:= rsum -    4950 * cos(2*DD-2*M-Mprime) *ee^2;
  rsum:= rsum +    4130 * cos(2*DD+Mprime-2*F);
  rsum:= rsum -    3958 * cos(4*DD-M-Mprime) *ee;
  rsum:= rsum +    3258 * cos(3*DD-Mprime);
  rsum:= rsum +    2616 * cos(2*DD+M+Mprime) *ee;
  rsum:= rsum -    1897 * cos(4*DD-M-2*Mprime) *ee;
  rsum:= rsum -    2117 * cos(2*M-Mprime) * ee^2;
  rsum:= rsum +    2354 * cos(2*DD+2*M-Mprime) *ee^2;
  rsum:= rsum -    1423 * cos(4*DD+Mprime);
  rsum:= rsum -    1117 * cos(4*Mprime);
  rsum:= rsum -    1571 * cos(4*DD-M) *ee;
  rsum:= rsum -    1739 * cos(DD-2*Mprime);
  rsum:= rsum -    4421 * cos(2*Mprime-2*F);
  rsum:= rsum +    1165 * cos(2*M+Mprime) *ee^2;
  rsum:= rsum +    8752 * cos(2*DD-Mprime-2*F);

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

  ts:= calcTS(lista);  // calcola tempo siderale
  Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST
  Hd:=simp_angle(→HMS(Hd));
  H:= Hd/15; // in HMS

  DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale
  SAD:= (90 + DA); // Semiarco Diurno
  // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno
  
  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)
  // s:= 358473400/d;  // semidiametro (in " arcsec)

  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 Calante"; END;
  IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;
  IF ph>100 AND ph<170 THEN fase:="Gibbosa Crescente"; 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
  moonlist:= {Lprime, l, b, lapp, α, δ, az, h, Mprime, ω, Ω, d, DD, P, s, H, DA, SAD,ph, fase, eta, ill};
RETURN moonlist;
END;

moon(lista)
BEGIN
  LOCAL ε, ts, JD, JDE, segno, Lprime, b;
  LOCAL l, Ω, ω, Mprime, lapp, α, δ;
  LOCAL az, h, sorg, tram, culm, Hd, DA, SAD;
  LOCAL ph, fase, eta, DD, d, s, ill;
  LOCAL dep, temp, listamoon, listamoontran, listamoonothers, ch;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
  LOCAL listaARh0, transito, Hor;

  ε:= epsilon(lista);
  ts:= calcTS(lista); //
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);
  dep:= deltaEpsilonPsi(lista);
  Lprime:= moonlist(1); l:= moonlist(2);
  b:= moonlist(3);
  α:= moonlist(5); δ:= moonlist(6); 
  az:= moonlist(7); h:= moonlist(8); 
  ω:= moonlist(10); Ω:= moonlist(11); 
  lapp:= moonlist(4);  Mprime:= moonlist(9);
  d:= moonlist(12); DD:= moonlist(13);
  P:= moonlist(14); s:= moonlist(15);
  Hor:= moonlist(16); Hd:= 15*moonlist(16);
  DA:= moonlist(17); SAD:= moonlist(18);
  ph:= moonlist(19); fase:= moonlist(20);

  eta:= moonlist(21); ill:= moonlist(22);
  segno:= zodiaco(l);

  // transit routine
  h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)
  // h0:= 0°7′30″; // medium value: 0.125°
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= moonAlphaDelta(makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  culm:= transito(1);
  sorg:= transito(2);
  tram:= transito(3);
  // end transit

  PRINT;
  PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);
  PRINT("L Long media " + trunc_angle(Lprime));
  PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("β Latitudine " + trunc_angle(b));
  PRINT("λ Longitudine apparente "+ trunc_angle(lapp));
  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("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");
  PRINT("Semidiametro "+  trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,3) + " km");
  PRINT(" ");
WAIT();
PRINT();
  PRINT("ω Perigeo medio " + trunc_angle(ω));
  PRINT("Ω Nodo ascendente " + trunc_angle(Ω));
  PRINT("M anomalia vera " + trunc_angle(Mprime));
  PRINT("D Elongazione " + trunc_angle(DD));
  PRINT(" ");
  PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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("Julian day of Ephemeris " + JDE);
  PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));
  PRINT(" ");
WAIT();
PRINT();
  PRINT("Transita a: " + trunc_angle(culm MOD 24 ) + " 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("ε obliquità ecl. " + ε);
  PRINT("Δψ Nutazione long "+ dep(1));
  PRINT("Δε Nutazione obl. " + dep(2));
  PRINT(" ");
  PRINT("π Parallasse " + trunc_angle(P));
  PRINT("Fase " + trunc_angle(ph) + " " + fase);
  PRINT("Età " + eta);
  PRINT("Parte illuminata " + ROUND(ill,5) + " " + ROUND(100*ill,2)+"%");
  listamoon:= {→HMS(l), →HMS(b), ROUND(HMS→(d),6), →HMS(α), →HMS(δ), →HMS(az), →HMS(h)};
  listamoonothers:= {→HMS(Lprime), →HMS(lapp), →HMS(Mprime),→HMS(ω),→HMS(Ω), →HMS(P), 
→HMS(s), →HMS(Hor), →HMS(α*15), ph, eta};
  listamoontran:= →HMS({culm, sorg, tram});
  Moon_Effemerids:= listamoon;
  Moon_Transit:= listamoontran;
  Moon_Others:= listamoonothers;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;
END;

moonAlphaDelta(lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= moonCalc(lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
END;

moonAh(lista)
// lista like dateshortlist
BEGIN
  LOCAL lha,alpha,delta,phi, ah, alfa;
  LOCAL temp,longitude, Azimuth, Altitude, ts;
  temp:= sunAlphaDelta(lista);
  alpha:=15*HMS→(temp(1));
  alfa:= →HMS(temp(1)); // α in HMS
  delta:= temp(2);
  phi:=lat;
  ts:= calcTS(lista);
  longitude:= - HMS→(long); // - if long -E
  lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);
  // lha:= 15*HMS→(ts(2)-alfa);  // 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;

planetCalc(nameplanet, lista)
BEGIN
  LOCAL ω,Ω, l, v, r, x, y, z;
  LOCAL Ls, Ms, b, xs, ys, b_ec;
  LOCAL x1, y1, z1, l_ec, 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, lapp, dep;
  LOCAL listaplanet, listaplanettran, listaplanetothers, ch;

  JD:= julianDay(lista);
  ε:= epsilon(lista);
  dep:= deltaEpsilonPsi(lista);

  // calc of Sun (as basis)
  sunCalc(lista);
  xs:= sunlist(13); ys:= sunlist(14);
  Ls:= sunlist(1); // Longitudine media
  Ms:= sunlist(9); // 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 eclittic heliocentric
  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 (dist Sun)
  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); // distanza Terra
  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

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

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

  ts:= calcTS(lista);  // 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

  planetlist:= {L, l, b, lapp, α,δ, az, h, M, Ω, ω, v, x,y, z, x1, y1, z1, 
    r, R, H, zc, dz, DA, SAD, ec, in, a, θ, l_ec, b_ec, Ls, Ms, xs, ys, nameplanet};
  RETURN planetlist;
END;

planetDisplay(nameplanet, lista)
BEGIN
  LOCAL ω,Ω, l, v, r;
  LOCAL x, y, z, x1, y1, z1;
  LOCAL Ls, Ms, b, xs, ys, l_ec, b_ec;
  LOCAL ec, in, a, δ, α, ts, az, h;
  LOCAL DA, SAD, tram, sorg, culm, dayY;
  LOCAL as, ac, zc, dz, θ, mov, segno, ε;
  LOCAL transito, JD, JDE, temp, Hor, Hd;
  LOCAL listaplanet, listaplanettran, listaplanetothers, ch;
  LOCAL alfa1, alfa2, alfa3, delta1, delta2, delta3, h0;
  LOCAL ee, magni, m0, aa, bb,cc, listamagn, ill;
  LOCAL d0,dapp;
  LOCAL listaARh0, lapp;

  ts:= calcTS(lista);
  JD:= julianDay(lista);
  JDE:= julianEphemerisDay(lista);

  L:= planetlist(1); l:= planetlist(2);
  b:= planetlist(3); lapp:= planetlist(4);
  α:= planetlist(5); δ:= planetlist(6);
  az:= planetlist(7); h:= planetlist(8);
  M:= planetlist(9); Ω:= planetlist(10);
  ω:= planetlist(11); v:= planetlist(12);
  x := planetlist(13);  y:= planetlist(14);  z:= planetlist(15);
  x1:= planetlist(16); y1:= planetlist(17); z1:= planetlist(18);
  r:= planetlist(19);   R:= planetlist(20);
  Hor:= planetlist(21); Hd:= 15*planetlist(21);
  zc:= planetlist(22); dz:= planetlist(23);
  DA:= planetlist(24); SAD:= planetlist(25);
  ec:= planetlist(26); in:= planetlist(27);
  a:= planetlist(28); θ:= planetlist(29);
  l_ec:= planetlist(30); b_ec:= planetlist(31);
  Ls:= planetlist(32); Ms:= planetlist(33);  // Sun values
  xs:= planetlist(34); ys:= planetlist(35);
  // nameplanet:= planetlist(36);

  segno:= zodiaco(l);

  ee:= abs(l - l_ec); // E (phase) to calc phase and magnitude
  ill:= 0.5*(1+cos(ee));
  CASE
  IF nameplanet=="Mercurius" THEN listamagn:= {6.7, -0.21,3.8,-3.4,2}; END;
  IF nameplanet=="Venus"  THEN listamagn:= {16.7, -0.21, 3.8, -3.4, 2};END;
  IF nameplanet=="Mars"  THEN listamagn:= {9.4, -1.36, 1.5,0, 0}; END;
  IF nameplanet=="Jupiter"  THEN listamagn:= {196.9, -9, 1.48, 0,0}; END;
  IF nameplanet=="Saturn"  THEN listamagn:= {165.5, -8.7, 1.7, 0,0}; END;
  IF nameplanet== "Uranus" THEN listamagn:= {71.4, -7, 0,0,0}; END;
  IF nameplanet=="Neptunus"  THEN listamagn:= {68.3, -7, 0,0,0}; END;
  IF nameplanet=="Pluto"  THEN listamagn:= {8.3, -1.5, 0, 0, 0}; END;
  DEFAULT
  END; // case
  d0:= listamagn(1);
  m0:= listamagn(2);
  aa:= listamagn(3); bb:= listamagn(4); cc:= listamagn(5);
  dapp:= d0/R;  // apparent diameter
  magni:= m0 + 5*LOG(r*R)+aa*ee/100+bb*(ee/100)^2+cc*(ee/100)^3;

  P:= →HMS(HMS→(0°0′8.794″)/R); // Parallasse

  // transit routine
  h0:= -0°34′; // 'standard' altitude of a planet
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D-1,0,0,0}));
  alfa1:= temp(1); delta1:= temp(2);
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D,0,0,0}));
  alfa2:= temp(1); delta2:= temp(2);
  temp:= planetAlphaDelta(nameplanet, makeDateListShort({Y,m,D+1,0,0,0}));
  alfa3:= temp(1); delta3:= temp(2);

  listaARh0:= {alfa1, alfa2, alfa3, delta1, delta2, delta3, h0};
  transito:= transit(lista, listaARh0);
  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(" ");
  PRINT("L Long media " + trunc_angle(L));
  PRINT("λ Longitudine " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");
  PRINT("β Latitudine " + trunc_angle(b));
  PRINT("Δ Distanza " + ROUND(R,4) + " (UA) from Earth ");
  PRINT("α Asc. retta "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));
  PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));
  PRINT(" ");
  PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));
  PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(simp_angle(v)));
  PRINT("r Raggio vettore " + ROUND(r,4) + " (UA) from Sun "); // distanza eliocentrica
WAIT();
PRINT();
  PRINT("ω Perielio " + trunc_angle(ω));
  PRINT("Ω Nodo " + trunc_angle(Ω));
  PRINT("M Anomalia vera " + trunc_angle(M));
  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(" ");
  PRINT("Ang. orario "+ trunc_angle(Hor) + " " + 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);
WAIT();
PRINT();
  PRINT("Diff. ascensionale " + trunc_angle(DA));
  PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));
  PRINT(" ");
  PRINT("π Parallasse " + ROUND(P, 5));
  PRINT("Diametro apparente "+ ROUND(dapp, 5));
  PRINT("Magnitudo " + ROUND(magni, 5));
  PRINT("Frazione illuninata " + ROUND(ill*100, 2)+"%");
  PRINT("Fase " + →HMS(ee));
  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);
  listaplanet:= {l,b, R, α, δ, az, h};
  listaplanettran:= {culm, sorg, tram};
  listaplanetothers:= {L, lapp, l_ec,b_ec, M, ω,Ω, v, x,y,z, x1,y1,z1, r, Hor, Hd, SAD, DA, mov, dapp, magni};
  Planet_Effemerids:= listaplanet;
  Planet_Transit:= listaplanettran;
  Planet_Others:= listaplanetothers;

  FREEZE();
  WAIT();
  ch:= MSGBOX("Another calculation?", 1);
  IF (ch<>0)  THEN changeData(); END;

END;

planetAlphaDelta(nameplanet, lista)
// lista like dateshortlist
BEGIN
LOCAL alpha,delta,temp, adelta;
temp:= planetCalc(nameplanet, lista);
alpha:=temp(5);
delta:=temp(6);
adelta:={alpha,delta};
RETURN adelta;
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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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
  N:= 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;



RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-02-2015 11:10 PM

New version:
corrected some bugs, added Δα and Δδ correction for daytime parallax, exported a function to calc precession (given RA -DEG- and declination), "graphic" restyling of printed values...

Enjoy!

EDIT: corrected a but for age of Moon, added "epatta" for Moon (old calculation for Easter date...)


RE: Astronomy: Effemeridi (Ephemeris) - Thomas_Sch - 07-03-2015 06:27 AM

(07-01-2015 09:48 PM)salvomic Wrote:  ...in the meantime this version introduce other values for planets: parallax, illuminated fraction, magnitudo, phase.

EDIT: Thomas, try this with an option to save coordinates at the starting, the program should save them in Long_Lat var in the program variables and load them when used again...
I need a better safeguard, however, then only control if the size is ...2 ;-)
Hello Salvo,
many thanks for your message & your efforts! I'll try it as soon as possible.
I've seen your other questions about checking the list, and the discussion about type(), very interesting! Maybe the special characters could be checked? Don't know... Range of values? (-90 >= latitude <= +90; -180 >= longitude <= +180)?


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-03-2015 06:53 AM

(07-03-2015 06:27 AM)Thomas_Sch Wrote:  Hello Salvo,
many thanks for your message & your efforts! I'll try it as soon as possible.
I've seen your other questions about checking the list, and the discussion about type(), very interesting! Maybe the special characters could be checked? Don't know... Range of values? (-90 >= latitude <= +90; -180 >= longitude <= +180)?

I'm investigating, Thomas,
for now is enough to check for size(): if there are two (or more elements is a list...); why the user should change that variable? ;-)

Do you think is important save only coordinates or also time zone (and / or DST)?
Time and date I think no, better to leave it as is: the user can change them easily...

Salvo


RE: Astronomy: Effemeridi (Ephemeris) - Thomas_Sch - 07-03-2015 07:05 AM

(07-03-2015 06:53 AM)salvomic Wrote:  I'm investigating, Thomas,
for now is enough to check for size(): if there are two (or more elements is a list...); why the user should change that variable? ;-)

Do you think is important save only coordinates or also time zone (and / or DST)?
Time and date I think no, better to leave it as is: the user can change them easily...

Salvo
You are right, time zone should also be saved, (also there is a rough correlation between longitude and time zone).
Maybe an option or checkbox ("Save?" yes / no ) could avoid uncertainties ..


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-03-2015 02:31 PM

(07-03-2015 07:05 AM)Thomas_Sch Wrote:  You are right, time zone should also be saved, (also there is a rough correlation between longitude and time zone).
Maybe an option or checkbox ("Save?" yes / no ) could avoid uncertainties ..

ok, now the program saves long, lat, time zone, DST and deltaT (like the app by Marcel, in the same way, so we can use both program and app together).

Exported also some useful function (Ascendant, making lists, Julian Day, epsilon, precession, epatta,Sun and Moon alpha/delta, Sun and Moon Ah...)

Salvo

File's date: July 5, 2015


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-06-2015 03:01 PM

New version:
rewritten some parts: now calculations for Sun, Moon, Planets are a bit slow, but more precise (rewritten the routine for nutation, deltapsi, deltaepsilon, more precise but require a loop of 63 linear combination of parameters, used also in others part of program, for precision)...
More precise transit routine for Sun and Planets and finally also the Moon (thanks Marcel!).
Reviewed Ascendant calculation (more precise) with indication of Descendant...
More precise latitude for Sun.
Corrected a but in Moon phase.

Salvo

Date last version: 2015, July, 8


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-09-2015 12:30 PM

New version:
further improved Sun precision (about 100 parameters, precision about 0.01").
Exported functions Sun_Coordinates, Moon_Coordinates, Planets_Coordinates.
They gives however only few important coordinates (lambda, beta, alfa, delta, R, parallax); to get more use Effemeridi that gives (for Sun, Moon, Planets) other data: anomaly, node, perigee/apogee, sidereal time, azimuth, height, cartesian coordinates, age of Moon..., and transit (exported also as independent function).

Date: 2015, July, 9


RE: Astronomy: Effemeridi (Ephemeris) - Marcel - 07-09-2015 01:50 PM

Hi Salvomic!

All that DATA (anomaly, node,...) is a great addition to my program.

Thank you.
Marcel


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-09-2015 04:17 PM

(07-09-2015 01:50 PM)Marcel Wrote:  Hi Salvomic!

All that DATA (anomaly, node,...) is a great addition to my program.

Thank you.
Marcel

thank you Marcel,
I think that the best experience for our users is to use both app and program: they give the same data for the principal aspects (Sun, Moon, Planets, time, most important parameters), then different approach (graphic and detailed). They are becoming "complementary"...
And we must always thank J. Meeus for the great input and his precise calculations!
Eventually, the Prime is a great machine to do Astronomy tools, and our programs demonstrate that.

Salvo


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-10-2015 08:15 PM

hi,
in this new version I added a function dedicated to the Earth (ellipsoid).
Earth(datelist);
With it the user can input a latitude (or use that already saved with the program) -and (or not) heigh above the sea level- and calculate some geodetic parameters: geocentric latitude, flattening, equatorial and meridian radius, eccentricity of the meridian, rho•sin(phi) and rho•cos(phi), difference between geographic and geocentric latitude, radius of parallel of latitude, rotational angular velocity, radius of curvature of the meridian...
Giving lon/lat for two points the program returns the angular and linear distance (with an error or ±50 meters!).
Then the program output actual heliocentric longitude and latitude for the Earth.

Examples:
• to input φ, h (height above sea level): ρ, ρ·sin(φ'), ρ·sin(φ') :: useful data to get diurnal parallax, eclipses, occultations)
• to input φ: Rp, Rm, λ°, φ° (radius of parallel of latitude, radius of curvature, degree of longitude or latitude length
• to input φ1, φ2: s, d (distance from two points given their coordinates: both linear and angular)

Enjoy!

Salvo M.

Date: 2015, July, 13
Version: 2.1


RE: Astronomy: Effemeridi (Ephemeris) - Marcel - 07-12-2015 12:38 PM

Hi Salvomic!
Great to see that this chapter on the Earth globe is in your program!
The book is covered from the start to the end..

If the probe New Horizons interest you, I sugest to you the download of the program "NASA's Eyes". This is a beta application from the JPL Lab. If you go with the probe New Horizons, you will have a "Live" option and a "Preview" option. This "Preview" option will show you what the probe will do tomorrow and the day after. Is you look at the right of the screen, you will see what is camera will do... This is really a great animation... I think Astro Lab 8 or may be 9 will show this view!

Marcel


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-12-2015 08:57 PM

(07-12-2015 12:38 PM)Marcel Wrote:  Hi Salvomic!
Great to see that this chapter on the Earth globe is in your program!
The book is covered from the start to the end..
yes, Marcel, the book is covered all! and it deserve: it's a great book!
Quote:If the probe New Horizons interest you, I sugest to you the download of the program "NASA's Eyes". This is a beta application from the JPL Lab. If you go with the probe New Horizons, you will have a "Live" option and a "Preview" option. This "Preview" option will show you what the probe will do tomorrow and the day after. Is you look at the right of the screen, you will see what is camera will do... This is really a great animation...
yes, I'm interested enough, that's a historic event with Pluto, we'll know many new things about this fascinating planet.
I'll download the program tomorrow.
Quote:I think Astro Lab 8 or may be 9 will show this view!

Marcel

sure! :-)
Thanks always for your great app!
Salvo


RE: Astronomy: Effemeridi (Ephemeris) - salvomic - 07-14-2015 10:17 AM

New version:
implemented a new method to calculate planets coordinates (for more precise result -max ±30', mostly ±2'- use the correspondent functions in Astro Lab 4 that uses another precise method to calculates values);
exported a function to calculate angular separation for two bodies, giving right ascension and declination for both (d = angular separation from the Earth, p = relative separation);
exported also the function to interpolate from three values, giving the "n" factor, useful to get also angular separation (and many other purpose) from a table of RA and dec;
exported Anomaly (mean, eccentric, true) and Moon elongation;
translated other parts;
corrected printing of perihelion and added "argument of perihelion";
corrected (and exported) the tool to get "eccentric anomaly" via Kepler equation...

I must also write a guide: soon.

Enjoy!
Salvo

Version: 2.3
Date: 2015, July, 18