太陽の高さの計算 - wetchのブログ

参考文献
仕様
標準モジュール：main
標準モジュール：数学関数
クラスモジュール：LocalTime
クラスモジュール：測地座標系
クラスモジュール：地平座標系
クラスモジュール：時角赤緯座標系
クラスモジュール：赤道座標系
クラスモジュール：黄道座標系
考え事

参考文献

仕様

角度関係はすべてラジアン単位で表す．
経度は東経をプラスとし， $0\sim 2\pi$ の間の値をとる．
緯度は北緯をプラスとし， $-\pi/2\sim+\pi/2$ の間の値をとる．

標準モジュール：main

year , month , day : 年月日
hour , minute , second : 時分秒
時差 : 時間単位で入力。日本なら+9．

Function 方位角(経度 As Double, 緯度 As Double, year As Long, month As Long, day As Long, hour As Long, minute As Long, second As Long, Optional 時差 As Long = 9)
    方位角 = main(経度, 緯度, year, month, day, hour, minute, second, 時差).方位角
End Function

Function 高度(経度 As Double, 緯度 As Double, year As Long, month As Long, day As Long, hour As Long, minute As Long, second As Long, Optional 時差 As Long = 9)
    Dim result As Double
    result = main(経度, 緯度, year, month, day, hour, minute, second, 時差).高度
    
    '大気差補正
    高度 = result + WorksheetFunction.Radians(0.0167) / Tan(result + WorksheetFunction.Radians(8.6) / (WorksheetFunction.Degrees(result) + 4.4))
End Function

Private Function main(経度 As Double, 緯度 As Double, year As Long, month As Long, day As Long, hour As Long, minute As Long, second As Long, 時差 As Long) As 地平座標系
    Dim 観測点 As 測地座標系
    Set 観測点 = New 測地座標系
    観測点.経度 = 経度
    観測点.緯度 = 緯度
    
    Dim t As LocalTime
    Set t = New LocalTime
    t.Set年月日時分秒 year, month, day, hour, minute, second
    t.Set時差 時差
     
    Set main = 太陽座標(t).赤道座標系へ変換(傾角(t)) _
                          .時角赤緯座標系へ変換(恒星時(t, 経度)) _
                          .地平座標系へ変換(観測点)
End Function

Private Function 太陽座標(nw As localtime) As 黄道座標系
    Dim t As Double
    t = nw.経過ユリウス年
    
    Dim 太陽黄経 As Double
    太陽黄経 = WorksheetFunction.Radians(280.4603 + 360.00769 * t _
            + (1.9146 - 0.00005 * t) * Sin(WorksheetFunction.Radians(357.538 + 359.991 * t)) _
            + 0.02 * Sin(WorksheetFunction.Radians(355.05 + 719.981 * t)) _
            + 0.0048 * Sin(WorksheetFunction.Radians(234.95 + 19.341 * t)) _
            + 0.002 * Sin(WorksheetFunction.Radians(247.1 + 329.64 * t)) _
            + 0.0018 * Sin(WorksheetFunction.Radians(297.8 + 4452.67 * t)) _
            + 0.0018 * Sin(WorksheetFunction.Radians(251.3 + 0.2 * t)) _
            + 0.0015 * Sin(WorksheetFunction.Radians(343.2 + 450.37 * t)) _
            + 0.0013 * Sin(WorksheetFunction.Radians(81.4 + 225.18 * t)) _
            + 0.0008 * Sin(WorksheetFunction.Radians(132.5 + 659.29 * t)) _
            + 0.0007 * Sin(WorksheetFunction.Radians(153.3 + 90.38 * t)) _
            + 0.0007 * Sin(WorksheetFunction.Radians(206.8 + 30.35 * t)) _
            + 0.0006 * Sin(WorksheetFunction.Radians(29.8 + 337.18 * t)) _
            + 0.0005 * Sin(WorksheetFunction.Radians(207.4 + 1.5 * t)) _
            + 0.0005 * Sin(WorksheetFunction.Radians(291.2 + 22.81 * t)) _
            + 0.0004 * Sin(WorksheetFunction.Radians(234.9 + 315.56 * t)) _
            + 0.0004 * Sin(WorksheetFunction.Radians(157.3 + 299.3 * t)) _
            + 0.0004 * Sin(WorksheetFunction.Radians(21.1 + 720.02 * t)) _
            + 0.0003 * Sin(WorksheetFunction.Radians(352.5 + 1079.97 * t)) _
            + 0.0003 * Sin(WorksheetFunction.Radians(329.7 + 44.43 * t)))
            
    Dim result As 黄道座標系
    Set result = New 黄道座標系
    result.黄経 = 太陽黄経
    result.黄緯 = 0
    Set 太陽座標 = result
