' Linear Approximation Progressive Leap Rule Functions
' by Dr. Irv Bromberg, University of Toronto, Canada
' <http://www.sym454.org/leap/>

' Version 1.3(4), revised Gregorian 2010.07.09 Changed LANEY offset from +0.9 to +0.4,
' and changed LANSY offset from +3 to +0.4.
' Version 1.2(3), revised Gregorian 2006.07.07 Expanded comments, defined syntax, minor improvements,
' revised arithmetic to compute integer parts of new year moments separately from the fractional part,
' for improved floating point precision.
' (previous version was 1.1(2), revised Gregorian 2006.01.05)

' Implements a fixed arithmetic leap rule with a progressively changing calendar mean year,
' based on a starting fixed calendar mean year until the begin-progressive-era year number,
' then a progressively changing calendar mean year until the end of the progressive era,
' after which the end-of-progressive-era fixed calendar mean year will apply.

' Typically, each era lasts about 10 millennia, for further information explaining why this
' is so see "The Lengths of the Seasons" at <http://www.sym454.org/seasons/>.

' These functions are designed to work in terms of the Reingold & Dershowitz "rata die",
' the number of days elapsed since the Gregorian epoch, with time of day expressed as a
' fraction of day elapsed since civil midnight, as described in their book
' "Calendrical Calculations:  The Millennium Edition" (abbreviated hereinafter as CCME)
' [see their web site at <http://www.calendarists.com/>].

' If you prefer to work in terms of the Julian Day, use the MomentToJulianDay() function
' or its converse JulianDayToMoment(), both shown below, or use a similar function to
' convert the rata die to an ordinal day count relative to any desired epoch.

' Further discussion about calendrical calculations and such interconversions is given in
' the Symmetry454 Arithmetic PDF that is available at <http://www.sym454.org/symmetry/>.

' The LinApprox( ) function returns the New Year Moment for the specified year using the
' specified linear approximation leap rule.  For normal calendrical calculations a higher-
' level function would be called, such as NewYearDay( ) for a leap day calendar,
' or NewYearDayLW( ) for a leap week calendar.  Once the New Year Day is known, the user
' can employ it to start the year of any desired solar calendar structure.

' The leap status of a year is given by the isLeapYear( ) function for a leap day calendar,
' or the isLeapYearLW( ) function for a leap week calendar.

' ApproxRule string constants to be passed as a parameter to functions that require it
Public Const LANEY As String = "LANEY"
Public Const LANSY As String = "LANSY"
Public Const LASEY As String = "LASEY"
Public Const LASSY As String = "LASSY"
Public Const LABY As String = "LABY"

' Calling syntax for the Linear Approximation Progressive Leap Rule functions herein:

' If only the New Year Moment is required, for example to compare to a mean astronomical algorithm:
' MyNewYearMoment = NewYearMoment( TheYear, ApproxRule )

' An even lower level degree of control of the New Year Moment is possible by calling the LinApprox( )
' function directly, passing the parameters that it needs.  This is normally handled for the user by
' calling the higher-level NewYearDay( ) or NewYearDayLW( ) function:
' MyNewYearMoment = LinApprox(TheYear, YearBeginProgressive, YearEndProgressive, StartFraction, EndFraction, Offset)
' where YearBeginProgressive is the first year number of the progressive leap rule era,
' YearEndProgressive is the year number at the end of the progressive leap rule era,
' StartFraction is the fixed calendar mean year before the progressive leap rule era,
' EndFraction is the fixed calendar mean year after the progressive leap rule era,
' and Offset is the number of days to add the epoch to account for the first New Year Moment
' not landing exactly at the moment of the epoch.  The user can empirically adjust the Offset parameter
' to obtain the desired average long-term alignment of an equinox or solstice moment relative to a
' chosen date on the target calendar.

' Note that in all cases the calendar mean years are expressed as the fraction of a day in excess
' of 365 days so that the calculations can be carried out to 3 more significant figures of precision
' than would be possible if the fixed integer value 365 were included with the fraction.  In other words,
' the LinApprox( ) function calculates the integer part separately from the fractional part, and then
' sums them to compute the New Year Moment.

