проблемы с поиском набора восходов и времени прохождения в высоких широтах с использованием алгоритмов Миуса

Я реализую алгоритм Жана Меуса из его книги «Астрономические алгоритмы», 2-е изд. для времени восхода, захода и прохождения в Swift. Кажется, он дает очень хорошие результаты (разница не более 1 минуты) по сравнению, например, с timeanddate.com, в том числе на высоких широтах (выше 66 градусов). Однако иногда это не срабатывает.

Я проследил это до обстоятельств, когда а) солнце восходит в определенный день UTC, но не заходит, или наоборот, б) когда солнце восходит дважды в данный день UTC. Конкретный пример последнего происходит 16 апреля 2020 года в Лонгйире (Шпицберген), Норвегия. И, наконец, c) когда время нарастания «перемещается» назад с «после» на «до полуночи», если смотреть на него в последовательные дни.

Приблизительное время подъема, установления и прохождения, соответственно, m1, m2 и m0, может быть достаточно легко рассчитано. Затем их необходимо уточнить с помощью значения дельта m. Это может/может/должно быть сделано итеративно.

Присмотревшись, я обнаружил, что для событий, описанных выше, значения m1 и m2, скорректированные на дельту m, как описано на стр. 103, обычно превышают 1. Миус утверждает, что значения m «должны быть между 0 и 1. Если одно или если они находятся за пределами этого диапазона, прибавьте или вычтите 1 дюйм. Чуть дальше, в примечании 1 в конце главы 15, есть дразнящее примечание, в котором говорится, что «если время установления [..] необходимо в местном времени, расчет должен выполняться с использованием [...] m2 = 1,12113, что больше 1.

Это заставляет меня подозревать — я не астроном, как вы уже догадались, — что значения m1 или m2 больше 1, вероятно, могут помочь мне рассчитать время нарастания в день, когда это происходит дважды, а также правильное время нарастания на день, когда солнце не заходит (и наоборот).

Это привело меня на github, где я нашел некоторый код JavaScript (и не только) от onekiloparsec , который сравнивает время нарастания со временем прохождения, и если первое позже второго, время нарастания принимается за предыдущий день. Точно так же, если установленное время предшествует времени транзита, оно считается следующим днем.

Я также просмотрел «Практическую астрономию с помощью вашего калькулятора или электронной таблицы» Питера Даффет-Смита и Джонатана Цварта, но не нашел там ответа. Это дало очень полезную информацию, которую Меус не предоставляет, а именно, что знак результата уравнения 15.1 Миуса, когда результат по абсолютной величине больше 1, позволяет определить, находится ли звезда постоянно внизу (поскольку H0 >1) или над горизонтом (cos H0 < -1).

Было бы здорово получить объяснение или ссылку и некоторые подробности о том, как интерпретировать результаты для m1 и m2, как, к сожалению, в книге Миуса, при описании алгоритмов до уровня, на котором неастроном может выполнить множество вычислений, оставляет меня со странным вопросом.

Ниже приведен (Swift) код, который я написал для итеративного уточнения значений m0, m1 и m2.

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

Для справки: я пилот, а не астроном, и в прошлом мне приходилось изучать методы навигации в полярных районах, где магнитный компас ненадежен. Таким образом, мы должны были знать, как использовать солнце в качестве компаса — найдя его истинный пеленг, а затем направив на него нос самолета, мы могли установить курс. Другой метод заключался в использовании астрокомпаса, установленного на самолете. Это сработало, если вы не забыли взять с собой альманах или таблицы!!

Конечно, я знаю, что в наши дни GPS... Это в основном для развлечения и в качестве резервной копии. Расчет истинного пеленга солнца, я думаю, не будет проблемой, так как я уже могу рассчитать RA и склонение солнца (используя адаптацию Meeus VSOP87). Но было бы хорошо знать заранее, действительно ли солнце взойдет, когда я планирую навести на него астрокомпас. И причина, по которой я ищу ответ в высоких широтах, заключается в том, что, в конце концов, его можно использовать «на севере».