Ich implementiere den Algorithmus von Jean Meeus aus seinem Buch "Astronomical Algorithms", 2nd ed. für Anstiegs-, Untergangs- und Laufzeiten in Swift. Es scheint sehr gute Ergebnisse zu liefern (nicht mehr als 1 Minute Unterschied), wenn man es beispielsweise mit timeanddate.com vergleicht, auch in hohen Breiten (über 66 Grad). Manchmal funktioniert es jedoch nicht.
Ich habe dies auf Umstände zurückgeführt, wo a) die Sonne an einem bestimmten UTC-Tag aufgeht, aber nicht untergeht, oder umgekehrt, b) wo die Sonne zweimal an einem bestimmten UTC-Tag aufgeht. Ein konkretes Beispiel für Letzteres ereignet sich am 16. April 2020 in Longyearbyen (Spitzbergen), Norwegen. Und schließlich c) wenn sich eine Anstiegszeit bei Betrachtung an aufeinanderfolgenden Tagen von nach auf vor Mitternacht zurückbewegt.
Ungefähre Anstiegs-, Setz- und Laufzeiten, bzw. m1, m2 und m0, können ziemlich einfach berechnet werden. Sie müssen dann mit einem Delta-m-Wert verfeinert werden. Dies kann/darf/sollte iterativ erfolgen.
Bei näherer Betrachtung stelle ich fest, dass für die oben beschriebenen Ereignisse die Werte für m1 und m2, wenn sie wie auf Seite 103 beschrieben um Delta m angepasst werden, typischerweise 1 überschreiten. Meeus gibt an, dass die m-Werte „zwischen 0 und 1 liegen sollten mehr von ihnen außerhalb dieses Bereichs liegen, addieren oder subtrahieren Sie 1". Etwas weiter, in Anmerkung 1 am Ende von Kapitel 15, gibt es einen verlockenden Hinweis, der besagt: „Wenn die Zeit der Einstellung [..] in Ortszeit benötigt wird, sollte die Berechnung mit [...] m2 = durchgeführt werden 1,12113, was größer als 1 ist.
Dies lässt mich vermuten – ich bin kein Astronom, wie Sie erraten haben – dass Werte von m1 oder m2 größer als 1 mir wahrscheinlich helfen können, die Anstiegszeiten an einem Tag zu berechnen, an dem dies zweimal vorkommt, und auch die Anstiegszeit an a zu korrigieren Tag, an dem die Sonne nicht untergeht (und umgekehrt).
Dies führte mich zu github, wo ich JavaScript-Code (und mehr) von onekiloparsec fand, der die Anstiegszeit mit der Transitzeit vergleicht, und wenn die erste später als die zweite ist, nimmt die Anstiegszeit den Vortag an. Wenn die eingestellte Zeit vor der Transitzeit liegt, wird sie in ähnlicher Weise als am folgenden Tag angenommen.
Ich habe mir auch "Practical Astronomy with your Calculator or Spreadsheet" von Peter Duffett-Smith und Jonathan Zwart angesehen, aber dort keine Antwort gefunden. Es lieferte eine sehr nützliche Information, die Meeus nicht liefert, nämlich dass das Vorzeichen des Ergebnisses von Meeus' Gleichung 15.1, wenn das Ergebnis im Absolutwert größer als 1 ist, es ermöglicht, zu unterscheiden, ob ein Stern dauerhaft unter (cos H0 >1) oder über dem Horizont (cos H0 < -1).
Es wäre großartig, eine Erklärung oder eine Referenz und einige Details zur Interpretation der Ergebnisse für m1 und m2 zu erhalten, wie leider Meeus 'Buch, während die Algorithmen auf einem Niveau beschrieben werden, auf dem ein Nicht-Astronom viele der Berechnungen implementieren kann. lässt mich mit der seltsamen Frage zurück.
Unten ist der (Swift-)Code, den ich geschrieben habe, um die Werte für m0, m1 und m2 iterativ zu verfeinern.
func iterateForPreciseM(m_fractionOfDay:Double, desiredTime:DesiredTime, calculationMode:CalculationMode, debugmsg:String = "") -> (Double?, CalculationQualityLevel) { //inline function. calculation mode allows to specify if rise/set or transit is to be calculated.
//returns refined fraction of day and an indicator of result quality. Quality "good" means it was calculated with no more than 3 passes. Quality "problematic" signals that more than 3 passes were required before deltaM reached "convergence" limit, but in less than 20 loops. If more than 20 passes, quality is set to "bad" to indicate failure to converge. I arrived at those values (3 and 20) arbitrarily.
//desired time is used to specify whether we are calculating transit, or rise & set time.
//calculationMode specifies whether we are calculating civil twilight times or sun rise & set times. NB rawValue feature used.
var m_fractionOfDay = m_fractionOfDay //shadow the passed-in value as I will need to modify it.
var loopCount = 1 , maxAcceptableLoopCount = 3 , maxLoopCount = 20 //arbitrary count limit for the loop.
let deltamLimit = 0.0001 ///0.0001 is arbitrary. So far 2020-05-07 I observe that it very often is a little above this, but on the second iteration becomes infinitesimal
repeat {
var small_theta0_degrees = GAST0_degrees + 360.985647 * m_fractionOfDay ///small theta0 is "sidereal time at Greenwich", per AA Meeus, top of p103. Don't know what the difference between that and Greenwich Apparent Sidereal Time means. Perhaps sidereal time at the observer's location, since that enters into the calculation of m-fractionOfDay ? Or, more likely, as in AA Chap 12 p87, small_theta_0 is defined as sidereal time at Greenwich for a specific instant 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 of -1 to +1 range - \(debugmsg)") }
}
/* Right ascension always lies in the range 0 to 360 degrees, and continuously increases with an increase in date. However when it reaches 360 degrees, which happens once a year at the spring (or i believe more accurately at the vernal) equinox, it "wraps around" to 0.
Per Wikipedia, "https://en.wikipedia.org/wiki/Right_ascension", RA is customarily measured in hours, minutes and seconds, ranging from 0 to 24. Interestingly, the article also states that SHA is the 24h-complement of RA.
Meeus' interpolation formula (eq.3.3) needs to be adjusted to handle this wrapping, (though this is not stated explicitly in AA - I discovered it during tracing). This means some of the RA values will need to be increased by 360 degrees.
*/
///copy the original RA values for the 3-day range obtained previously - since we are in an inline function which gets called multiple times and loops as well, we cannot modify the original values. I could modify them when I first calculate them - which happens outside this inline function, but doing this there makes it less obvious what I need to do.
let rightAscensionDegreesDay0 = rightAscensionDegrees[0]
var rightAscensionDegreesDay1 = rightAscensionDegrees[1]
var rightAscensionDegreesDay2 = rightAscensionDegrees[2]
//now adjust them if right ascension increases through 360 degrees during the 3 days for which we are interpolating.
if rightAscensionDegreesDay1 < rightAscensionDegreesDay0 { //for the case ra[2]=1.6 ra[1]=0.7 ra[0]=359.8
rightAscensionDegreesDay1 += 360
rightAscensionDegreesDay2 += 360 // now ra[2]=361.6, ra[1]=360.7, ra[0] unchanged 359.8
} // falling through to next check won't cause further modification to ra[] values.
if rightAscensionDegreesDay2 < rightAscensionDegreesDay1 { //for the case ra[2]=0.7 ra[1]=359.8 ra[0]=358.9
rightAscensionDegreesDay2 += 360 // now ra[2]= 360.7, ra[1] and ra[0] unchanged.
}
let a1 = rightAscensionDegreesDay1 - rightAscensionDegreesDay0
let b1 = rightAscensionDegreesDay2 - rightAscensionDegreesDay1
let c1 = b1 - a1
let alpha_degrees :Double = normalizedDegrees360(degrees: rightAscensionDegrees[1] + (n/2.0) * (a1 + b1 + n * c1 )) //need to normalize as some cases of wrapping at the equinox may cause alpha to go slightly above 360.
//interpolate declination using eq.3.3
/* Declination FOR THE SUN ranges from +23.4x to -23.4x degrees. It rises above 0 at the spring equinox, peaks at summer solstice, then descends through 0 at the fall equinox, bottoms out at winter solstice and rises again.
Tests reveal that Meeus' interpolation formula correctly handles inflection points at the solstices as well as passage from negative to positive and vice-versa, without requiring adaptation as was the case for right ascension.
*/
let a2 = declinationDegrees[1] - declinationDegrees[0]
let b2 = declinationDegrees[2] - declinationDegrees[1]
let c2 = b2 - a2
let delta_degrees :Double = declinationDegrees[1] + (n/2.0) * (a2 + b2 + n * c2 )
//calculate H - this is the LHA
var H_degrees = small_theta0_degrees - observerLongitudeDegrees - alpha_degrees
//Bring H (LHA) back into the -180 to +180 range - Per Meeus Chap 15 p103
H_degrees = normalizedDegreesPlusMinus180(angleDegrees: H_degrees)
//calculate the deltaM, for either transit or for rise/set
var deltam:Double = 0
var sin_h:Double = 0; var altitude_degrees:Double = 0 //for tracing, define outside the switch. Otwz both can be defined inside switch, not needed outside.
switch desiredTime {
case .transit:
//deltaM for transit chap 15 p103
deltam = -H_degrees / 360
case .riseSet:
//calculate Sun's altitude
///AA eq. 13.6
sin_h = sin(radians(degrees: observerLatitudeDegrees)) * sin(radians(degrees: delta_degrees)) + cos(radians(degrees: observerLatitudeDegrees)) * cos(radians(degrees: delta_degrees)) * cos(radians(degrees: H_degrees))
if abs(sin_h) > 1 {
// FIXME: asin may return NaN if abs(sin_h) is greater than 1. For now I will let this happen. Should find a way to handle this situation.
}
altitude_degrees = degrees(radians:asin(sin_h))
// deltaM for rise and set Chap 15 p 103
let geometricAltitudeOfCelestialBodyCenter_degrees = calculationMode.rawValue
deltam = ( altitude_degrees - geometricAltitudeOfCelestialBodyCenter_degrees ) / (360.0 * cos(radians(degrees: delta_degrees)) * cos(radians(degrees: observerLatitudeDegrees)) * sin(radians(degrees: H_degrees)) )
// FIXME: guard against division by 0 - everywhere in this class! If the observer latitude is 90N/S, div by 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 {
// abnormally high loop count at exit - m:\(m_fractionOfDay)
break
}
if loopCount > maxLoopCount { ///for debugging purposes only.
// maxLoopCount EXCEEDED
break
}
loopCount += 1
} while true
if loopCount > maxLoopCount {
return (m_fractionOfDay, CalculationQualityLevel.bad)
}
if loopCount > maxAcceptableLoopCount {
return (m_fractionOfDay, CalculationQualityLevel.problematic)
}
return (m_fractionOfDay, CalculationQualityLevel.good)
} ///end inline func
Als kleiner Hintergrund: Ich bin Pilot, kein Astronom, und musste in der Vergangenheit Navigationstechniken in Polargebieten lernen, in denen der Magnetkompass unzuverlässig ist. Also mussten wir wissen, wie man die Sonne als Kompass benutzt – indem wir ihre wahre Peilung fanden und dann die Nase des Flugzeugs darauf richteten, konnten wir einen Kurs bestimmen. Die andere Technik beinhaltete die Verwendung eines im Flugzeug installierten Astrokompass. Hat geklappt, vorausgesetzt man hat nicht vergessen einen Almanach oder Tische mitzunehmen!!
Natürlich weiß ich heutzutage, GPS ... Das ist hauptsächlich zum Spaß und als Backup. Die wahre Peilungsberechnung der Sonne wird meiner Meinung nach kein Problem sein, da ich bereits die RA und die Deklination der Sonne berechnen kann (unter Verwendung von Meeus 'Adaption von VSOP87). Aber es wäre gut, vorher zu wissen, ob die Sonne tatsächlich aufgehen wird, wenn ich vorhabe, einen Astrokompass darauf zu richten. Und der Grund, warum ich in hohen Breiten nach der Antwort suche, ist, dass es schließlich "im Norden" verwendbar sein sollte.
Beantwortung meiner eigenen Frage hier. Ich habe ziemlich viel Zeit damit verbracht, dies zu untersuchen, und das ist, was ich herausgefunden habe.
Zur Interpretation der Anstiegs- und Untergangswerte: Wenn m1 oder m2 negativ werden, habe ich nur festgestellt, dass dies bedeutet, dass der Anstieg oder Untergang am vorherigen Zulu-Tag auftritt. Und wenn die Werte über 1 gehen, ist/sind es am folgenden Zulu-Tag aufgetreten.
Ich war überrascht, dass der Algorithmus für einige Daten und Orte eine Anstiegszeit angibt, die später als die eingestellte Zeit ist. Das schien zunächst zweifelhaft, scheint aber mittlerweile durchaus plausibel zu sein, wie andere Quellen bestätigen. Ich bin auch zu dem Schluss gekommen, dass dies in keiner Weise bedeutet, dass der Sonnenaufgang dann an einem anderen Datum als dem Sonnenuntergang liegen könnte - in Zulu-Begriffen beziehen sich Ergebnisse zwischen 0 und 1 auf das Datum, für das die Berechnung durchgeführt wird. unabhängig von ihrer Reihenfolge.
In Bezug auf die Wiederholung der Berechnung (von m1 oder m2), bis deltaM klein genug wird, habe ich festgestellt, dass, wenn die Wiederholungen 20 überschreiten, dies ein sehr guter Hinweis darauf ist, dass es an diesem Tag tatsächlich keinen Anstieg (oder Abfall) gibt. Dies scheint auch der Realität zu entsprechen, da es tatsächlich Tage gibt, an denen es nur einen Sonnenaufgang, aber keinen Sonnenuntergang gibt (oder umgekehrt). Dieser Aspekt der von mir geposteten Berechnungsimplementierung scheint also auch korrekt zu sein.
Weiter: Die Ergebnisse des Meeus-Algorithmus zeigen, dass sich die Abfolge von Auf- und Untergang im Laufe der Zeit umkehrt, weil ein Sonnenuntergang „fehlt“. Zum Beispiel könnten Sie Rise-Set, Rise-Set, Rise-(NO set) haben, und dann ändert sich die Reihenfolge in Set-Rise, Set-Rise usw. Soweit ich das beurteilen kann, entspricht dies auch der Realität. Und der Tag, an dem es keinen Sonnenuntergang gibt, kann erkannt werden, wie im vorherigen Absatz beschrieben.
Im wirklichen Leben gibt es jedoch Situationen, in denen an einem Zulu-Tag die Sonne aufgeht, lange aufgeht, dann untergeht und dann wieder aufgeht. Während Meeus' Buch nicht sagt, wie man die Zeit des zweiten Sonnenuntergangs bestimmt, kann ich zumindest feststellen, dass es passiert, weil ich für diesen Tag eine Aufgangs-Untergangs-Sequenz habe (ohne fehlenden Aufgang oder Untergang) und am nächsten Tag wird es eine Set-Rise-Sequenz sein.
Dieser besondere Umstand scheint ziemlich selten zu sein, und die Arbeit von Meeus bleibt sehr beeindruckend, da sie es mir ermöglicht hat, sie ohne Vorkenntnisse in Astronomie umzusetzen. Und es hat sicherlich meine Neugier geweckt – ich habe mir ein Teleskop gekauft, weil ich viel Zeit mit seinem Buch verbracht habe!
PS. In dem von mir geposteten Beispielcode war ein Fehler, den ich jetzt entfernt habe. Der neue Wert von m1 (oder m2), der durch Addition von deltaM erhalten wird, sollte nicht „renormiert“ werden, um zwischen 0 und 1 zu liegen, da dies zu völlig falschen Ergebnissen führt.
Planetenmacher
Ugo
Huy Pham
Ugo
genannt2voyage