30 เมษายน 2559

สูตรการคำนวณตำแหน่งดาวเคราะห์แบบดาราศาสตร์ : พุธ

สูตรการคำนวณตำแหน่งดาวเคราะห์แบบดาราศาสตร์ : พุธ
ฟังก์ชั่น mercury คือฟังก์ชั่น VBA บน Excel ใช้ในการหาค่าตำแหน่งดาวพุธแบบดาราศาสตร์สากล หรือระบบสายนะวิธี มีตัวแปรที่รับเข้าดังต่อไปนี้
day = วันที่
month = เดือน
year = ปี ค.ศ.
hr = ชั่วโมง
mn = นาที
และฟังก์ชั่นเสริมที่ใช้ประกอบการคำนวณคือ
Modulo คือ การหารแบบเอาเศษ
Atan2 คือ ค่าอาร์กแทนเจนต์ หรือแทนเจนต์ผกผัน ของพิกัด x และ y ที่ระบุ
Asin คือ ค่า invers ของ sin
Acos คือ ค่า invers ของ cos

สูตรการคำนวณมีดังนี้
Function mercury(day, month, year, hr, mn As Double) As Double
Dim dj2000, days, Ma, Ta, L, RV As Double
dj2000 = 367 * year - Int(7 * (year + Int((month + 9) / 12)) / 4) + Int(275 * month / 9) + day - 730531.5 + ((hr + (mn / 60)) / 24)
days = dj2000 + 864.5
'องค์ประกอบของดาวเคราะห์
Ma = modulo((0.071425034 * days + 5.487728637 - 1.351827319), 6.283185307) 'Mean anomaly
Ta = Ma + (2 * 0.2056324 - 0.2056324 ^ 3 / 4 + 5 / 96 * 0.2056324 ^ 5) * Sin(Ma) + (5 * 0.2056324 ^ 2 / 4 - 11 / 24 * 0.2056324 ^ 4) * Sin(2 * Ma) + (13 * 0.2056324 ^ 3 / 12 - 43 / 64 * 0.2056324 ^ 5) * Sin(3 * Ma) + 103 / 96 * 0.2056324 ^ 4 * Sin(4 * Ma) + 1097 / 960 * 0.2056324 ^ 5 * Sin(5 * Ma) 'True anomaly
L = modulo(Ta + 1.351827319, 6.283185307) 'Longitude Planet
RV = 0.3870978 * (1 - 0.2056324 ^ 2) / (1 + 0.2056324 * Cos(Ta)) 'Radius Vector
'องค์ประกอบของโลก
Dim Mae, Tae, Le, RVe, Hlong, Hlat, Dis, lambda, beta, alpha, delta As Double
Mae = modulo(0.017201609 * days + 5.731722874 - 1.795100806, 6.283185307)  'Mean anomaly
Tae = Mae + (2 * 0.0166967 - 0.0166967 ^ 3 / 4) * Sin(Mae) + 5 * 0.0166967 ^ 2 / 4 * Sin(2 * Mae) + 13 * 0.0166967 ^ 3 / 12 * Sin(3 * Mae) 'True anomaly
Le = modulo(Tae + 1.795100806, 6.283185307) 'Longitude Earth
RVe = 1 * (1 - 0.0166967 ^ 2) / (1 + 0.0166967 * Cos(Tae)) 'Radius Vector
Hlong = Atan2(Cos(L - 0.843585695), Sin(L - 0.843585695) * Cos(0.122261536)) + 0.843585695 'Helio long
Hlat = Asin(Sin(L - 0.843585695) * Sin(0.122261536)) 'Helio lat (phi)
Dis = RV * Cos(Hlat) 'dist from Sun (au)
lambda = modulo(Atn(Dis * Sin(Le - Hlong) / (RVe - Dis * Cos(Le - Hlong))) + 3.141592654 + Le, 6.283185307)
mercury = lambda * (180 / 3.141592654)
End Function

Function modulo(a, b As Double) As Double
modulo = a / b
modulo = (modulo - Int(modulo)) * b
End Function
Function Atan2(ByVal x As Double, ByVal y As Double) As Double
    If y > 0 Then
      If x >= y Then
        Atan2 = Atn(y / x)
      ElseIf x <= -y Then
        Atan2 = Atn(y / x) + 3.141592654
      Else
        Atan2 = 3.141592654 / 2 - Atn(x / y)
      End If
    Else
      If x >= -y Then
        Atan2 = Atn(y / x)
      ElseIf x <= y Then
        Atan2 = Atn(y / x) - 3.141592654
      Else
        Atan2 = -Atn(x / y) - 3.141592654 / 2
      End If
    End If
End Function
Function Asin(x As Double) As Double ' ค่า invers ของ sin
    Asin = Atn(x / (Sqr(1 - x * x)))
End Function
Function Acos(x As Double) As Double ' ค่า invers ของ cos
    Acos = Atn((Sqr(1 - x * x)) / x)
End Function
ตัวอย่างการคำนวณ 1 พฤษภาคม 2016 เวลา 12.00 น.(กรีนิช)
รูปแบบสูตร mercury(1,5,2016,12,0)
ค่าที่ได้คือ 52.67609395 องศา
(เปรียบเทียบกับค่าที่ได้จากโมดูล swiss คือ 53.265609951987 ต่างกันอยู่ 0.589515998 องศา)
ในกรณีที่จะทำเป็นค่าตำแหน่งดาวในระบบนิรายนะวิธี ให้นำค่าอายนางศมาลบออกจากค่าที่คำนวณได้ กรณีนี้ใช้ค่าอายนางศแบบลาหิรีที่คำนวณได้ในวันเวลาตามตัวอย่าง คือ 24.0852731799126
ตำแหน่งดาวพุธ(นิรายนะวิธี)  52.67609395 - 24.0852731799126 = 28.59082077
1.ราศี 28.59082077 หาร 30 ลัพธ์เป็นราศี = 0
2.องศา จำนวนเต็มของเศษจากการหาร = 28
3.ลิปดา ทศนิยมที่เหลือคูณด้วย 60 = 0.59082077 X 60 = 35.4492462 = 35
4.ฟิลิปดา ทศนิยมที่เหลือในข้อ 3. คูณด้วย 60 = 0.4492462 X 60 = 26.954772 = 26
สรุป ตำแหน่งดาวพุธ(นิรายนะ) = ราศีเมษ 28 องศา 35 ลิปดา 26 ฟิลิปดา (เวลาประเทศไทย 19.00 น.)
เปรียบเทียบผลกับโมดูล swiss คือ ราศีเมษ 29 องศา 10 ลิปดา 49 ฟิลิปดา ต่างกันอยู่ 35 ลิปดา 22 ฟิลิปดา

ธีรพร  เพชรกำแพง
30 เมษายน 2559
(ต้องการขอรับไฟล์งาน Excel ติดต่อได้ที่ tepar2009@gmail.com)