[VB.NET] [Algorithm] 위도, 경도 값을 이용하여 일출 / 일몰 시각 계산하기

#Region "Math functions"
        Private Function degtorad(deg)
            Return Math.PI * deg / 180
        End Function

        Private Function radtodeg(rad)
            Return 180 * rad / Math.PI
        End Function

        Private Function sin(deg)
            Return Math.Sin(degtorad(deg))
        End Function

        Private Function cos(deg)
            Return Math.Cos(degtorad(deg))
        End Function

        Private Function tan(deg)
            Return Math.Tan(degtorad(deg))
        End Function

        Private Function asin(x)
            Return radtodeg(Math.Asin(x))
        End Function

        Private Function acos(x)
            Return radtodeg(Math.Acos(x))
        End Function
  
        Private Function atan(x)
            Return radtodeg(Math.Atan(x))
        End Function

        Private Function modpos(x, m)
            Return ((x Mod m) + m) Mod m
        End Function
#End Region

        ' 위도(lat)와 경도(lng)를 나타내는 상수를 선언하는 부분입니다.
        ' 위도와 경도는 지구 상의 특정 위치를 나타내는 좌표입니다.
        ' 위도는 적도로부터 북쪽 또는 남쪽으로 얼마나 떨어져 있는지를 의미하며,
        ' 경도는 그리니치 천문대로부터 동쪽 또는 서쪽으로 얼마나 떨어져 있는지를 의미합니다.

        Private Const lat As String = "[위도 값]"
        Private Const lng As String = "[경도 값]"

#Region "zenith variables"
        'official': 90.833333,
        '90deg 50'
        'civil': 96,
        'nautical': 102,
        'astronomical': 108