End Function

Private Function 太陽距離(localtime As LocalTime) As Double 
    '[天文単位]
    Dim t As Double, q As Double
    t = localtime.経過ユリウス年
    
    q = (0.007256 - 0.0000002 * t) * Sin(WorksheetFunction.Radians(267.54 + 359.991 * t)) _
        + 0.000091 * Sin(WorksheetFunction.Radians(265.1 + 719.98 * t)) _
        + 0.00003 * Sin(WorksheetFunction.Radians(90)) _
        + 0.000013 * Sin(WorksheetFunction.Radians(27.8 + 4452.67 * t)) _
        + 0.000007 * Sin(WorksheetFunction.Radians(254 + 450.4 * t)) _
        + 0.000007 * Sin(WorksheetFunction.Radians(156 + 329.6 * t))
            
    太陽距離 = 10 ^ q
End Function

Private Function 傾角(localtime As LocalTime) As Double
    Dim t As Double
    t = localtime.経過ユリウス年
    傾角 = WorksheetFunction.Radians(23.439291 - 0.000130042 * t)
End Function

Private Function 恒星時(localtime As LocalTime, 経度 As Double) As Double
    '恒星時 = 時角 + 赤経
    Dim t As Double, Theta_G As Double, Theta_0 As Double
    
    t = localtime.経過ユリウス年
    Theta_G = WorksheetFunction.Radians(100.4606 + 360.007700536 * t + 0.00000003879 * t ^ 2) 'グリニッジ恒星時
    Theta_0 = Theta_G _
                + 2 * Pi * CDbl(localtime.時分秒) _
                - 2 * Pi * CDbl(localtime.時差) '時差のあるところの地方時を修正したグリニッジ恒星時
    
    恒星時 = Modulo(Theta_0 + 経度, 2 * Pi, 0)
End Function

標準モジュール：数学関数

Function Pi() As Double
    '面倒なので関数化．
    Pi = WorksheetFunction.Pi
End Function

Function Modulo(x As Double, max As Double, Optional min As Double = 0) As Double
    '2πを法とすることが多いので関数にした．
    Dim result As Double, dx As Double
    result = x
    dx = max - min
    
    Do Until result >= min
        result = result + dx
    Loop
    
    Do Until result < max
        result = result - dx
    Loop
    
    Modulo = result
End Function

'簡単ベクトル回転
Function Rx(x As Double, y As Double, z As Double, theta As Double) As Double()
    'x軸周りの回転
    Dim result(0 To 2) As Double
    result(0) = x
    result(1) = y * Cos(theta) - z * Sin(theta)
    result(2) = y * Sin(theta) + z * Cos(theta)
    Rx = result
End Function

Function Ry(x As Double, y As Double, z As Double, theta As Double) As Double()
    Dim result(0 To 2) As Double
    result(0) = x * Cos(theta) + z * Sin(theta)
    result(1) = y
    result(2) = -x * Sin(theta) + z * Cos(theta)
    Ry = result
End Function

Function 方向余弦から極座標への変換(x As Double, y As Double, z As Double) As Double()
'方向余弦(x,y,z)を経度緯度系に変換．
'    x = Cos(経度) * Cos(緯度)
'    y = Sin(経度) * Cos(緯度)
'    z = Sin(緯度)
'    y/x = tan(経度)
    
    Dim 経度 As Double, 緯度 As Double
    経度 = Modulo(WorksheetFunction.Atan2(x, y), 2 * Pi, 0)
    緯度 = WorksheetFunction.Asin(z)
    
    Dim result(0 To 1) As Double
    result(0) = 経度
    result(1) = 緯度
    方向余弦から極座標への変換 = result
End Function

クラスモジュール：LocalTime

現地時刻を保持

Private 年月日時分秒_ As Date
Private 時差_ As Date

Property Get 時分秒() As Date
    時分秒 = Modulo(CDbl(年月日時分秒_), 1, 0)
End Property

Sub Set年月日時分秒(y As Long, mo As Long, d As Long, h As Long, mi As Long, s As Long)
    年月日時分秒_ = DateSerial(y, mo, d) + TimeSerial(h, mi, s)
End Sub

