Astronomia

Problemi nel trovare i tempi di salita e di transito alle alte latitudini utilizzando gli algoritmi di Meeus

Problemi nel trovare i tempi di salita e di transito alle alte latitudini utilizzando gli algoritmi di Meeus

Sto implementando l'algoritmo di Jean Meeus dal suo libro "Astronomical Algorithms", 2a ed. per tempi di salita, tramonto e transito in Swift. Sembra fornire ottimi risultati (non più di 1 minuto di differenza) se confrontato ad esempio con timeanddate.com, anche alle alte latitudini (sopra i 66 gradi). Tuttavia, a volte non funziona.

L'ho ricondotto a circostanze in cui a) il sole sorge in un particolare giorno UTC ma non tramonta, o viceversa, b) dove il sole sorge due volte in un determinato giorno UTC. Un esempio specifico di quest'ultimo si verifica il 16 aprile 2020 a Longyearbyen (Svalbard) Norvegia. E infine c) quando un tempo di levata "sposta" indietro da dopo a prima della mezzanotte se guardato in giorni consecutivi.

I tempi approssimativi di salita, presa e transito, rispettivamente m1, m2 e m0, possono essere calcolati abbastanza facilmente. Devono quindi essere raffinati di un valore delta m. Questo può/potrebbe/dovrebbe essere fatto in modo iterativo.

Guardando più da vicino, trovo che per gli eventi sopra descritti, i valori per m1 e m2, quando regolati per delta m come descritto a pagina 103, in genere superano 1. Meeus afferma che i valori m "dovrebbero essere compresi tra 0 e 1. Se uno o più di loro sono al di fuori di questo intervallo, aggiungere o sottrarre 1". Un po' più avanti, nella Nota 1 alla fine del Capitolo 15, c'è una nota allettante che dice "se il tempo di impostazione [... ] è necessario nell'ora locale, il calcolo deve essere eseguito utilizzando [... ] m2 = 1,12113, che è maggiore di 1.

Questo mi fa sospettare - non sono un astronomo, come avrai intuito - che valori di m1 o m2 maggiori di 1 possano probabilmente aiutarmi a calcolare i tempi di salita in un giorno in cui ciò accade due volte, e anche a correggere il tempo di salita su un giorno in cui il sole non tramonta (e viceversa).

Questo mi ha portato a github, dove ho trovato del codice JavaScript (e altro) di onekiloparsec, che confronta il tempo di salita con il tempo di transito e, se il primo è successivo al secondo, impiega il tempo di salita per essere il giorno precedente. Allo stesso modo, se l'orario impostato è precedente all'orario di transito, viene considerato il giorno successivo.

Ho anche guardato "Astronomia pratica con la tua calcolatrice o foglio di calcolo" di Peter Duffett-Smith e Jonathan Zwart, ma non ho trovato una risposta lì. Ha fornito un'informazione molto utile che Meeus non fornisce, ovvero che il segno del risultato dell'equazione di Meeus 15.1, quando il risultato in valore assoluto è maggiore di 1, permette di distinguere se una stella è permanentemente al di sotto (cos H0 >1) o sopra l'orizzonte (cos H0 < -1).

Sarebbe bello avere una spiegazione o un riferimento e alcuni dettagli su come interpretare i risultati per m1 e m2, come sfortunatamente il libro di Meeus, mentre descrive gli algoritmi a un livello in cui un non astronomo può implementare molti dei calcoli, mi lascia con la domanda strana.

Di seguito è riportato il codice (Swift) che ho scritto per perfezionare in modo iterativo i valori per m0, m1 e m2.

