[分享]彗星轨道根数算法文件（β版）

Attribute VB_Name = "TrackEquation"
'彗星轨道根数算法文件
Option Explicit
Private Const pi = 3.14159265358979   '给定π的值

Public Function ComputePosition(year As Long, month As Long, day As Double, q As Double, e As Double, w As Double, n As Double, i As Double, cType As Long) As Single

'定义数据库中轨道参数相应的变量
'定义名称，年，月，日，q，e，ω，Ω，i， H，G
'Dim name As String
'定义变量
Dim epsilon  As Double   '定义黄赤交角变量(角度)
Dim longradii As Double '定义长半径变量
epsilon = 23.439375     '给定2000.0的黄赤交角ε常量值
epsilon = epsilon * pi / 180
i = i * pi / 180    '将角度转换为弧度 (弧度＝角度×pi/180)
                   '(角度＝弧度×180/pi)
n = n * pi / 180
w = w * pi / 180
'在sin之类的函数计算中，sin(变量)一式中的变量应该为弧度值
'所有角度值应先转换为弧度值后进行计算

       '定义6个高斯常数A,B,C,a,b,c
       Dim aa As Double, bb As Double, cc As Double    '定义A,B,C变量
       Dim a As Double, b As Double, c As Double       '定义a,b,c变量
       longradii = q / (1 - e)         '计算长半径数值
       '---A--------------------------------------------------
       aa = Atn(-Cos(n) / Cos(i) / Sin(n)) '计算A的值(弧度)
       aa = aa * 180 / pi                  '弧度转换为角度
       '---B--------------------------------------------------
       bb = Atn(Sin(n) * Cos(epsilon) / (Cos(i) * Cos(epsilon) * Cos(n) - Sin(i) * Sin(epsilon))) '计算B的值(弧度)
       bb = bb * 180 / pi                  '弧度转换为角度
       '---C--------------------------------------------------
       cc = Atn(Sin(n) * Sin(epsilon) / (Cos(i) * Cos(n) * Sin(epsilon) + Sin(i) * Cos(epsilon)))
       cc = cc * 180 / pi                  '弧度转换为角度
       '_________________              ____________________
       '                 验算      开始
       Dim sina As Double, sinb As Double, sinc As Double '定义sina,sinb,sinc
           sina = Cos(n) / Sin(aa * pi / 180)
           sinb = Sin(n) * Cos(epsilon) / Sin(bb * pi / 180)
           sinc = Sin(n) * Sin(epsilon) / Sin(cc * pi / 180)
       '如果sa*sa+sb*sb+sc*sc-2的绝对值大于0.00001,说明数据计算有问题
           If Abs(sina * sina + sinb * sinb + sinc * sinc - 2) > 0.0000001 Then ComputePosition = Abs(sina * sina + sinb * sinb + sinc * sinc - 2)
       '_________________验算      完成____________________'验算公式二和三可以不必使用，通常用第一式即可
       
       '定义6个辅助量Px,Py,Pz,Qx,Qy,Qz
       Dim px As Double, py As Double, pz As Double
       Dim qx As Double, qy As Double, qz As Double
           px = sina * Sin(aa * pi / 180 + w)
           py = sinb * Sin(bb * pi / 180 + w)
           pz = sinc * Sin(cc * pi / 180 + w)
           qx = sina * Cos(aa * pi / 180 + w)
           qy = sinb * Cos(bb * pi / 180 + w)
           qz = sinc * Cos(cc * pi / 180 + w)
       '----验算-----
       If Abs(px * px + py * py + pz * pz - 1) > 0.0000001 Then ComputePosition = Abs(qx * qx + qy * qy + qz * qz - 1)
       '___________________________________________________
       
       Dim cosfi As Double     '定义cosφ变量
       cosfi = Sqr(1 - e * e) '计算出cosφ的数值
       '定义6个大写坐标Ax,Ay,Az,Bx,By,Bz
       Dim ax As Double, ay As Double, az As Double
       Dim bx As Double, by As Double, bz As Double
           longradii = q / (1 - e)     '计算出轨道半长径
           '下面分别计算Ax,Ay,Az,Bx,By,Bz的数值
           ax = longradii * px
           ay = longradii * py
           az = longradii * pz
           bx = longradii * cosfi * qx
           by = longradii * cosfi * qy
           bz = longradii * cosfi * qz
       Dim GS As Double    '定义天体平均运动速度 n 变量 其值为sqr(GS/longradii的3次方)
       '计算出的GS的值是 x.xxxxxx度/日
       GS = Sqr(0.000295912208 / (longradii * longradii * longradii)) * 180 / pi
