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

Some routines and data are still in Italian, sorry.

Salvo

The code (I'll post here the new versions for now):
Code:
 data(); calcN(); calcTS(); sunCalc(); planetCalc(); planets(); sun(); moon(); Mercurius(); Venus(); Mars(); Jupiter(); Saturn(); Uranus(); Neptunus(); Pluto(); ascendent(); zodiaco(); trunc_angle(); atan2(); thisday:= 2000.0101; lstd:=0; gmt:= 12; long:= 0; lat:=  0; m:= 1; dst:= 0; utc:=0; sunlist:= MAKELIST(X,X,1,25); planet_const:= MAKELIST(X,X,1,10); ε:= →HMS(23.4392911);  // approx 23.45 EXPORT Effemeridi() // from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015 // Credits: from a book, Calcule astronomique pour amateurs, by Serge Bouiges // Thanks Didier Lachieze, Marcel, Eried BEGIN   LOCAL ch;   HAngle:=1;   data(); // input data and calc Sun   CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc N", "Calc TS", "Ascendent", "Quit");   CASE   IF ch==1 THEN sun(); END;   IF ch==2 THEN moon(); END;   IF ch==3 THEN planets(); END;   IF ch==4 THEN calcN(); END;   IF ch==5 THEN calcTS(); END;   IF ch==6 THEN ascendent(); END;   IF ch==7 THEN HAngle:=0; RETURN; END;   DEFAULT   END; // case END; data() BEGIN   LOCAL dd, mm, yy, hh, mn, ss, loct;   yy:= IP(Date);   mm:= IP(FP(Date)*100);   dd:= FP(Date*100)*100;   hh:= IP(HMS→(Time));   mn:= IP(FP(HMS→(Time))*60);   ss:= IP(FP(FP(FP(HMS→(Time))*60))*60);   loct:= →HMS(hh+mn/60+ss/3600);   INPUT({ {m,[0],{15,15,1}}, {D,[0],{50,15,1}}, {Y,[0],{80,15,1}},   {lstd,[0],{65,30,2}},   {dst,0,{40,15,3}}, {utc,[0],{80,15,3}},    {long,[0],{20,25,4}},    {lat,[0],{70,25,4}} },   "Data: Use Shift+9 for H°M′S″",   {"Month:","Day :","Year :","Local Time (24):", "DST", "GMT",    "Long (-E):","Lat (+N)"},   {"Enter month", "Enter day", "Enter year (1901-2099)",    "Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC lag",   "Longitude", "Latitude"},    {1,1,1901,0, 0,0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″});   // Ragusa (-14°43′41″, 36°55′15″) -    // Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)   D:=EVAL(D); m:=EVAL(m); Y:=EVAL(Y);   lstd:= lstd;    // Greenwich Mean Time   gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)   thisday:= Y+m/100+D/10000;   sunCalc(); // calculate Sun data END; planets() BEGIN   LOCAL ch, nameplanet;   INPUT({{ch,{"Mercurius","Venus","Mars", "Jupiter", "Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},   "Choose planet",    "planet: ", "Choose the planet to calculate an press OK", 0,5 );   CASE   IF ch==1 THEN nameplanet:="Mercurius"; planetCalc(nameplanet); END;   IF ch==2 THEN nameplanet:="Venus"; planetCalc(nameplanet); END;   IF ch==3 THEN nameplanet:="Mars"; planetCalc(nameplanet); END;   IF ch==4 THEN nameplanet:="Jupiter"; planetCalc(nameplanet); END;   IF ch==5 THEN nameplanet:="Saturn"; planetCalc(nameplanet); END;   IF ch==6 THEN nameplanet:="Uranus"; planetCalc(nameplanet); END;   IF ch==7 THEN nameplanet:="Neptunus"; planetCalc(nameplanet); END;   IF ch==8 THEN nameplanet:="Pluto"; planetCalc(nameplanet); END;   DEFAULT   END; // case END; calcN() BEGIN   N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT   RETURN N; END; calcTS() BEGIN   LOCAL TS, TSd, N0, TSL;   calcN();   N0:= DDAYS(1901.0101, thisday)+1+0; // N at 0 GMT (for TS)   TS:= ((23750.3 + 236.555362*(N0+gmt/24)) / 3600) MOD 24; // Sideral Time (GMT) in seconds   TS:= →HMS(TS); // in HMS   TSL:= (((23750.3 + 236.555362*N+lstd/24) / 3600) MOD 24); // local ST   TSL:= (TSL + gmt - 4*(HMS→(long))/60) MOD 24;   TSL:= →HMS(TSL); // in HMS   TSd:= →HMS((98.965 + 0.985647342*N0) MOD 360); // ST at 0 GMT in degrees   RETURN {TS, TSL, TSd, gmt}; END; ascendent() BEGIN   LOCAL ASC, ts, segno;   ts:= calcTS();   ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));   IF ASC<0 THEN ASC:= ASC+360; END; IF ASC>360 THEN ASC:=ASC-360; END;   segno:= zodiaco(ASC); RETURN {→HMS(ASC), segno}; END; zodiaco(astro) BEGIN   LOCAL segno;   CASE  // Segni e case   IF astro>=0   AND astro <30  THEN segno:="Ariete"; END;   IF astro>=30  AND astro <60  THEN segno:="Toro"; END;   IF astro>=60  AND astro <90  THEN segno:="Gemelli"; END;   IF astro>=90  AND astro <120 THEN segno:="Cancro"; END;   IF astro>=120 AND astro <150 THEN segno:="Leone"; END;   IF astro>=150 AND astro <180 THEN segno:="Vergine"; END;   IF astro>=180 AND astro <210 THEN segno:="Bilancia"; END;   IF astro>=210 AND astro <240 THEN segno:="Scorpione"; END;   IF astro>=240 AND astro <270 THEN segno:="Sagittario"; END;   IF astro>=270 AND astro <300 THEN segno:="Capricorno"; END;   IF astro>=300 AND astro <330 THEN segno:="Acquario"; END;   IF astro>=330 AND astro <360 THEN segno:="Pesci"; END;   END; // case   RETURN segno; END; atan2(y,x) BEGIN // atan2 returns coordinates in correct angle for all four quadrants (DEG)     CASE     IF x==0 AND y==0 then return "undefined"; end;         // x=0, y=0       IF x>0 then return atan(y/x); end;                     // x>0     IF x<0 AND y>=0 then return atan(y/x)+180; end;         // x<0, y>=0     IF x==0 AND y>0 then return 90; end;                 // x=0, y>0     IF x<0 AND y<=0 then return atan(y/x)-180; end;         // x<0, y<0     IF x==0 AND y<0 then return -90; end;                // x=0, y<0    END;   RETURN;  END; trunc_angle(ang) // trunc decimal from a DEG form angle BEGIN →HMS((IP(HMS→(ang)*3600)/3600)); END; sunCalc() BEGIN   LOCAL ω, l, v, r, x, y, dz;   LOCAL α, δ, dayY, H, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc;   calcN();   L:= →HMS((278.965+0.985647342*N)) MOD 360; // long media   ω:= →HMS((281.235 + 0.0000469*N)) MOD 360; // long del perielio   M:= (L - ω); // anomalia vera   v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio   l:= v + ω; // Longitudine Sole   lt:= l + 180; // Longitudine eliocentrica della Terra   r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)   x:= r * cos(l); y:= r * sin(l);   dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper   δ:= →HMS((asin(sin(ε)*sin(l))) );   // α:= →HMS((acos(cos(l)/cos(δ))/15) MOD 24); // ascensione retta   α:= (acos(cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta   IF l>180 THEN α:= α +180; END;   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale   // alt:= →HMS(90 - z); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza   // as:= cos(δ)*sin(Hd)/sin(dz); // calc az   // ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   // IF (ac<0) THEN az:= (360 + az) MOD 360; END;   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   IF (az<0) THEN az:= 360+az; END;   E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   // SAD:= →HMS(acos(-tan(lat)*tan(δ))); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)   sunlist:= {L,ω,l, M, α,δ, az, h, x,y,sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA, SAD, E, lt};   RETURN sunlist; END; sun() BEGIN   LOCAL ω, l, v, r, x, y, dz;   LOCAL α, δ, dayY, H, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, segno;   L:= sunlist(1); ω:= sunlist(2);    l:= sunlist(3); M:= sunlist(4);    α:= sunlist(5); δ:= sunlist(6);    az:= sunlist(7); h:= sunlist(8);    x:= sunlist(9); y:= sunlist(10);   sorg:= sunlist(11); tram:= sunlist(12); culm:= sunlist(13);   r:= sunlist(14); H:= sunlist(15); Hd:= sunlist(16);   zc:= sunlist(17); dz:= sunlist(18);    ts:= calcTS(); // no 19   DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("");   PRINT("L longitudine media " + trunc_angle(L));   PRINT("ω long del perielio " + trunc_angle(ω));   PRINT("M anomalia vera " + trunc_angle(M));   PRINT("λ longitudine Sole " + trunc_angle(l) + " " + segno);   PRINT("r Distanza (UA) " + r);   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT(" ");   PRINT("Ang. orario "+ H + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS locale " + trunc_angle(ts(2)));   PRINT("Equazione del tempo (min) " + trunc_angle(E));   PRINT("Dif asc." + DA + " SAD " + trunc_angle(SAD));   PRINT("");   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT("Long elioc. Terra "+ trunc_angle(lt));   // RETURN {L, ω, M, v, l, r, x, y};   RETURN [[L, ω], [l, M], [α,δ], [az, h], [x,y], [sorg, tram], [culm, r], [H,ts(1)]]; END; moon() BEGIN   LOCAL ω,Ω, l, v, r, x, y, i;   LOCAL Ls, Ms, LongM, LatM, P, d, s;   LOCAL α,δ, ts, Hd, zc, dz, h, az;   LOCAL as, ac, DA, SAD, DD, culm, sorg, tram;   LOCAL b, dayY, xs, ys, ph, fase, eta;   LOCAL dd, hh, mm, segno;   i:= 5.13333333; // inclinazione 5°8' sull'eclittica   calcN();   // Sun   Ls:= sunlist(1); // Longitudine media   Ms:= sunlist(4); // Anomalia   // end Sun   L:= (33.231 + 13.17639653*N) MOD 360; // long. media   Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente   M:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera   ω:= (L - M) MOD 360; // perigeo medio   DD:= L - Ls; F:= L-Ω; // for the variations   LongM:= L + 6.28875 * sin(M); // calculus (h/ and above) of Moon longitude   LongM:= LongM + 0.2136 * sin(2*M); // equazione del centro   LongM:= LongM + 0.6583 * sin(2*DD); // variazione   LongM:= LongM - 0.1856 * sin(Ms); // equazione annuale   LongM:= LongM + 1.2740 * sin(2*DD-M); // evezione   LongM:= LongM - 0.1143 * sin(2*F); // riduzione all'eclittica   LongM:= LongM + 0.0588 * sin(2*DD-2*M);   LongM:= LongM + 0.0572 * sin(2*DD-M-Ms);   LongM:= LongM + 0.0533 * sin(2*DD+M);   LongM:= LongM + 0.0459 * sin(2*DD-Ms);   LongM:= LongM + 0.0410 * sin(M-Ms);   LongM:= LongM - 0.0305 * sin(M+Ms);   LongM:= LongM - 0.0348 * sin(DD);   LongM:= LongM MOD 360;   LatM := 5.1280 * sin(F); // termine principale   LatM := LatM + 0.2806 * sin(M+F);   LatM := LatM + 0.2777 * sin(M-F); // equazione del centro   LatM := LatM + 0.1732 * sin(2*DD-F); // grande ineguaglianza   LatM := LatM + 0.0554 * sin(2*DD-M+F); // evezione   LatM := LatM + 0.0462 * sin(2*DD-M-F); // evezione (2)   LatM := LatM + 0.0326 * sin(2*DD+F);   LatM := LatM MOD 360; b:= LatM; dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(LongM))) ); // declinazione (coord locali)   // α:= →HMS((acos(cos(b)*cos(LongM)/cos(δ))/15) MOD 24); // ascensione retta   α:= (acos(cos(b)*cos(IFTE(LongM>180, LongM-180, LongM))/cos(δ))); // ascensione retta   IF LongM>180 THEN α:= α +180; END;   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)   as:= cos(δ)*sin(Hd)/sin(dz); // calc az   ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   // IF (ac<0) THEN az:= ((360 + az) MOD 360); END;   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse   // Raggio quadratico medio della Terra: 6373,044737   d:= 1/sin(P); // distanza Terra-Luna in raggi terrestri   s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)   ph:= abs(LongM - Ls); // Phase   CASE   // IF ph == 0 THEN fase:="New Moon"; END;   IF ph>280 AND ph<350 THEN fase:="Decrescente"; END;   IF ph>=260 AND ph<=280 THEN fase:="Last Quart"; END;   IF ph>190 AND ph<260 THEN fase:="Gibbosa Cal."; END;   IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;   IF ph>100 AND ph<170 THEN fase:="Gibbosa Cresc."; END;   IF ph>=80 AND ph<=100 THEN fase:="1st Quart"; END;   IF ph>10 AND ph<80 THEN fase:="Crescente"; END;   DEFAULT fase:="New Moon";   END; // case   // età della Luna   dd:= IP(HMS→(ph/15));   hh:= IP(FP(HMS→(ph/15))*60);   IF (hh>24) THEN dd:= dd+1; hh:= hh MOD 24; END;   mm:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);   eta:= →HMS(dd+hh/60+mm/3600); // Age of the Moon   segno:= zodiaco(LongM);   PRINT;   PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L));   PRINT("M anomalia vera " + trunc_angle(M));   PRINT("λ Longitudine " + trunc_angle(LongM) + " " + segno);   PRINT("β Latitudine " + trunc_angle(LatM));   PRINT("ω Perigeo medio " + trunc_angle(ω));   PRINT("Ω Nodo ascendente " + trunc_angle(Ω));   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("α Ascensione retta (gradi) " + trunc_angle((α*15) MOD 360));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS (locale) " + trunc_angle(ts(2)));   PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));   PRINT("P Parallasse " + trunc_angle(P));   PRINT("d Distanza Terra - Luna "+ ROUND(d*6373.044737,3) + " km");   PRINT("Semidiametro "+  trunc_angle(→HMS(s)) + " :: " + ROUND(s*6373.044737,3) + " km");   PRINT(" ");   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");   PRINT("Età " + eta);   RETURN [[L, M],[LongM, LatM],[ω,Ω],[α,δ], [az, h], [P,0],[d,s], [sorg, tram], [culm, ts(1)],[ph, eta] ]; END; planetCalc(nameplanet) BEGIN   LOCAL ω,Ω, l, v, r, x, y, z;   LOCAL Ls, Ms, b, xs, ys, b_ec;   LOCAL x1, y1, z1, R, l_ec, H, Hd;   LOCAL ec, in, a, δ, α, ts, az, h;   LOCAL DA, SAD, tram, sorg, culm, dayY;   LOCAL as, ac, zc, dz, θ, mov, segno;   // calc of Sun (as basis)   xs:= sunlist(9); ys:= sunlist(10);   Ls:= sunlist(1); // Longitudine media   Ms:= sunlist(4); // Anomalia   // end Sun  EVAL(EXPR(nameplanet+"()")); // call planet-name function   L:= planet_const(1); // Longitudine media   ω:= planet_const(2); // perielio   Ω:= planet_const(3); // nodo   M:= planet_const(4); // anomalia   v:= planet_const(5);    ec:= planet_const(6); // eccentricità   in:= planet_const(7); // inclinazione   a:= planet_const(8);  // semiasse maggiore   θ:= planet_const(9);  // θs to calc if retrogade   b_ec:= asin(sin(in)*sin(v+ω-Ω)); // argomento di latitudine del perielio   l_ec:= acos(cos(v+ω-Ω)/cos(b_ec));  // long from eclittic   IF b_ec<0 THEN l_ec:= 360 - l_ec; END; // argomento secondario   l_ec:= →HMS((l_ec + Ω)) MOD 360; // aggiunge Nodo e riporta a 2PI   r:= a*(1-ec^2)/(1+ec*cos(v)); // 9.524899/(1+ 0.0559*cos(v)); raggio vettore   x1:= r*cos(b_ec)*cos(l_ec); y1:= r*cos(b_ec)*sin(l_ec); z1:=r*sin(b_ec); // cartesian   x:= x1+xs; y:= y1+ys; z:=z1; // cohordinates cartesian->spheric   R:= sqrt(x^2+y^2+z^2);    b:=asin(z/R); // latitudine  β   // l:=→HMS(atan(y/x)); // cart->spheric longtitudine λ   l:= atan2(y,x);   IF l<0 THEN l:= l+360; END; IF l>360 THEN l:=l-360; END;   dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper   δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(l))) ); // declinazione (coord locali)   α:= (acos(cos(b)*cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta IF l>180 THEN α:= α +180; END;   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)   as:= cos(δ)*sin(Hd)/sin(dz); // calc az   ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   DA:= →HMS(asin(tan(lat)*tan(δ))); // Differenza ascensionale   SAD:= →HMS(90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) MOD 24 ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   IF abs(l-220.4)>θ THEN mov:="diretto"; ELSE mov:="retrogrado"; END;   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi di " + nameplanet + " a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L));   PRINT("ω Perielio " + trunc_angle(ω));   PRINT("Ω Nodo " + trunc_angle(Ω));   PRINT("M Anomalia vera " + trunc_angle(M));   PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));   PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(v));   PRINT("r Raggio vettore " + r);   PRINT("Cartes. elioc.: x' " + ROUND(x1,3) + " y' " + ROUND(y1,3) + " z' " + ROUND(z1,3));   PRINT("Cartes. geoc.: x " + ROUND(x,3) + " y " + ROUND(y,3) + " z " + ROUND(z,3));   PRINT("Spher.: λ Longitudine " + trunc_angle(l) + " " + segno); // spheric   PRINT("Spher.: β Latitudine " + trunc_angle(b));   PRINT("Spher.: R Distanza " + ROUND(R,3));   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS (locale) " + trunc_angle(ts(2)));   PRINT("Diff. ascensionale " + trunc_angle(DA));   PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT("Movimento " + mov);   RETURN [[L,ω,Ω ],[M,l_ec,b_ec], [v,r,0], [x1,y1,z1],[x,y,z],[l,b,R], [α,δ,H], [az,h,ts(1)],[sorg, tram, culm]]; END; Mercurius() BEGIN   // Costanti per Mercurio   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.205615; // eccentricità   in:= 7.003; // inclinazione   a:= 0.387098; // semiasse maggiore   θ:= 35.6;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   ω:= →HMS((75.913 + 0.00004253*N) MOD 360); // perielio   Ω:= →HMS((47.157 + 0.00003244*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 23.4372*sin(M)+2.9810*sin(2*M)+0.5396*sin(3*M)+0.1098*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Venus() BEGIN   // Costanti per Venere   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.006816; // eccentricità   in:= 3.393636; // inclinazione   a:= 0.72333; // semiasse maggiore   θ:= 13.0;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((206.758 + 1.60216873*N) MOD 360); // Longitudine media   ω:= →HMS((130.154 + 0.00003757*N) MOD 360); // perielio   Ω:= →HMS((75.797 + 0.00002503*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 0.7811*sin(M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Mars() BEGIN   // Costanti per Marte   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.093309; // eccentricità   in:= 1.850303; // inclinazione   a:= 1.523678; // semiasse maggiore   θ:= 16.8;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((124.765 + 0.52407108*N) MOD 360); // Longitudine media   ω:= →HMS((334.251 + 0.00005038*N) MOD 360); // perielio   Ω:= →HMS((48.794 + 0.00002127*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 10.608*sin(M)+0.6216*sin(2*M)+0.0504*sin(3*M)+0.0047*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Jupiter() BEGIN   // Costanti per Giove   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.048335; // eccentricità   in:= 1.308736; // inclinazione   a:= 5.202561; // semiasse maggiore   θ:= 54.43;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((268.350 + 0.08312942*N) MOD 360); // Longitudine media   ω:= →HMS((12.737 + 0.00004408*N) MOD 360); // perielio   Ω:= →HMS((99.453 + 0.00002767*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5372*sin(M)+0.1662*sin(2*M)+0.007*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Saturn() BEGIN   // Costanti per Saturno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.055892; // eccentricità   in:= 2.492519; // inclinazione   a:= 9.554747; // semiasse maggiore   θ:= 65.53;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((278.774 + 0.03349788*N) MOD 360); // Longitudine media   ω:= →HMS((91.117 + 0.00005362*N) MOD 360); // perielio   Ω:= →HMS((112.79 + 0.00002391*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 6.4023*sin(M)+0.2235*sin(2*M)+0.011*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Uranus() BEGIN   // Costanti per Urano   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.046344; // eccentricità   in:= 0.772464; // inclinazione   a:= 19.21814; // semiasse maggiore   θ:= 73.92;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((248.487 + 0.01176902*N) MOD 360); // Longitudine media   ω:= →HMS((171.563 + 0.00004064*N) MOD 360); // perielio   Ω:= →HMS((73.482 + 0.00001365*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5.3092*sin(M)+0.1537*sin(2*M)+0.0062*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Neptunus() BEGIN   // Costanti per Nettuno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.008997; // eccentricità   in:= 1.779242; // inclinazione   a:= 30.10957; // semiasse maggiore   θ:= 77.63;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((86.652 + 0.00602015*N) MOD 360); // Longitudine media   ω:= →HMS((46.742 + 0.000039*N) MOD 360); // perielio   Ω:= →HMS((130.693 + 0.00003009*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 1.031*sin(M)+0.0058*sin(2*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Pluto() BEGIN   // Costanti per Plutone   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.250236; // eccentricità   in:= 17.17047; // inclinazione   a:= 39.438712; // semiasse maggiore   θ:= 79.41;  // θs to calc if retrogade   calcN();   // costanti planetarie   // Long, perielio, nodo, non precisi, considerati costanti   // L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   M:= →HMS((230.67 + 0.00397943*N) MOD 360); // Anomalia vera (calcolata direttamente)   Ω:= →HMS((109.06 + 0.00003823*N) MOD 360); // nodo   ω:= →HMS((114.27 + Ω) MOD 360); // perielio   L:= →HMS((M + ω) MOD 360); // L (non sicuro)      v:= M + 28.45*sin(M)+4.38*sin(2*M)+0.97*sin(3*M)+0.24*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-20-2015, 07:18 PM
Post: #2
 Marcel Member Posts: 177 Joined: Mar 2014
RE: Astronomy...
Hi Salvomic
We don't use the same technique for calculating the planetary position. I use the VSOP technique (Variation séculaire des orbites planétaires). In this technique, we use polynomials and series of DATA to obtain the heliocentric position of a planet...
I complete today two new functions for the location of Pluto. This is the year of Pluto because in 25 days, the probe New Horizon will arrive at Pluto. For this reason, I will post my astronomy app's just before that moment.
Good lock in your project, this is fun and we can use both programs.

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

Marcel

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

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

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

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-20-2015, 10:29 PM (This post was last modified: 06-20-2015 10:33 PM by Marcel.)
Post: #4
 Marcel Member Posts: 177 Joined: Mar 2014
RE: Astronomy...
Hi,
I send you my app's with Pluto, so you will be able to compare with your program. If you have Mathematica or WolframAlpha, check for precision. It is very good.
Marcel

Attached File(s)
06-21-2015, 06:38 AM (This post was last modified: 06-21-2015 06:31 PM by salvomic.)
Post: #5
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Astronomy...
(06-20-2015 10:29 PM)Marcel Wrote:  Hi,
I send you my app's with Pluto, so you will be able to compare with your program. If you have Mathematica or WolframAlpha, check for precision. It is very good.
Marcel

thank Marcel,
I'm trying it!

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

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

Salvo

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

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-22-2015, 12:42 PM
Post: #6
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Astronomy...
some improvements (calculation of times, sidereal times, julian day, epsilon, data for Pluto...)

Code:
 data(); calcN(); calcTS(); julianDay(); tempo(); sunCalc(); planetCalc(); planets(); sun(); moon(); Mercurius(); Venus(); Mars(); Jupiter(); Saturn(); Uranus(); Neptunus(); Pluto(); ascendent(); zodiaco(); trunc_angle(); atan2(); thisday:= 2000.0101; lstd:=0; gmt:= 12; long:= 0; lat:=  0; m:= 1; dst:= 0; utc:=0; deltaT:= 68; datetimelist:= MAKELIST(X,X,1,6); sunlist:= MAKELIST(X,X,1,25); planet_const:= MAKELIST(X,X,1,10); ε:= →HMS(23.4392911);  // approx 23.45 (23°26'21") EXPORT Effemeridi() // from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015 // Credits: from a book, Calcule astronomique pour amateurs, by Serge Bouiges // Thanks Didier Lachieze, Marcel, Eried BEGIN   LOCAL ch;   HAngle:=1;   data(); // input data and calc Sun   CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "Ascendent", "Quit");   CASE   IF ch==1 THEN sun(); END;   IF ch==2 THEN moon(); END;   IF ch==3 THEN planets(); END;   IF ch==4 THEN tempo(); END;   IF ch==5 THEN ascendent(); END;   IF ch==6 THEN HAngle:=0; RETURN; END;   DEFAULT   END; // case END; data() BEGIN   LOCAL dd, mm, yy, hh, mn, ss, loct;   LOCAL t, jd;   yy:= IP(Date);   mm:= IP(FP(Date)*100);   dd:= FP(Date*100)*100;   hh:= IP(HMS→(Time));   mn:= IP(FP(HMS→(Time))*60);   ss:= IP(FP(FP(FP(HMS→(Time))*60))*60);   loct:= →HMS(hh+mn/60+ss/3600);   INPUT({ {m,[0],{15,15,1}}, {D,[0],{50,15,1}}, {Y,[0],{80,15,1}},   {lstd,[0],{65,30,2}},   {dst,0,{40,15,3}}, {utc,[0],{80,15,3}},    {long,[0],{20,25,4}},    {lat,[0],{70,25,4}},   {deltaT, [0], {20,25,5}} },   "Data: Use Shift+9 for H°M′S″",   {"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",    "Long (-E):","Lat (+N)", "delta T"},   {"Enter month", "Enter day", "Enter year (1901-2099)",    "Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",   "Longitude", "Latitude", "Delta T (sec.)"},    {1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});   // Ragusa (-14°43′41″, 36°55′15″) -    // Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)  //  D:=EVAL(D); m:=EVAL(m); Y:=EVAL(Y);   datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };   // Greenwich Mean Time   gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)   thisday:= Y+m/100+D/10000; jd:= julianDay(); t:=(jd(1)-2451545)/36525; ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; ε:= →HMS(ε);   sunCalc(); // calculate Sun data END; planets() BEGIN   LOCAL ch, nameplanet;   INPUT({{ch,{"Mercurius","Venus","Mars", "Jupiter", "Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},   "Choose planet",    "planet: ", "Choose the planet to calculate an press OK", 0,5 );   CASE   IF ch==1 THEN nameplanet:="Mercurius"; planetCalc(nameplanet); END;   IF ch==2 THEN nameplanet:="Venus"; planetCalc(nameplanet); END;   IF ch==3 THEN nameplanet:="Mars"; planetCalc(nameplanet); END;   IF ch==4 THEN nameplanet:="Jupiter"; planetCalc(nameplanet); END;   IF ch==5 THEN nameplanet:="Saturn"; planetCalc(nameplanet); END;   IF ch==6 THEN nameplanet:="Uranus"; planetCalc(nameplanet); END;   IF ch==7 THEN nameplanet:="Neptunus"; planetCalc(nameplanet); END;   IF ch==8 THEN nameplanet:="Pluto"; planetCalc(nameplanet); END;   DEFAULT   END; // case END; calcN() BEGIN   N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT   RETURN N; END; calcTS() BEGIN   LOCAL TS, TSd, N0, TSL, TSapp;   calcN();   N0:= DDAYS(1901.0101, thisday)+1+0; // N at 0 GMT (for TS)   TS:=  ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds   TS:= →HMS(TS); // in HMS    TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST   TSL:= →HMS(TSL); // in HMS   TSapp:= (TS + gmt) MOD 24; // apparent ST (ST at Greenwitch at local our)   TSapp:= →HMS(TSapp); // in HMS   TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees   RETURN {trunc_angle(TS), trunc_angle(TSL), trunc_angle(TSapp), trunc_angle(TSd), gmt}; END; julianDay() BEGIN   LOCAL yy, mm, dd,dde, hh, mn, ss;   LOCAL aa,bb, jd, jde;   yy:= Y; mm:= m;   hh:= datetimelist(4);   mn:= datetimelist(5);   ss:= datetimelist(6);   dd:= D + hh/24 + mn/(60*24) + ss/(60*60*24);   dde:= D + hh/24 + mn/(60*24) + (ss+deltaT)/(60*60*24);   IF (m=1 OR m=2) THEN yy:=Y-1; mm:=m+12;END;   IF yy*100+mm+dd>=158225 THEN     aa:=IP(ABS(yy/100));     bb:=2-aa+IP(aa/4);   ELSE bb:=0; END;   jd:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dd+bb-1524.5;   jde:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dde+bb-1524.5;   RETURN {jd, jde}; END; tempo() BEGIN   LOCAL n, jd, TSid;   n:= calcN();   jd:= julianDay();   TSid:= calcTS;   PRINT;   PRINT("Date " + m + " " + D + " " + Y);   PRINT("Julian Day " + jd(1));   PRINT("Julian Day Effem. " + jd(2));   PRINT("N (from 1901jan1) " + n);   PRINT("Tempo Siderale 0h " + trunc_angle(TSid(1)));   PRINT("Tempo Siderale locale " + trunc_angle(TSid(2)));   PRINT("Tempo Siderale apparente " + trunc_angle(TSid(3)));   PRINT("Tempo Siderale gradi " + trunc_angle(TSid(4)));   PRINT("Time " + lstd + " GMT " + gmt);   RETURN({ROUND(jd(1),5), ROUND(jd(2),5), ROUND(n, 5), trunc_angle(TSid(1)), trunc_angle(TSid(2)),      trunc_angle(TSid(3)), trunc_angle(TSid(4)), gmt }); END; ascendent() BEGIN   LOCAL ASC, ts, segno;   ts:= calcTS();   ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));   IF ASC<0 THEN ASC:= ASC+360; END; IF ASC>360 THEN ASC:=ASC-360; END;   segno:= zodiaco(ASC); RETURN {→HMS(ASC), segno}; END; zodiaco(astro) BEGIN   LOCAL segno;   CASE  // Segni e case   IF astro>=0   AND astro <30  THEN segno:="Ariete"; END;   IF astro>=30  AND astro <60  THEN segno:="Toro"; END;   IF astro>=60  AND astro <90  THEN segno:="Gemelli"; END;   IF astro>=90  AND astro <120 THEN segno:="Cancro"; END;   IF astro>=120 AND astro <150 THEN segno:="Leone"; END;   IF astro>=150 AND astro <180 THEN segno:="Vergine"; END;   IF astro>=180 AND astro <210 THEN segno:="Bilancia"; END;   IF astro>=210 AND astro <240 THEN segno:="Scorpione"; END;   IF astro>=240 AND astro <270 THEN segno:="Sagittario"; END;   IF astro>=270 AND astro <300 THEN segno:="Capricorno"; END;   IF astro>=300 AND astro <330 THEN segno:="Acquario"; END;   IF astro>=330 AND astro <360 THEN segno:="Pesci"; END;   END; // case   RETURN segno; END; atan2(y,x) BEGIN // atan2 returns coordinates in correct angle for all four quadrants (DEG)     CASE     IF x==0 AND y==0 then return "undefined"; end;         // x=0, y=0       IF x>0 then return atan(y/x); end;                     // x>0     IF x<0 AND y>=0 then return atan(y/x)+180; end;         // x<0, y>=0     IF x==0 AND y>0 then return 90; end;                 // x=0, y>0     IF x<0 AND y<=0 then return atan(y/x)-180; end;         // x<0, y<0     IF x==0 AND y<0 then return -90; end;                // x=0, y<0    END;   RETURN;  END; trunc_angle(ang) // trunc decimal from a DEG form angle BEGIN →HMS((IP(HMS→(ang)*3600)/3600)); END; sunCalc() BEGIN   LOCAL ω, l, v, r, x, y, dz;   LOCAL α, δ, dayY, H, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, h0, H0, temp;   h0:= -0°50′; // 'standard' altitude of Sun   calcN();   julianDay();   L:= →HMS((278.965+0.985647342*N)) MOD 360; // long media   ω:= →HMS((281.235 + 0.0000469*N)) MOD 360; // long del perielio   M:= (L - ω); // anomalia vera   v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio   l:= v + ω; // Longitudine Sole   lt:= l + 180; // Longitudine eliocentrica della Terra   r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)   x:= r * cos(l); y:= r * sin(l);   dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper   δ:= →HMS((asin(sin(ε)*sin(l))) );   // α:= →HMS((acos(cos(l)/cos(δ))/15) MOD 24); // ascensione retta   α:= (acos(cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta   IF l>180 THEN α:= α +180; END;   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi // temp:= sin(h0)-sin(lat)*sin(δ)/cos(lat)*cos(δ); // IF temp<0 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END; // H0:= acos(temp); // MSGBOX("H0: " + H0);   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale   // alt:= →HMS(90 - z); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza   // as:= cos(δ)*sin(Hd)/sin(dz); // calc az   // ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   // IF (ac<0) THEN az:= (360 + az) MOD 360; END;   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   IF (az<0) THEN az:= 360+az; END;   E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   // SAD:= →HMS(acos(-tan(lat)*tan(δ))); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset) //temp:= ((α*15 + long - ts(3)*15) / 360); // IF temp<0 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END; // culm:= (temp - Hd/360) *24; // sorg:= (culm - H0/360) *24; // tram:= (culm + H0/360) *24;   sunlist:= {L,ω,l, M, α,δ, az, h, x,y,sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA, SAD, E, lt};   RETURN sunlist; END; sun() BEGIN   LOCAL ω, l, v, r, x, y, dz;   LOCAL α, δ, dayY, H, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, segno, jd;   L:= sunlist(1); ω:= sunlist(2);    l:= sunlist(3); M:= sunlist(4);    α:= sunlist(5); δ:= sunlist(6);    az:= sunlist(7); h:= sunlist(8);    x:= sunlist(9); y:= sunlist(10);   sorg:= sunlist(11); tram:= sunlist(12); culm:= sunlist(13);   r:= sunlist(14); H:= sunlist(15); Hd:= sunlist(16);   zc:= sunlist(17); dz:= sunlist(18);    ts:= calcTS(); //   jd:= julianDay();   DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("at Long " + long + "Lat " + lat);   PRINT("");   PRINT("L longitudine media " + trunc_angle(L));   PRINT("ω long del perielio " + trunc_angle(ω));   PRINT("M anomalia vera " + trunc_angle(M));   PRINT("λ longitudine Sole " + trunc_angle(l) + " " + segno);   PRINT("r Distanza (UA) " + r);   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT(" ");   PRINT("Ang. orario "+ H + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS locale " + trunc_angle(ts(2)));   PRINT("TS apparente " + trunc_angle(ts(3)));   PRINT("Equazione del tempo (min) " + trunc_angle(E));   PRINT("Julian day " + jd(1));   PRINT("Dif asc." + DA + " SAD " + trunc_angle(SAD));   PRINT("");   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT("Long elioc. Terra "+ trunc_angle(lt));   PRINT("ε epsilon " + ε);   // RETURN {L, ω, M, v, l, r, x, y};   RETURN [[L, ω], [l, M], [α,δ], [az, h], [x,y], [sorg, tram], [culm, r], [H,ts(1)], [lt, ε]]; END; moon() BEGIN   LOCAL ω,Ω, l, v, r, x, y, i;   LOCAL Ls, Ms, LongM, LatM, P, d, s;   LOCAL α,δ, ts, Hd, zc, dz, h, az;   LOCAL as, ac, DA, SAD, DD, culm, sorg, tram;   LOCAL b, dayY, xs, ys, ph, fase, eta;   LOCAL dd, hh, mm, segno;   i:= 5.13333333; // inclinazione 5°8' sull'eclittica   calcN();   // Sun   Ls:= sunlist(1); // Longitudine media   Ms:= sunlist(4); // Anomalia   // end Sun   L:= (33.231 + 13.17639653*N) MOD 360; // long. media   Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente   M:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera   ω:= (L - M) MOD 360; // perigeo medio   DD:= L - Ls; F:= L-Ω; // for the variations   LongM:= L + 6.28875 * sin(M); // calculus (h/ and above) of Moon longitude   LongM:= LongM + 0.2136 * sin(2*M); // equazione del centro   LongM:= LongM + 0.6583 * sin(2*DD); // variazione   LongM:= LongM - 0.1856 * sin(Ms); // equazione annuale   LongM:= LongM + 1.2740 * sin(2*DD-M); // evezione   LongM:= LongM - 0.1143 * sin(2*F); // riduzione all'eclittica   LongM:= LongM + 0.0588 * sin(2*DD-2*M);   LongM:= LongM + 0.0572 * sin(2*DD-M-Ms);   LongM:= LongM + 0.0533 * sin(2*DD+M);   LongM:= LongM + 0.0459 * sin(2*DD-Ms);   LongM:= LongM + 0.0410 * sin(M-Ms);   LongM:= LongM - 0.0305 * sin(M+Ms);   LongM:= LongM - 0.0348 * sin(DD);   LongM:= LongM MOD 360;   LatM := 5.1280 * sin(F); // termine principale   LatM := LatM + 0.2806 * sin(M+F);   LatM := LatM + 0.2777 * sin(M-F); // equazione del centro   LatM := LatM + 0.1732 * sin(2*DD-F); // grande ineguaglianza   LatM := LatM + 0.0554 * sin(2*DD-M+F); // evezione   LatM := LatM + 0.0462 * sin(2*DD-M-F); // evezione (2)   LatM := LatM + 0.0326 * sin(2*DD+F);   LatM := LatM MOD 360; b:= LatM; dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(LongM))) ); // declinazione (coord locali)   // α:= →HMS((acos(cos(b)*cos(LongM)/cos(δ))/15) MOD 24); // ascensione retta   α:= (acos(cos(b)*cos(IFTE(LongM>180, LongM-180, LongM))/cos(δ))); // ascensione retta   IF LongM>180 THEN α:= α +180; END;   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)   as:= cos(δ)*sin(Hd)/sin(dz); // calc az   ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   // IF (ac<0) THEN az:= ((360 + az) MOD 360); END;   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse   // Raggio quadratico medio della Terra: 6373,044737   d:= 1/sin(P); // distanza Terra-Luna in raggi terrestri   s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)   ph:= abs(LongM - Ls); // Phase   CASE   // IF ph == 0 THEN fase:="New Moon"; END;   IF ph>280 AND ph<350 THEN fase:="Decrescente"; END;   IF ph>=260 AND ph<=280 THEN fase:="Last Quart"; END;   IF ph>190 AND ph<260 THEN fase:="Gibbosa Cal."; END;   IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;   IF ph>100 AND ph<170 THEN fase:="Gibbosa Cresc."; END;   IF ph>=80 AND ph<=100 THEN fase:="1st Quart"; END;   IF ph>10 AND ph<80 THEN fase:="Crescente"; END;   DEFAULT fase:="New Moon";   END; // case   // età della Luna   dd:= IP(HMS→(ph/15));   hh:= IP(FP(HMS→(ph/15))*60);   IF (hh>24) THEN dd:= dd+1; hh:= hh MOD 24; END;   mm:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);   eta:= →HMS(dd+hh/60+mm/3600); // Age of the Moon   segno:= zodiaco(LongM);   PRINT;   PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L));   PRINT("M anomalia vera " + trunc_angle(M));   PRINT("λ Longitudine " + trunc_angle(LongM) + " " + segno);   PRINT("β Latitudine " + trunc_angle(LatM));   PRINT("ω Perigeo medio " + trunc_angle(ω));   PRINT("Ω Nodo ascendente " + trunc_angle(Ω));   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("α Ascensione retta (gradi) " + trunc_angle((α*15) MOD 360));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS (locale) " + trunc_angle(ts(2)));   PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));   PRINT("P Parallasse " + trunc_angle(P));   PRINT("d Distanza Terra - Luna "+ ROUND(d*6373.044737,3) + " km");   PRINT("Semidiametro "+  trunc_angle(→HMS(s)) + " :: " + ROUND(s*6373.044737,3) + " km");   PRINT(" ");   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");   PRINT("Età " + eta);   RETURN [[L, M],[LongM, LatM],[ω,Ω],[α,δ], [az, h], [P,0],[d,s], [sorg, tram], [culm, ts(1)],[ph, eta] ]; END; planetCalc(nameplanet) BEGIN   LOCAL ω,Ω, l, v, r, x, y, z;   LOCAL Ls, Ms, b, xs, ys, b_ec;   LOCAL x1, y1, z1, R, l_ec, H, Hd;   LOCAL ec, in, a, δ, α, ts, az, h;   LOCAL DA, SAD, tram, sorg, culm, dayY;   LOCAL as, ac, zc, dz, θ, mov, segno;   // calc of Sun (as basis)   xs:= sunlist(9); ys:= sunlist(10);   Ls:= sunlist(1); // Longitudine media   Ms:= sunlist(4); // Anomalia   // end Sun  EVAL(EXPR(nameplanet+"()")); // call planet-name function   L:= planet_const(1); // Longitudine media   ω:= planet_const(2); // perielio   Ω:= planet_const(3); // nodo   M:= planet_const(4); // anomalia   v:= planet_const(5);    ec:= planet_const(6); // eccentricità   in:= planet_const(7); // inclinazione   a:= planet_const(8);  // semiasse maggiore   θ:= planet_const(9);  // θs to calc if retrogade   b_ec:= asin(sin(in)*sin(v+ω-Ω)); // argomento di latitudine del perielio   l_ec:= acos(cos(v+ω-Ω)/cos(b_ec));  // long from eclittic   IF b_ec<0 THEN l_ec:= 360 - l_ec; END; // argomento secondario   l_ec:= →HMS((l_ec + Ω)) MOD 360; // aggiunge Nodo e riporta a 2PI   r:= a*(1-ec^2)/(1+ec*cos(v)); // 9.524899/(1+ 0.0559*cos(v)); raggio vettore   x1:= r*cos(b_ec)*cos(l_ec); y1:= r*cos(b_ec)*sin(l_ec); z1:=r*sin(b_ec); // cartesian   x:= x1+xs; y:= y1+ys; z:=z1; // cohordinates cartesian->spheric   R:= sqrt(x^2+y^2+z^2);    b:=asin(z/R); // latitudine  β   // l:=→HMS(atan(y/x)); // cart->spheric longtitudine λ   l:= atan2(y,x);   IF l<0 THEN l:= l+360; END; IF l>360 THEN l:=l-360; END;   dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper   δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(l))) ); // declinazione (coord locali)   α:= (acos(cos(b)*cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta IF l>180 THEN α:= α +180; END;   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)   as:= cos(δ)*sin(Hd)/sin(dz); // calc az   ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   DA:= →HMS(asin(tan(lat)*tan(δ))); // Differenza ascensionale   SAD:= →HMS(90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) MOD 24 ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   IF abs(l-220.4)>θ THEN mov:="diretto"; ELSE mov:="retrogrado"; END;   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi di " + nameplanet + " a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L));   PRINT("ω Perielio " + trunc_angle(ω));   PRINT("Ω Nodo " + trunc_angle(Ω));   PRINT("M Anomalia vera " + trunc_angle(M));   PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));   PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(v));   PRINT("r Raggio vettore " + r);   PRINT("Cartes. elioc.: x' " + ROUND(x1,3) + " y' " + ROUND(y1,3) + " z' " + ROUND(z1,3));   PRINT("Cartes. geoc.: x " + ROUND(x,3) + " y " + ROUND(y,3) + " z " + ROUND(z,3));   PRINT("Spher.: λ Longitudine " + trunc_angle(l) + " " + segno); // spheric   PRINT("Spher.: β Latitudine " + trunc_angle(b));   PRINT("Spher.: R Distanza " + ROUND(R,3));   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS (locale) " + trunc_angle(ts(2)));   PRINT("Diff. ascensionale " + trunc_angle(DA));   PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT("Movimento " + mov);   RETURN [[L,ω,Ω ],[M,l_ec,b_ec], [v,r,0], [x1,y1,z1],[x,y,z],[l,b,R], [α,δ,H], [az,h,ts(1)],[sorg, tram, culm]]; END; Mercurius() BEGIN   // Costanti per Mercurio   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.205615; // eccentricità   in:= 7.003; // inclinazione   a:= 0.387098; // semiasse maggiore   θ:= 35.6;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   ω:= →HMS((75.913 + 0.00004253*N) MOD 360); // perielio   Ω:= →HMS((47.157 + 0.00003244*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 23.4372*sin(M)+2.9810*sin(2*M)+0.5396*sin(3*M)+0.1098*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Venus() BEGIN   // Costanti per Venere   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.006816; // eccentricità   in:= 3.393636; // inclinazione   a:= 0.72333; // semiasse maggiore   θ:= 13.0;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((206.758 + 1.60216873*N) MOD 360); // Longitudine media   ω:= →HMS((130.154 + 0.00003757*N) MOD 360); // perielio   Ω:= →HMS((75.797 + 0.00002503*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 0.7811*sin(M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Mars() BEGIN   // Costanti per Marte   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.093309; // eccentricità   in:= 1.850303; // inclinazione   a:= 1.523678; // semiasse maggiore   θ:= 16.8;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((124.765 + 0.52407108*N) MOD 360); // Longitudine media   ω:= →HMS((334.251 + 0.00005038*N) MOD 360); // perielio   Ω:= →HMS((48.794 + 0.00002127*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 10.608*sin(M)+0.6216*sin(2*M)+0.0504*sin(3*M)+0.0047*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Jupiter() BEGIN   // Costanti per Giove   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.048335; // eccentricità   in:= 1.308736; // inclinazione   a:= 5.202561; // semiasse maggiore   θ:= 54.43;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((268.350 + 0.08312942*N) MOD 360); // Longitudine media   ω:= →HMS((12.737 + 0.00004408*N) MOD 360); // perielio   Ω:= →HMS((99.453 + 0.00002767*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5372*sin(M)+0.1662*sin(2*M)+0.007*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Saturn() BEGIN   // Costanti per Saturno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.055892; // eccentricità   in:= 2.492519; // inclinazione   a:= 9.554747; // semiasse maggiore   θ:= 65.53;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((278.774 + 0.03349788*N) MOD 360); // Longitudine media   ω:= →HMS((91.117 + 0.00005362*N) MOD 360); // perielio   Ω:= →HMS((112.79 + 0.00002391*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 6.4023*sin(M)+0.2235*sin(2*M)+0.011*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Uranus() BEGIN   // Costanti per Urano   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.046344; // eccentricità   in:= 0.772464; // inclinazione   a:= 19.21814; // semiasse maggiore   θ:= 73.92;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((248.487 + 0.01176902*N) MOD 360); // Longitudine media   ω:= →HMS((171.563 + 0.00004064*N) MOD 360); // perielio   Ω:= →HMS((73.482 + 0.00001365*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5.3092*sin(M)+0.1537*sin(2*M)+0.0062*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Neptunus() BEGIN   // Costanti per Nettuno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.008997; // eccentricità   in:= 1.779242; // inclinazione   a:= 30.10957; // semiasse maggiore   θ:= 77.63;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((86.652 + 0.00602015*N) MOD 360); // Longitudine media   ω:= →HMS((46.742 + 0.000039*N) MOD 360); // perielio   Ω:= →HMS((130.693 + 0.00003009*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 1.031*sin(M)+0.0058*sin(2*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Pluto() BEGIN   // Costanti per Plutone   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.250236; // eccentricità   in:= 17.17047; // inclinazione   a:= 39.438712; // semiasse maggiore   θ:= 79.41;  // θs to calc if retrogade   calcN();   // costanti planetarie   // Long, perielio, nodo, non precisi, considerati costanti   // L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   M:= →HMS((230.67 + 0.00397943*N) MOD 360); // Anomalia vera (calcolata direttamente)   Ω:= →HMS((109.06 + 0.00003823*N) MOD 360); // nodo   ω:= →HMS((114.27 + Ω) MOD 360); // perielio   L:= →HMS((M + ω) MOD 360); // L (non sicuro)      v:= M + 28.45*sin(M)+4.38*sin(2*M)+0.97*sin(3*M)+0.24*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-22-2015, 06:26 PM
Post: #7
 Marcel Member Posts: 177 Joined: Mar 2014
RE: Astronomy...
Hi Salvomic
We can calculate Delta T with some analytical way but this is not always precise. In fact we know the value of Delta T after the year. For instance, we will know the value for the year 2015, in 2016! This is a correction or a difference between two times scales.

In my app's, many data are in variable associate with the App. You can see these var if you click on the "var" tab in the conn kit. At this point this is more complex than your technique because many calculations most been done with these data in a FOR statement. If you check for the book of Jean Meeus on the internet, you will learn how to use that data.
Marcel
06-22-2015, 06:34 PM (This post was last modified: 06-22-2015 06:37 PM by salvomic.)
Post: #8
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Astronomy...
(06-22-2015 06:26 PM)Marcel Wrote:  Hi Salvomic
We can calculate Delta T with some analytical way but this is not always precise. In fact we know the value of Delta T after the year. For instance, we will know the value for the year 2015, in 2016! This is a correction or a difference between two times scales.

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

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

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

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

thank you
Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-22-2015, 06:48 PM
Post: #9
 Marcel Member Posts: 177 Joined: Mar 2014
RE: Astronomy...
Hi,
One of the chapter in the book describe a way to solve the Kepler equation (Chapter 29). You will find in this chapter a program to do it.
Marcel
06-22-2015, 07:29 PM
Post: #10
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Astronomy...
(06-22-2015 06:48 PM)Marcel Wrote:  Hi,
One of the chapter in the book describe a way to solve the Kepler equation (Chapter 29). You will find in this chapter a program to do it.
Marcel

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-23-2015, 03:16 PM (This post was last modified: 06-23-2015 10:25 PM by salvomic.)
Post: #11
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Astronomy...
New version, greater precision (especially for the Moon) and rationalization...

Code:
 data(); calcN(); calcTS(); julianDay(); julianDayAtZero(); tempo(); precession(); parametri(); sunCalc(); planetCalc(); planets(); sun(); moon(); Mercurius(); Venus(); Mars(); Jupiter(); Saturn(); Uranus(); Neptunus(); Pluto(); ascendent(); zodiaco(); trunc_angle(); simp_angle(); atan2(); thisday:= 2000.0101; lstd:=0; gmt:= 12; long:= 0; lat:=  0; m:= 1; dst:= 0; utc:=0; deltaT:= 68; // 2015 datetimelist:= MAKELIST(X,X,1,6); sunlist:= MAKELIST(X,X,1,25); planet_const:= MAKELIST(X,X,1,10); ε:= →HMS(23.4392911);  // intial value, approx 23.45 (23°26'21") EXPORT Effemeridi() // from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015 // Credits: Serge Bouiges, Calcule astronomique pour amateurs // Jean Meeus, Astronomical Algorithms // Thanks Didier Lachieze, Marcel Pelletier, Eried BEGIN   LOCAL ch;   HAngle:=1;   data(); // input data and calc Sun   CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendent", "Quit");   CASE   IF ch==1 THEN sun(); END;   IF ch==2 THEN moon(); END;   IF ch==3 THEN planets(); END;   IF ch==4 THEN tempo(); END;   IF ch==5 THEN parametri(); END;   IF ch==6 THEN ascendent(); END;   IF ch==7 THEN HAngle:=0; RETURN; END;   DEFAULT   END; // case END; data() BEGIN   LOCAL dd, mm, yy, hh, mn, ss, loct;   LOCAL t, jd;   yy:= IP(Date);   mm:= IP(FP(Date)*100);   dd:= FP(Date*100)*100;   hh:= IP(HMS→(Time));   mn:= IP(FP(HMS→(Time))*60);   ss:= IP(FP(FP(FP(HMS→(Time))*60))*60);   loct:= →HMS(hh+mn/60+ss/3600);   INPUT({ {m,[0],{15,15,1}}, {D,[0],{50,15,1}}, {Y,[0],{80,15,1}},   {lstd,[0],{65,30,2}},   {dst,0,{40,15,3}}, {utc,[0],{80,15,3}},    {long,[0],{20,25,4}},    {lat,[0],{70,25,4}},   {deltaT, [0], {20,25,5}} },   "Data: Use Shift+9 for H°M′S″",   {"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",    "Long (-E):","Lat (+N)", "delta T"},   {"Enter month", "Enter day", "Enter year (1901-2099)",    "Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",   "Longitude", "Latitude", "Delta T (sec.)"},    {1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});   // Ragusa (-14°43′41″, 36°55′15″) -    // Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)  //  D:=EVAL(D); m:=EVAL(m); Y:=EVAL(Y);   datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };   // Greenwich Mean Time   gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)   thisday:= Y+m/100+D/10000;   jd:= julianDay();   t:=(jd(1)-2451545)/36525;   U:= t/100;    // ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)   ε:= 23°26′21.448″ - 1°18′00.93″*U-1.55*U^2+1999.25*U^3-51.38*U^4-249.67*U^5;   ε:= ε-39.05*U^6+7.12*U^7+27.87*U^8+5.79*U^9+2.45*U^10; // Laskar, Meeus   ε:= →HMS(ε);   sunCalc(); // calculate Sun data END; planets() BEGIN   LOCAL ch, nameplanet;   INPUT({{ch, nameplanet:={"Mercurius","Venus","Mars", "Jupiter",    "Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},   "Choose planet",    "planet: ", "Choose the planet to calculate an press OK", 0,5 );   planetCalc(nameplanet(ch)); END; calcN() BEGIN   N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT   RETURN N; END; calcTS() BEGIN   LOCAL T, TS, TSd, N0, TSL, TSapp;   LOCAL θ0, θ00, jd, jd0;   calcN();   jd:= julianDay();   jd0:= julianDayAtZero();   T:= (jd(1)-2451545.0)/36525;   θ0:= 280.46061837+360.98564736629*(jd(1)-2451545.0)+0.000387933*T^2-T^3/38710000; // DEG at 0h   θ0:= θ0 MOD 360;   θ00:= 280.46061837+360.98564736629*(jd0(1)-2451545.0)+0.000387933*T^2-T^3/38710000; // DEG at 0h   θ00:= θ00 MOD 360;    N0:= DDAYS(1901.0101, thisday)+1+0; // N at 0 GMT (for TS)   TS:=  ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds   TS:= →HMS(TS); // in HMS    TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST   TSL:= →HMS(TSL); // in HMS   TSapp:= (TS + gmt) MOD 24; // apparent ST (ST at Greenwitch at local our)   TSapp:= →HMS(TSapp); // in HMS   TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees   RETURN {TS, TSL, TSapp, TSd, θ0, θ00, T}; END; julianDay() BEGIN   LOCAL yy, mm, dd,dde, hh, mn, ss;   LOCAL aa,bb, jd, jde;   yy:= Y; mm:= m;   hh:= datetimelist(4);   mn:= datetimelist(5);   ss:= datetimelist(6);   dd:= D + hh/24 + mn/(60*24) + ss/(60*60*24);   dde:= D + hh/24 + mn/(60*24) + (ss+deltaT)/(60*60*24);   IF (m=1 OR m=2) THEN yy:=Y-1; mm:=m+12;END;   IF yy*100+mm+dd>=158225 THEN   aa:=IP(ABS(yy/100));   bb:=2-aa+IP(aa/4);   ELSE bb:=0; END;   jd:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dd+bb-1524.5;   jde:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dde+bb-1524.5;   RETURN {jd, jde}; END; julianDayAtZero() BEGIN   LOCAL yy, mm, dd,dde, hh, mn, ss;   LOCAL aa,bb, jd, jde;   yy:= Y; mm:= m;   hh:= 0;   mn:= 0;   ss:= 0;   dd:= D;   dde:= D + (ss+deltaT)/(60*60*24);   IF (m=1 OR m=2) THEN yy:=Y-1; mm:=m+12;END;   IF yy*100+mm+dd>=158225 THEN   aa:=IP(ABS(yy/100));   bb:=2-aa+IP(aa/4);   ELSE bb:=0; END;   jd:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dd+bb-1524.5;   jde:=IP(365.25*(yy+4716))+IP(30.6001*(mm+1))+dde+bb-1524.5;   RETURN {jd, jde}; END; tempo() BEGIN   LOCAL n, jd, TSid;   n:= calcN();   jd:= julianDay();   TSid:= calcTS();   PRINT;   PRINT("Date " + m + " " + D + " " + Y);   PRINT("Julian Day " + jd(1));   PRINT("Julian Day Effem. " + jd(2));   PRINT("N (from 1901jan1) " + n);   PRINT("T (from JD) " + TSid(7));   PRINT("Tempo Siderale 0h " + trunc_angle(TSid(1)));   PRINT("Tempo Siderale locale " + trunc_angle(TSid(2)));   PRINT("Tempo Siderale apparente " + trunc_angle(TSid(3)));   PRINT("Tempo Siderale (DEG) " + trunc_angle(TSid(4)));   PRINT("θ0 Mean ST 0h " + trunc_angle(TSid(5)/15));   PRINT("θ0 Mean ST 0h (DEG) " + trunc_angle(TSid(5)));   PRINT("θ00 MST 00h GMT " + trunc_angle(TSid(6)/15));   PRINT("θ00 MST 00h GMT (DEG) " + trunc_angle(TSid(6)));   PRINT("Time " + lstd + " GMT " + gmt);   RETURN({ROUND(jd(1),5), ROUND(jd(2),5), ROUND(n, 5), trunc_angle(TSid(1)), trunc_angle(TSid(2)),      trunc_angle(TSid(3)), trunc_angle(TSid(4)), trunc_angle(TSid(5)) }); END; deltaEpsilonPsi() BEGIN   // Nutazione Δψ (longitudine), Δε (obliquità)   LOCAL L1, deltaPsi, deltaEpsilon, Ω, jd;   jd:= julianDayAtZero();   T:= (jd(2)-2451545)/36525; // use JDE   Ω:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)   L:= (→HMS(280.4665)+36000.76983*T+→HMS(0.0003032*T^2)) MOD 360; // Long media Sole   L1:= →HMS(218.3164591) + →HMS(481267.88134236)*T-0.0013268*T^2; //  Long media Luna   L1:= L1+(T^3)/538841-(T^4)/65194000;   L1:= L1 MOD 360;   deltaPsi:= -0°0′17.20″*sin(Ω)-0°0′1.32″*sin(2*L)-0°0′0.23″*sin(2*L1)+0°0′0.21″*sin(2*Ω);   deltaEpsilon:= 0°0′9.2″*cos(Ω)+0°0′0.57″*cos(2*L)+0°0′0.10″*cos(2*L1)-0°0′0.09″*cos(2*Ω);   RETURN({deltaPsi, deltaEpsilon});    END; precession(alfa, delta) BEGIN   // precessione dati asce retta e declinazione   LOCAL m,n,n1, deltaAlfa, deltaDelta, Tsec;   Tsec:= IP((Y-2000)/100)+1; // T century from 2000.0   m:= 0°0′3.07496″+0°0′0.00186″*Tsec;   n:= 0°0′1.33621″+0°0′0.00057″*Tsec; // sec   n1:= 0°0′20.0431″+0°0′0.0085″*Tsec; // ° ' "   deltaAlfa:= m+n*sin(alfa)*tan(delta); // Δα   deltaDelta:= n1*cos(alfa); // Δδ   RETURN({deltaAlfa, deltaDelta}); END; parametri() BEGIN   LOCAL dep, prec, α, δ;   dep:= deltaEpsilonPsi();   α:= sunlist(5); δ:= sunlist(6);   prec:= precession(α, δ);   // Nutazione Δψ (longitudine), Δε (obliquità), ε (inclinazione eclittica)...   PRINT;   PRINT("ε Inclinazione eclittica "+ ε);   PRINT("Δψ Nutazione longitudine "+ dep(1));   PRINT("Δε Nutazione longitudine "+ dep(2));   PRINT(" ");   PRINT("Δα Precessione AR "+ prec(1));   PRINT("Δδ Precessione dec "+ prec(2));   RETURN({ε, dep(1), dep(2), prec(1), prec(2)}); // ε, Δψ, Δε, Δα, Δδ END; ascendent() BEGIN   LOCAL ASC, ts, segno;   ts:= calcTS();   ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));   IF ASC<0 THEN ASC:= ASC+360; END; IF ASC>360 THEN ASC:=ASC-360; END;   segno:= zodiaco(ASC); RETURN {→HMS(ASC), segno}; END; zodiaco(astro) BEGIN   LOCAL segno;   CASE  // Segni e case   IF astro>=0   AND astro <30  THEN segno:="Ariete"; END;   IF astro>=30  AND astro <60  THEN segno:="Toro"; END;   IF astro>=60  AND astro <90  THEN segno:="Gemelli"; END;   IF astro>=90  AND astro <120 THEN segno:="Cancro"; END;   IF astro>=120 AND astro <150 THEN segno:="Leone"; END;   IF astro>=150 AND astro <180 THEN segno:="Vergine"; END;   IF astro>=180 AND astro <210 THEN segno:="Bilancia"; END;   IF astro>=210 AND astro <240 THEN segno:="Scorpione"; END;   IF astro>=240 AND astro <270 THEN segno:="Sagittario"; END;   IF astro>=270 AND astro <300 THEN segno:="Capricorno"; END;   IF astro>=300 AND astro <330 THEN segno:="Acquario"; END;   IF astro>=330 AND astro <360 THEN segno:="Pesci"; END;   END; // case   RETURN segno; END; atan2(y,x) BEGIN // atan2 returns coordinates in correct angle for all four quadrants (DEG)     CASE     IF x==0 AND y==0 then return "undefined"; end;         // x=0, y=0       IF x>0 then return atan(y/x); end;                     // x>0     IF x<0 AND y>=0 then return atan(y/x)+180; end;         // x<0, y>=0     IF x==0 AND y>0 then return 90; end;                 // x=0, y>0     IF x<0 AND y<=0 then return atan(y/x)-180; end;         // x<0, y<0     IF x==0 AND y<0 then return -90; end;                // x=0, y<0    END;   RETURN;  END; trunc_angle(ang) // trunc decimal from a DEG form angle BEGIN →HMS((IP(HMS→(ang)*3600)/3600)); END; simp_angle(ang) BEGIN LOCAL angle; angle:=360*FP(ang/360); IF angle<0 THEN angle:=angle+360 END; RETURN angle; END; kepler(ec, M)  // needs eccentricity and true anomaly M BEGIN   LOCAL M1, j, u;   HAngle:=0; // RAD   M1:= HMS→(M)*PI/180; // anomalia RAD   F:= SIGN(M1); M1:= ABS(M1)/(2*PI);   M1:=(M1-IP(M1))*2*PI*F;   IF M1<0 THEN M1:= M1+2*PI; END;   F:= 1;    IF M1>PI THEN F:= -1; END;   IF M1>PI THEN M1:=2*PI-M; END;   u:= PI/2; D:=PI/4;   FOR j FROM 1 TO 33 DO // 33 iterations (Meeus, Roger Sinnot)   M2:= u-ec*sin(u); u:= u+D*SIGN(M1-M2); D:= D/2;   END; // for   u:= u * F; // anomalia eccentrica    HAngle:=1;   u:= u*180/PI;   RETURN u; END; sunCalc() BEGIN   LOCAL ω, l, v, r, x, y, dz;   LOCAL α, δ, dayY, H, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, h0, H0, temp;   LOCAL ec, a, u, M1, M2, j, t;   LOCAL m0, m1, m2, jd, c, lapp, Ω;   calcN();   jd:= julianDay();   a:= 1.00000023; // semiasse maggiore   t:= N/36525;   T:= (jd(2)-2451545)/36525; // use JDE   // ec:= 0.01675062-0.0000418*t-0.000000137*t^2; // eccentricità (RAD)   ec:= 0.016708617-0.000042037*T-0.0000001236*T^2; // eccentricità Meeus   // ec:= ec*180/PI; // DEG   h0:= -0°50′; // 'standard' altitude of Sun   // L:= →HMS(278.9654+0.985647342*N + 0.0003*t^2) MOD 360; // long media   L:= (→HMS(280.4665)+36000.76983*T+→HMS(0.0003032*T^2)) MOD 360; // Long media Sole Meeus   ω:= →HMS(281.238 + 0.000047067*N + 0.00045*t^2) MOD 360; // long del perielio   // M:= (L - ω); // anomalia vera   //IF M<0 THEN M:= (360+M); END; // potrebbe essere negativa   Ω:= (→HMS(125.04)-→HMS(1934.136)*T) MOD 360; // Nutazione e aberrazione (nodo)   M:= →HMS(357.52910)+→HMS(35999.05030)*T+→HMS(0.0001559)*T^2-→HMS(0.00000048)*T^3; // anom vera Meeus   M:= M MOD 360;   u:=atan2(sin(M),(cos(M)-ec)); // only for small eccentricity   c:= (→HMS(1.914600)-→HMS(0.004817)*T-→HMS(0.000016)*T^2)*sin(M);  // equation of center   c:= c+(→HMS(0.019993)-→HMS(0.000101)*T)*sin(2*M)+(→HMS(0.000290))*sin(3*M);   // v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio (Bouiges)   // v:= acos((cos(u)-ec)/(1-ec*cos(u))); // kepler   v:= M + c; // true anomaly, Meeus   // l:= (v + ω) MOD 360; // Longitudine Sole   l:= (L+c) MOD 360; // Ө Longitudine Sole, Meeus   lt:= l + 180; // Longitudine eliocentrica della Terra   lapp:= l - →HMS(0.00569)-→HMS(0.00478)*sin(Ω); // Longitudine apparente (ref true equinox)   lapp:= simp_angle(lapp);   // r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)   r:= (1.000001018*(1-ec^2))/(1+ec*cos(v)); // Distanza Sole-Terra (UA) Meeus   x:= r * cos(l); y:= r * sin(l);   dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper   δ:= →HMS((asin(sin(ε)*sin(l))) );   // α:= (acos(cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta   // IF l>180 THEN α:= α +180; END;   α:= →HMS(atan2(cos(ε)*sin(l),cos(l))); // Ascensione retta, Meeus (b never > 1.2", l=Ө)   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale   // alt:= →HMS(90 - z); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd)); // Azimuth   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   IF (az<0) THEN az:= 360+az; END;   E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   //   SAD:= →HMS(acos(-tan(lat)*tan(δ))); // Semiarco Diurno  culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione, transito  sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)  tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset) // temp:= (sin(h0)-sin(lat)*sin(δ))/(cos(lat)*cos(δ)); //IF temp<-1 THEN temp:= temp +1; END; IF temp>1 THEN temp:=temp-1; END; //H0:= acos(temp); // m0:= ((α*15 + long - ts(6)) / 360); //IF m0<0 THEN m0:= m0 +1; END; IF m0>1 THEN m0:=m0-1; END; // culm:= m0*24; //temp:= ts(3)*15+360.985647*m0; // H1:= temp - long - α*15; //        IF temp>360 THEN temp:= temp MOD 360; END; //culm:= (m0 - H/360) * 24 ; // sorg:= (culm - H0/360) *24; // tram:= (culm + H0/360) *24;   sunlist:= {L,ω,l, M, α,δ, az, h, x,y,     sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA,     SAD, E, lt, v, c, Ω, lapp, ec};   RETURN sunlist; END; sun() BEGIN   LOCAL ω, l, v, r, x, y, dz;   LOCAL α, δ, dayY, H, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, segno, jd, v, c;   LOCAL Ω, lapp, ec;   L:= sunlist(1); ω:= sunlist(2);    l:= sunlist(3); M:= sunlist(4);   Ω:= sunlist(26); lapp:=sunlist(27);    α:= sunlist(5); δ:= sunlist(6);    az:= sunlist(7); h:= sunlist(8);    x:= sunlist(9); y:= sunlist(10);   sorg:= sunlist(11); tram:= sunlist(12); culm:= sunlist(13);   r:= sunlist(14); H:= sunlist(15); Hd:= sunlist(16);   zc:= sunlist(17); dz:= sunlist(18);   v:= sunlist(24); c:= sunlist(25);    ts:= calcTS(); //   jd:= julianDay();   DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);   ec:= sunlist(28);   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("at Long " + long + " Lat " + lat);   PRINT("");   PRINT("L longitudine media " + trunc_angle(L));   PRINT("ω long del perielio " + trunc_angle(ω));   PRINT("M anomalia " + trunc_angle(M));   PRINT("v anomalia vera " + v);   PRINT("c equazione del centro " + c);   PRINT(" ");   PRINT("Ө Longitudine Sole " + trunc_angle(l) + " " + segno);   PRINT("λ longitudine apparente " + trunc_angle(lapp));   PRINT("Ω aberr, nutaz nodo "+ Ω);   PRINT("r distanza (UA) " + r);   PRINT(" ");   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT(" ");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS locale " + trunc_angle(ts(2)));   PRINT("TS apparente " + trunc_angle(ts(3)));   PRINT("θ0 Mean ST 0h " + trunc_angle(ts(5)/15));   PRINT("Equazione del tempo (min) " + trunc_angle(E));   PRINT("Julian day " + jd(1));   PRINT("Dif asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));   PRINT("");   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT(" ");   PRINT("Long elioc. Terra "+ trunc_angle(lt));   PRINT("ε epsilon " + ε);   PRINT("ec eccentricità "+ ec);   // RETURN {L, ω, M, v, l, r, x, y};   RETURN [[L, ω], [l, lapp], [M,v], [Ω, c], [α,δ], [az, h], [x,y],    [sorg, tram], [culm, r], [H,ts(1)], [lt, 0], [ε, ec]]; END; moon() BEGIN   LOCAL ω,Ω, l, v, r, x, y, i;   LOCAL Ls, Ms, LongM, LatM, P, d, s;   LOCAL α,δ, ts, Hd, zc, dz, h, az;   LOCAL as, ac, DA, SAD, DD, culm, sorg, tram;   LOCAL b, dayY, xs, ys, ph, fase, eta;   LOCAL d1, h1, m1, segno, jd, rr;   LOCAL lapp, dep, ill, ee, a1,a2,a3;   i:= 5.13333333; // inclinazione 5°8' sull'eclittica   calcN();   jd:= julianDay();   T:= (jd(2)-2451545)/36525; // use JDE   dep:= deltaEpsilonPsi();   // Sun   Ls:= sunlist(1); // Longitudine media   Ms:= sunlist(4); // Anomalia   // end Sun   //L:= (33.231 + 13.17639653*N) MOD 360; // long. media   L:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+(T^3)/538841-(T^4)/65194000);   //  Long media Luna Meeus   // Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente   Ω:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)   Ω:= simp_angle(Ω);   // M:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera   M:= 134.9634114+477198.8676313*T+0.0089970*T^2+(T^3)/69699-(T^4)/14712000; // anomalia media   M:= simp_angle(M);   ω:= (L - M) MOD 360; // perigeo medio   // DD:= L - Ls; F:= L-Ω; // for the variations   DD:= 297.8502042+445267.1115168*T-0.0016300*T^2+(T^3)/545868-(T^4)/113065000; // mean elongation   DD:= simp_angle(DD);   F:= 93.2720993+483202.0175273*T-0.0034029*T^2-(T^3)/3526000+(T^4)/863310000; // dist from asc. node   F:= simp_angle(F);   ee:= 1-0.002516*T-0.0000074*T^2;  // multiply those w/ Ms *ee   a1:= →HMS(119.75)+→HMS(131.849)*T; // action of Venus   a2:= →HMS(53.09)+→HMS(479264.290)*T; // action of Jupiter   a3:= →HMS(313.45)+→HMS(481266.484)*T;    LongM:= L     + 6.288774 * sin(M); // calc Moon longitude   LongM:= LongM + 1.274027 * sin(2*DD-M); // evezione   LongM:= LongM + 0.658314 * sin(2*DD); // variazione   LongM:= LongM + 0.213618 * sin(2*M); // equazione del centro   LongM:= LongM - 0.185116 * sin(Ms) *ee; // equazione annuale   LongM:= LongM - 0.114332 * sin(2*F); // riduzione all'eclittica   LongM:= LongM + 0.058793 * sin(2*DD-2*M);   LongM:= LongM + 0.057066 * sin(2*DD-M-Ms) *ee;   LongM:= LongM + 0.053322 * sin(2*DD+M);   LongM:= LongM + 0.045758 * sin(2*DD-Ms) *ee;   LongM:= LongM - 0.040923 * sin(Ms-M) *ee;   LongM:= LongM - 0.034720 * sin(DD);   LongM:= LongM - 0.030383 * sin(M+Ms) *ee;   LongM:= LongM + 0.015327 * sin(2*DD-2*F);   LongM:= LongM - 0.012528 * sin(M+2*F);   LongM:= LongM + 0.010980 * sin(M-2*F);   LongM:= LongM + 0.010980 * sin(M-2*F);   LongM:= LongM + 0.010675 * sin(4*DD-M);   LongM:= LongM + 0.010034 * sin(3*M);   LongM:= LongM + 0.008548 * sin(4*DD-2*M);   LongM:= LongM - 0.007888 * sin(2*DD+Ms-M) *ee;   LongM:= LongM - 0.006766 * sin(2*DD+Ms) *ee;   LongM:= LongM - 0.005163 * sin(DD-M);   LongM:= LongM + 0.004987 * sin(DD+Ms) *ee;   LongM:= LongM + 0.004036 * sin(2*DD-Ms+M);   LongM:= LongM + 0.003994 * sin(2*DD+2*M);   LongM:= LongM + 0.003861 * sin(4*DD);   LongM:= LongM + 0.003665 * sin(2*DD-3*M);   LongM:= LongM - 0.002689 * sin(Ms-2*M) *ee;   LongM:= LongM - 0.002602 * sin(2*DD-M+2*F);   LongM:= LongM + 0.002390 * sin(2*DD-Ms-2*M) *ee;   LongM:= LongM - 0.002348 * sin(DD+M);   LongM:= LongM + 0.002236 * sin(2*DD-2*Ms) *ee^2;   LongM:= LongM - 0.002120 * sin(Ms+2*M) *ee;   LongM:= LongM - 0.002069 * sin(2*Ms) *ee^2;   LongM:= LongM + 0.002048 * sin(2*DD-2*Ms-M) *ee^2;   LongM:= LongM - 0.001773 * sin(2*DD+M-2*F);   LongM:= LongM - 0.001595 * sin(2*DD+2*F);   LongM:= LongM + 0.001215 * sin(4*DD-Ms-M) *ee;   LongM:= LongM - 0.001110 * sin(2*M+2*F);   LongM:= LongM - 0.000892 * sin(3*DD-M);   LongM:= LongM - 0.000810 * sin(2*DD+Ms+M) *ee;   LongM:= LongM + 0.000759 * sin(4*DD-Ms-2*M) *ee;   LongM:= LongM - 0.000713 * sin(2*Ms-M) * ee^2;   LongM:= LongM - 0.000700 * sin(2*DD+2*Ms-M) *ee^2;   LongM:= LongM + 0.000691 * sin(2*DD+Ms-2*M) *ee;   LongM:= LongM + 0.000596 * sin(2*DD-Ms-2*F) *ee;   LongM:= LongM + 0.000549 * sin(4*DD+M);   LongM:= LongM + 0.000537 * sin(4*M);   LongM:= LongM + 0.000520 * sin(4*DD-Ms) *ee;   LongM:= LongM - 0.000487 * sin(DD-2*M);   LongM:= LongM - 0.000399 * sin(2*DD+Ms-2*F) *ee;   LongM:= LongM - 0.000381 * sin(2*M-2*F);   LongM:= LongM + 0.000351 * sin(DD+Ms+M) *ee;   LongM:= LongM - 0.000340 * sin(3*DD-2*M);   LongM:= LongM + 0.000330 * sin(4*DD-3*M);   LongM:= LongM + 0.000327 * sin(2*DD-Ms+2*M) *ee;   LongM:= LongM - 0.000323 * sin(2*Ms+M) *ee^2;   LongM:= LongM + 0.000299 * sin(DD+Ms-M) *ee;   LongM:= LongM + 0.000294 * sin(2*DD+3*M);   LongM:= LongM + (3958*sin(a1) + 1962*sin(L-F)+318*sin(a2))/1000000;   LongM:= simp_angle(LongM);   lapp:= LongM + dep(1); // longitudine apparente (+ nutazione long, Δψ)   LatM :=        5.128122 * sin(F); // termine principale   LatM := LatM + 0.280026 * sin(M+F);   LatM := LatM + 0.277693 * sin(M-F); // equazione del centro   LatM := LatM + 0.173237 * sin(2*DD-F); // grande ineguaglianza   LatM := LatM + 0.055413 * sin(2*DD-M+F); // evezione   LatM := LatM + 0.046271 * sin(2*DD-M-F); // evezione (2)   LatM := LatM + 0.032573 * sin(2*DD+F);   LatM := LatM + 0.017198 * sin(2*M+F);   LatM := LatM + 0.009266 * sin(2*D+M-F);   LatM := LatM + 0.008822 * sin(2*M-F);   LatM := LatM + 0.008216 * sin(2*DD-Ms-F) *ee;   LatM := LatM + 0.004324 * sin(2*DD-2*M-F);   LatM := LatM + 0.004200 * sin(2*DD+M-F);   LatM := LatM - 0.003359 * sin(2*DD+Ms-F) *ee;   LatM := LatM + 0.002463 * sin(2*DD-Ms-M+F) *ee;   LatM := LatM + 0.002211 * sin(2*DD-Ms+F) *ee;   LatM := LatM + 0.002065 * sin(2*DD-Ms-M-F) *ee;   LatM := LatM - 0.001870 * sin(Ms-M-F) *ee;   LatM := LatM + 0.001828 * sin(4*DD-M-F);   LatM := LatM - 0.001794 * sin(Ms-F) *ee;   LatM := LatM - 0.001749 * sin(3*F);   LatM := LatM - 0.001565 * sin(Ms-M+F) *ee;   LatM := LatM - 0.001491 * sin(DD+F);   LatM := LatM - 0.001475 * sin(Ms+M+F) *ee;   LatM := LatM - 0.001410 * sin(Ms+M-F) *ee;   LatM := LatM - 0.001344 * sin(Ms-F) *ee;   LatM := LatM - 0.001335 * sin(DD-F);   LatM := LatM + 0.001107 * sin(3*M+F);   LatM := LatM + 0.001021 * sin(4*DD-F);   LatM := LatM + 0.000833 * sin(4*DD-M+F);   LatM := LatM + 0.000777 * sin(M-3*F);   LatM := LatM + 0.000671 * sin(4*DD-2*M+F);   LatM := LatM + 0.000607 * sin(2*DD-3*F);   LatM := LatM + 0.000596 * sin(2*DD-F);   LatM := LatM + 0.000491 * sin(2*DD-M+F);   LatM := LatM - 0.000451 * sin(2*DD-M-F);   LatM := LatM + 0.000439 * sin(2*DD+F);   LatM := LatM + 0.000422 * sin(2*M-F);   LatM := LatM + 0.000421 * sin(2*DD+M-F);   LatM := LatM - 0.000366 * sin(2*M-F);   LatM := LatM - 0.000351 * sin(2*DD+Ms-F) *ee;   LatM := LatM + 0.000331 * sin(4*DD+F);   LatM := LatM + 0.000315 * sin(2*DD-Ms+M+F) *ee;   LatM := LatM + 0.000302 * sin(2*DD-2*Ms-F) *ee^2;   LatM := LatM - 0.000283 * sin(M+3*F);   LatM := LatM - 0.000229 * sin(2*DD+Ms+M-F) *ee;   LatM := LatM + 0.000223 * sin(DD+Ms-F) *ee;   LatM := LatM + 0.000223 * sin(DD+Ms+F) *ee;   LatM := LatM - 0.000220 * sin(Ms-2*M-F) *ee;   LatM := LatM - 0.000220 * sin(2*DD+Ms-M-F) *ee;   LatM := LatM - 0.000185 * sin(DD+M+F);   LatM := LatM + 0.000181 * sin(2*DD-Ms-2*M-F) *ee;   LatM := LatM - 0.000177 * sin(Ms+2*M+F);   LatM := LatM + 0.000176 * sin(4*DD-2*M-F);   LatM := LatM + 0.000166 * sin(4*DD-Ms-M+F) *ee;   LatM := LatM - 0.000164 * sin(DD+M-F);   LatM := LatM + 0.000132 * sin(4*DD+M-F);   LatM := LatM - 0.000119 * sin(DD-M-F);   LatM := LatM + 0.000115 * sin(4*DD-Ms-F) *ee;   LatM := LatM + 0.000107 * sin(2*DD-2*Ms+F) *ee^2;   LatM := LatM + (0-2235*sin(L)+382*sin(a3)+175*sin(a1-F)     +175*sin(a1+F)+127*sin(L-M)-115*sin(L+M)) / 1000000;   // correzioni per la distanza dalla Terra (Δ)   rr:=0;   rr:= rr -20.905355 * cos(M); // coseni   rr:= rr - 3.699111 * cos(2*DD-M); // evezione   rr:= rr - 2.955968 * cos(2*DD); // variazione   rr:= rr - 0.569925 * cos(2*M); // equazione del centro   rr:= rr + 0.048888 * cos(Ms) *ee; // equazione annuale   rr:= rr - 0.003149 * cos(2*F); // riduzione all'eclittica   rr:= rr + 0.246158 * cos(2*DD-2*M);   rr:= rr - 0.152138 * cos(2*DD-M-Ms) *ee;   rr:= rr - 0.170733 * cos(2*DD+M);   rr:= rr - 0.204586 * cos(2*DD-Ms) *ee;   rr:= rr - 0.129620 * cos(Ms-M);   rr:= rr + 0.108743 * cos(DD);   rr:= rr + 0.104755 * cos(M+Ms) *ee;   rr:= rr + 0.010321 * cos(2*DD-2*F);   rr:= rr - 0.000000 * cos(M+2*F); // null   rr:= rr + 0.079661 * cos(M-2*F);   rr:= rr - 0.034782 * cos(4*DD-M);   rr:= rr - 0.023210 * cos(3*M);   rr:= rr - 0.021636 * cos(4*DD-2*M);   rr:= rr + 0.024208 * cos(2*DD+Ms-M) *ee;   rr:= rr + 0.030824 * cos(2*DD+Ms) *ee;   rr:= rr - 0.008379 * cos(DD-M);   rr:= rr - 0.016675 * cos(DD+Ms) *ee;   rr:= rr - 0.012831 * cos(2*DD-Ms+M);   rr:= rr - 0.010445 * cos(2*DD+2*M);   rr:= rr - 0.011650 * cos(4*DD);   rr:= rr + 0.014403 * cos(2*DD-3*M);   rr:= rr - 0.007003 * cos(Ms-2*M) *ee;   rr:= rr - 0.000000 * cos(2*DD-M+2*F); // null   rr:= rr + 0.010056 * cos(2*DD-Ms-2*M) * ee;   rr:= rr + 0.006322 * cos(DD+M);   rr:= rr - 0.009884 * cos(2*DD-2*Ms) * ee;   rr:= rr + 0.005751 * cos(Ms+2*M) *ee;   rr:= rr - 0.000000 * cos(2*Ms) *ee^2; // null   rr:= rr - 0.004950 * cos(2*DD-2*Ms-M) *ee^2;   rr:= rr + 0.004130 * cos(2*DD+M-2*F);   rr:= rr - 0.000000 * cos(2*DD+2*F); // null   rr:= rr - 0.003958 * cos(4*DD-Ms-M) *ee;   rr:= rr - 0.000000 * cos(2*M+2*F); // null   rr:= rr + 0.003258 * cos(3*DD-M);   rr:= rr + 0.002616 * cos(2*DD+Ms+M) *ee;   rr:= rr - 0.001897 * cos(4*DD-Ms-2*M) *ee;   rr:= rr - 0.002117 * cos(2*Ms-M) * ee^2;   rr:= rr + 0.002354 * cos(2*DD+2*Ms-M) *ee^2;   rr:= rr + 0.000000 * cos(2*DD+Ms-2*M) *ee; // null   rr:= rr + 0.000000 * cos(2*DD-Ms-2*F) *ee; // null   rr:= rr - 0.001423 * cos(4*DD+M);   rr:= rr - 0.001117 * cos(4*M);   rr:= rr - 0.001571 * cos(4*DD-Ms) *ee;   rr:= rr - 0.001739 * cos(DD-2*M);   rr:= rr - 0.000000 * cos(2*DD+Ms-2*F) *ee; // null   rr:= rr - 0.004421 * cos(2*M-2*F);   rr:= rr + 0.000000 * cos(DD+Ms+M) *ee; // null   rr:= rr - 0.000000 * cos(3*DD-2*M); // null   rr:= rr + 0.000000 * cos(4*DD-3*M); // null   rr:= rr + 0.000000 * cos(2*DD-Ms+2*M) *ee; // null   rr:= rr + 0.001165 * cos(2*Ms+M) *ee^2;   rr:= rr + 0.000000 * cos(DD+Ms-M) *ee; // null   rr:= rr + 0.000000 * cos(2*DD+3*M); // null   rr:= rr + 0.008752 * cos(2*DD-M-2*F); // no sin   b:= LatM; // Latitudine   // dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(LongM))) ); // declinazione (coord locali)   //  α:= (acos(cos(b)*cos(IFTE(LongM>180, LongM-180, LongM))/cos(δ))); // ascensione retta   //IF LongM>180 THEN α:= α +180; END;   α:= atan2(sin(LongM)*cos(ε)-tan(b)*sin(ε), cos(LongM));   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)   as:= cos(δ)*sin(Hd)/sin(dz); // calc az   ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   // IF (ac<0) THEN az:= ((360 + az) MOD 360); END;   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   // Raggio quadratico medio della Terra: 6373,044737   d:= 385000.56 + rr*1000; // Δ distanza Terra-Luna in km (=1/sin(P)),rr in 10^-3   //P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse   P:= asin(6378.14/d); // Parallasse   s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)   ph:= abs(LongM - Ls); // Phase   CASE   // IF ph == 0 THEN fase:="New Moon"; END;   IF ph>280 AND ph<350 THEN fase:="Decrescente"; END;   IF ph>=260 AND ph<=280 THEN fase:="Last Quart"; END;   IF ph>190 AND ph<260 THEN fase:="Gibbosa Cal."; END;   IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;   IF ph>100 AND ph<170 THEN fase:="Gibbosa Cresc."; END;   IF ph>=80 AND ph<=100 THEN fase:="1st Quart"; END;   IF ph>10 AND ph<80 THEN fase:="Crescente"; END;   DEFAULT fase:="New Moon";   END; // case   // età della Luna   d1:= IP(HMS→(ph/15));   h1:= IP(FP(HMS→(ph/15))*60);   IF (h1>24) THEN d1:= d1+1; h1:= h1 MOD 24; END;   m1:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);   eta:= →HMS(d1+h1/60+m1/3600); // Age of the Moon   segno:= zodiaco(LongM);   ill:= (1+cos(ph))/2; // illuminated fraction   PRINT;   PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L));   PRINT("M anomalia vera " + trunc_angle(M));   PRINT("Ө Longitudine " + trunc_angle(LongM) + " " + segno);   PRINT("β Latitudine " + trunc_angle(LatM));   PRINT("λ Longitudine apparente "+ trunc_angle(lapp));   PRINT("Δψ Nutazione long "+ dep(1));   PRINT("Δε Nutazione obl. " + dep(2));   PRINT("ε inclinazione ecl. " + ε);   PRINT("ω Perigeo medio " + trunc_angle(ω));   PRINT("Ω Nodo ascendente " + trunc_angle(Ω));   PRINT(" ");   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("α Ascensione retta (gradi) " + trunc_angle((α*15) MOD 360));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS (locale) " + trunc_angle(ts(2)));   PRINT("TS (apparente) " + trunc_angle(ts(3)));   PRINT("θ0 Mean ST 0h " + trunc_angle(ts(5)/15));   PRINT("Julian day " + jd(1));   PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));   PRINT("P Parallasse " + trunc_angle(P));   PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");   PRINT("Semidiametro "+  trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,3) + " km");   PRINT(" ");   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");   PRINT("Età " + eta);   PRINT("Parte illuminata " + ill);   RETURN [[L, M],[LongM, LatM],[lapp,P], [ω,Ω],[α,δ], [az, h], [d,s], [sorg, tram], [culm, ts(1)],[ph, eta] ]; END; planetCalc(nameplanet) BEGIN   LOCAL ω,Ω, l, v, r, x, y, z;   LOCAL Ls, Ms, b, xs, ys, b_ec;   LOCAL x1, y1, z1, R, l_ec, H, Hd;   LOCAL ec, in, a, δ, α, ts, az, h;   LOCAL DA, SAD, tram, sorg, culm, dayY;   LOCAL as, ac, zc, dz, θ, mov, segno;   // calc of Sun (as basis)   xs:= sunlist(9); ys:= sunlist(10);   Ls:= sunlist(1); // Longitudine media   Ms:= sunlist(4); // Anomalia   // end Sun  EVAL(EXPR(nameplanet+"()")); // call planet-name function   L:= planet_const(1); // Longitudine media   ω:= planet_const(2); // perielio   Ω:= planet_const(3); // nodo   M:= planet_const(4); // anomalia   v:= planet_const(5);    ec:= planet_const(6); // eccentricità   in:= planet_const(7); // inclinazione   a:= planet_const(8);  // semiasse maggiore   θ:= planet_const(9);  // θs to calc if retrogade   b_ec:= asin(sin(in)*sin(v+ω-Ω)); // argomento di latitudine del perielio   l_ec:= acos(cos(v+ω-Ω)/cos(b_ec));  // long from eclittic   IF b_ec<0 THEN l_ec:= 360 - l_ec; END; // argomento secondario   l_ec:= →HMS((l_ec + Ω)) MOD 360; // aggiunge Nodo e riporta a 2PI   r:= a*(1-ec^2)/(1+ec*cos(v)); // 9.524899/(1+ 0.0559*cos(v)); raggio vettore   x1:= r*cos(b_ec)*cos(l_ec); y1:= r*cos(b_ec)*sin(l_ec); z1:=r*sin(b_ec); // cartesian   x:= x1+xs; y:= y1+ys; z:=z1; // cohordinates cartesian->spheric   R:= sqrt(x^2+y^2+z^2);    b:=asin(z/R); // latitudine  β   IF b> 360 THEN b:= b MOD 360; END;   // l:=→HMS(atan(y/x)); // cart->spheric longtitudine λ   l:= atan2(y,x);   IF l<0 THEN l:= l+360; END; IF l>360 THEN l:=l-360; END;   dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper   δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(l))) ); // declinazione (coord locali)   α:= (acos(cos(b)*cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta   IF l>180 THEN α:= α +180; END;   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)   as:= cos(δ)*sin(Hd)/sin(dz); // calc az   ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   DA:= →HMS(asin(tan(lat)*tan(δ))); // Differenza ascensionale   SAD:= →HMS(90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) MOD 24 ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   IF abs(l-220.4)>θ THEN mov:="diretto"; ELSE mov:="retrogrado"; END;   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi di " + nameplanet + " a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L));   PRINT("ω Perielio " + trunc_angle(ω));   PRINT("Ω Nodo " + trunc_angle(Ω));   PRINT("M Anomalia vera " + trunc_angle(M));   PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));   PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(v));   PRINT("r Raggio vettore " + r);   PRINT("Cartes. elioc.: x' " + ROUND(x1,3) + " y' " + ROUND(y1,3) + " z' " + ROUND(z1,3));   PRINT("Cartes. geoc.: x " + ROUND(x,3) + " y " + ROUND(y,3) + " z " + ROUND(z,3));   PRINT("Spher.: λ Longitudine " + trunc_angle(l) + " " + segno); // spheric   PRINT("Spher.: β Latitudine " + trunc_angle(b));   PRINT("Spher.: R Distanza " + ROUND(R,3));   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS (locale) " + trunc_angle(ts(2)));   PRINT("Diff. ascensionale " + trunc_angle(DA));   PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));   PRINT("Culmina a: " + trunc_angle(culm) + " UTC");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC");   PRINT("Movimento " + mov);   RETURN [[L,ω,Ω ],[M,l_ec,b_ec], [v,r,0], [x1,y1,z1],[x,y,z],[l,b,R], [α,δ,H], [az,h,ts(1)],[sorg, tram, culm]]; END; Mercurius() BEGIN   // Costanti per Mercurio   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.205615; // eccentricità   in:= 7.003; // inclinazione   a:= 0.387098; // semiasse maggiore   θ:= 35.6;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   ω:= →HMS((75.913 + 0.00004253*N) MOD 360); // perielio   Ω:= →HMS((47.157 + 0.00003244*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 23.4372*sin(M)+2.9810*sin(2*M)+0.5396*sin(3*M)+0.1098*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Venus() BEGIN   // Costanti per Venere   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.006816; // eccentricità   in:= 3.393636; // inclinazione   a:= 0.72333; // semiasse maggiore   θ:= 13.0;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((206.758 + 1.60216873*N) MOD 360); // Longitudine media   ω:= →HMS((130.154 + 0.00003757*N) MOD 360); // perielio   Ω:= →HMS((75.797 + 0.00002503*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 0.7811*sin(M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Mars() BEGIN   // Costanti per Marte   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.093309; // eccentricità   in:= 1.850303; // inclinazione   a:= 1.523678; // semiasse maggiore   θ:= 16.8;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((124.765 + 0.52407108*N) MOD 360); // Longitudine media   ω:= →HMS((334.251 + 0.00005038*N) MOD 360); // perielio   Ω:= →HMS((48.794 + 0.00002127*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 10.608*sin(M)+0.6216*sin(2*M)+0.0504*sin(3*M)+0.0047*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Jupiter() BEGIN   // Costanti per Giove   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.048335; // eccentricità   in:= 1.308736; // inclinazione   a:= 5.202561; // semiasse maggiore   θ:= 54.43;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((268.350 + 0.08312942*N) MOD 360); // Longitudine media   ω:= →HMS((12.737 + 0.00004408*N) MOD 360); // perielio   Ω:= →HMS((99.453 + 0.00002767*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5372*sin(M)+0.1662*sin(2*M)+0.007*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Saturn() BEGIN   // Costanti per Saturno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.055892; // eccentricità   in:= 2.492519; // inclinazione   a:= 9.554747; // semiasse maggiore   θ:= 65.53;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((278.774 + 0.03349788*N) MOD 360); // Longitudine media   ω:= →HMS((91.117 + 0.00005362*N) MOD 360); // perielio   Ω:= →HMS((112.79 + 0.00002391*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 6.4023*sin(M)+0.2235*sin(2*M)+0.011*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Uranus() BEGIN   // Costanti per Urano   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.046344; // eccentricità   in:= 0.772464; // inclinazione   a:= 19.21814; // semiasse maggiore   θ:= 73.92;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((248.487 + 0.01176902*N) MOD 360); // Longitudine media   ω:= →HMS((171.563 + 0.00004064*N) MOD 360); // perielio   Ω:= →HMS((73.482 + 0.00001365*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5.3092*sin(M)+0.1537*sin(2*M)+0.0062*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Neptunus() BEGIN   // Costanti per Nettuno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.008997; // eccentricità   in:= 1.779242; // inclinazione   a:= 30.10957; // semiasse maggiore   θ:= 77.63;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((86.652 + 0.00602015*N) MOD 360); // Longitudine media   ω:= →HMS((46.742 + 0.000039*N) MOD 360); // perielio   Ω:= →HMS((130.693 + 0.00003009*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 1.031*sin(M)+0.0058*sin(2*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Pluto() BEGIN   // Costanti per Plutone   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.250236; // eccentricità   in:= 17.17047; // inclinazione   a:= 39.438712; // semiasse maggiore   θ:= 79.41;  // θs to calc if retrogade   calcN();   // costanti planetarie   // Long, perielio, nodo, non precisi, considerati costanti   // L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   M:= →HMS((230.67 + 0.00397943*N) MOD 360); // Anomalia vera (calcolata direttamente)   Ω:= →HMS((109.06 + 0.00003823*N) MOD 360); // nodo   ω:= →HMS((114.27 + Ω) MOD 360); // perielio   L:= →HMS((M + ω) MOD 360); // L (non sicuro)      v:= M + 28.45*sin(M)+4.38*sin(2*M)+0.97*sin(3*M)+0.24*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-25-2015, 05:35 PM (This post was last modified: 06-27-2015 05:17 PM by salvomic.)
Post: #12
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Astronomy...
Marcel (and Manuel and others),
please, help me to find this problem now: I'm using for Sun and Moon the Meeus formulae mostly... If I calc Sun I get longitude almost equal to your program, if I calculate Moon I get the same results only if I use your app "at 0h time" (like {2015, 6, 25.0}) and not for the "justnow" time...
I'm trying using formulas like in your program, but including parameters (for sin and so on) into program, not in a variable...
Every numbers seems to be ok, but...
At the real local hour my program returns for Moon now some degrees of difference in longitude (with the first version was almost ok) if your app is used with "justnow" parameter; the same or so, if your is calculated at 0h...
So, the program could be in Julian day (but it seems to be ok) or in a global variable (I'm trying to put if off the program, as far as I can)...

Salvo

Code:
 data(); calcN(); calcTS(); meanSideralTime(); apparentSideralTime(); julianDay(); julianEphemerisDay(); julianDayAtZero(); makeDateList(); makeDateListShort(); deltaEpsilonPsi(); epsilon(); transformEclipticalToEquatorial(); transformEquatorialToHorizontal(); transformEclipticalToHorizontal(); tempo(); precession(); parametri(); transit(); sunCalc(); planetCalc(); planets(); sun(); sunAh(); moon(); Mercurius(); Venus(); Mars(); Jupiter(); Saturn(); Uranus(); Neptunus(); Pluto(); ascendent(); zodiaco(); trunc_angle(); simp_angle(); atan2(); interpolation(); thisday:= 2000.0101; lstd:=0; gmt:= 12; long:= 0; lat:=  0; m:= 1; dst:= 0; utc:=0; deltaT:= 68; // 2015 datetimelist:= MAKELIST(X,X,1,6); dateshortlist:= MAKELIST(X,X,1,3); dateshortlistgmt:= MAKELIST(X,X,1,3); sunlist:= MAKELIST(X,X,1,25); planet_const:= MAKELIST(X,X,1,10); // ε:= epsilon();  // approx 23.45 (23°26'21") EXPORT Effemeridi() // from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015 // Credits: Serge Bouiges, Calcule astronomique pour amateurs // Jean Meeus, Astronomical Algorithms // Thanks Didier Lachieze, Marcel Pelletier, Eried BEGIN   LOCAL ch, ε;   HAngle:=1;   data(); // input data and calc Sun   ε:=epsilon();   CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendent", "Quit");   CASE   IF ch==1 THEN sunCalc(); sun(); END;   IF ch==2 THEN moon(); END;   IF ch==3 THEN planets(); END;   IF ch==4 THEN tempo(); END;   IF ch==5 THEN parametri(); END;   IF ch==6 THEN ascendent(); END;   IF ch==7 THEN HAngle:=0; RETURN; END;   DEFAULT sunCalc(); sun();   END; // case END; data() BEGIN   LOCAL dd, mm, yy, hh, mn, ss, loct;   LOCAL t, JD, JDE, ε, dep, eps;   yy:= IP(Date);   mm:= IP(FP(Date)*100);   dd:= FP(Date*100)*100;   hh:= IP(HMS→(Time));   mn:= IP(FP(HMS→(Time))*60);   ss:= IP(FP(FP(FP(HMS→(Time))*60))*60);   loct:= →HMS(hh+mn/60+ss/3600);   INPUT({ {m,[0],{15,15,1}}, {D,[0],{50,15,1}}, {Y,[0],{80,15,1}},   {lstd,[0],{65,30,2}},   {dst,0,{40,15,3}}, {utc,[0],{80,15,3}},    {long,[0],{20,25,4}},    {lat,[0],{70,25,4}},   {deltaT, [0], {20,25,5}} },   "Data: Use Shift+9 for H°M′S″",   {"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",    "Long (-E):","Lat (+N)", "delta T"},   {"Enter month", "Enter day", "Enter year (1901-2099)",    "Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",   "Longitude", "Latitude", "Delta T (sec.)"},    {1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});   // Ragusa (-14°43′41″, 36°55′15″) -    // Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)   gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)   thisday:= Y+m/100+D/10000;   datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };   dateshortlist:= {Y, m, D+HMS→(lstd)/24};   dateshortlistgmt:= {Y, m, D+HMS→(gmt)/24};   // Greenwich Mean Time   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   T:=(JD-2451545)/36525;   U:= T/100;    ε:=epsilon(); // (formula IAU J2000.0)   sunCalc(); END; epsilon() BEGIN LOCAL ε, ε0, dep, eps, JD, JDE, U;   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   T:=(JD-2451545)/36525;   U:= T/100;   // ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)   ε0:= 23°26′21.448″ - 1°18′00.93″*U-1.55*U^2+1999.25*U^3-51.38*U^4-249.67*U^5;   ε0:= ε0-39.05*U^6+7.12*U^7+27.87*U^8+5.79*U^9+2.45*U^10; // Laskar, Meeus   ε0:= →HMS(ε0);   dep:= deltaEpsilonPsi();   eps:= dep(2);   ε:= ε0+eps; // true obliquity ε   RETURN ε; END; deltaEpsilonPsi() BEGIN   // Nutazione Δψ (longitudine), Δε (obliquità)   LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;   LOCAL DD;   JD:= julianDayAtZero(Y,m,D);   JDE:= julianEphemerisDay(dateshortlist);   T:= (JDE-2451545)/36525; // use JDE   DD:=297.85036+445267.111480*T-0.0019142*T^2+(T^3)/189474;  // elongazione   DD:=FP(DD/360)*360;   Ωmoon:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)   Lsun:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Sole   Lmoon:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);    //  Long media Luna   deltaPsi:= (-0°0′17.20″)*sin(Ωmoon)-(0°0′1.32″)*sin(2*Lsun)+(-0°0′0.23″)*sin(2*Lmoon)+(0°0′0.21″)*sin(2*Ωmoon);   deltaEpsilon:= (0°0′9.2″)*cos(Ωmoon)+(0°0′0.57″)*cos(2*Lsun)+(0°0′0.10″)*cos(2*Lmoon)-(0°0′0.09″)*cos(2*Ωmoon);   RETURN({deltaPsi, deltaEpsilon});    END; precession(alfa, delta) BEGIN   // precessione dati asce retta e declinazione   LOCAL m,n,n1, deltaAlfa, deltaDelta, Tsec;   Tsec:= IP((Y-2000)/100)+1; // T century from 2000.0   m:= 0°0′3.07496″+0°0′0.00186″*Tsec;   n:= 0°0′1.33621″+0°0′0.00057″*Tsec; // sec   n1:= 0°0′20.0431″+0°0′0.0085″*Tsec; // ° ' "   deltaAlfa:= m+n*sin(alfa)*tan(delta); // Δα   deltaDelta:= n1*cos(alfa); // Δδ   RETURN({deltaAlfa, deltaDelta}); END; planets() BEGIN   LOCAL ch, nameplanet;   INPUT({{ch, nameplanet:={"Mercurius","Venus","Mars", "Jupiter",    "Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},   "Choose planet",    "planet: ", "Choose the planet to calculate an press OK", 0,5 );   planetCalc(nameplanet(ch)); END; transformEclipticalToEquatorial(lambda, beta) BEGIN   // trasnform λ and β (long, lat) into α and δ (Asc retta, decl)   LOCAL ε,Epsilon,Lambda, Beta, alpha, delta;   Lambda:=HMS→(lambda);   Beta:=HMS→(beta);   Epsilon:= epsilon();   ε:= epsilon();   alpha:=atan2(sin(Lambda)*cos(Epsilon)-tan(Beta)*sin(Epsilon),cos(Lambda));   IF alpha>=360 THEN alpha:=alpha-360 END;   IF alpha<0 THEN alpha:=alpha+360 END;   delta:=asin(sin(Beta)*cos(Epsilon)+cos(Beta)*sin(Epsilon)*sin(Lambda));   alpha:=→HMS(alpha/15);   delta:=→HMS(delta);   RETURN ({alpha,delta}); END; transformEquatorialToHorizontal(Lha, Phi, Delta) BEGIN   LOCAL azimuth,altitude,lha,phi,delta;   lha:=HMS→(Lha);   phi:=HMS→(Phi);   // latitude   delta:=HMS→(Delta);   azimuth:=atan2(sin(lha),cos(lha)*sin(phi)-tan(delta)*cos(phi));   altitude:=asin(sin(phi)*sin(delta)+cos(phi)*cos(delta)*cos(lha));   azimuth:=→HMS(azimuth);   altitude:=→HMS(altitude);   RETURN({azimuth,altitude}); END; transformEclipticalToHorizontal(Lambda, Beta, Lha) BEGIN   LOCAL ε;   // transform λ and β (long, lat) into azimuth and height   LOCAL lambda, beta, epsilon,phi, lha;   LOCAL alpha,delta,azimuth,altitude;   lambda:=HMS→(Lambda); // longitudine astro   beta:=HMS→(Beta); // latitudine astro   epsilon:= epsilon();   ε:= epsilon();   lha:=HMS→(Lha);  // angolo orario   phi:=HMS→(lat); // lat geo   alpha:=atan2(SIN(lambda)*COS(epsilon)-TAN(beta)*SIN(epsilon),COS(lambda));   delta:=asin(sin(beta)*cos(epsilon)+cos(beta)*sin(epsilon)*sin(lambda));   azimuth:=atan2(SIN(lha),COS(lha)*SIN(phi)-TAN(delta)*COS(phi));   altitude:=ASIN(SIN(phi)*SIN(delta)+COS(phi)*COS(delta)*COS(lha));   azimuth:=→HMS(azimuth);   altitude:=→HMS(altitude);   RETURN ({azimuth,altitude}); END; calcN() BEGIN   N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT   RETURN N; END; calcTS() BEGIN   LOCAL T, TS, TSd, N0, TSL, TSapp, ε;   LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;   ε:= epsilon;   calcN();   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   JD0:= julianDayAtZero(Y,m,D);   dep:= deltaEpsilonPsi();   psi:=dep(1);   eps:=dep(2);   T:= (JD-2451545.0)/36525;   θ0:= meanSideralTime(dateshortlist); // MST DEG at 0h   θ0GMT:= meanSideralTime({Y,m,D}); // MST DEG at 0h at Greenwitch   // TS:=  ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds   // TS:= →HMS(TS); // in HMS   TS:= θ0GMT;   // TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST   // TSL:= →HMS(TSL); // in HMS   TSapp:= apparentSideralTime(dateshortlist); // apparent ST (ST at Greenwitch at local our) TSL:= apparentSideralTime(dateshortlist) - 1*dst;   TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees   RETURN {TS, TSL, TSapp, θ0, T}; END; meanSideralTime(lista) // lista like dateshortlist BEGIN   LOCAL T,MST, JD, ε;   ε:= epsilon;   JD:= julianDay(lista);   T:=(julianDay(lista)-2451545.0)/36525;   MST:=280.46061837+360.98564736629*(JD-2451545.0)+0.000387933*T^2-(T^3)/38710000;   MST:=FP(MST/360)*360; // θ0   IF MST<0 THEN MST:=MST+360 END;   MST:=→HMS(MST/15);   RETURN MST; END; apparentSideralTime(lista) // lista like dateshortlist BEGIN   LOCAL MST,correction, deltaPsi, AST, dep, ε;   ε:= epsilon;   MST:= meanSideralTime(lista);   dep:= deltaEpsilonPsi();   deltaPsi:= dep(1);   correction:=(HMS→(deltaPsi)*3600*COS(HMS→(ε)))/15;   AST:=→HMS(HMS→(MST)+(correction/3600));   RETURN AST; END; julianDay(lista) // lista like dateshortlist BEGIN   LOCAL y,m,d,a,b, JD;   y:=lista(1);   m:=lista(2);   d:=lista(3);   IF (m=1 OR m=2) THEN   y:=y-1;   m:=m+12   END;   IF y*100+m+d>=158225 THEN   a:=IP(ABS(y/100));   b:=2-a+IP(a/4);   ELSE   b:=0   END;   JD:=IP(365.25*(y+4716))+IP(30.6001*(m+1))+d+b-1524.5;   RETURN JD; END; julianEphemerisDay(lista) BEGIN   LOCAL JDE;   JDE:=julianDay(makeDateListShort(makeDateList(lista)+{0,0,0,0,0,deltaT}));   RETURN JDE; END; julianDayAtZero(Y,m,D) BEGIN   // JD at 0h GMT   LOCAL y,mm,d,a,b, JD;   y:=Y;   mm:=m;   d:=dateshortlist(3);   IF (mm=1 OR mm=2) THEN   y:=y-1;   mm:=mm+12   END;   IF y*100+mm+d>=158225 THEN   a:=IP(ABS(y/100));   b:=2-a+IP(a/4);   ELSE   b:=0   END;   JD:=IP(365.25*(y+4716))+IP(30.6001*(mm+1))+d+b-1524.5;   RETURN JD; END; makeDateList(lista) // lista {Y m D.d} (dateshortlist) BEGIN LOCAL y,o,m,d,h,s,t; y:=lista(1); o:=lista(2); d:=IP(lista(3)); t:=FP(lista(3)); h:=IP(t*24); t:=FP(t*24); m:=IP(t*60); t:=FP(t*60); s:=IP(t*60); lista:={y,o,d,h,m,s}; RETURN lista; END; makeDateListShort(lista) BEGIN // lista {Y m D h mn, s} LOCAL y,m,d, listacorta; y:=lista(1); m:=lista(2); d:=lista(3)+lista(4)/24+lista(5)/1440+lista(6)/86400; listacorta:={y,m,d}; RETURN listacorta; END; tempo() BEGIN   LOCAL n, JD, JDE, TSid, dep, psi, eps;   LOCAL ε;   ε:= epsilon();   n:= calcN();   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   TSid:= calcTS();   dep:= deltaEpsilonPsi();   psi:=dep(1);   eps:=dep(2);   PRINT;   PRINT("Date " + m + " " + D + " " + Y);   PRINT("Julian Day " + JD);   PRINT("Julian Day Effem. " + JDE);   PRINT("N (from 1901jan1) " + n);   PRINT("T (from JD) " + TSid(5));   PRINT(" ");   PRINT("Tempo Siderale 0h " + →HMS(TSid(1)));   PRINT("Tempo Siderale locale " + →HMS(TSid(2)));   PRINT("Tempo Siderale apparente " + →HMS(TSid(3)));   PRINT("θ0 Mean ST GMT " + →HMS(TSid(4)));   PRINT("Tempo Siderale (DEG) " + trunc_angle(TSid(1)*15));   PRINT("θ0 Mean ST GMT (DEG) " + trunc_angle(TSid(4)*15));   PRINT(" ");   PRINT("Δψ Nutazione (longit.) "+ psi);   PRINT("Δε Nutazione (obl.) "+ eps);   PRINT("ε Obliquità Ecl. "+ ε);   PRINT(" ");   PRINT("Time " + lstd + " GMT " + gmt);   RETURN({ROUND(JD,5), ROUND(JDE,5), trunc_angle(TSid(1)), trunc_angle(TSid(2)),      trunc_angle(TSid(3)), trunc_angle(TSid(4)), trunc_angle(TSid(5)) }); END; transit(alfa, delta, h, h0) BEGIN // α, δ (AR, dec), height, h0 = standard height of body   LOCAL temp, H0, tsideral, Hor, m0, m1, m2;   LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;   α:= alfa*15; // DEG   δ:= delta;   θ0:= apparentSideralTime(makeDateListShort({Y,m,D,0,0,0})) *15; // TS in DEG   ts:= calcTS();   // θ0:= ts(2)*15; // TS=θ00 GMT   temp:= (sin(h0)-sin(lat)*sin(δ))/(cos(lat)*cos(δ));   IF temp<-1 OR temp>1 THEN culm:=0; END; // above horizon?   H0:= acos(temp);   IF H0<-180 THEN H0:=H0+180; END; IF H0>180 THEN H0:= H0-180; END; // it should be from -180 to 180     Hor:=15*HMS→(ts(2)-α);  // Hour angle, local ST vs apparent ST   Hor:=simp_angle(→HMS(Hor));   // Hor:= simp_angle(tsideral - long - α);     m0:= HMS→(α + long - θ0) / 360;   IF m0<0 THEN m0:= m0 +1; END; IF m0>1 THEN m0:=m0-1; END;   tsideral:= simp_angle(θ0 + 360.985647*m0); // TS transit DEG   deltaM:= -HMS→(H/360);   culm:= (m0 + deltaM)*24;  // transit   m1:= m0-H0/360;   IF m1<0 THEN m1:= m1 +1; END; IF m1>1 THEN m1:=m1-1; END;   tsideral:= simp_angle(θ0 + 360.985647*m1); // TS rising   deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));   sorg:= ((m1 + deltaM)*24) MOD 24; // rising   m2:= m0+H0/360;   IF m2<0 THEN m2:= m2 +1; END; IF m2>1 THEN m2:=m2-1; END;   tsideral:= simp_angle(θ0 + 360.985647*m2); // TS setting   deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));   tram:= ((m2 + deltaM) * 24) MOD 24; // setting RETURN ({culm, sorg, tram}); END; parametri() BEGIN   LOCAL dep, prec, α, δ, ε;   dep:= deltaEpsilonPsi();   α:= sunlist(5); δ:= sunlist(6);   prec:= precession(α, δ);   ε:= epsilon();   // Nutazione Δψ (longitudine), Δε (obliquità), ε (inclinazione eclittica)...   PRINT;   PRINT("ε Obliquità eclittica "+ ε);   PRINT("Δψ Nutazione longitudine "+ dep(1));   PRINT("Δε Nutazione longitudine "+ dep(2));   PRINT(" ");   PRINT("Δα Precessione AR "+ prec(1));   PRINT("Δδ Precessione dec "+ prec(2));   RETURN({ε, dep(1), dep(2), prec(1), prec(2)}); // ε, Δψ, Δε, Δα, Δδ END; ascendent() BEGIN   LOCAL ASC, ts, segno, ε;   ε:= epsilon;   ts:= calcTS();   ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));   ASC:= simp_angle(ASC);   segno:= zodiaco(ASC); RETURN {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno}; END; zodiaco(astro) BEGIN   LOCAL segno;   CASE  // Segni e case   IF astro>=0   AND astro <30  THEN segno:="Ariete"; END;   IF astro>=30  AND astro <60  THEN segno:="Toro"; END;   IF astro>=60  AND astro <90  THEN segno:="Gemelli"; END;   IF astro>=90  AND astro <120 THEN segno:="Cancro"; END;   IF astro>=120 AND astro <150 THEN segno:="Leone"; END;   IF astro>=150 AND astro <180 THEN segno:="Vergine"; END;   IF astro>=180 AND astro <210 THEN segno:="Bilancia"; END;   IF astro>=210 AND astro <240 THEN segno:="Scorpione"; END;   IF astro>=240 AND astro <270 THEN segno:="Sagittario"; END;   IF astro>=270 AND astro <300 THEN segno:="Capricorno"; END;   IF astro>=300 AND astro <330 THEN segno:="Acquario"; END;   IF astro>=330 AND astro <360 THEN segno:="Pesci"; END;   END; // case   RETURN segno; END; atan2(y,x) BEGIN // atan2 returns coordinates in correct angle for all four quadrants (DEG)     CASE     IF x==0 AND y==0 then return "undefined"; end;         // x=0, y=0       IF x>0 then return atan(y/x); end;                     // x>0     IF x<0 AND y>=0 then return atan(y/x)+180; end;         // x<0, y>=0     IF x==0 AND y>0 then return 90; end;                 // x=0, y>0     IF x<0 AND y<=0 then return atan(y/x)-180; end;         // x<0, y<0     IF x==0 AND y<0 then return -90; end;                // x=0, y<0    END;   RETURN;  END; trunc_angle(ang) // trunc decimal from a DEG form angle BEGIN →HMS((IP(HMS→(ang)*3600)/3600)); END; simp_angle(ang) BEGIN LOCAL angle; angle:=360*FP(ang/360); IF angle<0 THEN angle:=angle+360 END; RETURN angle; END; kepler(ec, M)  // needs eccentricity and true anomaly M BEGIN   LOCAL M1, j, u;   HAngle:=0; // RAD   M1:= HMS→(M)*PI/180; // anomalia RAD   F:= SIGN(M1); M1:= ABS(M1)/(2*PI);   M1:=(M1-IP(M1))*2*PI*F;   IF M1<0 THEN M1:= M1+2*PI; END;   F:= 1;    IF M1>PI THEN F:= -1; END;   IF M1>PI THEN M1:=2*PI-M; END;   u:= PI/2; D:=PI/4;   FOR j FROM 1 TO 33 DO // 33 iterations (Meeus, Roger Sinnot)   M2:= u-ec*sin(u); u:= u+D*SIGN(M1-M2); D:= D/2;   END; // for   u:= u * F; // anomalia eccentrica    HAngle:=1;   u:= u*180/PI;   RETURN u; END; interpolation(lista,n) BEGIN LOCAL a,b,c,R,dif1,dif2; dif1:=ΔLIST(lista); dif2:=ΔLIST(dif1); a:=dif1(1); b:=dif1(2); c:=dif2(1); R:=lista(2)+(n/2)*(a+b+n*c); RETURN R; END; sunCalc() BEGIN   LOCAL ω, l, v, r, x, y, dz, ε;   LOCAL α, δ, dayY, H, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, h0, H0, temp;   LOCAL ec, a, u, M1, M2, j, t, tau;   LOCAL JD, JDE, c, lapp, Ω, θ0, transito;   LOCAL dep, deltapsi, longitude;   ε:= epsilon;   N:= calcN();   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   a:= 1.00000023; // semiasse maggiore   dep:= deltaEpsilonPsi();   deltapsi:= dep(1);   t:= N/36525;   T:=   (JD-2451545.0)/36525;    tau:= (JDE-2451545.0)/365250; // millenni  - JDE   ec:= 0.016708617-0.000042037*T-0.0000001236*T^2; // eccentricità Meeus   // ec:= ec*180/PI; // DEG   // L:= →HMS(278.9654+0.985647342*N + 0.0003*t^2) MOD 360; // long media Bouiges   // L:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Meeus   L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);   ω:= →HMS(281.238 + 0.000047067*N + 0.00045*t^2) MOD 360; // long del perielio   // M:= (L - ω); // anomalia vera   //IF M<0 THEN M:= (360+M); END; // potrebbe essere negativa   Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)   M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-0.00000048*T^3); // anom vera Meeus   u:=atan2(sin(M),(cos(M)-ec)); // only for small eccentricity   c:= (1.914600-0.004817*T-0.000016*T^2)*sin(M);  // equation of center   c:= c+(0.019993-0.000101*T)*sin(2*M)+(0.000290)*sin(3*M);   // v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio (Bouiges)   // v:= acos((cos(u)-ec)/(1-ec*cos(u))); // kepler   v:= M + c; // true anomaly, Meeus   // l:= (v + ω) MOD 360; // Longitudine Sole   l:= simp_angle(L+c); // Ө Longitudine Sole, Meeus   lt:= l + 180; // Longitudine eliocentrica della Terra   lapp:= simp_angle(l - 0.00569-0.00478*sin(Ω)); // Longitudine apparente (ref true equinox)   // r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)   r:= (1.000001018*(1-ec^2))/(1+ec*cos(v)); // Distanza Sole-Terra (UA) Meeus   x:= r * cos(l); y:= r * sin(l);   δ:= →HMS((asin(sin(ε)*sin(l))) );   α:= →HMS(atan2(cos(ε)*sin(l),cos(l))); // Ascensione retta, Meeus (b never > 1.2", l=Ө)   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale  // H:= (ts(2) - α) MOD 24; // Angolo orario  // Hd:= →HMS(15*H); // Angolo orario in gradi   // longitude:=-HMS→(long); // - if long -E   Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST   IF Hd<0 THEN Hd:=Hd+360 END;   IF Hd>=360 THEN Hd:=Hd-360 END;   Hd:=→HMS(Hd);   H:= Hd/15; // in HMS      temp:= sunAh(dateshortlist);   az:= temp(1);   h:=  temp(2);   // E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)   E:= 4*( L - 0.0057183 - α*15 + deltapsi*cos(ε) );  // α in DEG   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   //   SAD:= →HMS(acos(-tan(lat)*tan(δ))); // Semiarco Diurno      // culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione, transito (Bouiges)    // sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)    // tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)   // transit routine   h0:= -0°50′; // 'standard' altitude of Sun   transito:= transit(α, δ, h, h0);   culm:= transito(1);   sorg:= transito(2);   tram:= transito(3);   // end transit   sunlist:= {L,ω,l, M, α,δ, az, h, x,y,              sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA,              SAD, E, lt, v, c, Ω, lapp, ec};   RETURN sunlist; END; sun() BEGIN   LOCAL ω, l, v, r, x, y, dz, ε;   LOCAL α, δ, dayY, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, segno, JD, JDE, v, c;   LOCAL Ω, lapp, ec;   ε:= epsilon;   L:= sunlist(1); ω:= sunlist(2);    l:= sunlist(3); M:= sunlist(4);   Ω:= sunlist(26); lapp:=sunlist(27);    α:= sunlist(5); δ:= sunlist(6);    az:= sunlist(7); h:= sunlist(8);    x:= sunlist(9); y:= sunlist(10);   sorg:= sunlist(11); tram:= sunlist(12); culm:= sunlist(13);   r:= sunlist(14); H:= sunlist(15); Hd:= sunlist(16);   zc:= sunlist(17); dz:= sunlist(18);   v:= sunlist(24); c:= sunlist(25);    ts:= calcTS(); //   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);   ec:= sunlist(28);   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("at Long " + long + " Lat " + lat);   PRINT("");   PRINT("L longitudine media " + trunc_angle(L));   PRINT("ω long del perielio " + trunc_angle(ω));   PRINT("M anomalia media " + trunc_angle(M));   PRINT("v anomalia vera " + trunc_angle(v));   PRINT("c equazione del centro " + c);   PRINT(" ");   PRINT("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");   PRINT("λ longitudine apparente " + trunc_angle(lapp));   PRINT("Ω aberr, nutaz nodo "+ Ω);   PRINT("r distanza (UA) " + r);   PRINT(" ");   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("α Asc.r. (DEG) " + trunc_angle(α*15));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT(" ");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS locale " + trunc_angle(ts(2)));   PRINT("TS apparente " + ts(3));   PRINT("θ0 Mean ST GMT " + ts(4));   PRINT("Equazione del tempo (min) " + trunc_angle(E));   PRINT("Julian day " + JD);   PRINT("Dif asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));   PRINT("");   PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");   PRINT(" ");   PRINT("Long elioc. Terra "+ trunc_angle(lt));   PRINT("ε obliquità ecl. " + ε);   PRINT("ec eccentricità "+ ec);   // RETURN {L, ω, M, v, l, r, x, y};   RETURN [[L, ω], [l, lapp], [M,v], [Ω, c], [α,δ], [az, h], [x,y],    [sorg, tram], [culm, r], [H,ts(1)], [lt, 0], [ε, ec]]; END; sunAh(lista) // lista like dateshortlist BEGIN   LOCAL lha,alpha,delta,phi, ah, α;   LOCAL temp,longitude, Azimuth, Altitude, ts;   alpha:=15*HMS→(sunlist(5));   α:= sunlist(5); // α in HMS   delta:=sunlist(6);   phi:=lat;   ts:= calcTS();   // longitude:= HMS→(long); // - if long -E   // lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);   lha:= 15*HMS→(ts(2)-α);  // local ST vs apparent ST;   IF lha<0 THEN lha:=lha+360 END;   IF lha>=360 THEN lha:=lha-360 END;   lha:=→HMS(lha);   temp:= transformEquatorialToHorizontal(lha, phi, delta);   Azimuth:= HMS→(temp(1))+180;   IF Azimuth>360 THEN Azimuth:=Azimuth-360 END;   Azimuth:= →HMS(Azimuth);   Altitude:= temp(2);   ah:= {Azimuth,Altitude};   RETURN ah; END; moon() BEGIN   LOCAL ω,Ω, L1, M1, l, lsum, bsum, rsum;    LOCAL v, r, x, y, dep, ε;   LOCAL b, P, d, s, temp;   LOCAL α,δ, ts, Hd, zc, dz, h, az;   LOCAL i, as, ac, DA, SAD, DD;   LOCAL culm, sorg, tram, h0, transito;   LOCAL xs, ys, ph, fase, eta;   LOCAL d1, h1, m1, segno, JD, JDE, rr;   LOCAL lapp, ill, ee, a1,a2,a3, longitude;   i:= 5.13333333; // inclinazione 5°8' sull'eclittica   ε:= epsilon();   N:= calcN();   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   T:= (JDE-2451545.0)/36525; // use JDE   dep:= deltaEpsilonPsi();   // Sun   L:= simp_angle(280.4665+36000.76983*T+0.0003032*T^2); // Long media Sole Meeus   M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+T^3/24490000); // mean anomaly Sun      // end Sun   // L1:= (33.231 + 13.17639653*N) MOD 360; // long. media   L1:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);   //  Long media Luna Meeus   // Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente   Ω:= simp_angle(125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000); // Nodo Luna (Moon)   // M1:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera   M1:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media   ω:= (L1 - M1) MOD 360; // perigeo medio   // DD:= L1 - L; F:= L1-Ω; // for the variations   DD:= simp_angle(297.8502042+445267.1115168*T-0.0016300*T^2+T^3/545868-T^4/113065000); // mean elongation   F:= simp_angle(93.2720993+483202.0175273*T-0.0034029*T^2-T^3/3526000+T^4/863310000); // dist from asc. node   // ee:= 1-0.002516*T-0.0000074*T^2;  // multiply those w/ M *ee   a1:= simp_angle(119.75+131.849*T); // action of Venus   a2:= simp_angle(53.09+479264.290*T); // action of Jupiter   a3:= simp_angle(313.45+481266.484*T);    lsum:=0;  // da sommare a L1   lsum:= lsum + 6288774 * sin(M1); // calc Moon longitude   lsum:= lsum + 1274027 * sin(2*DD-M1); // evezione   lsum:= lsum +  658314 * sin(2*DD); // variazione   lsum:= lsum +  213618 * sin(2*M1); // equazione del centro   lsum:= lsum -  185116 * sin(M) *ee; // equazione annuale   lsum:= lsum -  114332 * sin(2*F); // riduzione all'eclittica   lsum:= lsum +   58793 * sin(2*DD-2*M1);   lsum:= lsum +   57066 * sin(2*DD-M1-M) *ee;   lsum:= lsum +   53322 * sin(2*DD+M1);   lsum:= lsum +   45758 * sin(2*DD-M) *ee;   lsum:= lsum -   40923 * sin(M-M1) *ee;   lsum:= lsum -   34720 * sin(DD);   lsum:= lsum -   30383 * sin(M1+M) *ee;   lsum:= lsum +   15327 * sin(2*DD-2*F);   lsum:= lsum -   12528 * sin(M1+2*F);   lsum:= lsum +   10980 * sin(M1-2*F);   lsum:= lsum +   10675 * sin(4*DD-M1);   lsum:= lsum +   10034 * sin(3*M1);   lsum:= lsum +    8548 * sin(4*DD-2*M1);   lsum:= lsum -    7888 * sin(2*DD+M-M1) *ee;   lsum:= lsum -    6766 * sin(2*DD+M) *ee;   lsum:= lsum -    5163 * sin(DD-M1);   lsum:= lsum +    4987 * sin(DD+M) *ee;   lsum:= lsum +    4036 * sin(2*DD-M+M1) *ee;   lsum:= lsum +    3994 * sin(2*DD+2*M1);   lsum:= lsum +    3861 * sin(4*DD);   lsum:= lsum +    3665 * sin(2*DD-3*M1);   lsum:= lsum -    2689 * sin(M-2*M1) *ee;   lsum:= lsum -    2602 * sin(2*DD-M1+2*F);   lsum:= lsum +    2390 * sin(2*DD-M-2*M1) *ee;   lsum:= lsum -    2348 * sin(DD+M1);   lsum:= lsum +    2236 * sin(2*DD-2*M) *ee^2;   lsum:= lsum -    2120 * sin(M+2*M1) *ee;   lsum:= lsum -    2069 * sin(2*M) *ee^2;   lsum:= lsum +    2048 * sin(2*DD-2*M-M1) *ee^2;   lsum:= lsum -    1773 * sin(2*DD+M1-2*F);   lsum:= lsum -    1595 * sin(2*DD+2*F);   lsum:= lsum +    1215 * sin(4*DD-M-M1) *ee;   lsum:= lsum -    1110 * sin(2*M1+2*F);   lsum:= lsum -     892 * sin(3*DD-M1);   lsum:= lsum -     810 * sin(2*DD+M+M1) *ee;   lsum:= lsum +     759 * sin(4*DD-M-2*M1) *ee;   lsum:= lsum -     713 * sin(2*M-M1) * ee^2;   lsum:= lsum -     700 * sin(2*DD+2*M-M1) *ee^2;   lsum:= lsum +     691 * sin(2*DD+M-2*M1) *ee;   lsum:= lsum +     596 * sin(2*DD-M-2*F) *ee;   lsum:= lsum +     549 * sin(4*DD+M1);   lsum:= lsum +     537 * sin(4*M1);   lsum:= lsum +     520 * sin(4*DD-M) *ee;   lsum:= lsum -     487 * sin(DD-2*M1);   lsum:= lsum -     399 * sin(2*DD+M-2*F) *ee;   lsum:= lsum -     381 * sin(2*M1-2*F);   lsum:= lsum +     351 * sin(DD+M+M1) *ee;   lsum:= lsum -     340 * sin(3*DD-2*M1);   lsum:= lsum +     330 * sin(4*DD-3*M1);   lsum:= lsum +     327 * sin(2*DD-M+2*M1) *ee;   lsum:= lsum -     323 * sin(2*M+M1) *ee^2;   lsum:= lsum +     299 * sin(DD+M-M1) *ee;   lsum:= lsum +     294 * sin(2*DD+3*M1);   lsum:= lsum + 3958*sin(a1)+1962*sin(L1-F)+318*sin(a2);   l:= (simp_angle(L1 + lsum/1000000)); // λ divided by 1 milion   lapp:= →HMS(l + dep(1)); // longitudine apparente (+ nutazione long, Δψ)   // fattori per latitudine   bsum:=        5128122 * sin(F); // termine principale latitudine   bsum:= bsum +  280026 * sin(M1+F);   bsum:= bsum +  277693 * sin(M1-F); // equazione del centro   bsum:= bsum +  173237 * sin(2*DD-F); // grande ineguaglianza   bsum:= bsum +   55413 * sin(2*DD-M1+F); // evezione   bsum:= bsum +   46271 * sin(2*DD-M1-F); // evezione (2)   bsum:= bsum +   32573 * sin(2*DD+F);   bsum:= bsum +   17198 * sin(2*M1+F);   bsum:= bsum +    9266 * sin(2*D+M1-F);   bsum:= bsum +    8822 * sin(2*M1-F);   bsum:= bsum +    8216 * sin(2*DD-M-F) *ee;   bsum:= bsum +    4324 * sin(2*DD-2*M1-F);   bsum:= bsum +    4200 * sin(2*DD+M1+F);   bsum:= bsum -    3359 * sin(2*DD+M-F) *ee;   bsum:= bsum +    2463 * sin(2*DD-M-M1+F) *ee;   bsum:= bsum +    2211 * sin(2*DD-M+F) *ee;   bsum:= bsum +    2065 * sin(2*DD-M-M1-F) *ee;   bsum:= bsum -    1870 * sin(M-M1-F) *ee;   bsum:= bsum +    1828 * sin(4*DD-M1-F);   bsum:= bsum -    1794 * sin(M+F) *ee;   bsum:= bsum -    1749 * sin(3*F);   bsum:= bsum -    1565 * sin(M-M1+F) *ee;   bsum:= bsum -    1491 * sin(DD+F);   bsum:= bsum -    1475 * sin(M+M1+F) *ee;   bsum:= bsum -    1410 * sin(M+M1-F) *ee;   bsum:= bsum -    1344 * sin(M-F) *ee;   bsum:= bsum -    1335 * sin(DD-F);   bsum:= bsum +    1107 * sin(3*M1+F);   bsum:= bsum +    1021 * sin(4*DD-F);   bsum:= bsum +     833 * sin(4*DD-M1+F);   bsum:= bsum +     777 * sin(M1-3*F);   bsum:= bsum +     671 * sin(4*DD-2*M1+F);   bsum:= bsum +     607 * sin(2*DD-3*F);   bsum:= bsum +     596 * sin(2*DD+2*M1-F);   bsum:= bsum +     491 * sin(2*DD-M+M1-F) *ee;   bsum:= bsum -     451 * sin(2*DD-2*M1+F);   bsum:= bsum +     439 * sin(3*M1-F);   bsum:= bsum +     422 * sin(2*DD+2*M1+F);   bsum:= bsum +     421 * sin(2*DD-3*M1-F);   bsum:= bsum -     366 * sin(2*DD+M-M1+F) *ee;   bsum:= bsum -     351 * sin(2*DD+M+F) *ee;   bsum:= bsum +     331 * sin(4*DD+F);   bsum:= bsum +     315 * sin(2*DD-M+M1+F) *ee;   bsum:= bsum +     302 * sin(2*DD-2*M-F) *ee^2;   bsum:= bsum -     283 * sin(M1+3*F);   bsum:= bsum -     229 * sin(2*DD+M+M1-F) *ee;   bsum:= bsum +     223 * sin(DD+M-F) *ee;   bsum:= bsum +     223 * sin(DD+M+F) *ee;   bsum:= bsum -     220 * sin(M-2*M1-F) *ee;   bsum:= bsum -     220 * sin(2*DD+M-M1-F) *ee;   bsum:= bsum -     185 * sin(DD+M1+F);   bsum:= bsum +     181 * sin(2*DD-M-2*M1-F) *ee;   bsum:= bsum -     177 * sin(M+2*M1+F) *ee;   bsum:= bsum +     176 * sin(4*DD-2*M1-F);   bsum:= bsum +     166 * sin(4*DD-M-M1-F) *ee;   bsum:= bsum -     164 * sin(DD+M1-F);   bsum:= bsum +     132 * sin(4*DD+M1-F);   bsum:= bsum -     119 * sin(DD-M1-F);   bsum:= bsum +     115 * sin(4*DD-M-F) *ee;   bsum:= bsum +     107 * sin(2*DD-2*M+F) *ee^2;   bsum:= bsum-2235*sin(L1)+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(L1-M1)-115*sin(L1+M1);   b := →HMS(bsum/1000000); // β div by 1 milion   // correzioni per la distanza dalla Terra (Δ)   rsum:=0;   rsum:= rsum -20905355 * cos(M1); // coseni   rsum:= rsum - 3699111 * cos(2*DD-M1); // evezione   rsum:= rsum - 2955968 * cos(2*DD); // variazione   rsum:= rsum -  569925 * cos(2*M1); // equazione del centro   rsum:= rsum +   48888 * cos(M) *ee; // equazione annuale   rsum:= rsum -    3149 * cos(2*F); // riduzione all'eclittica   rsum:= rsum +  246158 * cos(2*DD-2*M1);   rsum:= rsum -  152138 * cos(2*DD-M1-M) *ee;   rsum:= rsum -  170733 * cos(2*DD+M1);   rsum:= rsum -  204586 * cos(2*DD-M) *ee;   rsum:= rsum -  129620 * cos(M-M1) *ee;   rsum:= rsum +  108743 * cos(DD);   rsum:= rsum +  104755 * cos(M1+M) *ee;   rsum:= rsum +   10321 * cos(2*DD-2*F);   rsum:= rsum +   79661 * cos(M1-2*F);   rsum:= rsum -   34782 * cos(4*DD-M1);   rsum:= rsum -   23210 * cos(3*M1);   rsum:= rsum -   21636 * cos(4*DD-2*M1);   rsum:= rsum +   24208 * cos(2*DD+M-M1) *ee;   rsum:= rsum +   30824 * cos(2*DD+M) *ee;   rsum:= rsum -    8379 * cos(DD-M1);   rsum:= rsum -   16675 * cos(DD+M) *ee;   rsum:= rsum -   12831 * cos(2*DD-M+M1) *ee;   rsum:= rsum -   10445 * cos(2*DD+2*M1);   rsum:= rsum -   11650 * cos(4*DD);   rsum:= rsum +   14403 * cos(2*DD-3*M1);   rsum:= rsum -    7003 * cos(M-2*M1) *ee;   rsum:= rsum +   10056 * cos(2*DD-M-2*M1) * ee;   rsum:= rsum +    6322 * cos(DD+M1);   rsum:= rsum -    9884 * cos(2*DD-2*M) * ee^2;   rsum:= rsum +    5751 * cos(M+2*M1) *ee;   rsum:= rsum -    4950 * cos(2*DD-2*M-M1) *ee^2;   rsum:= rsum +    4130 * cos(2*DD+M1-2*F);   rsum:= rsum -    3958 * cos(4*DD-M-M1) *ee;   rsum:= rsum +    3258 * cos(3*DD-M1);   rsum:= rsum +    2616 * cos(2*DD+M+M1) *ee;   rsum:= rsum -    1897 * cos(4*DD-M-2*M1) *ee;   rsum:= rsum -    2117 * cos(2*M-M1) * ee^2;   rsum:= rsum +    2354 * cos(2*DD+2*M-M1) *ee^2;   rsum:= rsum -    1423 * cos(4*DD+M1);   rsum:= rsum -    1117 * cos(4*M1);   rsum:= rsum -    1571 * cos(4*DD-M) *ee;   rsum:= rsum -    1739 * cos(DD-2*M1);   rsum:= rsum -    4421 * cos(2*M1-2*F);   rsum:= rsum +    1165 * cos(2*M+M1) *ee^2;   rsum:= rsum +    8752 * cos(2*DD-M1-2*F);   rsum:= rsum/1000;  // div by 1000     temp:= transformEclipticalToEquatorial(l,b);  // transform into AR, decl   α:= temp(1);   δ:= temp(2);   ts:= calcTS();  // calcola tempo siderale   //H:= (ts(2) - α) MOD 24; // Angolo orario   //    // H:= asin(-tan(lat)*tan(δ));   //Hd:= →HMS(15*H); // Angolo orario in gradi   longitude:=-HMS→(long); // - if long -E   Hd:=15*HMS→(apparentSideralTime(dateshortlist))-longitude-α;   IF Hd<0 THEN Hd:=Hd+360 END;   IF Hd>=360 THEN Hd:=Hd-360 END;   Hd:=→HMS(Hd);   H:= Hd/15;   temp:= transformEquatorialToHorizontal(Hd, lat, δ);     //temp:= azimuthFromLBPhi(l, b, Hd);   az:= temp(1);   h:= temp(2);   // Raggio quadratico medio della Terra: 6373,044737   d:= (385000.56 + rsum); // Δ distanza Terra-Luna in km (=1/sin(P)),rr in 10^-3   //P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse   P:= asin(6378.14/d); // Parallasse   s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno     // culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione (Bouiges)     // sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)     // tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   // transit routine   h0:= 0.7275*P -0°34′; // 'standard' altitude of Moon (depend on parallax)   // h0:= 0°7′30″; // medium value: 0.125°   transito:= transit(α, δ, h, h0);   culm:= transito(1);   sorg:= transito(2);   tram:= transito(3);   // end transit   ph:= abs(l - L); // Phase  Moon-Sun   CASE   // IF ph == 0 THEN fase:="New Moon"; END;   IF ph>280 AND ph<350 THEN fase:="Decrescente"; END;   IF ph>=260 AND ph<=280 THEN fase:="Last Quart"; END;   IF ph>190 AND ph<260 THEN fase:="Gibbosa Cal."; END;   IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;   IF ph>100 AND ph<170 THEN fase:="Gibbosa Cresc."; END;   IF ph>=80 AND ph<=100 THEN fase:="1st Quart"; END;   IF ph>10 AND ph<80 THEN fase:="Crescente"; END;   DEFAULT fase:="New Moon";   END; // case   // età della Luna   d1:= IP(HMS→(ph/15));   h1:= IP(FP(HMS→(ph/15))*60);   IF (h1>24) THEN d1:= d1+1; h1:= h1 MOD 24; END;   m1:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);   eta:= →HMS(d1+h1/60+m1/3600); // Age of the Moon   segno:= zodiaco(l);   ill:= (1+cos(ph))/2; // illuminated fraction   PRINT;   PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L1));   PRINT("M anomalia vera " + trunc_angle(M1));   PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");   PRINT("β Latitudine " + trunc_angle(b));   PRINT("λ Longitudine apparente "+ trunc_angle(lapp));   PRINT("Δψ Nutazione long "+ dep(1));   PRINT("Δε Nutazione obl. " + dep(2));   PRINT("ε inclinazione ecl. " + ε);   PRINT("ω Perigeo medio " + trunc_angle(ω));   PRINT("Ω Nodo ascendente " + trunc_angle(Ω));   PRINT(" ");   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("α Ascensione retta (gradi) " + trunc_angle((α*15) MOD 360));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS (locale) " + trunc_angle(ts(2)));   PRINT("TS (apparente) " + trunc_angle(ts(3)));   PRINT("θ0 Mean ST GMT " + trunc_angle(ts(5)));   PRINT("Julian day " + JD);   PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));   PRINT("P Parallasse " + trunc_angle(P));   PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");   PRINT("Semidiametro "+  trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,3) + " km");   PRINT(" ");   PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");   PRINT(" ");   PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");   PRINT("Età " + eta);   PRINT("Parte illuminata " + ill);   RETURN [[L1, M1],[l, b],[lapp,P], [ω,Ω],[α,δ], [az, h],      [ROUND(HMS→(d),6),s], [sorg, tram], [culm, ts(3)],[ph, eta] ]; END; planetCalc(nameplanet) BEGIN   LOCAL ω,Ω, l, v, r, x, y, z;   LOCAL Ls, Ms, b, xs, ys, b_ec;   LOCAL x1, y1, z1, R, l_ec, H, Hd;   LOCAL ec, in, a, δ, α, ts, az, h;   LOCAL DA, SAD, tram, sorg, culm, dayY;   LOCAL as, ac, zc, dz, θ, mov, segno, ε;   // calc of Sun (as basis)   xs:= sunlist(9); ys:= sunlist(10);   Ls:= sunlist(1); // Longitudine media   Ms:= sunlist(4); // Anomalia   // end Sun  EVAL(EXPR(nameplanet+"()")); // call planet-name function   L:= planet_const(1); // Longitudine media   ω:= planet_const(2); // perielio   Ω:= planet_const(3); // nodo   M:= planet_const(4); // anomalia   v:= planet_const(5);    ec:= planet_const(6); // eccentricità   in:= planet_const(7); // inclinazione   a:= planet_const(8);  // semiasse maggiore   θ:= planet_const(9);  // θs to calc if retrogade   b_ec:= asin(sin(in)*sin(v+ω-Ω)); // argomento di latitudine del perielio   l_ec:= acos(cos(v+ω-Ω)/cos(b_ec));  // long from eclittic   IF b_ec<0 THEN l_ec:= 360 - l_ec; END; // argomento secondario   l_ec:= →HMS((l_ec + Ω)) MOD 360; // aggiunge Nodo e riporta a 2PI   r:= a*(1-ec^2)/(1+ec*cos(v)); // 9.524899/(1+ 0.0559*cos(v)); raggio vettore   x1:= r*cos(b_ec)*cos(l_ec); y1:= r*cos(b_ec)*sin(l_ec); z1:=r*sin(b_ec); // cartesian   x:= x1+xs; y:= y1+ys; z:=z1; // cohordinates cartesian->spheric   R:= sqrt(x^2+y^2+z^2);    b:=asin(z/R); // latitudine  β   IF b> 360 THEN b:= b MOD 360; END;   // l:=→HMS(atan(y/x)); // cart->spheric longtitudine λ   l:= atan2(y,x);   IF l<0 THEN l:= l+360; END; IF l>360 THEN l:=l-360; END;   dayY:= DDAYS(IP(thisday)+1/100+1/10000, thisday); // calc day number   // δ:= →HMS((ε*sin((360*(284+dayY))/365)) MOD 360); // declinazione, formula di Cooper   δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(l))) ); // declinazione (coord locali)   α:= (acos(cos(b)*cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta   IF l>180 THEN α:= α +180; END;   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)   as:= cos(δ)*sin(Hd)/sin(dz); // calc az   ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   DA:= →HMS(asin(tan(lat)*tan(δ))); // Differenza ascensionale   SAD:= →HMS(90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno   culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) MOD 24 ; // Culminazione   sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   IF abs(l-220.4)>θ THEN mov:="diretto"; ELSE mov:="retrogrado"; END;   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi di " + nameplanet + " a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L));   PRINT("ω Perielio " + trunc_angle(ω));   PRINT("Ω Nodo " + trunc_angle(Ω));   PRINT("M Anomalia vera " + trunc_angle(M));   PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));   PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(v));   PRINT("r Raggio vettore " + r);   PRINT("Cartes. elioc.: x' " + ROUND(x1,3) + " y' " + ROUND(y1,3) + " z' " + ROUND(z1,3));   PRINT("Cartes. geoc.: x " + ROUND(x,3) + " y " + ROUND(y,3) + " z " + ROUND(z,3));   PRINT("Spher.: λ Longitudine " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");   PRINT("Spher.: β Latitudine " + trunc_angle(b));   PRINT("Spher.: R Distanza " + ROUND(R,3));   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS (locale) " + trunc_angle(ts(2)));   PRINT("Diff. ascensionale " + trunc_angle(DA));   PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));   PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");   PRINT("Movimento " + mov);   RETURN [[L,ω,Ω ],[M,l_ec,b_ec], [v,r,0], [x1,y1,z1],[x,y,z],[l,b,R], [α,δ,H], [az,h,ts(1)],[sorg, tram, culm]]; END; Mercurius() BEGIN   // Costanti per Mercurio   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.205615; // eccentricità   in:= 7.003; // inclinazione   a:= 0.387098; // semiasse maggiore   θ:= 35.6;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   ω:= →HMS((75.913 + 0.00004253*N) MOD 360); // perielio   Ω:= →HMS((47.157 + 0.00003244*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 23.4372*sin(M)+2.9810*sin(2*M)+0.5396*sin(3*M)+0.1098*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Venus() BEGIN   // Costanti per Venere   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.006816; // eccentricità   in:= 3.393636; // inclinazione   a:= 0.72333; // semiasse maggiore   θ:= 13.0;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((206.758 + 1.60216873*N) MOD 360); // Longitudine media   ω:= →HMS((130.154 + 0.00003757*N) MOD 360); // perielio   Ω:= →HMS((75.797 + 0.00002503*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 0.7811*sin(M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Mars() BEGIN   // Costanti per Marte   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.093309; // eccentricità   in:= 1.850303; // inclinazione   a:= 1.523678; // semiasse maggiore   θ:= 16.8;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((124.765 + 0.52407108*N) MOD 360); // Longitudine media   ω:= →HMS((334.251 + 0.00005038*N) MOD 360); // perielio   Ω:= →HMS((48.794 + 0.00002127*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 10.608*sin(M)+0.6216*sin(2*M)+0.0504*sin(3*M)+0.0047*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Jupiter() BEGIN   // Costanti per Giove   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.048335; // eccentricità   in:= 1.308736; // inclinazione   a:= 5.202561; // semiasse maggiore   θ:= 54.43;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((268.350 + 0.08312942*N) MOD 360); // Longitudine media   ω:= →HMS((12.737 + 0.00004408*N) MOD 360); // perielio   Ω:= →HMS((99.453 + 0.00002767*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5372*sin(M)+0.1662*sin(2*M)+0.007*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Saturn() BEGIN   // Costanti per Saturno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.055892; // eccentricità   in:= 2.492519; // inclinazione   a:= 9.554747; // semiasse maggiore   θ:= 65.53;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((278.774 + 0.03349788*N) MOD 360); // Longitudine media   ω:= →HMS((91.117 + 0.00005362*N) MOD 360); // perielio   Ω:= →HMS((112.79 + 0.00002391*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 6.4023*sin(M)+0.2235*sin(2*M)+0.011*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Uranus() BEGIN   // Costanti per Urano   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.046344; // eccentricità   in:= 0.772464; // inclinazione   a:= 19.21814; // semiasse maggiore   θ:= 73.92;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((248.487 + 0.01176902*N) MOD 360); // Longitudine media   ω:= →HMS((171.563 + 0.00004064*N) MOD 360); // perielio   Ω:= →HMS((73.482 + 0.00001365*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5.3092*sin(M)+0.1537*sin(2*M)+0.0062*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Neptunus() BEGIN   // Costanti per Nettuno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.008997; // eccentricità   in:= 1.779242; // inclinazione   a:= 30.10957; // semiasse maggiore   θ:= 77.63;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((86.652 + 0.00602015*N) MOD 360); // Longitudine media   ω:= →HMS((46.742 + 0.000039*N) MOD 360); // perielio   Ω:= →HMS((130.693 + 0.00003009*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 1.031*sin(M)+0.0058*sin(2*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Pluto() BEGIN   // Costanti per Plutone   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.250236; // eccentricità   in:= 17.17047; // inclinazione   a:= 39.438712; // semiasse maggiore   θ:= 79.41;  // θs to calc if retrogade   calcN();   // costanti planetarie   // Long, perielio, nodo, non precisi, considerati costanti   // L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   M:= →HMS((230.67 + 0.00397943*N) MOD 360); // Anomalia vera (calcolata direttamente)   Ω:= →HMS((109.06 + 0.00003823*N) MOD 360); // nodo   ω:= →HMS((114.27 + Ω) MOD 360); // perielio   L:= →HMS((M + ω) MOD 360); // L (non sicuro)      v:= M + 28.45*sin(M)+4.38*sin(2*M)+0.97*sin(3*M)+0.24*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
06-27-2015, 05:51 PM (This post was last modified: 06-28-2015 02:45 PM by salvomic.)
Post: #13
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: Astronomy...
new version, better handling of Sun and Moon RA, dec and Az,h...
Now should work almost well also the routine for transit, rising and setting of Sun, Moon, planets.

When the program will be complete, I'll put the actual version in the Software Library of the Forum.

Salvo

Code:
 data(); calcN(); calcTS(); meanSideralTime(); apparentSideralTime(); julianDay(); julianEphemerisDay(); julianDayAtZero(); makeDateList(); makeDateListShort(); deltaEpsilonPsi(); epsilon(); transformEclipticalToEquatorial(); transformEquatorialToHorizontal(); transformEclipticalToHorizontal(); tempo(); precession(); parametri(); transit(); sunCalc(); planetCalc(); planets(); sun(); sunAh(); moon(); Mercurius(); Venus(); Mars(); Jupiter(); Saturn(); Uranus(); Neptunus(); Pluto(); ascendent(); zodiaco(); trunc_angle(); simp_angle(); atan2(); interpolation(); thisday:= 2000.0101; lstd:=0; gmt:= 12; long:= 0; lat:=  0; m:= 1; dst:= 0; utc:=0; deltaT:= 68; // 2015 datetimelist:= MAKELIST(X,X,1,6); dateshortlist:= MAKELIST(X,X,1,3); dateshortlistgmt:= MAKELIST(X,X,1,3); sunlist:= MAKELIST(X,X,1,25); planet_const:= MAKELIST(X,X,1,10); // ε:= epsilon();  // approx 23.45 (23°26'21") EXPORT Effemeridi() // from "Astrolabio ed Effemeridi", by Salvo Micciché, salvomic, 2015 // Credits: Serge Bouiges, Calcule astronomique pour amateurs // Jean Meeus, Astronomical Algorithms // Thanks Didier Lachieze, Marcel Pelletier, Eried BEGIN   LOCAL ch, ε;   HAngle:=1;   data(); // input data and calc Sun   ε:=epsilon();   CHOOSE(ch, "Effemeridi", "Sun", "Moon", "Planets", "Calc TS JD N", "ε, Δψ, Δε...", "Ascendent", "Quit");   CASE   IF ch==1 THEN sunCalc(); sun(); END;   IF ch==2 THEN moon(); END;   IF ch==3 THEN planets(); END;   IF ch==4 THEN tempo(); END;   IF ch==5 THEN parametri(); END;   IF ch==6 THEN ascendent(); END;   IF ch==7 THEN HAngle:=0; RETURN; END;   DEFAULT sunCalc(); sun();   END; // case END; data() BEGIN   LOCAL dd, mm, yy, hh, mn, ss, loct;   LOCAL t, JD, JDE, ε, dep, eps;   yy:= IP(Date);   mm:= IP(FP(Date)*100);   dd:= FP(Date*100)*100;   hh:= IP(HMS→(Time));   mn:= IP(FP(HMS→(Time))*60);   ss:= IP(FP(FP(FP(HMS→(Time))*60))*60);   loct:= →HMS(hh+mn/60+ss/3600);   INPUT({ {m,[0],{15,15,1}}, {D,[0],{50,15,1}}, {Y,[0],{80,15,1}},   {lstd,[0],{65,30,2}},   {dst,0,{40,15,3}}, {utc,[0],{80,15,3}},    {long,[0],{20,25,4}},    {lat,[0],{70,25,4}},   {deltaT, [0], {20,25,5}} },   "Data: Use Shift+9 for H°M′S″",   {"Month:","Day :","Year :","Local Time (24):", "DST", "UTC",    "Long (-E):","Lat (+N)", "delta T"},   {"Enter month", "Enter day", "Enter year (1901-2099)",    "Enter local time (decimal or H°M′S″)", "Daylight Save Time?", "UTC time zone",   "Longitude", "Latitude", "Delta T (sec.)"},    {1,1,1901,0, 0,0, 0, 0, 0},{mm,dd,yy,loct, 0, 1, -14°43′41″, 36°55′15″, 68});   // Ragusa (-14°43′41″, 36°55′15″) -    // Marina di Ragusa (-14°33′15″, 36°47′06″) - Scicli (-14°42'22″, 36°47' 22″)   gmt:= lstd-utc-(1*dst); // UTC (DST, daylight save time)   thisday:= Y+m/100+D/10000;   datetimelist:= {Y, m, D, IP(HMS→(lstd)), IP(FP(HMS→(lstd))*60), ROUND(FP(HMS→(lstd)*60)*60,3) };   dateshortlist:= {Y, m, D+HMS→(lstd)/24};   dateshortlistgmt:= {Y, m, D+HMS→(gmt)/24};   // Greenwich Mean Time   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   T:=(JD-2451545)/36525;   U:= T/100;    ε:=epsilon(); // (formula IAU J2000.0)   sunCalc(); END; epsilon() BEGIN LOCAL ε, ε0, dep, eps, JD, JDE, U;   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   T:=(JD-2451545)/36525;   U:= T/100;   // ε:=23.4392911111-0.013004166667*t-0.000001638889*t^2+0.000000503611*t^3; (formula IAU J2000.0)   ε0:= 23°26′21.448″ - 1°18′00.93″*U-1.55*U^2+1999.25*U^3-51.38*U^4-249.67*U^5;   ε0:= ε0-39.05*U^6+7.12*U^7+27.87*U^8+5.79*U^9+2.45*U^10; // Laskar, Meeus   ε0:= →HMS(ε0);   dep:= deltaEpsilonPsi();   eps:= dep(2);   ε:= ε0+eps; // true obliquity ε   RETURN ε; END; deltaEpsilonPsi() BEGIN   // Nutazione Δψ (longitudine), Δε (obliquità)   LOCAL Lmoon, Lsun, deltaPsi, deltaEpsilon, Ωmoon, JD, JDE;   LOCAL DD;   JD:= julianDayAtZero(Y,m,D);   JDE:= julianEphemerisDay(dateshortlist);   T:= (JDE-2451545)/36525; // use JDE   DD:=297.85036+445267.111480*T-0.0019142*T^2+(T^3)/189474;  // elongazione   DD:=FP(DD/360)*360;   Ωmoon:= 125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000; // Nodo Luna (Moon)   Lsun:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Sole   Lmoon:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);    //  Long media Luna   deltaPsi:= (-0°0′17.20″)*sin(Ωmoon)-(0°0′1.32″)*sin(2*Lsun)+(-0°0′0.23″)*sin(2*Lmoon)+(0°0′0.21″)*sin(2*Ωmoon);   deltaEpsilon:= (0°0′9.2″)*cos(Ωmoon)+(0°0′0.57″)*cos(2*Lsun)+(0°0′0.10″)*cos(2*Lmoon)-(0°0′0.09″)*cos(2*Ωmoon);   RETURN({deltaPsi, deltaEpsilon});    END; precession(alfa, delta) BEGIN   // precessione dati asce retta e declinazione   LOCAL m,n,n1, deltaAlfa, deltaDelta, Tsec;   Tsec:= IP((Y-2000)/100)+1; // T century from 2000.0   m:= 0°0′3.07496″+0°0′0.00186″*Tsec;   n:= 0°0′1.33621″+0°0′0.00057″*Tsec; // sec   n1:= 0°0′20.0431″+0°0′0.0085″*Tsec; // ° ' "   deltaAlfa:= m+n*sin(alfa)*tan(delta); // Δα   deltaDelta:= n1*cos(alfa); // Δδ   RETURN({deltaAlfa, deltaDelta}); END; planets() BEGIN   LOCAL ch, nameplanet;   INPUT({{ch, nameplanet:={"Mercurius","Venus","Mars", "Jupiter",    "Saturn", "Uranus", "Neptunus", "Pluto"}, {20,30,1}}},   "Choose planet",    "planet: ", "Choose the planet to calculate an press OK", 0,5 );   planetCalc(nameplanet(ch)); END; transformEclipticalToEquatorial(lambda, beta) BEGIN   // trasnform λ and β (long, lat) into α and δ (Asc retta, decl)   LOCAL ε,Epsilon,Lambda, Beta, alpha, delta;   Lambda:=HMS→(lambda);   Beta:=HMS→(beta);   Epsilon:= epsilon();   ε:= epsilon();   alpha:=atan2(sin(Lambda)*cos(Epsilon)-tan(Beta)*sin(Epsilon),cos(Lambda));   IF alpha>=360 THEN alpha:=alpha-360 END;   IF alpha<0 THEN alpha:=alpha+360 END;   delta:=asin(sin(Beta)*cos(Epsilon)+cos(Beta)*sin(Epsilon)*sin(Lambda));   alpha:=→HMS(alpha/15);   delta:=→HMS(delta);   RETURN ({alpha,delta}); END; transformEquatorialToHorizontal(Lha, Phi, Delta) BEGIN   LOCAL azimuth,altitude,lha,phi,delta;   lha:=HMS→(Lha);   phi:=HMS→(Phi);   // latitude   delta:=HMS→(Delta);   azimuth:=atan2(sin(lha),cos(lha)*sin(phi)-tan(delta)*cos(phi));   altitude:=asin(sin(phi)*sin(delta)+cos(phi)*cos(delta)*cos(lha));   azimuth:=→HMS(azimuth);   altitude:=→HMS(altitude);   RETURN({azimuth,altitude}); END; transformEclipticalToHorizontal(Lambda, Beta, Lha) BEGIN   LOCAL ε;   // transform λ and β (long, lat) into azimuth and height   LOCAL lambda, beta, epsilon,phi, lha;   LOCAL alpha,delta,azimuth,altitude;   lambda:=HMS→(Lambda); // longitudine astro   beta:=HMS→(Beta); // latitudine astro   epsilon:= epsilon();   ε:= epsilon();   lha:=HMS→(Lha);  // angolo orario   phi:=HMS→(lat); // lat geo   alpha:=atan2(SIN(lambda)*COS(epsilon)-TAN(beta)*SIN(epsilon),COS(lambda));   delta:=asin(sin(beta)*cos(epsilon)+cos(beta)*sin(epsilon)*sin(lambda));   azimuth:=atan2(SIN(lha),COS(lha)*SIN(phi)-TAN(delta)*COS(phi));   altitude:=ASIN(SIN(phi)*SIN(delta)+COS(phi)*COS(delta)*COS(lha));   azimuth:=→HMS(azimuth);   altitude:=→HMS(altitude);   RETURN ({azimuth,altitude}); END; calcN() BEGIN   N:= DDAYS(1901.0101, thisday)+1+gmt/24; // N to GMT   RETURN N; END; calcTS() BEGIN   LOCAL T, TS, TSd, N0, TSL, TSapp, ε;   LOCAL θ0, θ0GMT, JD, JDE, JD0, dep, psi, eps;   ε:= epsilon;   calcN();   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   JD0:= julianDayAtZero(Y,m,D);   dep:= deltaEpsilonPsi();   psi:=dep(1);   eps:=dep(2);   T:= (JD-2451545.0)/36525;   θ0:= meanSideralTime(dateshortlist); // MST DEG at 0h   θ0GMT:= meanSideralTime({Y,m,D}); // MST DEG at 0h at Greenwitch   // TS:=  ((23750.3 + 236.555362*N) / 3600) MOD 24; // Sideral Time (GMT) in seconds   // TS:= →HMS(TS); // in HMS   TS:= θ0GMT;   // TSL:= (TS + gmt - 4*(HMS→(long))/60) MOD 24; // local ST   // TSL:= →HMS(TSL); // in HMS   TSapp:= apparentSideralTime(dateshortlist); // apparent ST (ST at Greenwitch at local our)   TSL:= apparentSideralTime(dateshortlist) - 1*dst;   TSd:= →HMS(TS*15) MOD 360; // ST at 0 GMT in degrees   RETURN {TS, TSL, TSapp, θ0, T}; END; meanSideralTime(lista) // lista like dateshortlist BEGIN   LOCAL T,MST, JD, ε;   ε:= epsilon;   JD:= julianDay(lista);   T:=(julianDay(lista)-2451545.0)/36525;   MST:=280.46061837+360.98564736629*(JD-2451545.0)+0.000387933*T^2-(T^3)/38710000;   MST:=FP(MST/360)*360; // θ0   IF MST<0 THEN MST:=MST+360 END;   MST:=→HMS(MST/15);   RETURN MST; END; apparentSideralTime(lista) // lista like dateshortlist BEGIN   LOCAL MST,correction, deltaPsi, AST, dep, ε;   ε:= epsilon;   MST:= meanSideralTime(lista);   dep:= deltaEpsilonPsi();   deltaPsi:= dep(1);   correction:=(HMS→(deltaPsi)*3600*COS(HMS→(ε)))/15;   AST:=→HMS(HMS→(MST)+(correction/3600));   RETURN AST; END; julianDay(lista) // lista like dateshortlist BEGIN   LOCAL y,m,d,a,b, JD;   y:=lista(1);   m:=lista(2);   d:=lista(3);   IF (m=1 OR m=2) THEN   y:=y-1;   m:=m+12   END;   IF y*100+m+d>=158225 THEN   a:=IP(ABS(y/100));   b:=2-a+IP(a/4);   ELSE   b:=0   END;   JD:=IP(365.25*(y+4716))+IP(30.6001*(m+1))+d+b-1524.5;   RETURN JD; END; julianEphemerisDay(lista) BEGIN   LOCAL JDE;   JDE:=julianDay(makeDateListShort(makeDateList(lista)+{0,0,0,0,0,deltaT}));   RETURN JDE; END; julianDayAtZero(Y,m,D) BEGIN   // JD at 0h GMT   LOCAL y,mm,d,a,b, JD;   y:=Y;   mm:=m;   d:=dateshortlist(3);   IF (mm=1 OR mm=2) THEN   y:=y-1;   mm:=mm+12   END;   IF y*100+mm+d>=158225 THEN   a:=IP(ABS(y/100));   b:=2-a+IP(a/4);   ELSE   b:=0   END;   JD:=IP(365.25*(y+4716))+IP(30.6001*(mm+1))+d+b-1524.5;   RETURN JD; END; makeDateList(lista) // lista {Y m D.d} (dateshortlist) BEGIN LOCAL y,o,m,d,h,s,t; y:=lista(1); o:=lista(2); d:=IP(lista(3)); t:=FP(lista(3)); h:=IP(t*24); t:=FP(t*24); m:=IP(t*60); t:=FP(t*60); s:=IP(t*60); lista:={y,o,d,h,m,s}; RETURN lista; END; makeDateListShort(lista) BEGIN // lista {Y m D h mn, s} LOCAL y,m,d, listacorta; y:=lista(1); m:=lista(2); d:=lista(3)+lista(4)/24+lista(5)/1440+lista(6)/86400; listacorta:={y,m,d}; RETURN listacorta; END; tempo() BEGIN   LOCAL n, JD, JDE, TSid, dep, psi, eps;   LOCAL ε;   ε:= epsilon();   n:= calcN();   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   TSid:= calcTS();   dep:= deltaEpsilonPsi();   psi:=dep(1);   eps:=dep(2);   PRINT;   PRINT("Date " + m + " " + D + " " + Y);   PRINT("Julian Day " + JD);   PRINT("Julian Day Effem. " + JDE);   PRINT("N (from 1901jan1) " + n);   PRINT("T (from JD) " + TSid(5));   PRINT(" ");   PRINT("Tempo Siderale 0h " + →HMS(TSid(1)));   PRINT("Tempo Siderale locale " + →HMS(TSid(2)));   PRINT("Tempo Siderale apparente " + →HMS(TSid(3)));   PRINT("θ0 Mean ST GMT " + →HMS(TSid(4)));   PRINT("Tempo Siderale (DEG) " + trunc_angle(TSid(1)*15));   PRINT("θ0 Mean ST GMT (DEG) " + trunc_angle(TSid(4)*15));   PRINT(" ");   PRINT("Δψ Nutazione (longit.) "+ psi);   PRINT("Δε Nutazione (obl.) "+ eps);   PRINT("ε Obliquità Ecl. "+ ε);   PRINT(" ");   PRINT("Time " + lstd + " GMT " + gmt);   RETURN({ROUND(JD,5), ROUND(JDE,5), trunc_angle(TSid(1)), trunc_angle(TSid(2)),      trunc_angle(TSid(3)), trunc_angle(TSid(4)), trunc_angle(TSid(5)) }); END; transit(alfa, delta, h0) BEGIN // α, δ (AR, dec), height, h0 = standard height of body   LOCAL temp, H0, tsideral, Hor, m0, m1, m2, h;   LOCAL α, δ, θ0,ts, culm, sorg, tram, deltaM;   α:= alfa; // AR   δ:= delta;   temp:= calcTS();   // θ0:= apparentSideralTime(makeDateListShort({Y,m,D,0,0,0})) *15; // TS in DEG   θ0 := temp(1) * 15;   temp:= (sin(h0)-sin(lat)*sin(δ))/(cos(lat)*cos(δ));   IF temp<-1 OR temp>1 THEN culm:=0; END; // above horizon?   H0:= acos(temp);   //IF H0<-180 THEN H0:=H0+180; END; IF H0>180 THEN H0:= H0-180; END; // it should be from -180 to 180     m0:= HMS→(α*15 + long - θ0) / 360;   IF m0<0 THEN m0:= m0 +1; END; IF m0>1 THEN m0:=m0-1; END;   tsideral:= simp_angle(θ0 + 360.985647*m0); // TS transit DEG   Hor:= HMS→(tsideral - long -α*15);  // Hour angle, local ST vs apparent ST   Hor:=simp_angle(→HMS(Hor));   deltaM:= -HMS→(Hor/360);   culm:= (m0 + deltaM)*24;  // transit   m1:= m0-H0/360;   IF m1<0 THEN m1:= m1 +1; END; IF m1>1 THEN m1:=m1-1; END;   tsideral:= simp_angle(θ0 + 360.985647*m1); // TS rising   Hor:= HMS→(tsideral - long - α*15);   Hor:=simp_angle(→HMS(Hor));   h:=asin(sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hor));   deltaM:= HMS→((h-h0)/(360*cos(δ)*cos(lat)*sin(Hor)));   sorg:= ((m1 + deltaM)*24); // rising   m2:= m0+H0/360;   IF m2<0 THEN m2:= m2 +1; END; IF m2>1 THEN m2:=m2-1; END;   tsideral:= simp_angle(θ0 + 360.985647*m2); // TS setting   Hor:= HMS→(tsideral - long - α*15);    Hor:=simp_angle(→HMS(Hor));   h:=asin(sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hor));   deltaM:= (HMS→(h)-HMS→(h0))/(360*cos(δ)*cos(lat)*sin(Hor));   tram:= ((m2 + deltaM) * 24); // setting RETURN ({culm, sorg, tram}); END; parametri() BEGIN   LOCAL dep, prec, α, δ, ε;   dep:= deltaEpsilonPsi();   α:= sunlist(5); δ:= sunlist(6);   prec:= precession(α, δ);   ε:= epsilon();   // Nutazione Δψ (longitudine), Δε (obliquità), ε (inclinazione eclittica)...   PRINT;   PRINT("ε Obliquità eclittica "+ ε);   PRINT("Δψ Nutazione longitudine "+ dep(1));   PRINT("Δε Nutazione longitudine "+ dep(2));   PRINT(" ");   PRINT("Δα Precessione AR "+ prec(1));   PRINT("Δδ Precessione dec "+ prec(2));   RETURN({ε, dep(1), dep(2), prec(1), prec(2)}); // ε, Δψ, Δε, Δα, Δδ END; ascendent() BEGIN   LOCAL ASC, ts, segno, ε;   ε:= epsilon;   ts:= calcTS();   ASC:= atan2(-cos(ts(2)), sin(ts(2))*cos(ε)+tan(lat)*sin(ε));   ASC:= simp_angle(ASC);   segno:= zodiaco(ASC); RETURN {trunc_angle(→HMS(ASC)), trunc_angle(→HMS(ASC) MOD 30), segno}; END; zodiaco(astro) BEGIN   LOCAL segno;   CASE  // Segni e case   IF astro>=0   AND astro <30  THEN segno:="Ariete"; END;   IF astro>=30  AND astro <60  THEN segno:="Toro"; END;   IF astro>=60  AND astro <90  THEN segno:="Gemelli"; END;   IF astro>=90  AND astro <120 THEN segno:="Cancro"; END;   IF astro>=120 AND astro <150 THEN segno:="Leone"; END;   IF astro>=150 AND astro <180 THEN segno:="Vergine"; END;   IF astro>=180 AND astro <210 THEN segno:="Bilancia"; END;   IF astro>=210 AND astro <240 THEN segno:="Scorpione"; END;   IF astro>=240 AND astro <270 THEN segno:="Sagittario"; END;   IF astro>=270 AND astro <300 THEN segno:="Capricorno"; END;   IF astro>=300 AND astro <330 THEN segno:="Acquario"; END;   IF astro>=330 AND astro <360 THEN segno:="Pesci"; END;   END; // case   RETURN segno; END; atan2(y,x) BEGIN // atan2 returns coordinates in correct angle for all four quadrants (DEG)     CASE     IF x==0 AND y==0 then return "undefined"; end;         // x=0, y=0       IF x>0 then return atan(y/x); end;                     // x>0     IF x<0 AND y>=0 then return atan(y/x)+180; end;         // x<0, y>=0     IF x==0 AND y>0 then return 90; end;                 // x=0, y>0     IF x<0 AND y<=0 then return atan(y/x)-180; end;         // x<0, y<0     IF x==0 AND y<0 then return -90; end;                // x=0, y<0    END;   RETURN;  END; trunc_angle(ang) // trunc decimal from a DEG form angle BEGIN →HMS((IP(HMS→(ang)*3600)/3600)); END; simp_angle(ang) BEGIN LOCAL angle; angle:=360*FP(ang/360); IF angle<0 THEN angle:=angle+360 END; RETURN angle; END; kepler(ec, M)  // needs eccentricity and true anomaly M BEGIN   LOCAL M1, j, u;   HAngle:=0; // RAD   M1:= HMS→(M)*PI/180; // anomalia RAD   F:= SIGN(M1); M1:= ABS(M1)/(2*PI);   M1:=(M1-IP(M1))*2*PI*F;   IF M1<0 THEN M1:= M1+2*PI; END;   F:= 1;    IF M1>PI THEN F:= -1; END;   IF M1>PI THEN M1:=2*PI-M; END;   u:= PI/2; D:=PI/4;   FOR j FROM 1 TO 33 DO // 33 iterations (Meeus, Roger Sinnot)   M2:= u-ec*sin(u); u:= u+D*SIGN(M1-M2); D:= D/2;   END; // for   u:= u * F; // anomalia eccentrica    HAngle:=1;   u:= u*180/PI;   RETURN u; END; interpolation(lista,n) BEGIN LOCAL a,b,c,R,dif1,dif2; dif1:=ΔLIST(lista); dif2:=ΔLIST(dif1); a:=dif1(1); b:=dif1(2); c:=dif2(1); R:=lista(2)+(n/2)*(a+b+n*c); RETURN R; END; sunCalc() BEGIN   LOCAL ω, l, v, r, x, y, dz, ε;   LOCAL α, δ, dayY, H, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, h0, H0, temp;   LOCAL ec, a, u, M1, M2, j, t, tau;   LOCAL JD, JDE, c, lapp, Ω, θ0, transito;   LOCAL dep, deltapsi, longitude;   ε:= epsilon;   N:= calcN();   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   a:= 1.00000023; // semiasse maggiore   dep:= deltaEpsilonPsi();   deltapsi:= dep(1);   t:= N/36525;   T:=   (JD-2451545.0)/36525;    tau:= (JDE-2451545.0)/365250; // millenni  - JDE   ec:= 0.016708617-0.000042037*T-0.0000001236*T^2; // eccentricità Meeus   // ec:= ec*180/PI; // DEG   // L:= →HMS(278.9654+0.985647342*N + 0.0003*t^2) MOD 360; // long media Bouiges   // L:= simp_angle(280.46645+36000.76983*T+0.0003032*T^2); // Long media Meeus   L:= simp_angle(280.4664567+360007.6982779*tau+0.03032028*tau^2+(tau^3)/49931-(tau^4)/15299-(tau^5)/1988000);   ω:= →HMS(281.238 + 0.000047067*N + 0.00045*t^2) MOD 360; // long del perielio   // M:= (L - ω); // anomalia vera   //IF M<0 THEN M:= (360+M); END; // potrebbe essere negativa   Ω:= simp_angle(125.04-1934.136*T); // Nutazione e aberrazione (nodo)   M:= simp_angle(357.52910+35999.05030*T+0.0001559*T^2-0.00000048*T^3); // anom vera Meeus   u:=atan2(sin(M),(cos(M)-ec)); // only for small eccentricity   c:= (1.914600-0.004817*T-0.000016*T^2)*sin(M);  // equation of center   c:= c+(0.019993-0.000101*T)*sin(2*M)+(0.000290)*sin(3*M);   // v:= M + 1.91934 * sin(M)+0.02*sin(2*M); // calcolo intermedio (Bouiges)   // v:= acos((cos(u)-ec)/(1-ec*cos(u))); // kepler   v:= M + c; // true anomaly, Meeus   // l:= (v + ω) MOD 360; // Longitudine Sole   l:= simp_angle(L+c); // Ө Longitudine Sole, Meeus   lt:= l + 180; // Longitudine eliocentrica della Terra   lapp:= simp_angle(l - 0.00569-0.00478*sin(Ω)); // Longitudine apparente (ref true equinox)   // r:= ((0.999721)/(1+0.01675*cos(v))); // Distanza Sole-Terra (UA)   r:= (1.000001018*(1-ec^2))/(1+ec*cos(v)); // Distanza Sole-Terra (UA) Meeus   x:= r * cos(l); y:= r * sin(l);   δ:= →HMS((asin(sin(ε)*sin(l))) );   α:= →HMS(atan2(cos(ε)*sin(l),cos(l))); // Ascensione retta, Meeus (b never > 1.2", l=Ө)   α:= →HMS(α/15) MOD 24; // α in hms   ts:= calcTS();  // calcola tempo siderale  // H:= (ts(2) - α) MOD 24; // Angolo orario  // Hd:= →HMS(15*H); // Angolo orario in gradi   Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST   Hd:=simp_angle(→HMS(Hd));   H:= Hd/15; // in HMS      temp:= sunAh(dateshortlist);   az:= temp(1);  // Azimuth   h:=  temp(2);  // height   // E:= →HMS((460*sin(M)-592*sin(2*L))/60); // equazione del tempo in min (460s, 592s)   E:= 4*( L - 0.0057183 - α*15 + deltapsi*cos(ε) );  // α in DEG   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   //   SAD:= →HMS(acos(-tan(lat)*tan(δ))); // Semiarco Diurno      // culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione, transito (Bouiges)    // sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)    // tram:= (culm + SAD/15) MOD 24; // Tramonto (Sunset)   // transit routine   h0:= -0°50′; // 'standard' altitude of Sun   transito:= transit(α, δ, h0);   culm:= transito(1);   sorg:= transito(2);   tram:= transito(3);   // end transit   sunlist:= {L,ω,l, M, α,δ, az, h, x,y,              sorg, tram, culm, r, H, Hd, zc, dz, ts(1), DA,              SAD, E, lt, v, c, Ω, lapp, ec};   RETURN sunlist; END; sun() BEGIN   LOCAL ω, l, v, r, x, y, dz, ε;   LOCAL α, δ, dayY, ts, az, z, h;   LOCAL SAD, DA, sorg, tram, culm, Hd, lt;   LOCAL as, ac, zc, segno, JD, JDE, v, c;   LOCAL Ω, lapp, ec;   ε:= epsilon;   L:= sunlist(1); ω:= sunlist(2);    l:= sunlist(3); M:= sunlist(4);   Ω:= sunlist(26); lapp:=sunlist(27);    α:= sunlist(5); δ:= sunlist(6);    az:= sunlist(7); h:= sunlist(8);    x:= sunlist(9); y:= sunlist(10);   sorg:= sunlist(11); tram:= sunlist(12); culm:= sunlist(13);   r:= sunlist(14); H:= sunlist(15); Hd:= sunlist(16);   zc:= sunlist(17); dz:= sunlist(18);   v:= sunlist(24); c:= sunlist(25);    ts:= calcTS(); //   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   DA:= sunlist(20); SAD:= sunlist(21); E:= sunlist(22); lt:= sunlist(23);   ec:= sunlist(28);   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi del Sole a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("at Long " + long + " Lat " + lat);   PRINT("");   PRINT("L longitudine media " + trunc_angle(L));   PRINT("ω long del perielio " + trunc_angle(ω));   PRINT("M anomalia media " + trunc_angle(M));   PRINT("v anomalia vera " + trunc_angle(v));   PRINT("c equazione del centro " + c);   PRINT(" ");   PRINT("Ө Longitudine Sole " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");   PRINT("λ longitudine apparente " + trunc_angle(lapp));   PRINT("Ω aberr, nutaz nodo "+ Ω);   PRINT("r distanza (UA) " + r);   PRINT(" ");   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("α Asc.r. (DEG) " + trunc_angle(α*15));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT(" ");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS locale " + trunc_angle(ts(2)));   PRINT("TS apparente " + ts(3));   PRINT("θ0 Mean ST GMT " + ts(4));   PRINT("Equazione del tempo (min) " + trunc_angle(E));   PRINT("Julian day " + JD);   PRINT("Dif asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));   PRINT("");   PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");   PRINT(" ");   PRINT("Long elioc. Terra "+ trunc_angle(lt));   PRINT("ε obliquità ecl. " + ε);   PRINT("ec eccentricità "+ ec);   // RETURN {L, ω, M, v, l, r, x, y};   RETURN [[L, ω], [l, lapp], [M,v], [Ω, c], [α,δ], [az, h], [x,y],    [sorg, tram], [culm, r], [H,ts(1)], [lt, 0], [ε, ec]]; END; sunAh(lista) // lista like dateshortlist BEGIN   LOCAL lha,alpha,delta,phi, ah, α;   LOCAL temp,longitude, Azimuth, Altitude, ts;   alpha:=15*HMS→(sunlist(5));   α:= sunlist(5); // α in HMS   delta:=sunlist(6);   phi:=lat;   ts:= calcTS();   // longitude:= HMS→(long); // - if long -E   // lha:=(15*HMS→(apparentSideralTime(lista))-longitude-alpha);   lha:= 15*HMS→(ts(2)-α);  // local ST vs apparent ST;   // IF lha<0 THEN lha:=lha+360 END;   // IF lha>=360 THEN lha:=lha-360 END;   lha:=simp_angle(→HMS(lha));   temp:= transformEquatorialToHorizontal(lha, phi, delta);   Azimuth:= HMS→(temp(1))+180;   IF Azimuth>360 THEN Azimuth:=Azimuth-360 END;   Azimuth:= →HMS(Azimuth);   Altitude:= temp(2);   ah:= {Azimuth,Altitude};   RETURN ah; END; moon() BEGIN   LOCAL ω,Ω, L1, M1, l, lsum, bsum, rsum;    LOCAL v, r, x, y, dep, ε;   LOCAL b, P, d, s, temp;   LOCAL α,δ, ts, Hd, zc, dz, h, az;   LOCAL i, as, ac, DA, SAD, DD;   LOCAL culm, sorg, tram, h0, transito;   LOCAL xs, ys, ph, fase, eta;   LOCAL d1, h1, m1, segno, JD, JDE, rr;   LOCAL lapp, ill, ee, a1,a2,a3, longitude;   i:= 5.13333333; // inclinazione 5°8' sull'eclittica   ε:= epsilon();   N:= calcN();   JD:= julianDay(dateshortlist);   JDE:= julianEphemerisDay(dateshortlist);   T:= (JDE-2451545.0)/36525; // use JDE   dep:= deltaEpsilonPsi();   // Sun   L:= simp_angle(280.4665+36000.76983*T+0.0003032*T^2); // Long media Sole Meeus   M:= simp_angle(357.5291092+35999.0502909*T-0.0001536*T^2+T^3/24490000); // mean anomaly Sun      // end Sun   // L1:= (33.231 + 13.17639653*N) MOD 360; // long. media   L1:= simp_angle(218.3164591+481267.88134236*T-0.0013268*T^2+T^3/538841-T^4/65194000);   //  Long media Luna Meeus   // Ω:= (239.882 - 0.052953922*N) MOD 360; // Long nodo ascendente   Ω:= simp_angle(125.04452-1934.136261*T-0.0020708*T^2+(T^3)/450000); // Nodo Luna (Moon)   // M1:= (18.294 + 13.06499245*N) MOD 360; // anomalia vera   M1:= simp_angle(134.9634114+477198.8676313*T+0.0089970*T^2+T^3/69699-T^4/14712000); // anomalia media   ω:= (L1 - M1) MOD 360; // perigeo medio   // DD:= L1 - L; F:= L1-Ω; // for the variations   DD:= simp_angle(297.8502042+445267.1115168*T-0.0016300*T^2+T^3/545868-T^4/113065000); // mean elongation   F:= simp_angle(93.2720993+483202.0175273*T-0.0034029*T^2-T^3/3526000+T^4/863310000); // dist from asc. node   // ee:= 1-0.002516*T-0.0000074*T^2;  // multiply those w/ M *ee   a1:= simp_angle(119.75+131.849*T); // action of Venus   a2:= simp_angle(53.09+479264.290*T); // action of Jupiter   a3:= simp_angle(313.45+481266.484*T);    lsum:=0;  // da sommare a L1   lsum:= lsum + 6288774 * sin(M1); // calc Moon longitude   lsum:= lsum + 1274027 * sin(2*DD-M1); // evezione   lsum:= lsum +  658314 * sin(2*DD); // variazione   lsum:= lsum +  213618 * sin(2*M1); // equazione del centro   lsum:= lsum -  185116 * sin(M) *ee; // equazione annuale   lsum:= lsum -  114332 * sin(2*F); // riduzione all'eclittica   lsum:= lsum +   58793 * sin(2*DD-2*M1);   lsum:= lsum +   57066 * sin(2*DD-M1-M) *ee;   lsum:= lsum +   53322 * sin(2*DD+M1);   lsum:= lsum +   45758 * sin(2*DD-M) *ee;   lsum:= lsum -   40923 * sin(M-M1) *ee;   lsum:= lsum -   34720 * sin(DD);   lsum:= lsum -   30383 * sin(M1+M) *ee;   lsum:= lsum +   15327 * sin(2*DD-2*F);   lsum:= lsum -   12528 * sin(M1+2*F);   lsum:= lsum +   10980 * sin(M1-2*F);   lsum:= lsum +   10675 * sin(4*DD-M1);   lsum:= lsum +   10034 * sin(3*M1);   lsum:= lsum +    8548 * sin(4*DD-2*M1);   lsum:= lsum -    7888 * sin(2*DD+M-M1) *ee;   lsum:= lsum -    6766 * sin(2*DD+M) *ee;   lsum:= lsum -    5163 * sin(DD-M1);   lsum:= lsum +    4987 * sin(DD+M) *ee;   lsum:= lsum +    4036 * sin(2*DD-M+M1) *ee;   lsum:= lsum +    3994 * sin(2*DD+2*M1);   lsum:= lsum +    3861 * sin(4*DD);   lsum:= lsum +    3665 * sin(2*DD-3*M1);   lsum:= lsum -    2689 * sin(M-2*M1) *ee;   lsum:= lsum -    2602 * sin(2*DD-M1+2*F);   lsum:= lsum +    2390 * sin(2*DD-M-2*M1) *ee;   lsum:= lsum -    2348 * sin(DD+M1);   lsum:= lsum +    2236 * sin(2*DD-2*M) *ee^2;   lsum:= lsum -    2120 * sin(M+2*M1) *ee;   lsum:= lsum -    2069 * sin(2*M) *ee^2;   lsum:= lsum +    2048 * sin(2*DD-2*M-M1) *ee^2;   lsum:= lsum -    1773 * sin(2*DD+M1-2*F);   lsum:= lsum -    1595 * sin(2*DD+2*F);   lsum:= lsum +    1215 * sin(4*DD-M-M1) *ee;   lsum:= lsum -    1110 * sin(2*M1+2*F);   lsum:= lsum -     892 * sin(3*DD-M1);   lsum:= lsum -     810 * sin(2*DD+M+M1) *ee;   lsum:= lsum +     759 * sin(4*DD-M-2*M1) *ee;   lsum:= lsum -     713 * sin(2*M-M1) * ee^2;   lsum:= lsum -     700 * sin(2*DD+2*M-M1) *ee^2;   lsum:= lsum +     691 * sin(2*DD+M-2*M1) *ee;   lsum:= lsum +     596 * sin(2*DD-M-2*F) *ee;   lsum:= lsum +     549 * sin(4*DD+M1);   lsum:= lsum +     537 * sin(4*M1);   lsum:= lsum +     520 * sin(4*DD-M) *ee;   lsum:= lsum -     487 * sin(DD-2*M1);   lsum:= lsum -     399 * sin(2*DD+M-2*F) *ee;   lsum:= lsum -     381 * sin(2*M1-2*F);   lsum:= lsum +     351 * sin(DD+M+M1) *ee;   lsum:= lsum -     340 * sin(3*DD-2*M1);   lsum:= lsum +     330 * sin(4*DD-3*M1);   lsum:= lsum +     327 * sin(2*DD-M+2*M1) *ee;   lsum:= lsum -     323 * sin(2*M+M1) *ee^2;   lsum:= lsum +     299 * sin(DD+M-M1) *ee;   lsum:= lsum +     294 * sin(2*DD+3*M1);   lsum:= lsum + 3958*sin(a1)+1962*sin(L1-F)+318*sin(a2);   l:= (simp_angle(L1 + lsum/1000000)); // λ divided by 1 milion   lapp:= →HMS(l + dep(1)); // longitudine apparente (+ nutazione long, Δψ)   // fattori per latitudine   bsum:=        5128122 * sin(F); // termine principale latitudine   bsum:= bsum +  280026 * sin(M1+F);   bsum:= bsum +  277693 * sin(M1-F); // equazione del centro   bsum:= bsum +  173237 * sin(2*DD-F); // grande ineguaglianza   bsum:= bsum +   55413 * sin(2*DD-M1+F); // evezione   bsum:= bsum +   46271 * sin(2*DD-M1-F); // evezione (2)   bsum:= bsum +   32573 * sin(2*DD+F);   bsum:= bsum +   17198 * sin(2*M1+F);   bsum:= bsum +    9266 * sin(2*D+M1-F);   bsum:= bsum +    8822 * sin(2*M1-F);   bsum:= bsum +    8216 * sin(2*DD-M-F) *ee;   bsum:= bsum +    4324 * sin(2*DD-2*M1-F);   bsum:= bsum +    4200 * sin(2*DD+M1+F);   bsum:= bsum -    3359 * sin(2*DD+M-F) *ee;   bsum:= bsum +    2463 * sin(2*DD-M-M1+F) *ee;   bsum:= bsum +    2211 * sin(2*DD-M+F) *ee;   bsum:= bsum +    2065 * sin(2*DD-M-M1-F) *ee;   bsum:= bsum -    1870 * sin(M-M1-F) *ee;   bsum:= bsum +    1828 * sin(4*DD-M1-F);   bsum:= bsum -    1794 * sin(M+F) *ee;   bsum:= bsum -    1749 * sin(3*F);   bsum:= bsum -    1565 * sin(M-M1+F) *ee;   bsum:= bsum -    1491 * sin(DD+F);   bsum:= bsum -    1475 * sin(M+M1+F) *ee;   bsum:= bsum -    1410 * sin(M+M1-F) *ee;   bsum:= bsum -    1344 * sin(M-F) *ee;   bsum:= bsum -    1335 * sin(DD-F);   bsum:= bsum +    1107 * sin(3*M1+F);   bsum:= bsum +    1021 * sin(4*DD-F);   bsum:= bsum +     833 * sin(4*DD-M1+F);   bsum:= bsum +     777 * sin(M1-3*F);   bsum:= bsum +     671 * sin(4*DD-2*M1+F);   bsum:= bsum +     607 * sin(2*DD-3*F);   bsum:= bsum +     596 * sin(2*DD+2*M1-F);   bsum:= bsum +     491 * sin(2*DD-M+M1-F) *ee;   bsum:= bsum -     451 * sin(2*DD-2*M1+F);   bsum:= bsum +     439 * sin(3*M1-F);   bsum:= bsum +     422 * sin(2*DD+2*M1+F);   bsum:= bsum +     421 * sin(2*DD-3*M1-F);   bsum:= bsum -     366 * sin(2*DD+M-M1+F) *ee;   bsum:= bsum -     351 * sin(2*DD+M+F) *ee;   bsum:= bsum +     331 * sin(4*DD+F);   bsum:= bsum +     315 * sin(2*DD-M+M1+F) *ee;   bsum:= bsum +     302 * sin(2*DD-2*M-F) *ee^2;   bsum:= bsum -     283 * sin(M1+3*F);   bsum:= bsum -     229 * sin(2*DD+M+M1-F) *ee;   bsum:= bsum +     223 * sin(DD+M-F) *ee;   bsum:= bsum +     223 * sin(DD+M+F) *ee;   bsum:= bsum -     220 * sin(M-2*M1-F) *ee;   bsum:= bsum -     220 * sin(2*DD+M-M1-F) *ee;   bsum:= bsum -     185 * sin(DD+M1+F);   bsum:= bsum +     181 * sin(2*DD-M-2*M1-F) *ee;   bsum:= bsum -     177 * sin(M+2*M1+F) *ee;   bsum:= bsum +     176 * sin(4*DD-2*M1-F);   bsum:= bsum +     166 * sin(4*DD-M-M1-F) *ee;   bsum:= bsum -     164 * sin(DD+M1-F);   bsum:= bsum +     132 * sin(4*DD+M1-F);   bsum:= bsum -     119 * sin(DD-M1-F);   bsum:= bsum +     115 * sin(4*DD-M-F) *ee;   bsum:= bsum +     107 * sin(2*DD-2*M+F) *ee^2;   bsum:= bsum-2235*sin(L1)+382*sin(a3)+175*sin(a1-F)+175*sin(a1+F)+127*sin(L1-M1)-115*sin(L1+M1);   b := →HMS(bsum/1000000); // β div by 1 milion   // correzioni per la distanza dalla Terra (Δ)   rsum:=0;   rsum:= rsum -20905355 * cos(M1); // coseni   rsum:= rsum - 3699111 * cos(2*DD-M1); // evezione   rsum:= rsum - 2955968 * cos(2*DD); // variazione   rsum:= rsum -  569925 * cos(2*M1); // equazione del centro   rsum:= rsum +   48888 * cos(M) *ee; // equazione annuale   rsum:= rsum -    3149 * cos(2*F); // riduzione all'eclittica   rsum:= rsum +  246158 * cos(2*DD-2*M1);   rsum:= rsum -  152138 * cos(2*DD-M1-M) *ee;   rsum:= rsum -  170733 * cos(2*DD+M1);   rsum:= rsum -  204586 * cos(2*DD-M) *ee;   rsum:= rsum -  129620 * cos(M-M1) *ee;   rsum:= rsum +  108743 * cos(DD);   rsum:= rsum +  104755 * cos(M1+M) *ee;   rsum:= rsum +   10321 * cos(2*DD-2*F);   rsum:= rsum +   79661 * cos(M1-2*F);   rsum:= rsum -   34782 * cos(4*DD-M1);   rsum:= rsum -   23210 * cos(3*M1);   rsum:= rsum -   21636 * cos(4*DD-2*M1);   rsum:= rsum +   24208 * cos(2*DD+M-M1) *ee;   rsum:= rsum +   30824 * cos(2*DD+M) *ee;   rsum:= rsum -    8379 * cos(DD-M1);   rsum:= rsum -   16675 * cos(DD+M) *ee;   rsum:= rsum -   12831 * cos(2*DD-M+M1) *ee;   rsum:= rsum -   10445 * cos(2*DD+2*M1);   rsum:= rsum -   11650 * cos(4*DD);   rsum:= rsum +   14403 * cos(2*DD-3*M1);   rsum:= rsum -    7003 * cos(M-2*M1) *ee;   rsum:= rsum +   10056 * cos(2*DD-M-2*M1) * ee;   rsum:= rsum +    6322 * cos(DD+M1);   rsum:= rsum -    9884 * cos(2*DD-2*M) * ee^2;   rsum:= rsum +    5751 * cos(M+2*M1) *ee;   rsum:= rsum -    4950 * cos(2*DD-2*M-M1) *ee^2;   rsum:= rsum +    4130 * cos(2*DD+M1-2*F);   rsum:= rsum -    3958 * cos(4*DD-M-M1) *ee;   rsum:= rsum +    3258 * cos(3*DD-M1);   rsum:= rsum +    2616 * cos(2*DD+M+M1) *ee;   rsum:= rsum -    1897 * cos(4*DD-M-2*M1) *ee;   rsum:= rsum -    2117 * cos(2*M-M1) * ee^2;   rsum:= rsum +    2354 * cos(2*DD+2*M-M1) *ee^2;   rsum:= rsum -    1423 * cos(4*DD+M1);   rsum:= rsum -    1117 * cos(4*M1);   rsum:= rsum -    1571 * cos(4*DD-M) *ee;   rsum:= rsum -    1739 * cos(DD-2*M1);   rsum:= rsum -    4421 * cos(2*M1-2*F);   rsum:= rsum +    1165 * cos(2*M+M1) *ee^2;   rsum:= rsum +    8752 * cos(2*DD-M1-2*F);   rsum:= rsum/1000;  // div by 1000     temp:= transformEclipticalToEquatorial(l,b);  // transform into AR, decl   α:= temp(1);   δ:= temp(2);   ts:= calcTS();  // calcola tempo siderale   Hd:=15*HMS→(ts(2)-α);  // local ST vs apparent ST   Hd:=simp_angle(→HMS(Hd));   H:= Hd/15; // in HMS   temp:= transformEquatorialToHorizontal(Hd, lat, δ);     az:= HMS→(temp(1))+180; // azimuth   IF az>360 THEN az:=az-360 END;   az:= →HMS(az);   h:= temp(2);  // height   // Raggio quadratico medio della Terra: 6373,044737   d:= (385000.56 + rsum); // Δ distanza Terra-Luna in km (=1/sin(P)),rr in 10^-3   //P:= →HMS((0.9508+0.0518*cos(M)) MOD 360); // parallasse   P:= asin(6378.14/d); // Parallasse   s:= asin((3/11)*sin(P)); // semidiametro (raggi terrestri)   DA:= asin(tan(lat)*tan(δ)); // Differenza ascensionale   SAD:= (90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno     // culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) ; // Culminazione (Bouiges)     // sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)     // tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   // transit routine   h0:= 0.7275*P - HMS→(0°34′); // 'standard' altitude of Moon (depend on parallax)   // h0:= 0°7′30″; // medium value: 0.125°   transito:= transit(α, δ, h0);   culm:= transito(1);   sorg:= transito(2);   tram:= transito(3);   // end transit   ph:= abs(l - L); // Phase  Moon-Sun   CASE   // IF ph == 0 THEN fase:="New Moon"; END;   IF ph>280 AND ph<350 THEN fase:="Decrescente"; END;   IF ph>=260 AND ph<=280 THEN fase:="Last Quart"; END;   IF ph>190 AND ph<260 THEN fase:="Gibbosa Cal."; END;   IF ph>=170 AND ph<=190 THEN fase:="Full Moon"; END;   IF ph>100 AND ph<170 THEN fase:="Gibbosa Cresc."; END;   IF ph>=80 AND ph<=100 THEN fase:="1st Quart"; END;   IF ph>10 AND ph<80 THEN fase:="Crescente"; END;   DEFAULT fase:="New Moon";   END; // case   // età della Luna   d1:= IP(HMS→(ph/15));   h1:= IP(FP(HMS→(ph/15))*60);   IF (h1>24) THEN d1:= d1+1; h1:= h1 MOD 24; END;   m1:= IP(FP(FP(FP(HMS→(ph/15))*60))*60);   eta:= →HMS(d1+h1/60+m1/3600); // Age of the Moon   segno:= zodiaco(l);   ill:= 1-(1+cos(ph))/2; // illuminated fraction   PRINT;   PRINT("Effemeridi della Luna a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L1));   PRINT("M anomalia vera " + trunc_angle(M1));   PRINT("Ө Longitudine " + trunc_angle(l) + " (" + trunc_angle(l MOD 30) + " " + segno + ")");   PRINT("β Latitudine " + trunc_angle(b));   PRINT("λ Longitudine apparente "+ trunc_angle(lapp));   PRINT("Δψ Nutazione long "+ dep(1));   PRINT("Δε Nutazione obl. " + dep(2));   PRINT("ε inclinazione ecl. " + ε);   PRINT("ω Perigeo medio " + trunc_angle(ω));   PRINT("Ω Nodo ascendente " + trunc_angle(Ω));   PRINT(" ");   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("α Ascensione retta (gradi) " + trunc_angle((α*15) MOD 360));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS locale " + trunc_angle(ts(2)));   PRINT("TS apparente " + ts(3));   PRINT("θ0 Mean ST GMT " + ts(4));   PRINT("Julian day " + JD);   PRINT("Diff. asc." + trunc_angle(DA) + " SAD " + trunc_angle(SAD));   PRINT("P Parallasse " + trunc_angle(P));   PRINT("Δ Distanza Terra-Luna "+ ROUND(d,4) + " km");   PRINT("Semidiametro "+  trunc_angle(→HMS(s)) + " :: " + ROUND(s*6378.14,3) + " km");   PRINT(" ");   PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");   PRINT(" ");   PRINT("Fase " + trunc_angle(ph) + " " + fase + " " + ROUND(100*→HMS(ph)/360,2)+"%");   PRINT("Età " + eta);   PRINT("Parte illuminata " + ill);   RETURN [[L1, M1],[l, b],[lapp,P], [ω,Ω],[α,δ], [az, h],      [ROUND(HMS→(d),6),s], [sorg, tram], [culm, ts(3)],[ph, eta] ]; END; planetCalc(nameplanet) BEGIN   LOCAL ω,Ω, l, v, r, x, y, z;   LOCAL Ls, Ms, b, xs, ys, b_ec;   LOCAL x1, y1, z1, R, l_ec, H, Hd;   LOCAL ec, in, a, δ, α, ts, az, h;   LOCAL DA, SAD, tram, sorg, culm, dayY;   LOCAL as, ac, zc, dz, θ, mov, segno, ε;   LOCAL transito, h0, JD, temp;   JD:= julianDay(dateshortlist);   // calc of Sun (as basis)   xs:= sunlist(9); ys:= sunlist(10);   Ls:= sunlist(1); // Longitudine media   Ms:= sunlist(4); // Anomalia   // end Sun  EVAL(EXPR(nameplanet+"()")); // call planet-name function   L:= planet_const(1); // Longitudine media   ω:= planet_const(2); // perielio   Ω:= planet_const(3); // nodo   M:= planet_const(4); // anomalia   v:= planet_const(5);    ec:= planet_const(6); // eccentricità   in:= planet_const(7); // inclinazione   a:= planet_const(8);  // semiasse maggiore   θ:= planet_const(9);  // θs to calc if retrogade   b_ec:= asin(sin(in)*sin(v+ω-Ω)); // argomento di latitudine del perielio   l_ec:= acos(cos(v+ω-Ω)/cos(b_ec));  // long from eclittic   IF b_ec<0 THEN l_ec:= 360 - l_ec; END; // argomento secondario   l_ec:= →HMS((l_ec + Ω)) MOD 360; // aggiunge Nodo e riporta a 2PI   r:= a*(1-ec^2)/(1+ec*cos(v)); // 9.524899/(1+ 0.0559*cos(v)); raggio vettore   x1:= r*cos(b_ec)*cos(l_ec); y1:= r*cos(b_ec)*sin(l_ec); z1:=r*sin(b_ec); // cartesian   x:= x1+xs; y:= y1+ys; z:=z1; // cohordinates cartesian->spheric   R:= sqrt(x^2+y^2+z^2);    b:=asin(z/R); // latitudine  β   IF b> 360 THEN b:= b MOD 360; END;   // l:=→HMS(atan(y/x)); // cart->spheric longtitudine λ   l:= atan2(y,x);   IF l<0 THEN l:= l+360; END; IF l>360 THEN l:=l-360; END;   IF l>180 THEN α:= α +180; END;   // δ:= →HMS((asin(cos(ε)*sin(b)+sin(ε)*cos(b)*sin(l))) ); // declinazione (coord locali)   // α:= (acos(cos(b)*cos(IFTE(l>180, l-180, l))/cos(δ))); // ascensione retta   //α:= →HMS(α/15) MOD 24; // α in hms   temp:= transformEclipticalToEquatorial(l,b);  // transform into AR, decl   α:= temp(1);   δ:= temp(2);   ts:= calcTS();  // calcola tempo siderale   H:= (ts(2) - α) MOD 24; // Angolo orario   Hd:= →HMS(15*H); // Angolo orario in gradi   zc:= sin(lat)*sin(δ)+cos(lat)*cos(δ)*cos(Hd); // calc z, alt   dz:= →HMS(acos(zc)) ; // distanza zenitale :: h:= →HMS(90 - dz); // Altezza   h:= →HMS(ATAN(zc/sqrt(1-zc^2))); // Altezza (height)   as:= cos(δ)*sin(Hd)/sin(dz); // calc az   ac:= (-cos(lat)*sin(δ)+sin(lat)*cos(δ)*cos(Hd))/sin(dz);   // az:= →HMS(ATAN(as/ac) ); // Azimuth   az:= atan2((-cos(δ) * sin(Hd)), cos(lat) *sin(δ)-sin(lat)*cos(δ)*cos(Hd));   IF az<0 THEN az:= az+360; END; IF az>360 THEN az:=az-360; END;   DA:= →HMS(asin(tan(lat)*tan(δ))); // Differenza ascensionale   SAD:= →HMS(90 + DA); // Semiarco Diurno   // SAD:= acos(-tan(lat)*tan(δ)); // Semiarco Diurno   // culm:= (IFTE(α<ts(1),24+α,α) - ts(1)) MOD 24 ; // Culminazione   // sorg:= (culm - SAD/15) MOD 24; // Alba (Dawn)   // tram:= (culm + SAD/15) MOD 24; // Tramonto (Set)   // transit routine   h0:= -0°34′; // 'standard' altitude of a planet   transito:= transit(α, δ, h0);   culm:= transito(1);   sorg:= transito(2);   tram:= transito(3);   // end transit   IF abs(l-220.4)>θ THEN mov:="diretto"; ELSE mov:="retrogrado"; END;   segno:= zodiaco(l);   PRINT;   PRINT("Effemeridi di " + nameplanet + " a " + m+" "+D+" "+Y+" "+ lstd);   PRINT("L Long media " + trunc_angle(L));   PRINT("ω Perielio " + trunc_angle(ω));   PRINT("Ω Nodo " + trunc_angle(Ω));   PRINT("M Anomalia vera " + trunc_angle(M));   PRINT("l_ec Long. eclittica " + trunc_angle(l_ec));   PRINT("b_ec Lat. eclittica " + trunc_angle(b_ec) + " v " + trunc_angle(v));   PRINT("r Raggio vettore " + r);   PRINT("Cartes. elioc.: x' " + ROUND(x1,3) + " y' " + ROUND(y1,3) + " z' " + ROUND(z1,3));   PRINT("Cartes. geoc.: x " + ROUND(x,3) + " y " + ROUND(y,3) + " z " + ROUND(z,3));   PRINT("Spher.: λ Longitudine " + trunc_angle(l) + " (" +trunc_angle(l MOD 30) + " " + segno + ")");   PRINT("Spher.: β Latitudine " + trunc_angle(b));   PRINT("Spher.: R Distanza " + ROUND(R,3));   PRINT("α Asc. r. "+ trunc_angle(α) + " δ Decl. "+ trunc_angle(δ));   PRINT("Az Azimuth " + trunc_angle(az) + " h Alt. " + trunc_angle(h));   PRINT("");   PRINT("Ang. orario "+ trunc_angle(H) + " " + trunc_angle(Hd));   PRINT("TS (UTC) " + trunc_angle(ts(1)));   PRINT("TS locale " + trunc_angle(ts(2)));   PRINT("TS apparente " + ts(3));   PRINT("θ0 Mean ST GMT " + ts(4));   PRINT("Julian day " + JD);   PRINT("Diff. ascensionale " + trunc_angle(DA));   PRINT("SAD semiarco diurno " + trunc_angle(SAD) + " :: " + trunc_angle(→HMS(SAD/15)));   PRINT(" ");   PRINT("Transita a: " + trunc_angle(culm) + " UTC (" + trunc_angle((culm+utc+dst) MOD 24) + " loc)");   PRINT("Sorge a: " + trunc_angle(sorg) + " UTC (" + trunc_angle((sorg+utc+dst) MOD 24) + " loc)");    PRINT("Tramonta a: " + trunc_angle(tram) + " UTC (" + trunc_angle((tram+utc+dst) MOD 24) + " loc)");   PRINT("Movimento " + mov);   RETURN [[L,ω,Ω ],[M,l_ec,b_ec], [v,r,0], [x1,y1,z1],[x,y,z],[l,b,R], [α,δ,H], [az,h,ts(1)],[sorg, tram, culm]]; END; Mercurius() BEGIN   // Costanti per Mercurio   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.205615; // eccentricità   in:= 7.003; // inclinazione   a:= 0.387098; // semiasse maggiore   θ:= 35.6;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   ω:= →HMS((75.913 + 0.00004253*N) MOD 360); // perielio   Ω:= →HMS((47.157 + 0.00003244*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 23.4372*sin(M)+2.9810*sin(2*M)+0.5396*sin(3*M)+0.1098*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Venus() BEGIN   // Costanti per Venere   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.006816; // eccentricità   in:= 3.393636; // inclinazione   a:= 0.72333; // semiasse maggiore   θ:= 13.0;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((206.758 + 1.60216873*N) MOD 360); // Longitudine media   ω:= →HMS((130.154 + 0.00003757*N) MOD 360); // perielio   Ω:= →HMS((75.797 + 0.00002503*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 0.7811*sin(M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Mars() BEGIN   // Costanti per Marte   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.093309; // eccentricità   in:= 1.850303; // inclinazione   a:= 1.523678; // semiasse maggiore   θ:= 16.8;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((124.765 + 0.52407108*N) MOD 360); // Longitudine media   ω:= →HMS((334.251 + 0.00005038*N) MOD 360); // perielio   Ω:= →HMS((48.794 + 0.00002127*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 10.608*sin(M)+0.6216*sin(2*M)+0.0504*sin(3*M)+0.0047*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Jupiter() BEGIN   // Costanti per Giove   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.048335; // eccentricità   in:= 1.308736; // inclinazione   a:= 5.202561; // semiasse maggiore   θ:= 54.43;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((268.350 + 0.08312942*N) MOD 360); // Longitudine media   ω:= →HMS((12.737 + 0.00004408*N) MOD 360); // perielio   Ω:= →HMS((99.453 + 0.00002767*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5372*sin(M)+0.1662*sin(2*M)+0.007*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Saturn() BEGIN   // Costanti per Saturno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.055892; // eccentricità   in:= 2.492519; // inclinazione   a:= 9.554747; // semiasse maggiore   θ:= 65.53;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((278.774 + 0.03349788*N) MOD 360); // Longitudine media   ω:= →HMS((91.117 + 0.00005362*N) MOD 360); // perielio   Ω:= →HMS((112.79 + 0.00002391*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 6.4023*sin(M)+0.2235*sin(2*M)+0.011*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Uranus() BEGIN   // Costanti per Urano   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.046344; // eccentricità   in:= 0.772464; // inclinazione   a:= 19.21814; // semiasse maggiore   θ:= 73.92;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((248.487 + 0.01176902*N) MOD 360); // Longitudine media   ω:= →HMS((171.563 + 0.00004064*N) MOD 360); // perielio   Ω:= →HMS((73.482 + 0.00001365*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 5.3092*sin(M)+0.1537*sin(2*M)+0.0062*sin(3*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Neptunus() BEGIN   // Costanti per Nettuno   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.008997; // eccentricità   in:= 1.779242; // inclinazione   a:= 30.10957; // semiasse maggiore   θ:= 77.63;  // θs to calc if retrogade   calcN();   // costanti planetarie   L:= →HMS((86.652 + 0.00602015*N) MOD 360); // Longitudine media   ω:= →HMS((46.742 + 0.000039*N) MOD 360); // perielio   Ω:= →HMS((130.693 + 0.00003009*N) MOD 360); // nodo   M:= (L - ω); // Anomalia vera   v:= M + 1.031*sin(M)+0.0058*sin(2*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END; Pluto() BEGIN   // Costanti per Plutone   LOCAL ω,Ω, L, M, v;   LOCAL ec, in, a, θ;   ec:= 0.250236; // eccentricità   in:= 17.17047; // inclinazione   a:= 39.438712; // semiasse maggiore   θ:= 79.41;  // θs to calc if retrogade   calcN();   // costanti planetarie   // Long, perielio, nodo, non precisi, considerati costanti   // L:= →HMS((229.851 + 4.09237703*N) MOD 360); // Longitudine media   M:= →HMS((230.67 + 0.00397943*N) MOD 360); // Anomalia vera (calcolata direttamente)   Ω:= →HMS((109.06 + 0.00003823*N) MOD 360); // nodo   ω:= →HMS((114.27 + Ω) MOD 360); // perielio   L:= →HMS((M + ω) MOD 360); // L (non sicuro)      v:= M + 28.45*sin(M)+4.38*sin(2*M)+0.97*sin(3*M)+0.24*sin(4*M);    planet_const:= {L, ω, Ω, M, v, ec, in, a, θ};   RETURN planet_const; END;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)