牛顿迭代法求“二分二至“精确时刻

zopey

由VSOP87理论计算出来的几何位置黄经，经过坐标转换，章动修正和光行差修正后，就可以得到比较准确的太阳地心视黄经。
"二分二至"是指春分、秋分、夏至、冬至。历法计算需要根据太阳地心视黄经反求出此时的时间，要做的事情就是求解行星运动轨道方程的根。
牛顿迭代法是一种在实数域和复数域上近似求解方程的方法，因此用到牛顿迭代法，每迭代一次，结果的有效数字将增加一倍。



节气与易卦.zip (54.31 KB, 下载次数: 146)

2019-1-18 14:56 上传

点击文件名下载附件

zopey

1. Function getjq_24(yy, mm)
   
2. W = (mm + (yy - 1999) * 24) * 15 * pai / 180
   
3. jd1 = earth_t(W - pai) * 365250 + 2451545#
   
4. 
   
5. Do
   
6.    jd0 = jd1
   
7.    stDegree = earth_L(jd0) - W
   
8.    stDegreep = (earth_L(jd0 + 0.000005) - earth_L(jd0 - 0.000005)) / 0.00001
   
9.    jd1 = jd0 - stDegree / stDegreep
   
10. Loop Until Abs(jd1 - jd0) < 0.0000001
    
11. 
    
12. getjq_24 = JD2GL(jd1 + 8 / 24 - deltatT(yy) / 86400)
    
13. End Function
    
14. 
    
15. Function getjq_64(yy, mm)
    
16. W = (mm + (yy - 1999) * 64) * 5.625 * pai / 180
    
17. jd1 = earth_t(W - pai) * 365250 + 2451545#
    
18. 
    
19. Do
    
20.    jd0 = jd1
    
21.    stDegree = earth_L(jd0) - W
    
22.    stDegreep = (earth_L(jd0 + 0.000005) - earth_L(jd0 - 0.000005)) / 0.00001
    
23.    jd1 = jd0 - stDegree / stDegreep
    
24. Loop Until Abs(jd1 - jd0) < 0.0000001
    
25. 
    
26. getjq_64 = JD2GL(jd1 + 8 / 24 - deltatT(yy) / 86400)
    
27. End Function

复制代码

zopey

‘求解农历初一

Function getjq_12(yy, mm)
W = (mm + (yy - 2000) * 12) * pai * 2

Do
   jd0 = jd1
   stDegree = moon_L(jd0) - earth_L(jd0) - W
   stDegreep = (moon_L(jd0 + 0.000005) - earth_L(jd0 + 0.000005) - moon_L(jd0 - 0.000005) + earth_L(jd0 - 0.000005)) / 0.00001
   jd1 = jd0 - stDegree / stDegreep
Loop Until Abs(jd1 - jd0) < 0.0000001

getjq_12 = JD2GL(jd1 + 8 / 24 - deltatT(yy) / 86400)
End Function



日月合塑.zip (84.94 KB, 下载次数: 94)

2019-1-18 17:08 上传

点击文件名下载附件

zopey

论坛搜索 曲线“交点”， 发现和法师前辈的 想法有点雷同（纯属巧合）。

http://club.excelhome.net/thread-805428-1-1.html

jinhaitiao16

谢谢楼主分享！！！！

zopey

逐年计算 节气的 具体日期：

节气日期变动.zip (116.21 KB, 下载次数: 117)

2019-2-26 15:52 上传

点击文件名下载附件

xq1234

Function JD2GL(JD As Double)         

    'J2000 = 2451545
    '地球自转修正项 需完善
    'JD = J2000
   
    '转换日期
    z = Int(JD + 0.5)
    F = JD + 0.5 - z
   
    a0 = Int((z - 1867216.25) / 36524.25)
    If z < 2299161 Then
       a = z
    Else
       a = z + 1 + a0 - Int(a0 / 4)
    End If
    B = a + 1524
    c = Int((B - 122.1) / 365.25)
    D = Int(365.25 * c)
    E = Int((B - D) / 30.6001)
   
    ''''
    d1 = B - D - Int(30.6001 * E) + F
   
    If E < 14 Then
      m1 = E - 1
    Else
      m1 = E - 13
    End If
   
    If m1 > 2 Then
      y1 = c - 4716
    Else
      y1 = c - 4715
    End If
   
    d2 = Int((d1 - Int(d1)) * 86400)
    hh1 = Int(d2 / 3600)
   
    mm1 = Int((d2 - hh1 * 3600) / 60)
    If mm1 < 10 Then mm1 = "0" & mm1
    If mm1 = 60 Then mm1 = 59
   
    ss1 = Int((d1 - Int(d1)) * 86400 - hh1 * 3600 - mm1 * 60)
    If ss1 < 10 Then ss1 = "0" & ss1
    If ss1 = 60 Then ss1 = 59
        
    JD2GL = CDate(y1 & "-" & m1 & "-" & Int(d1) & " " & hh1 & ":" & mm1 & ":" & ss1)
End Function

修改了一下，这样才是真正的“日期”数值。