Property Get 時差() As Date
    時差 = 時差_
End Property

Sub Set時差(t As Long)  '[時]で入力
    時差_ = TimeSerial(t, 0, 0)
End Sub

Property Get 世界時() As Date
    世界時 = 年月日時分秒_ - 時差_
End Property

Property Get 経過日数() As Date
    Dim j2000 As Date
    j2000 = DateSerial(2000, 1, 1) + TimeSerial(12, 0, 0)
    
    経過日数 = 世界時 - j2000
End Property

Property Get 経過ユリウス年() As Double
    Dim deltaT As Double
    deltaT = 0.707017544 * year(年月日時分秒_) - 1349.54386 '地球自転の遅れ[秒]
    経過ユリウス年 = (CDbl(経過日数) + deltaT / 86400) / 365.25
End Property

クラスモジュール：測地座標系

Private 経度_ As Double
Private 緯度_ As Double

Property Get 経度() As Double
    経度 = 経度_
End Property

Property Let 経度(l As Double)
    経度_ = Modulo(l, 2 * Pi, 0)
End Property

Property Get 緯度() As Double
    緯度 = 緯度_
End Property

Property Let 緯度(p As Double)
    緯度_ = Modulo(p, Pi / 2, -Pi / 2)
End Property

クラスモジュール：地平座標系

Private 方位角_ As Double '真北を0として，北から見て時計回り
Private 高度_ As Double '地平線を0，直上をpi/2とする．

Property Get 方位角() As Double
    方位角 = 方位角_
End Property

Property Let 方位角(l As Double)
    方位角_ = Modulo(l, 2 * Pi, 0)
End Property

Property Get 高度() As Double
    高度 = 高度_
End Property

Property Let 高度(p As Double)
    高度_ = Modulo(p, Pi / 2, -Pi / 2)
End Property

'方向余弦
Property Get x() As Double '北向きプラス
    x = Cos(方位角_) * Cos(高度_)
End Property

Property Get y() As Double '東向きプラス
    y = Sin(方位角_) * Cos(高度_)
End Property

Property Get z() As Double
    z = Sin(高度_)
End Property

Sub 方向余弦(x As Double, y As Double, z As Double)
    Dim result() As Double
    result = 方向余弦から極座標への変換(x, y, z)
    方位角_ = result(0)
    高度_ = result(1)
End Sub

Function 時角赤緯座標系へ変換(観測点 As 測地座標系) As 時角赤緯座標系
    Dim 回転後の方向余弦() As Double
    回転後の方向余弦 = Ry(-x, -y, z, Pi / 2 - 観測点.経度)
 
    Dim result As 時角赤緯座標系
    Set result = New 時角赤緯座標系
    result.方向余弦 回転後の方向余弦(0), 回転後の方向余弦(1), 回転後の方向余弦(2)
    Set 時角赤緯座標系へ変換 = result
End Function

クラスモジュール：時角赤緯座標系

地球の自転軸・赤道面を基準とした座標．地球に固定されている．

Private 時角_ As Double '天の赤道上で最も南に位置する点を0とし，北から見て時計回り
Private 赤緯_ As Double '天の赤道を0，北極星を+π/2とする

Property Get 時角() As Double
    時角 = 時角_
End Property

Property Let 時角(t As Double)
    時角_ = Modulo(t, 2 * Pi, 0)
End Property

Property Get 赤緯() As Double
    赤緯 = 赤緯_
End Property

Property Let 赤緯(d As Double)
    赤緯_ = Modulo(d, Pi / 2, -Pi / 2)
End Property

'方向余弦
Property Get x() As Double '南がプラス
    x = Cos(時角_) * Cos(赤緯_)
End Property

Property Get y() As Double '西がプラス
    y = Sin(時角_) * Cos(赤緯_)
End Property

Property Get z() As Double '天の赤道上で0，北極星で+π/2
    z = Sin(赤緯_)
End Property

Property Let 方向余弦(x As Double, y As Double, z As Double)
    Dim result() As Double
    result = 方向余弦から極座標への変換(x, y, z)
    時角_ = result(0)
    赤緯_ = result(1)
End Property

Function 赤道座標系へ変換(恒星時 As Double) As 赤道座標系
    Dim 赤経 As Double
    赤経 = 恒星時 - 時角_
    
    Dim result As 赤道座標系
    Set result = New 赤道座標系
    result.赤経 = 赤経
    result.赤緯 = 赤緯_
    Set 赤道座標系へ変換 = result