func iterateForPreciseM(m_fractionOfDay:Double, wantTime:DesiredTime, calcoloMode:CalculationMode, debugmsg:String = "") -> (Double?, CalculationQualityLevel) { //funzione in linea. la modalità di calcolo permette di specificare se si vuole calcolare salita/tramontata o transito. //restituisce una frazione di giorno raffinata e un indicatore della qualità del risultato. La qualità "buona" significa che è stata calcolata con non più di 3 passaggi. La qualità "problematica" segnala che sono stati necessari più di 3 passaggi prima che deltaM raggiunga il limite di "convergenza", ma in meno di 20 cicli. Se più di 20 passaggi, la qualità viene impostata su "cattiva" per indicare la mancata convergenza. Sono arrivato a quei valori (3 e 20) arbitrariamente. //il tempo desiderato viene utilizzato per specificare se stiamo calcolando il transito o il tempo di salita e di tramonto. //calculationMode specifica se stiamo calcolando i tempi civili del crepuscolo o il sorgere e il tramonto del sole. NB funzione rawValue utilizzata. var m_fractionOfDay = m_fractionOfDay //shadow il valore passato poiché dovrò modificarlo. var loopCount = 1 , maxAcceptableLoopCount = 3 , maxLoopCount = 20 //limite di conteggio arbitrario per il ciclo. let deltamLimit = 0,0001 ///0,0001 è arbitrario. Finora 07/05/2020 osservo che molto spesso è un po' al di sopra di questo, ma alla seconda iterazione diventa ripetitiva infinitesimale { var small_theta0_degrees = GAST0_degrees + 360.985647 * m_fractionOfDay ///small theta0 è "tempo siderale a Greenwich", per AA Meeus, parte superiore di p103. Non so quale sia la differenza tra questo e il tempo siderale apparente di Greenwich. Forse il tempo siderale nella posizione dell'osservatore, dal momento che entra nel calcolo di m-fractionOfDay ? O, più probabilmente, come in AA Cap 12 p87, small_theta_0 è definito come tempo siderale a Greenwich per uno specifico istante UT. small_theta0_degrees = normalizedDegrees360(degrees: small_theta0_degrees) let n = m_fractionOfDay + ( deltaTseconds / constants.SECS_IN_DAY ) if abs(n) > 1 { if verbose { wl(#function,#line," --**n (n) outside di -1 a +1 intervallo - (debugmsg)") } } /* L'ascensione retta si trova sempre nell'intervallo da 0 a 360 gradi e aumenta continuamente con l'aumento della data. Tuttavia, quando raggiunge i 360 gradi, cosa che accade una volta all'anno all'equinozio di primavera (o credo più precisamente all'equinozio di primavera), si "avvolge" a 0. Per Wikipedia, "https://en.wikipedia.org/wiki /Right_ascension", RA viene solitamente misurata in ore, minuti e secondi, che vanno da 0 a 24. È interessante notare che l'articolo afferma anche che SHA è il complemento 24 ore di RA. La formula di interpolazione di Meeus (eq.3.3) deve essere regolata per gestire questo wrapping (sebbene ciò non sia dichiarato esplicitamente in AA - l'ho scoperto durante il tracciamento). Ciò significa che alcuni dei valori di AR dovranno essere aumentati di 360 gradi. */ ///copia i valori RA originali per l'intervallo di 3 giorni ottenuto in precedenza - poiché siamo in una funzione in linea che viene chiamata più volte e anche in loop, non possiamo modificare i valori originali. Potrei modificarli quando li calcolo per la prima volta, il che accade al di fuori di questa funzione in linea, ma farlo rende meno ovvio ciò che devo fare. let rightAscensionDegreesDay0 = rightAscensionDegrees[0] var rightAscensionDegreesDay1 = rightAscensionDegrees[1] var rightAscensionDegreesDay2 = rightAscensionDegrees[2] //ora regolali se l'ascensione retta aumenta di 360 gradi durante i 3 giorni per i quali stiamo interpolando. if rightAscensionDegreesDay1 < rightAscensionDegreesDay0 { //per il caso ra[2]=1,6 ra[1]=0,7 ra[0]=359,8 rightAscensionDegreesDay1 += 360 rightAscensionDegreesDay2 += 360 // ora ra[2]=361.6, ra[1] =360.7, ra[0] invariato 359.8 } // passando al controllo successivo non causerà ulteriori modifiche ai valori ra[]. if rightAscensionDegreesDay2 < rightAscensionDegreesDay1 { //per il caso ra[2]=0.7 ra[1]=359.8 ra[0]=358.9 rightAscensionDegreesDay2 += 360 // ora ra[2]= 360,7, ra[1] e ra[0 ] invariato. } let a1 = rightAscensionDegreesDay1 - rightAscensionDegreesDay0 let b1 = rightAscensionDegreesDay2 - rightAscensionDegreesDay1 let c1 = b1 - a1 let alpha_degrees :Double = normalizedDegrees360(gradi: rightAscensionDegrees) *) + 1 (a/) /necessità di normalizzare poiché alcuni casi di avvolgimento all'equinozio possono far sì che l'alfa vada leggermente al di sopra di 360. //interpola la declinazione usando l'eq.3.3 /* La declinazione PER IL SOLE varia da +23.4x a -23,4x gradi. Sale sopra lo 0 all'equinozio di primavera, raggiunge il picco al solstizio d'estate, poi scende attraverso lo 0 all'equinozio d'autunno, raggiunge il minimo al solstizio d'inverno e risale. I test rivelano che la formula di interpolazione di Meeus gestisce correttamente i punti di flesso ai solstizi, nonché il passaggio da negativo a positivo e viceversa, senza richiedere adattamenti come nel caso dell'ascensione retta. */ let a2 = gradi declinazione[1] - gradi declinazione[0] let b2 = gradi declinazione[2] - gradi declinazione[1] let c2 = b2 - a2 let gradi_delta :Double = gradi declinazione[1] + (n/2.0) * (a2 + b2 + n * c2 ) //calcola H - questo è l'LHA var H_degrees = small_theta0_degrees - osservatoreLongitudeDegrees - alpha_degrees //Riporta H (LHA) nell'intervallo da -180 a +180 - Per Meeus Cap 15 p103 H_degrees = normalized(PlusMinus180) angleDegrees: H_degrees) //calcola il deltaM, sia per il transito che per l'ascesa/imposta var deltam:Double = 0 var sin_h:Double = 0; var height_degrees:Double = 0 //per la traccia, definire all'esterno dell'interruttore. Otwz entrambi possono essere definiti all'interno dell'interruttore, non necessari all'esterno. commutare l'ora desiderata { case .transit: //deltaM per il transito cap 15 p103 deltam = -H_degrees / 360 case .riseSet: //calcola l'altitudine del Sole ///AA eq. 13.6 sin_h = sin(radianti(gradi: ObserverLatitudeDegrees)) * sin(radianti(gradi: delta_degrees)) + cos(radianti(gradi: ObserverLatitudeDegrees)) * cos(radianti(gradi: delta_degrees)) * cos(radianti(gradi: H_degrees) )) if abs(sin_h) > 1 { // FIXME: asin può restituire NaN se abs(sin_h) è maggiore di 1. Per ora lascerò che accada. Dovrebbe trovare un modo per gestire questa situazione. } altitudine_gradi = gradi(radianti:asin(sin_h)) // deltaM per salire e impostare Cap 15 p 103 let geometricAltitudeOfCelestialBodyCenter_degrees = calcoloMode.rawValue deltam = (altitudine_gradi - geometricAltitudeOfCelestialBodyCenter_degrees)(360.0 *gradi_corpo celeste) (360.0 *radianti cos(radianti(gradi: ObserverLatitudeDegrees)) * sin(radianti(gradi: H_degrees)) ) // FIXME: guarda dalla divisione per 0 - ovunque in questa classe! Se la latitudine dell'osservatore è 90N/S, div per 0!!! } //endswitch m_fractionOfDay += deltam if m_fractionOfDay > 1.0 { wl(#function,#line,"!! --m_frac WENT ABOVE 1 = (debugmsg) -: (m_fractionOfDay) :- at loop #(loopCount) (calculationMode) (desiredTime)") } if m_fractionOfDay < 0.0 { wl(#function,#line,"!! --m_frac WENT BELOW 0 = (debugmsg) -: (m_fractionOfDay) :- at loop # (loopCount) (calculationMode) (desiredTime)") } if fabs(deltam) < deltamLimit { if loopCount > maxAcceptableLoopCount { // conteggio dei loop anormalmente alto all'uscita - m:(m_fractionOfDay) break } if loopCount > maxLoopCount { / //solo per scopi di debug. // maxLoopCount EXCEEDED break } loopCount += 1 } while true if loopCount > maxLoopCount { return (m_fractionOfDay, CalculationQualityLevel.bad) } if loopCount > maxAcceptableLoopCount { return (m_fractionOfDay, CalculationQualityLevel.problematic) } return Quality}, m_fraction ///fine funzione in linea

Per fare un po' di background, sono un pilota, non un astronomo, e in passato ho dovuto imparare tecniche di navigazione in aree polari dove la bussola magnetica è inaffidabile. Quindi dovevamo sapere come usare il sole come bussola: trovando la sua vera direzione e poi puntandoci il muso dell'aereo potevamo stabilire una direzione. L'altra tecnica prevedeva l'uso di un'astrobussola installata sull'aeromobile. Ha funzionato, a patto di non aver dimenticato di portare con sé un almanacco o dei tavoli!!

Certo, lo so, in questi giorni, il GPS... Questo è principalmente per il gusto di farlo, e come backup. Il calcolo del vero rilevamento del sole penso non sarà un problema poiché posso già calcolare l'AR del sole e la declinazione (usando l'adattamento di Meeus di VSOP87). Ma sarebbe bene sapere in anticipo se il sole sta per sorgere quando ho intenzione di puntargli un'astrobussola. E il motivo per cui cerco la risposta alle alte latitudini è che, dopotutto, dovrebbe essere utilizzabile "al nord".


