[转贴] 回2EASY同好：《彗星轨道根数算法文件》

今天有时间，看了你的程序，很好啊！“*”是我和你的不同想法。多多交流：）
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 epsilon(year As Long) As Double
  Dim T1,T2,T3,T as Double
  T1=23.45022944
  T2=0.0130125
  T3=0.0000050277
  Y_year=year-1900   '我从1900年算起
  T=Y_year/100         '以百年为单位
  epsilon=(T1+T2+T3) * T
End Function
这是95年8月用BC写的，现在的是VB的，你明白的，特提供给你参考参考。


** 太阳直角坐标变量，如果可以以天文年历中提供的数据（几年以上，最好是50年啦）来作个分析，
   找出其变动的规律（摄动也如此，很烦！），编出来的程序就不要老是查天文年历了。Skymap Pro怎么做的？不得而知！
   希望你编出个超级程序来，我为你祈祷，阿门……

看不明白！ :shock:  

有谁有计算彗星轨道参数的软件呢？

[quote:e7bd8266ba="1290982"]看不明白！ :shock:  

有谁有计算彗星轨道参数的软件呢？[/quote]SNP、SKYMAP都可以

SKYMAP都可以

怎样用？能举个例子吗？

它里面不是有彗星的资料么

VB的代码哦