End Function

Function 地平座標系へ変換(観測点 As 測地座標系) As 地平座標系
    Dim 回転後の方向余弦() As Double
    回転後の方向余弦 = Ry(-x, -y, z, Pi / 2 - 観測点.緯度)
    
    Dim result As 地平座標系
    Set result = New 地平座標系
    result.方向余弦 回転後の方向余弦(0), 回転後の方向余弦(1), 回転後の方向余弦(2)
    Set 地平座標系へ変換 = result
End Function

クラスモジュール：赤道座標系

地球の自転軸・赤道面を基準とした座標．天球上の天体に対し固定．太陽系外の天体を表すのに用いられる．

Private 赤経_ As Double '.春分点（うお座方向）を0とし，北極星から見て反時計回り
Private 赤緯_ As Double '天の赤道を0，北極星を+π/2とする

Property Get 赤経() As Double
    赤経 = 赤経_
End Property

Property Let 赤経(x As Double)
    赤経_ = Modulo(x, 2 * Pi, 0)
End Property

Property Get 赤緯() As Double
    赤緯 = 赤緯_
End Property

Property Let 赤緯(x As Double)
    赤緯_ = Modulo(x, Pi / 2, -Pi / 2)
End Property

'方向余弦
Property Get x() As Double '春分点が+1
    x = Cos(赤経_) * Cos(赤緯_)
End Property

Property Get y() As Double
    y = Sin(赤経_) * Cos(赤緯_)
End Property

Property Get z() As Double
    z = Sin(赤緯_)
End Property

Sub 方向余弦(x As Double, y As Double, z As Double)
    Dim result() As Double
    result = 方向余弦から極座標への変換(x, y, z)
    赤経_ = result(0)
    赤緯_ = result(1)
End Sub

Function 黄道座標系へ変換(傾角 As Double) As 黄道座標系
    Dim 回転後の方向余弦() As Double
    回転後の方向余弦 = Rx(x, y, z, -傾角)
    
    Dim result As 黄道座標系
    Set result = New 黄道座標系
    result.方向余弦 回転後の方向余弦(0), 回転後の方向余弦(1), 回転後の方向余弦(2)
    Set 黄道座標系へ変換 = result
End Function

Function 時角赤緯座標系へ変換(恒星時 As Double) As 時角赤緯座標系
    Dim 時角 As Double
    時角 = 恒星時 - 赤経_
    
    Dim result As 時角赤緯座標系
    Set result = New 時角赤緯座標系
    result.時角 = 時角
    result.赤緯 = 赤緯_
    Set 時角赤緯座標系へ変換 = result
End Function

クラスモジュール：黄道座標系

地球の公転面を基準とした座標．天球上の天体に対し固定．太陽系内の天体を表すのに用いられる．

Private 黄経_ As Double '春分点を0とし北から見て反時計回り
Private 黄緯_ As Double '黄道を0とし，黄道の北極をプラスとする

Property Get 黄経() As Double
    黄経 = 黄経_
End Property

Property Let 黄経(x As Double)
    黄経_ = Modulo(x, 2 * Pi, 0)
End Property

Property Get 黄緯() As Double
    黄緯 = 黄緯_
End Property

Property Let 黄緯(x As Double)
    黄緯_ = Modulo(x, Pi / 2, -Pi / 2)
End Property

'方向余弦
Property Get x() As Double '春分点が+1
    x = Cos(黄経_) * Cos(黄緯_)
End Property

Property Get y() As Double
    y = Sin(黄経_) * Cos(黄緯_)
End Property

Property Get z() As Double
    z = Sin(黄緯_)
End Property

Property Let 方向余弦(x As Double, y As Double, z As Double)
    Dim result() As Double
    result = 方向余弦から極座標への変換(x, y, z)
    黄径_ = result(0)
    黄緯_ = result(1)
End Property

Function 赤道座標系へ変換(傾角 As Double) As 赤道座標系
    Dim 回転後の方向余弦() As Double
    回転後の方向余弦 = Rx(x, y, z, 傾角)
   
    Dim result As 赤道座標系
    Set result = New 赤道座標系
    result.方向余弦 回転後の方向余弦(0), 回転後の方向余弦(1), 回転後の方向余弦(2)
    Set 赤道座標系へ変換 = result
End Function

考え事

各座標系毎にクラスなんか作らなくてもよかった．統一してVector3クラスにしたほうが却ってすっきりしたかも．