module1.bas

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

   M = (357.51716 + 0.985599987 * tdays) * 3600#
   l = (280.46435 + 0.9856091 * tdays) * 3600#
   lj = (34.40438 + 0.083086762 * tdays) * 3600#
   M = fmod(M, rev) / RAD2SEC
   l = fmod(l, rev) / RAD2SEC
   lj = fmod(lj, rev) / RAD2SEC
   E = M + ecc * Sin(M) + 0.5 * ecc * ecc * Sin(2# * M)
   secdiff = 0.001658 * Sin(E) + 0.00002073 * Sin(l - lj)
   tdtjd = tdb - secdiff / 86400#
End Function
Private Function julian_date(ByVal year As Long, ByVal month As Long, ByVal day As Long, ByVal hour As Double) As Double
   Dim jd12h As Long, mm As Long
   mm = clong(month - 14, 12)
   jd12h = day - 32075 + clong(1461 * (year + 4800 + mm), 4)
   jd12h = jd12h + clong((367 * (month - 2 - mm * 12)), 12)
   jd12h = jd12h - clong(3 * clong((year + 4900 + mm), 100), 4)
   julian_date = jd12h - 0.5 + hour / 24#
End Function
Private Function clong(ByVal x As Long, ByVal y As Long) As Long
   clong = Sgn(x) * Int(Abs(x) / y)
End Function
Private Function cal_date(ByVal tjd As Double, ByRef year As Long, ByRef month As Long, ByRef day As Long, ByRef hour As Double)
   Dim jd As Long, k As Long, M As Long, n As Long, djd As Double
   djd = tjd + 0.5
   jd = Int(djd)
   hour = (djd - Int(djd)) * 24#
   k = jd + 68569
   n = Int(4 * k / 146097)
   k = k - clong((146097 * n + 3), 4)
   M = clong(4000 * (k + 1), 1461001)
   k = k - clong(1461 * M, 4) + 31
   month = clong(80 * k, 2447)
   day = Int(k - clong(2447 * month, 80))
   k = clong(month, 11)
   month = Int(month + 2 - 12 * k)
   year = Int(100 * (n - 49) + M + k)
End Function

Private Function fmod(ByVal A As Double, ByVal b As Double) As Double
   fmod = Int(A / b)
   fmod = A - fmod * b
End Function
Private Function IsNormalYear(ByVal jdws As Double, ByRef vmoons() As Double) As Boolean
   i = 0: nmonth = 0
   jend = jdws + 0.5
   jend = Int(jend) + 0.5
   Do While (vmoons(i) < jend)
      nmonth = nmonth + 1
      i = i + 1: If i > UBound(vmoons) Then Exit Do
   Loop
   IsNormalYear = (nmonth = 12) ' Normal year
End Function
Private Function IsZhongQiInMonth(ByVal jstart As Double, ByVal jend As Double, ByRef vterms() As Double) As Boolean
    jstart = jstart - 0.5
    jstart = Int(jstart) + 0.5
    jend = jend - 0.5
    jend = Int(jend) + 0.5
    IsZhongQiInMonth = False
    For i = 1 To UBound(vterms) - LBound(vterms) Step 2
        IsZhongQiInMonth = (vterms(i) >= jstart And vterms(i) < jend)
        If IsZhongQiInMonth Then Exit Function
        If (vterms(i) >= jend) Then Exit Function
    Next i
End Function
Private Function TrimHour(ByRef vjds() As Double)
    For i = 0 To UBound(vjds)
        vjds(i) = Int(vjds(i) + 0.5)
    Next i
End Function
Private Function lunaryear(ByVal year As Long, ByRef vterms() As Double, ByRef lastnew As Double, ByRef lastmon As Double, ByRef vmoons() As Double, ByRef vmonth() As Double, ByRef nextnew As Double) As Double
       Dim jdpws As Double, jstart As Double, jend As Double
       Call initial_const
       If TypeName(pre_year) = "Long" And Val(pre_year) = year Then Exit Function
       pre_year = year
       Call solarterm(year, jdpws, vterms)
       ' Determine new moons since day after previous winter solstice
       jstart = jdpws + 0.5
       jstart = Int(jstart) + 0.5
       jend = julian_date(year, 12, 31, 23.999)
       Call mphases(jstart, jend, 0, vmoons)
       ' Determine the month numbers
       ReDim vmonth(UBound(vmoons))
       If (IsNormalYear(vterms(23), vmoons)) Then ' Normal year
           n = 12
           For i = 0 To UBound(vmoons)
                vmonth(i) = n
                If (n = 12) Then n = 0
                n = n + 1
           Next i
           If (n = 1) Then ' Involves next lunar year
              Dim jdnwsp As Double
              Dim vjdnterms() As Double
              Call solarterm(year + 1, jdnwsp, vjdnterms)
              Dim vnmoons() As Double
              jstart = jdnwsp + 0.5
              jstart = Int(jstart) + 0.5
              jend = julian_date(year + 1, 12, 31, 23.999)
              Call mphases(jstart, jend, 0, vnmoons)
              If (Not IsNormalYear(vjdnterms(23), vnmoons)) Then
                  'Check if a zhongqi falls inside the month
                 If (Not IsZhongQiInMonth(vnmoons(0), vnmoons(1), vjdnterms)) Then vmonth(UBound(vmonth)) = 11.5
              End If
           End If
       Else    'Leap year
              n = 11
              bleaped = False
              For i = 0 To UBound(vmoons) - LBound(vmoons) - 1
                ' Check if a zhongqi falls inside the month
                If (bleaped Or IsZhongQiInMonth(vmoons(i), vmoons(i + 1), vterms)) Then
                    n = n + 1: vmonth(i) = n
                Else
                   vmonth(i) = n + 0.5
                   bleaped = True
                End If
                If n = 12 Then n = 0
              Next i
              vmonth(UBound(vmonth)) = n + 1
        End If
    If (vmoons(0) < julian_date(year, 1, 1, 0#)) Then
        lastnew = Int(front(vmoons) + 0.5)
        lastmon = front(vmonth)
        Call Aerase(vmoons, begin(vmoons))
        Call Aerase(vmonth, begin(vmonth))
    Else 'Need to find the last new moon for previous year
        Dim vlastnew() As Double
        Call mphases(vmoons(0) - 35#, vmoons(0) - 25#, 0, vlastnew)
        Call TrimHour(vlastnew)
        lastnew = back(vlastnew)
        lastmon = 11#
    End If
    'Need to find the first new moon for next year
    Dim vnextnew() As Double
    Call mphases(back(vmoons) + 25#, back(vmoons) + 35#, 0, vnextnew)
    Call TrimHour(vnextnew)
    nextnew = back(vnextnew)
    ' Convert to whole day numbers
    Call TrimHour(vmoons)
    Call TrimHour(vterms)
    'Scan for leap month and return the calendar month
    If (Int(lastmon + 0.9) <> Int(lastmon)) Then 'lastmon is a leap month
        If (julian_date(year, 1, 1, 12#) < vmoons(0)) Then lunaryear = 0.5: Exit Function ' runs into new year
    End If
    Dim monnum As Long: monnum = 2
    For i = 0 To UBound(vmoons)
        If (Int(vmonth(i) + 0.9) <> Int(vmonth(i))) Then 'found leap month
            Dim jdfirst As Double
            jdfirst = julian_date(year, monnum, 1, 12#)
            Do While (monnum <= 12 And jdfirst < vmoons(i))
                monnum = monnum + 1: jdfirst = julian_date(year, monnum, 1, 12#)
            Loop
            If (monnum = 13) Then lunaryear = 12#: Exit Function
            ' See if leap month runs into next month
            tt = (i <> UBound(vmoons)): If tt Then tt = tt And jdfirst < vmoons(i + 1)
            If tt Then
               lunaryear = (monnum - 0.5): Exit Function  'Yes
            Else
               lunaryear = (monnum - 1#): Exit Function   'No
            End If
         End If
    Next i
    lunaryear = 0#
End Function
Private Function back(ar As Variant)
    back = ar(UBound(ar))
End Function
Private Function begin(ar As Variant)
    begin = LBound(ar)
End Function
Private Function Aerase(ar As Variant, ByVal i)
    n = UBound(ar) - 1
    For j = i To n
        ar(j) = ar(j + 1)
    Next j
    If n >= 0 Then ReDim Preserve ar(n)
End Function
Private Function push_back(ar As Variant, ByVal x)
   ReDim Preserve ar(UBound(ar) + 1): ar(UBound(ar)) = x
End Function
Private Function front(ar As Variant)
    front = ar(LBound(ar))
End Function

Private Function torad(ByVal x As Double) As Double
  torad = fmod(x, 360#) * 1.74532925199433E-02       ' normalize the angle
End Function

Private Function moonphasebylunation(ByVal lun As Long, ByVal phi As Long) As Double
  Dim k As Double
  k = lun + phi / 4#
  moonphasebylunation = moonphase(k, phi)
End Function

Private Function moonphase(ByVal k As Double, ByVal phi As Long) As Double
  Dim T As Double   ' time parameter, Julian Centuries since J2000
  Dim JDE As Double  ' Julian Ephemeris Day of phase event
  Dim E As Double  '                    /* Eccentricity anomaly
  Dim M As Double  '                    /* Sun's mean anomaly */
  Dim M1 As Double '                  /* Moon's mean anomaly */
  Dim f As Double  '                    /* Moon's argument of latitude */
  Dim O As Double  '                    /* Moon's longitude of ascenfing node */
  Dim A(15) As Double '                /* planetary arguments */
  Dim W As Double    '                  /* added correction for quarter phases */

  T = k / 1236.85                            ' /* (47.3) */

  '/* this is the first approximation.  all else is for style points! */
  JDE = 2451550.09765 + 29.530588853 * k + T * T * (0.0001337 + T * (-0.00000015 + 0.00000000073 * T))    '/* (47.1) */
 
  '/* these are correction parameters used below */
  E = 1# + T * (-0.002516 + -0.0000074 * T)                                    ' /* (45.6) */
  M = 2.5534 + 29.10535669 * k + T * T * (-0.0000218 + -0.00000011 * T)               '/* (47.4) */
  M1 = 201.5643 + 385.81693528 * k + T * T * (0.0107438 + T * (0.00001239 + -0.000000058 * T))           '/* (47.5) */
  f = 160.7108 + 390.67050274 * k + T * T * (-0.0016341 * T * (-0.00000227 + 0.000000011 * T))            '/* (47.6) */
  O = 124.7746 - 1.5637558 * k + T * T * (0.0020691 + 0.00000215 * T)                '/* (47.7) */

  '/* planetary arguments */
  A(0) = 0  '/* unused! */
  A(1) = 299.77 + 0.107408 * k - 0.009173 * T * T
  A(2) = 251.88 + 0.016321 * k
  A(3) = 251.83 + 26.651886 * k
  A(4) = 349.42 + 36.412478 * k
  A(5) = 84.66 + 18.206239 * k
  A(6) = 141.74 + 53.303771 * k
  A(7) = 207.14 + 2.453732 * k
  A(8) = 154.84 + 7.30686 * k
  A(9) = 34.52 + 27.261239 * k
  A(10) = 207.19 + 0.121824 * k
  A(11) = 291.34 + 1.844379 * k
  A(12) = 161.72 + 24.198154 * k
  A(13) = 239.56 + 25.513099 * k
  A(14) = 331.55 + 3.592518 * k

  '/* all of the above crap must be made into radians!!! */
  '/* except for E... */

  M = torad(M)
  M1 = torad(M1)
  f = torad(f)
  O = torad(O)

  '/* all those planetary arguments, too! */
  For i = 1 To 14
    A(i) = torad(A(i))
  Next i
  ' ok, we have all the parameters, let's apply them to the JDE.
  '  (remember the JDE?  this is a program about the JDE...)
  Select Case phi
       '/* a special case for each different phase.  NOTE!,
       ' i 'm not treating these in a 0123 order!!!  Payattention, there,  you!
    Case 0  '/* New Moon */
      JDE = JDE _
          - 0.4072 * Sin(M1) _