#End Region

        ' 천정각은 태양이 지평선과 이루는 각도를 말합니다.
        ' 천정각이 "90.833333" 으로 태양이 지평선에서 정확히 90.833333 도 높이에 있을 시점으로,
        ' 태양이 지평선 바로 위 정오에 가장 높이 떠 있을 시점을 의미합니다.
        Private Const zenith = 90.833333

        Private Function SolarCalculation(eventDate As Date, isRising As Boolean) As DateTimeOffset
            ' First calculate the day of the year.
            ' 먼저 올해의 날짜를 계산합니다.

            ' eventDate가 1 년 중 몇 번째 날인지를 계산하는 부분입니다.
            ' 계산한 결과는 변 N 에 저장됩니다. 태양의 위치를 계산하는 데 중요한 조건 중 하나입니다.
            Dim N = eventDate.DayOfYear

            ' Convert the longitude to hour value and calculate an approximate time.
            ' 경도 값을 시간 값으로 변환하고 대략적인 시간을 계산합니다.

            ' 경도를 시간 값으로 변환하고, 일출 또는 일몰의 대략적인 시간을 계산하는 부분입니다.
            ' 계산한 결과는 변수 t 에 저장합니다.
            ' 여기서 경도를 15 로 나누는 이유는 １ 시간  지구가 15 도 자전하기 때문입니다.
            ' 따라서 경도를 15 로 나누면 해당 위치의 시간대를 얻을 수 있습니다.
            Dim lngHour = lng / 15

            ' t 는 1 년 중 몇 번째 날인지 (N) 와 그 날의 특정 시간대를 합친 값입니다.
            ' 일출이면 6 (오전 6 시), 일몰이면 18 (오후 6 시) 을 반환합니다.
            Dim t = N + ((If(isRising, 6, 18) - lngHour) / 24)

            ' Calculate the Sun's mean anomaly.
            ' 태양의 평균 편차를 계산합니다.

            ' 계산한 결과는 M에 저장됩니다.
            ' 태양의 평균 편차는 태양이 지구 주위를 원형 궤도로 돌아가는 것처럼 가정했을 때의 위치를 뜻합니다.
            ' 실제로는 태양의 궤도는 완벽한 원형이 아니라 약간 타원형이므로, 결과 값은 실제 태양의 위치와 약간 차이가 있습니다.
            ' t 는 특정 날짜를 기준으로 한 시간입니다. 이 공식은 태양의 움직임을 근사화 한 것으로, 태양이 일 년 동안 하늘을 움직이는 경로를 모형화합니다.
            ' M 은 평균 경도를 나타냅니다. 0.9856 은 태양이 하루에 움직이는 평균 각도입니다. -3.289은 초기 평균 경도를 조정하는 상수 입니다.
            Dim M = (0.9856 * t) - 3.289

            ' Calculate the Sun's true longitude.
            ' 태양의 실제 경도를 계산합니다.

            ' 계산한 결과는 변수 L에 저장됩니다. 실제 경도는 태양의 실제 위치를 나타냅니다.
            ' 이 값은 태양의 평균 편차 (M) 와 태양의 궤도가 타원형인 것을 보정하는 항목들을 합친 것입니다.
            ' L 은 실제 경도를 나타냅니다.
            ' 1.916 * sin(M) 와 0.02 * sin(2 * M)은 태양의 궤도가 완전한 원형이 아니라 타원형이므로 발생하는 경도의 변화를 보정하기 위한 계산식입니다.
            ' 282.634 는 초기 실제 경도를 조정하는 상수입니다.
            Dim L = M + (1.916 * sin(M)) + (0.02 * sin(2 * M)) + 282.634

            ' NOTE: L potentially needs to be adjusted into the range [0, 360] by adding / subtracting 360.
            ' 참고: 경도는 잠재적으로 360 을 가산하거나 차감해서 [0, 360] 내의 범위로 조정해야 합니다.
            L = modpos(L, 360)

            ' Calculate the Sun's right ascension
            ' 태양의 적경을 계산합니다.

            ' 태양의 적경을 계산하는 부분입니다.
            ' 0.91764 는 태양의 궤도 기울기를 나타내는 상수 값입니다.
            ' 태양의 궤도가 지구의 적도에 대해 약 23.44 도로 기울어져 있음을 고려하여 도출된 값입니다.
            ' 기울기를 고려하지 않으면 태양의 위치를 정확하게 계산할 수 없습니다.
            Dim RA = atan(0.91764 * tan(L))

            ' RA potentially needs to be adjusted into the range [0, 360] by adding / subtracting 360.
            ' 적경은 잠재적으로 360 을 가산하거나 차감해서 [0, 360] 내의 범위로 조정해야 합니다.

            ' 적경을 0 도 이상 360 도 미만의 범위로 조정하는 부분입니다. 적경은 하늘에서 별이나 태양이 어디에 위치해 있는지를 나타내는 좌표 중 하나입니다.
            ' 적경은 보통 도 단위로 표현되며, 0 도에서 360 도 사이의 값으로 나타낼 수 있습니다.
            ' 해당 범위를 벗어나는 값은 혼동을 일으킬 수 있으므로, modpos 함수를 사용하여 적경을 이 범위 내로 조정합니다.
            ' 시계가 12 시를 넘어가면 다시 1 시로 돌아가는 것처럼, 적경도 360 도 이상이라면 이후부터는 0 도 부터 가산됩니다.
            RA = modpos(RA, 360)

            ' Right ascension value needs to be in the same quadrant as L
            ' 적경 값은 경도와 동일한 사분면에 있어야 합니다.

            ' 적경 (RA) 을 실제 경도 (L) 와 같은 사분면에 위치하도록 조정하는 부분입니다. 이를 이해하기 위해서는 먼저 사분면에 대한 이해가 필요합니다.
            ' 평면을 4 개의 영역으로 나눈 각 영역을 통틀어 사분면이라 합니다.
            ' 보통 원이나 원형의 도표에서 사용되며, 각각의 사분면은 90 도의 각도를 가집니다.
            ' 따라서 360 도를 4 로 나누면 90 도가 되는데, 하나의 사분면을 의미합니다.

            ' 실제 경도(L)와 적경(RA)이 속한 사분면을 계산하는 부분입니다.
            ' Math.Floor(L / 90) * 90는 실제 경도를 90으로 나눈 후 버림하여 정수 부분만 가져온 후 다시 90 을 곱합니다.
            ' 이를 통해 실제 경도가 속한 사분면의 시작점을 얻을 수 있습니다. 적경에 대해서도 동일한 연산을 수행합니다.
            Dim Lquad = (Math.Floor(L / 90)) * 90
            Dim RAquad = (Math.Floor(RA / 90)) * 90

            ' 적경을 실제 경도와 같은 사분면에 위치하도록 조정하는 부분입니다.
            ' Lquad - RAquad는 실제 경도와 적경이 속한 사분면 사이의 차이를 계산합니다.
            ' 계산한 차를 적경에 더함으로써 적경을 실제 경도와 같은 사분면에 위치하도록 조정합니다.
            ' 결과적으로 적경과 실제 경도가 같은 사분면에 위치하게 되므로, 태양의 위치를 보 정확하게 계산할 수 있습니다. 
            RA += (Lquad - RAquad)

            ' Right ascension value needs to be converted into hours.
            ' 적경 값을 시간으로 변환합니다.
            ' 적경은 원래 도 단위로 표현되지만, 여기서는 시간 단위로 변환하여 사용합니다.
            RA /= 15

            ' Calculate the Sun's declination.
            ' 태양의 적위를 계산합니다.
            ' 태양의 적위를 계산하는 부분입니다. 여기서 0.39782 는 태양의 궤도 기울기의 사인 (Sine) 값을 나타내는 상수 값이며,
            ' 태양의 궤도가 지구의 적도에 대해 약 23.44 도로 기울어져 있음을 고려한 것입니다.
            ' 태양의 기울기를 고려하지 않으면 위치를 정확하게 계산할 수 없습니다.
            Dim sinDec = 0.39782 * sin(L)
            Dim cosDec = cos(asin(sinDec))

            ' Calculate the Sun's local hour angle.
            ' 태양의 현지 시간 각도를 계산합니다.
            ' 만약 계산한 결과가 1보다 크거나 -1보다 작다면, 즉 계산 결과가 유효하지 않다면 Nothing (Null) 을 반환합니다.
            Dim cosH = (cos(zenith) - (sinDec * sin(lat))) / (cosDec * cos(lat))
            If cosH > 1 Or cosH < -1 Then Return Nothing

            ' Finish calculating local hour angle and convert into hours.
            ' 현지 시간 각도 계산을 마치고 시간으로 변환합니다.

            ' If(isRising, 360 - acos(cosH), acos(cosH)) 부분은 일출이면 360 - acos(cosH), 일몰이면 acos(cosH) 를 적용합니다.
            ' 그리고 그 결과를 15 로 나누어 시간으로 변환합니다.
            ' 변수 H 는 앞서 계산한 현지 시간 각도, RA 는 적경, t 는 일출 또는 일몰의 대략적인 시간을 의미합니다.
            Dim H = If(isRising, 360 - acos(cosH), acos(cosH)) / 15

            ' Calculate local mean time of rising / setting.
            ' 해가 뜨거나 지는 현지 평균 시간을 계산합니다.

            ' 태양의 위치를 계산하는 공식의 일부입니다.
            ' Tt 는 해가 뜨거나 지는 현지 평균 시간을 나타냅니다.
            ' H 는 현지 시간 각도, RA는 적경,
            ' t 는 일출 또는 일몰의 대략적인 시간을 나타냅니다.

            ' 0.06571 과 -6.622 는 태양의 움직임을 근사화하는 데 사용되는 상수 값입니다.
            ' 0.06571 은 태양의 움직임에 대한 일일 변화율을 나타내는 상수 값입니다. 이 값은 태양이 하루 동안 얼마나 움직이는지를 나타내는 값으로,
            ' 0.06571 이라는 값에 t 를 곱하면 특정 시간에 대한 태양의 움직임을 도출해 낼 수 있습니다.
            ' -6.622는 초기 시간 차를 의미하는 상수로 태양의 위치를 계산할 때 사용되는 시간 체계와 실제 시간 체계 사이의 차이를 보정하는 데 사용됩니다.
            ' 두 가지 상수는 태양의 움직임을 정확하게 모형화하는 데 필요한 값들로, 천문학자들이 오랜 기간 수집한 데이터를 바탕으로 결정된 값들입니다. 
            Dim Tt = H + RA - (0.06571 * t) - 6.622

            ' Adjust back to UTC.
            ' 협정 세계 시 기준으로 다시 조정합니다.
            Dim utc = Tt - lngHour

            ' 협정 세계 시를 [0, 24] 범위로 조정하는 부분입니다.
            ' NOTE: UTC potentially needs to be adjusted into the range [0, 24] by adding / subtracting 24.
            ' UTC는 잠재적으로 24 를 가산하거나 차감해서 [0, 24] 범위로 조정해야 합니다.
            utc = modpos(utc, 24)

            ' 협정 세계 시를 시간과 분으로 변환하는 부분입니다.
            Dim hours As Integer = Math.Floor(utc)
            Dim minutes As Integer = Math.Round(60 * (utc - hours))

            ' 분 값이 60 이상이라면 시간을 1 증가시키고 분을 0 으로 설정하는 부분입니다.
            ' 만약 시간이 24 이상이라면 날짜 (일자) 를 1 증가시키고 시간을 0 으로 설정하는 부분입니다.
            If minutes = 60 Then
                minutes = 0
                hours += 1
                If hours = 24 Then
                    hours = 0
                    eventDate.AddDays(1)
                End If
            End If

            ' 계산한 결과를 DateTimeOffset 타입으로 변환하고 반환하는 부분입니다. 만약 일출 시간을 계산했다면, 반환되는 시간에서 하루를 차감합니다.
            ' SolarCalculation 함수는 주어진 날짜에 대해 일출 또는 일몰 시간을 계산하는 역할을 합니다.
            Dim dateOut As New DateTimeOffset(eventDate.Year, eventDate.Month, eventDate.Day, hours, minutes, 0, New TimeSpan(0, 0, 0))
            Return If(isRising, dateOut.AddDays(-1), dateOut)
        End Function