Rispondendo alla mia domanda qui. Ho passato un bel po' di tempo a indagare su questo, questo è quello che mi è venuto in mente.

Per quanto riguarda l'interpretazione dei valori di salita e cala: se m1 o m2 diventano negativi, ho trovato solo che significa che la salita o il tramonto si verifica il giorno Zulu precedente. E se i valori superano 1, si è verificato il giorno Zulu successivo.

Sono rimasto sorpreso dal fatto che per alcune date e località, l'algoritmo dia un tempo di salita che è successivo al tempo impostato. All'inizio sembrava dubbio, ma ora sembra essere del tutto plausibile, come confermano altre fonti. Ho anche concluso che ciò non implica in alcun modo che l'alba possa essere in una data diversa dal tramonto - in termini Zulu, i risultati compresi tra 0 e 1 sono nella data per la quale viene effettuato il calcolo, indipendentemente dalla loro sequenza.

Per quanto riguarda la ripetizione del calcolo (di m1 o m2) fino a quando deltaM diventa abbastanza piccolo, ho scoperto che se le ripetizioni superano 20, questa è un'ottima indicazione che in realtà non c'è aumento (o set) in quel giorno. Anche questo sembra essere coerente con la realtà, in quanto ci sono infatti giorni in cui c'è solo l'alba ma non il tramonto (o viceversa). Quindi anche quell'aspetto dell'implementazione del calcolo che ho postato sembra corretto.