' For a leap day calendar, having 365 days in a non-leap year, 366 days in a leap year:
' MyNewYearDay = NewYearDay( TheYear, ApproxRule )
' for example MyNewYearDay = NewYearDay( 2006, LANSY )
' isLeap = isLeapYear( TheYear, ApproxRule )
' for example isLeap = isLeapYear( 2010, LASSY )

' For a leap week calendar, having 364 days in a non-leap year, 371 days in a leap year:
' MyNewYearDay = NewYearDayLW( TheYear, ApproxRule, StartOnWeekday )
' isLeap = isLeapYearLW( TheYear, ApproxRule, StartOnWeekday )
' where StartOnWeekday is an integer or enumeration from Sunday=0, Monday=1 ... Saturday=6:
Public Enum WeekDay
    Sunday
    Monday
    Tuesday
    Wednesday
    Thursday
    Friday
    Saturday
End Enum


' This source code has been written as a general purpose set of routines, for experimental evaluation.
' In a non-experimental permanent calendar implementation, the ApproxRule and StartOnWeekday can be
' made constants and used directly, eliminating much of the general purpose structure given herein.


' define some minor constants just to make the program code more self-evident

Public Const HoursPerDay As Integer = 24
Public Const MinutesPerDay As Integer = HoursPerDay * 60 ' 1440
Public Const SecondsPerDay As Long = 86400 ' can't show as expression MinutesPerDay * 60, gets overflow error in VB6

Public Const DaysPerWeek As Integer = 7
Public Const Days365 As Integer = 365

Public Const jdEPOCH As Double = -1721424.5    ' Noon on proleptic Gregorian date Nov 24, -4713
Public Const GregorianEPOCH As Double = 1#     ' rata die number 1 started on proleptic Gregorian January 1, 1 AD.


' define some general purpose auxiliary functions in case not built into the programming language

Public Function floor(ByVal x As Double) As Long

    floor = Int(x)

' The following is a more involved method for programming languages where INT() truncates
' instead of returning the greatest integer equal to or less than x
' Dim intx As Long
'
'    intx = Int(x)
'    If intx <= x Then floor = intx Else floor = intx - 1

End Function


Public Function modulus(ByVal x As Double, ByVal Y As Double) As Double
    
    modulus = x - Y * floor(x / Y)

' use this modulus function instead of built-in mod operator,
' to ensure that the following are true:
    ' x and y need not be integers
    ' if y>0 then x modulus y is >=0 for all values of x (even negative)
    ' (-x) modulus y = y - (x modulus y) for y>0
    ' for y<>0 and z<>0:  a = (x modulus y)  if and only if  az = (xz modulus yz)
    ' x - (x modulus y)  is always a multiple of y
    ' for y<>0:  0 <= sign(y) x (x modulus y) < absolute value of y

End Function


Public Function minimum(ParamArray Numbers() As Variant) As Double

    ' generic function to return the minimum numeric value
    ' from any number of parameters passed to this function
    ' (although the source code herein only calls it with two parameters)

    Dim v As Variant
    Dim ThisNumber As Double

    minimum = 1E+308
    
    For Each v In Numbers
        ThisNumber = CDbl(v)
        If ThisNumber < minimum Then minimum = ThisNumber
    Next v

End Function


' The following show examples to interconvert the rata die to ordinal numbering relative to
' an alternative epoch, in this case the epoch of the Julian Day Number:

Public Function MomentToJulianDay(ByVal tee As Double) As Double

    ' convert a CCME moment to a Julian Day number

    MomentToJulianDay = tee - jdEPOCH

End Function


Public Function JulianDayToMoment(ByVal JD As Double) As Double

    ' convert a Julian Day number to a CCME moment

    JulianDayToMoment = JD + jdEPOCH

End Function



