module1.bas

vb编写的中国日梭万年历,可查询年1583-99999年

          + 0.17241 * E * Sin(M) _
          + 0.01608 * Sin(2# * M1) _
          + 0.01039 * Sin(2# * f) _
          + 0.00739 * E * Sin(M1 - M) _
          - 0.00514 * E * Sin(M1 + M) _
          + 0.00208 * E * E * Sin(2# * M) _
          - 0.00111 * Sin(M1 - 2# * f) _
          - 0.00057 * Sin(M1 + 2# * f) _
          + 0.00056 * E * Sin(2# * M1 + M) _
          - 0.00042 * Sin(3# * M1) _
          + 0.00042 * E * Sin(M + 2# * f) _
          + 0.00038 * E * Sin(M - 2# * f) _
          - 0.00024 * E * Sin(2# * M1 - M) _
          - 0.00017 * Sin(O) - 0.00007 * Sin(M1 + 2# * M) _
          + 0.00004 * Sin(2# * M1 - 2# * f) _
          + 0.00004 * Sin(3# * M) _
          + 0.00003 * Sin(M1 + M - 2# * f) _
          + 0.00003 * Sin(2# * M1 + 2# * f) _
          - 0.00003 * Sin(M1 + M + 2# * f) _
          + 0.00003 * Sin(M1 - M + 2# * f) _
          - 0.00002 * Sin(M1 - M - 2# * f) _
          - 0.00002 * Sin(3# * M1 + M) _
          + 0.00002 * Sin(4# * M1)

    Case 2 ' /* Full Moon */
      JDE = JDE _
          - 0.40614 * Sin(M1) _
          + 0.17302 * E * Sin(M) _
          + 0.01614 * Sin(2# * M1) _
          + 0.01043 * Sin(2# * f) _
          + 0.00734 * E * Sin(M1 - M) _
          - 0.00515 * E * Sin(M1 + M) _
          + 0.00209 * E * E * Sin(2# * M) _
          - 0.00111 * Sin(M1 - 2# * f) _
          - 0.00057 * Sin(M1 + 2# * f) _
          + 0.00056 * E * Sin(2# * M1 + M) _
          - 0.00042 * Sin(3# * M1) _
          + 0.00042 * E * Sin(M + 2# * f) _
          + 0.00038 * E * Sin(M - 2# * f) _
          - 0.00024 * E * Sin(2# * M1 - M) _
          - 0.00017 * Sin(O) - 0.00007 * Sin(M1 + 2# * M) _
          + 0.00004 * Sin(2# * M1 - 2# * f) _
          + 0.00004 * Sin(3# * M) _
          + 0.00003 * Sin(M1 + M - 2# * f) _
          + 0.00003 * Sin(2# * M1 + 2# * f) _
          - 0.00003 * Sin(M1 + M + 2# * f) _
          + 0.00003 * Sin(M1 - M + 2# * f) _
          - 0.00002 * Sin(M1 - M - 2# * f) _
          - 0.00002 * Sin(3# * M1 + M) _
          + 0.00002 * Sin(4# * M1)

    Case 3 Or 1 ' /* Last Quarter  /* First Quarter */
      JDE = JDE _
          - 0.62801 * Sin(M1) + 0.17172 * E * Sin(M) _
          - 0.01183 * E * Sin(M1 + M) + 0.00862 * Sin(2# * M1) _
          + 0.00804 * Sin(2# * f) + 0.00454 * E * Sin(M1 - M) _
          + 0.00204 * E * E * Sin(2# * M) _
          - 0.0018 * Sin(M1 - 2# * f) _
          - 0.0007 * Sin(M1 + 2# * f) _
          - 0.0004 * Sin(3# * M1) _
          - 0.00034 * E * Sin(2# * M1 - M) _
          + 0.00032 * E * Sin(M + 2# * f) _
          + 0.00032 * E * Sin(M - 2# * f) _
          - 0.00028 * E * E * Sin(M1 + 2# * M) _
          + 0.00027 * E * Sin(2# * M1 + M) _
          - 0.00017 * Sin(O) - 0.00005 * Sin(M1 - M - 2# * f) _
          + 0.00004 * Sin(2# * M1 + 2# * f) _
          - 0.00004 * Sin(M1 + M + 2# * f) _
          + 0.00004 * Sin(M1 - 2# * M) _
          + 0.00003 * Sin(M1 + M - 2# * f) _
          + 0.00003 * Sin(3# * M) + 0.00002 * Sin(2# * M1 - 2# * f) _
          + 0.00002 * Sin(M1 - M + 2# * f) - 0.00002 * Sin(3# * M1 + M)
      W = 0.00306 _
        - 0.00038 * E * Cos(M) _
        + 0.00026 * Cos(M1) _
        - 0.00002 * Cos(M1 - M) _
        + 0.00002 * Cos(M1 + M) _
        + 0.00002 * Cos(2# * f)
      If (phi = 3) Then W = -W
      JDE = JDE + W
    Case Else
      Exit Function
      '/* now there are some final correction to everything */
  End Select
  JDE = JDE _
      + 0.000325 * Sin(A(1)) _
      + 0.000165 * Sin(A(2)) _
      + 0.000164 * Sin(A(3)) _
      + 0.000126 * Sin(A(4)) _
      + 0.00011 * Sin(A(5)) _
      + 0.000062 * Sin(A(6)) _
      + 0.00006 * Sin(A(7)) _
      + 0.000056 * Sin(A(8)) _
      + 0.000047 * Sin(A(9)) _
      + 0.000042 * Sin(A(10)) _
      + 0.00004 * Sin(A(11)) _
      + 0.000037 * Sin(A(12)) _
      + 0.000035 * Sin(A(13)) _
      + 0.000023 * Sin(A(14))
  moonphase = JDE
End Function
Private Function mphases(ByVal tstart As Double, ByVal tend As Double, ByVal phase As Long, ByRef vjdphases() As Double)
    Dim offs As Double
    If (tstart < 2425245) Then
        offs = (116# + 25# / 60#) / 360#
    Else
        offs = 120# / 360#
    End If
    'Find the lunation number of the first given phase after tstart
    Dim lun As Long: lun = CLng((tstart - offs - 2451550.09765 - phase * 7.375) / period)
    Dim jd1 As Double: jd1 = moonphasebylunation(lun, phase) + offs
    Do While (jd1 - tstart > 29)
        lun = lun - 1
        jd1 = moonphasebylunation(lun, phase) + offs
    Loop
    Do While (tstart > jd1)
        lun = lun + 1
        jd1 = moonphasebylunation(lun, phase) + offs
    Loop
    'Compute subsequent phases until after tend
    ReDim vjdphases(0): vjdphases(0) = jd1
    Do While (jd1 < tend - 29)
        lun = lun + 1
        jd1 = moonphasebylunation(lun, phase) + offs
        If (jd1 < tend) Then Call push_back(vjdphases, jd1)
    Loop
End Function
Private Function solarsystem(ByVal tjd As Double, ByVal body As Long, ByVal origin As Long, ByRef pos() As Double, ByRef vel() As Double) As Long
   solarsystem = 0
   Dim tlast As Double, sine As Double, cose As Double, tmass As Double, pbary(3) As Double, vbary(3) As Double
   Dim oblr As Double, qjd As Double, ras As Double, decs As Double, diss As Double, pos1(3) As Double, p(3, 3) As Double, dlon As Double, sinl As Double, cosl As Double, x As Double, y As Double, z As Double, xdot As Double, ydot As Double, zdot As Double, f As Double
   tlast = 0#
   If (tlast = 0#) Then
      oblr = obl * TWOPI / 360#
      sine = Sin(oblr)
      cose = Cos(oblr)
      tmass = 1#
      For i = 0 To 3
         tmass = tmass + 1# / pm(i)
      Next i
      tlast = 1#
   End If
   If ((body = 0) Or (body = 1) Or (body = 10)) Then
      For i = 0 To 2
         pos(i) = 0: vel(i) = 0#
      Next i
   ElseIf ((body = 2) Or (body = 3)) Then
      For i = 0 To 2
         qjd = tjd + CDbl((i - 1) * 0.1)
         Call sun_eph(qjd, ras, decs, diss)
         Call radec2vector(ras, decs, diss, pos1)
         Call precession(qjd, pos1, T0, pos)
         p(i, 0) = -pos(0)
         p(i, 1) = -pos(1)
         p(i, 2) = -pos(2)
      Next i
      For i = 0 To 2
         pos(i) = p(1, i)
         vel(i) = (p(2, i) - p(0, i)) / 0.2
      Next i
   Else
         solarsystem = 2: Exit Function
   End If
   If (origin = 0) Then
      If (Abs(tjd - tlast) >= 0.000001) Then
         For i = 0 To 2
            pbary(i) = 0: vbary(i) = 0#
         Next i
         For i = 0 To 3
            dlon = pl(i) + pn(i) * (tjd - T0)
            dlon = fmod(dlon, TWOPI)
            sinl = Sin(dlon)
            cosl = Cos(dlon)

            x = pa(i) * cosl
            y = pa(i) * sinl * cose
            z = pa(i) * sinl * sine
            xdot = -pa(i) * pn(i) * sinl
            ydot = pa(i) * pn(i) * cosl * cose
            zdot = pa(i) * pn(i) * cosl * sine

            f = 1# / (pm(i) * tmass)

            pbary(0) = pbary(0) + x * f
            pbary(1) = pbary(1) + y * f
            pbary(2) = pbary(2) + z * f
            vbary(0) = vbary(0) + xdot * f
            vbary(1) = vbary(1) + ydot * f
            vbary(2) = vbary(2) + zdot * f
         Next i

         tlast = tjd
      End If
      For i = 0 To 2
         pos(i) = pos(i) - pbary(i)
         vel(i) = vel(i) - vbary(i)
      Next i
   End If
End Function
Private Function sun_eph(ByVal jd As Double, ByRef ra As Double, ByRef dec As Double, ByRef dis As Double)
   Dim sum_lon As Double: sum_lon = 0#
   Dim sum_r As Double: sum_r = 0#
   Dim u As Double, arg As Double, lon As Double, lat As Double, T As Double, t2 As Double, emean As Double, sin_lon As Double
    u = (jd - T0) / 3652500#
'   Compute longitude and distance terms from the series.
 
   For i = 0 To 49
      arg = con(i)(2) + con(i)(3) * u
      sum_lon = sum_lon + con(i)(0) * Sin(arg)
      sum_r = sum_r + con(i)(1) * Cos(arg)
   Next i
'   Compute longitude, latitude, and distance referred to mean equinox and ecliptic of date.
   lon = 4.9353929 + 62833.196168 * u + factor * sum_lon
   lon = fmod(lon, TWOPI)
   If (lon < 0#) Then lon = lon + TWOPI
   lat = 0#
   dis = 1.0001026 + factor * sum_r

'   Compute mean obliquity of the ecliptic.
   T = u * 100#
   t2 = T * T
   emean = (0.001813 * t2 * T - 0.00059 * t2 - 46.815 * T + 84381.448) / RAD2SEC
   
'   Compute equatorial spherical coordinates referred to the mean equator   and equinox of date.
   sin_lon = Sin(lon)
   jj = (Cos(emean) * sin_lon)
   kk = Cos(lon)
   ra = Atn2(jj, kk) * RAD2DEG
   ra = fmod(ra, 360#)
   If (ra < 0#) Then ra = ra + 360#
   ra = ra / 15#
   dec = Asin(Sin(emean) * sin_lon) * RAD2DEG
End Function
Private Function Atn2(ByVal y As Double, ByVal x As Double) As Double
   If x = 0 Then
      Atn2 = Sgn(y) * PI / 2
   Else
      Atn2 = Atn(y / x)
      If x < 0 Then Atn2 = Atn2 + PI
   End If
End Function
Private Function Asin(ByVal x As Double) As Double
     If x = 1 Or x = -1 Then
          Asin = Sgn(x) * PI / 2
     Else
          Asin = Atn(x / Sqr(-x * x + 1))
     End If
End Function
Private Function timeangle(ByVal T As Double) As Double
    Const BARYC = 0
    Const HELIOC = 1
    Dim sposg(3) As Double, spos(3) As Double