Inoltre: i risultati dell'algoritmo di Meeus mostrano che la sequenza di sorgere e tramontare si inverte nel tempo, a causa di un tramonto "mancante". Ad esempio, potresti avere Rise-Set, Rise-Set, Rise-(NO set), e quindi l'ordine cambia in Set-Rise, Set-Rise ecc. Anche coerente con la realtà per quanto posso dire. E il giorno in cui non c'è il tramonto può essere rilevato, come descritto nel paragrafo precedente.

Tuttavia, esistono nella vita reale alcune circostanze in cui in un giorno Zulu, il sole sorge, rimane alto a lungo, poi tramonta e poi sorge di nuovo. Sebbene il libro di Meeus non dica come determinare l'ora del secondo tramonto, posso almeno rilevare che accade perché per quel giorno ho una sequenza di alzarsi e tramontare (senza alba o tramonto mancanti), e il prossimo giorno sarà una sequenza di alzate.

Questa particolare circostanza sembra essere piuttosto rara, e il lavoro di Meeus rimane molto impressionante in quanto mi ha permesso di implementarlo senza alcuna precedente conoscenza dell'astronomia. E ha sicuramente stuzzicato la mia curiosità: sono uscito e ho comprato un telescopio dopo aver passato del tempo con il suo libro!

Quindi sono contento di questo in questa fase, anche se potrei ancora guardare la libreria CSPICE a cui fa riferimento User21 in un commento a questa domanda.

PS. C'era un errore nel codice di esempio che ho pubblicato, che ora ho rimosso. Il nuovo valore di m1 (o m2) ottenuto aggiungendo il deltaM non dovrebbe essere "rinormalizzato" per cadere tra 0 e 1, in questo modo si ottengono risultati del tutto errati.