Public Function NewYearDay(ByVal TheYear As Long, ByVal ApproxRule As String) As Long

  ' for a leap DAY calendar, simply floor the New Year Moment

  NewYearDay = floor(NewYearMoment(TheYear, ApproxRule))
        
End Function



Public Function NewYearDayLW(ByVal TheYear As Long, ByVal ApproxRule As String, ByVal StartOnWeekday As WeekDay) As Long

   ' for a leap WEEK calendar, this function works like the CCME k-day-nearest function
   ' k-day-nearest functionality MUST have an integer to work with, otherwise returns erroneous results

   ' +3 offset required for k-day-nearest functionality
    NewYearDayLW = floor(NewYearMoment(TheYear, ApproxRule)) + 3#
    
   ' now adjust the New Year Day to the nearest k-day, depending on which weekday begins the calendar year
   ' set the following to the desired weekday number, where Sunday=0, Monday=1, ... Saturday=6.

   NewYearDayLW = NewYearDayLW - modulus(NewYearDayLW - StartOnWeekday, DaysPerWeek)
        
End Function



Public Function isLeapYear(ByVal TheYear As Long, ByVal ApproxRule As String) As Boolean

   ' Generalized for a leap DAY calendar using the specified Linear Approximation Progressive Leap Rule.
   ' Works by simply calculating the length of the year in days, can only be 365 or 366 days.
   ' Returns TRUE if it is a leap year (366 days), false if a non-leap year (365 days).

   isLeapYear = (NewYearDay(TheYear + 1, ApproxRule) - NewYearDay(TheYear, ApproxRule) > Days365)

End Function



Public Function isLeapYearLW(ByVal TheYear As Long, ByVal ApproxRule As String, ByVal StartOnWeekday As WeekDay) As Boolean

   ' Generalized for a leap WEEK calendar using the specified Linear Approximation Progressive Leap Rule.
   ' Works by simply calculating the length of the year in days, can only be 364 or 371 days.
   ' Returns TRUE if it is a leap year (371 days), false if a non-leap year (364 days).

   isLeapYearLW = (NewYearDayLW(TheYear + 1, ApproxRule, StartOnWeekday) - NewYearDayLW(TheYear, ApproxRule, StartOnWeekday) > 364)

End Function