'_____________________________________________________________________________
       
       Dim M As Double     '定义平近点角M变量 该变量为开普勒方程的一个常量
           '需要判断已知日期同计算日期之间的时间差
           M = days360("1998-11-21", "1999-1-12") * GS
           M = 7.71
'此处应计算两个日期之间的时间差，精确到小数点后5位
'_____________________________________________________________________________

       '通过渐进法(牛顿迭代法)计算开普勒方程并找出其解来
       Dim ee As Double    '定义偏近点角E变量
       Dim e0 As Double    '定义中间变量E0
       '将M转换为弧度值，如果M仍旧是角度值时，需要根据上面的计算而知
       M = M * pi / 180
       e0 = M '将M赋值给E0
count:      ee = M + e * Sin(e0)
               If ee <> e0 Then e0 = ee: GoTo count
               If ee = e0 Then ee = e0
       '变量E通过计算并已经给出数据
       'ee(数据是弧度数)
       '定义cosE，sinE，cosE-e
       Dim cose As Double, sine As Double, cosee As Double
           cose = Cos(ee)
           sine = Sin(ee)
           cosee = Cos(ee) - e
       Dim R As Double     '定义彗日距离，单位为A.U.
           R = longradii * (1 - e * cose)
       '定义日心赤道坐标(x',y',z')
       Dim lx As Double, ly As Double, lz As Double
           lx = ax * cosee + bx * sine
           ly = ay * cosee + by * sine
           lz = az * cosee + bz * sine
       '定义X,Y,Z太阳直角坐标变量 需要从天文年历中提取数据
       '计算公式
       Dim xx As Double, yy As Double, zz As Double
           xx = 0.35782
           yy = -0.84049
           zz = -0.3644
           
       '定义赤经,赤纬(α,δ)变量
       Dim alpha As Double, delta As Double
           alpha = Atn((ly + yy) / (lx + xx))
           delta = Atn((lz + zz) * Sin(alpha) / (ly + yy))
Select Case cType
   Case 1
     ComputePosition = alpha * 180 / pi    '角度值
   Case 2
     ComputePosition = delta * 180 / pi    '角度值
End Select
End Function


















Public Function days360(dt1 As Date, dt2 As Date) As Long   '计算两个日期间的天数
   
   Dim z1 As Long, z2 As Long
   Dim d1 As Long, d2 As Long
   Dim m1 As Long, m2 As Long
   Dim y1 As Long, y2 As Long
   
   d1 = day(dt1)
   m1 = month(dt1)
   y1 = year(dt1)
   
   d2 = day(dt2)
   m2 = month(dt2)
   y2 = year(dt2)
   
   If d1 = 31 Then
       z1 = 30
   Else
       z1 = d1
   End If
   
   If d2 = 31 And d1 >= 30 Then
       z2 = 30
   Else
       z2 = d2
   End If

   days360 = (y2 - y1) * 360 + (m2 - m1) * 30 + (z2 - z1)

End Function

来挖坟，怎么不编出来直接让人下载啊？

谁能把瑞士星历表的函数 说明汉化一下
http://www.astro.com/swisseph/swephprg.htm