PROGRAMMA PER IL CALCOLO DEL POLO DELL'ORBITA APPARENTE DI UNA COMETA

Il programma da me scritto per questo calcolo è stato sviluppato in IDL (Interactive Data Language), un linguaggio specifico per scopi scientifici, ma può agevolmente essere riscritto in qualsiasi altro linguaggio di programmazione.

Di seguito riportiamo l'algoritmo con le formule utilizzate (tratte dal libro "Astronomical formulae for calculators" di Jean Meeus, Willmann-Bell Inc., 1985).

1. Lettura delle effemeridi da un file di testo:
   Il programma deve leggere le coordinate dalle effemeridi pubblicate nel sito della NASA, aventi sempre il formato:

     Date (UT)       R.A. J2000  Dec.     Delta Deldot    r   Theta Beta Moon PsAng PsAMV TMag

   delle quali ci interessano unicamente la data e le coordinate equatoriali.
   
   
2. Calcolo del Giorno Giuliano dalla data e ora delle effemeridi
   La Data Giuliana è calendario assoluto usato dagli astronomi per il calcolo delle effemeridi ed è un numero progressivo che conta i giorni a partire dalle 12:00 di T.U. invece che dalle 00:00.
   L'algoritmo è il seguente:
   
   1. Se il mese è da Marzo a Dicembre allora: y=anno, m=mese, D.d=giorno e decimali
   2. Se il mese è Gennaio o Febbraio allora: y=anno-1, m=mese+12, D.d=giorno e decimali
   3. A = INT(y/100)
   4. B = 2 - A + INT(A/4)
   5. Data Giuliana: jd = INT(365.25*y) + INT(30.6001*(m+1)) + D.d + 1720994.5 + B
   Attenzione: questa formula vale solo dopo il 4 Ottobre 1582 (inizio del calendario gregoriano).
   
   
   
3. Calcolo del Tempo Siderale
   
   1. T = (jd_alle_00:00_di T.U. - 2415020) / 36525
   2. TS0 = .276919398 + 100.0021359*T + .000001075* (T^2)
   3. TS1 = TS0 + (ora.decimali/24)*1.002737908
   4. Tempo Siderale in ore: TS = (TS1 - INT(TS1)) * 24
   
   
   
4. Conversione delle effemeridi (alfa, delta) in coordinate altazimutali (lat=latitudine)
   
   1. Angolo Orario: hg = TS - alfa
   2. Convertire tutti gli angoli da gradi in radianti
   3. azy=sin(hg)
   4. azx=(cos(hg)*sin(lat)-tan(delta)*cos(lat))
   5. alt=(sin(lat)*sin(delta)+cos(lat)*cos(delta)*cos(hg))
   6. az=atan(azy/azx)+pi_greco/2*(azx lt 0)
   7. alt=asin(alt)
   8. Convertire da radianti in gradi tra 0 e 360
   
   
   
5. Calcolo del centro del cerchio osculatore istantaneo dell'orbita altazimutale
   
   1. Per ogni posizione "i" delle coordinate altazimutali (az,alt) siano:
      
      • p1r=[az(i-1),alt(i-1)]
      • p2r=[az(i),alt(i)]
      • p3r=[az(i+1),alt(i+1)]
      
      
   2. Convertire in xyz:
      • lon = [p1r(0), p2r(0), p3r(0)]
      • lat = [p1r(1), p2r(1), p3r(1)]
      • x = cos(lat) * cos(lon)
      • y = cos(lat)* sin(lon)
      • z = sin(lat)
      
      
   3. Calcolo della normale al piano contenente i tre punti:
      
      • aa=(y(1)-y(0))*(z(2)-z(0))-(y(2)-y(0))*(z(1)-z(0))
      • bb=-((x(1)-x(0))*(z(2)-z(0))-(x(2)-x(0))*(z(1)-z(0)))
      • cc=(x(1)-x(0))*(y(2)-y(0))-(x(2)-x(0))*(y(1)-y(0))
      • dd=aa*x(0)+bb*y(0)+cc*z(0)
      • m=sqrt(aa^2+bb^2+cc^2)
      • nx=aa/m ;normale al piano
      • ny=bb/m
      • nz=cc/m
      • if nz lt 0 then nx=-nx & ny=-ny & nz=-nz
      
      
   4. Calcolo delle coordinate del polo nord (nz>0) del cerchio per i tre punti:
      
      • alt_polo=asin(nz/(sqrt(nx^2+ny^2+nz^2)))
      • if nx eq 0 then azimut_polo=PI/4+PI/2*(ny lt 0) else azimut_polo=atan(ny/nx)+PI/2*(nx lt 0)
      
      
   5. Riconvertire (alt_polo, azimut_polo) in gradi.
   
   
   
6. Conversione delle posizioni del polo (alt,az) nuovamente in coordinate equatoriali
   
   1. alfay=sin(az)
   2. alfax=(cos(az)*sin(lat)+tan(alt)*cos(lat))
   3. if alfax ne 0 then alfa_polo=atan(alfay/alfax)+PI/2*(il minore tra "alfax" e "0") else alfa_polo=PI/4+PI/2*(il minore tra "alfay" e "0")
   4. delta_polo= - (asin(cos(lat)*cos(alt)*cos(az)-sin(lat)*sin(alt)))
   5. Convertire da radianti in gradi tra 0 e 360
   6. Convertire alfa_polo in ore da 0 a 24
   7. alfa_polo=(TS - alfa_polo) mod 24
   
   
   
7. Scrittura delle effemeridi del polo (alfa_polo, delta_polo) in un file di testo

INTRODUZIONE