Public Function NewYearMoment(ByVal TheYear As Long, ByVal ApproxRule As String) As Double

  ' returns the New Year Moment using the specified Linear Approximation rule

  Select Case ApproxRule

    Case LANEY  ' Linear Approximation to the Northward Equinoctial year
            
        Const LANEYstartFract As Double = 71# / 293#    ' 5h48m57s
        Const LANEYendFract As Double = 5# / HoursPerDay + 46# / MinutesPerDay + 22# / SecondsPerDay

        NewYearMoment = LinApprox(TheYear, 5500, 14000, LANEYstartFract, LANEYendFract, 0.4)

    Case LANSY  ' Linear Approximation to the North Solstitial year
            
        Const LANSYstartFract As Double = 94# / 389#    ' 5h47m58s
        Const LANSYendFract As Double = 5# / HoursPerDay + 46# / MinutesPerDay + 20# / SecondsPerDay

        NewYearMoment = LinApprox(TheYear, 10500, 18200, LANSYstartFract, LANSYendFract, 0.4)

    Case LASEY   ' Linear Approximation to the Southward Equinoctial year
            
        Const LASEYstartFract As Double = 5# / HoursPerDay + 50# / MinutesPerDay + 35# / SecondsPerDay
        Const LASEYendFract As Double = 5# / HoursPerDay + 47# / MinutesPerDay + 13# / SecondsPerDay

        NewYearMoment = LinApprox(TheYear, -3000, 5200, LASEYstartFract, LASEYendFract, 0.25)


    Case LASSY   ' Linear Approximation to the South Solstitial year
            
        Const LASSYstartFract As Double = 5# / HoursPerDay + 49# / MinutesPerDay + 47# / SecondsPerDay
        Const LASSYendFract As Double = 5# / HoursPerDay + 46# / MinutesPerDay + 45# / SecondsPerDay

        NewYearMoment = LinApprox(TheYear, 1400, 9250, LASSYstartFract, LASSYendFract, -1#)


    Case LABY   ' Linear Approximation to the Besselian year
            
        Const LABYstartFract As Double = 5# / HoursPerDay + 49# / MinutesPerDay + 44# / SecondsPerDay
        Const LABYendFract As Double = 5# / HoursPerDay + 46# / MinutesPerDay + 46# / SecondsPerDay

        NewYearMoment = LinApprox(TheYear, 1600, 9700, LABYstartFract, LABYendFract, -1.3)

End Function



Public Function LinApprox(ByVal TheYear As Long, ByVal YearBeginProgressive As Long, ByVal YearEndProgressive As Long, _
        ByVal StartFraction As Double, ByVal EndFraction As Double, _
        ByVal Offset As Double) As Double
        
    ' Returns the rata die New Year moment using the specified linear approximation leap cycle parameters.
    ' For a leap day calendar either FLOOR the moment or ROUND it to the nearest midnight, as desired.
    ' For a leap week calendar, start the year on the desired weekday closest to the moment.
    ' For example, if the leap week calendar starts every year on Monday then find the Monday that is closest
    ' to the New Year moment.
    
    ' Given the rata die of the target New Year Day and then calculating the rata die of the following New Year Day,
    ' the year is a leap year if the number of days that separate those New Year Days is greater than the number
    ' of days in a non-leap year.  This distributes leap years at intervals that are as uniformly spread as possible.

    Dim Slope As Double
    Dim Intercept As Double
    Dim YearsBefore As Long
    Dim SlopeYears As Long
    Dim YearsAfter As Long
    
    ' work with prior year, because we need the number of days from the epoch to the LinApprox New Year moment
    TheYear = TheYear - 1

    ' Set up parameters used for integrating the sloped linear approximation region:
    ' (Just for info, to convert Slope to seconds per century, multiply by 86400
    ' and then divide by the number of centuries from YearBeginProgressive to YearEndProgressive.)
    Slope = (EndFraction - StartFraction) / (YearEndProgressive - YearBeginProgressive)
    Intercept = StartFraction - Slope * YearBeginProgressive
    
    ' Adjust optional Offset value if necessary to obtain desired long-term mean equinox or solstice alignment
    LinApprox = GregorianEPOCH + Offset

    If TheYear > 0 Then     ' if the target year is after the calendar epoch then...
    
        If YearBeginProgressive > 0 Then   ' if the start of the sloped linear approximation region is after the epoch then...
        
            If TheYear <= YearBeginProgressive Then ' if we are before the sloped linear approximation region...
            
                ' just calculate the simple rectangular area of fixed mean year region from epoch to TheYear
                ' (calculate integer part separately from fractional part for higher floating point precision)
                ' LinApprox = LinApprox + TheYear * (Days365 + StartFraction)
                LinApprox = LinApprox + TheYear * Days365 + TheYear * StartFraction
            
            Else ' otherwise TheYear is after the sloped region YearBeginProgressive
            
                ' calculate the area of the rectangle from epoch to YearBeginProgressive
                ' (calculate integer part separately from fractional part for higher floating point precision)
                ' LinApprox = LinApprox + YearBeginProgressive * (Days365 + StartFraction)
                LinApprox = LinApprox + YearBeginProgressive * Days365 + YearBeginProgressive * StartFraction
                
                ' compute the number of years to include from the sloped region
                SlopeYears = minimum(TheYear, YearEndProgressive) - YearBeginProgressive
                
                ' add in the integral of the sloped region
                ' (note that this is a quadratic equation in terms of SlopeYears)
                ' (calculate integer parts separately from fractional part for higher floating point precision)
                ' LinApprox = LinApprox + SlopeYears * (Days365 + StartFraction + Slope * SlopeYears / 2#)
                LinApprox = LinApprox + SlopeYears * Days365 + SlopeYears * StartFraction + SlopeYears * SlopeYears * Slope / 2#
                
                ' below we will add in the final rectangular area from YearEndProgressive to TheYear, if applicable
                
            End If
            
        Else ' the start of the sloped linear approx region is before the epoch
        
            ' TheYear must be after YearBeginProgressive, because from above we know that it is after epoch.
            ' Integrate the area of the sloped region from the epoch to Minimum(TheYear, YearEndProgressive)
            SlopeYears = minimum(TheYear, YearEndProgressive)  ' don't subtract from YearBeginProgressive, just area since epoch
            ' (calculate integer parts separately from fractional part for higher floating point precision)
            ' LinApprox = LinApprox + SlopeYears * (Days365 + Intercept + Slope * SlopeYears / 2#)
            LinApprox = LinApprox + SlopeYears * Days365 + SlopeYears * Intercept + SlopeYears * SlopeYears * Slope / 2#
            
            ' below we will add in rectangle area from YearEndProgressive to TheYear, if applicable
        
        End If
        
        If TheYear > YearEndProgressive Then   ' if the year is beyond the end of the sloped linear approx region then...
        
            ' also add in rectangle area from YearEndProgressive to TheYear
            YearsAfter = TheYear - YearEndProgressive
            ' (calculate integer part separately from fractional part for higher floating point precision)
            ' LinApprox = LinApprox + YearsAfter * (Days365 + EndFraction)
            LinApprox = LinApprox + YearsAfter * Days365 + YearsAfter * EndFraction
        
        End If
        
    Else ' the year is prior to the calendar epoch
    
        If YearBeginProgressive > 0 Then   ' if the sloped region starts after the epoch then...
        
            ' YearBeginProgressive is after epoch but TheYear is before epoch
            ' so we only need the negative of the area of early mean year from TheYear to epoch
            ' (since TheYear is negative the area will also be negative)
            ' (calculate integer part separately from fractional part for higher floating point precision)
            ' LinApprox = LinApprox + TheYear * (Days365 + StartFraction)
            LinApprox = LinApprox + TheYear * Days365 + TheYear * StartFraction
            
        Else ' the sloped region starts prior to the epoch, now we check where relative to the target year
        
            If TheYear <= YearBeginProgressive Then    ' if the target year is prior to the sloped region
            
                ' assume that sloped region crosses epoch, we won't bother with sloped region entirely prior to epoch
                
                ' calculate the area of the rectangle from TheYear to YearBeginProgressive
                YearsBefore = TheYear - YearBeginProgressive
                ' (calculate integer part separately from fractional part for higher floating point precision)
                ' LinApprox = LinApprox + YearsBefore * (Days365 + StartFraction)
                LinApprox = LinApprox + YearsBefore * Days365 + YearsBefore * StartFraction
                
                ' next add in the area of the sloped region from YearBeginProgressive to epoch
                ' since YearBeginProgressive is negative (prior to epoch), this area is negative
                ' (calculate integer part separately from fractional part for higher floating point precision)
                ' LinApprox = LinApprox + YearBeginProgressive * (Days365 + Intercept + Slope * YearBeginProgressive / 2#)
                LinApprox = LinApprox + YearBeginProgressive * Days365 + YearBeginProgressive * Intercept + YearBeginProgressive * YearBeginProgressive * Slope / 2#
                
            Else ' otherwise the target year is after the start of the sloped region
                 ' and we know that the sloped region starts before the epoch
            
                ' just calculate the area of the sloped region from the target year to the epoch
                ' (we assume that the sloped region must continue beyond epoch)
                ' (calculate integer part separately from fractional part for higher floating point precision)
                ' LinApprox = LinApprox + TheYear * (Days365 + Intercept + Slope * TheYear / 2#)
                LinApprox = LinApprox + TheYear * Days365 + TheYear * Intercept + TheYear * TheYear * Slope / 2#
                
            End If  ' TheYear <= YearBeginProgressive
        
        End If  ' YearBeginProgressive > 0
        
    End If  ' TheYear > 0
    
End Function