VBA functions to calculate sun and moon positions

Question

I have the following VBA code (across three modules) to make UDFs to calculate sun & moon position data. The issue I'm facing is that they are very slow to run as I have over 6000 rows (over 10 years of sun/moon position data for multiple locations), and a sheet calculation takes ages when opening/closing/saving the workbook.

I have tried to make all the functions static but that does not seem to have helped, either because I didn't do it properly or that's not the way to speed it up.

Adivce will be appreciated to speed up the workbook.

Module Name: MoonFunctions

Function UT(ByVal time As Double, ByVal timezone As Integer) As Double
    
    UT = hour(time) - timezone + Minute(time) / 60 + Second(time) / 3600
    
End Function
    
Function JDay(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double) As Double

Dim UT
    UT = MoonFunctions.UT(time, timezone)

Dim INT1
    If Month(dte) > 2 Then
        INT1 = Year(dte) + 4716
    Else: INT1 = Year(dte) - 1 + 4716
    End If
    
Dim INT2
    If Month(dte) > 2 Then
        INT2 = Month(dte) + 1
    Else: INT2 = Month(dte) + 12 + 1
    End If

Dim INT3
    If Month(dte) > 2 Then
        INT3 = Year(dte)
    Else: INT3 = Year(dte) - 1
    End If


    JDay = Int(365.25 * INT1) _
        + Int(30.6001 * INT2) _
        + Day(dte) _
        + UT / 24 _
        + 2 _
        - Int(INT3 / 100) _
        + Int(Int(INT3 / 100) / 4) _
        - 1524.5
    
End Function

Function JCsince2000(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double) As Double  'number of Julian centuries since Jan 1, 2000, 12 UT

Dim JDay
    JDay = MoonFunctions.JDay(time, timezone, dte)
    JCsince2000 = (JDay - 2451545) / 36525
    
'(Number of Julian centuries since Jan 1, 2000, 12 UT)
'=(JD-2451545)/36525
End Function

Function XLMod(a, b) As Double
    XLMod = a - (b * Int(a / b))
    'XLMod = a - (b * (a \ b))

'Equivalent formula
'The expression a Mod b is equivalent to either of the following formulas:
'XLMod = a - (b * (a \ b))
'XLMod = a - (b * Fix(a / b))
'XLMod = a - (b * Int(a / b))

End Function

Function LST(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, Longitude As Double) As Double 'Mean local sidereal time 'looks like in degrees

Dim JDay
    JDay = MoonFunctions.JDay(time, timezone, dte)
    
Dim T
    T = JCsince2000(time, timezone, dte)
    
LST = XLMod((280.46061837 + 360.98564736629 * (JDay - 2451545) + 0.000387933 * T ^ 2 - T ^ 3 / 38710000 + Longitude), _
                360)
    

'LST (mean local sidereal time) (S2)    
End Function

Function MoonEclipticLongLat(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, LongLat As Integer, Optional RadDeg As Integer) As Double

Dim T
    T = JCsince2000(time, timezone, dte)
    
Dim meanMoonLong
    meanMoonLong = 0.60643382 + 1336.85522467 * T - 0.00000313 * T ^ 2 - Int(0.60643382 + 1336.85522467 * T - 0.00000313 * T ^ 2)
    
Dim L
    L = 2 * WorksheetFunction.Pi _
    * (0.374897 + 1325.55241 * T - Int(0.374897 + 1325.55241 * T))
    
Dim D
    D = 2 * WorksheetFunction.Pi _
    * (0.827361 + 1236.853086 * T - Int(0.827361 + 1236.853086 * T))
    
Dim F
    F = 2 * WorksheetFunction.Pi _
    * (0.259086 + 1342.227825 * T - Int(0.259086 + 1342.227825 * T))

Dim LS
    LS = 2 * WorksheetFunction.Pi _
    * (0.99312619 + 99.99735956 * T - 0.00000044 * T ^ 2 - Int(0.99312619 + 99.99735956 * T - 0.00000044 * T ^ 2))

Dim DL
    DL = 22640 * Sin(L) - 4586 * Sin(L - 2 * D) + 2370 * Sin(2 * D) + 769 * Sin(2 * L) - 668 * Sin(LS) - 412 * Sin(2 * F) - 212 * Sin(2 * L - 2 * D) - 206 * Sin(L + LS - 2 * D) + 192 * Sin(L + 2 * D) - 165 * Sin(LS - 2 * D) - 125 * Sin(D) - 110 * Sin(L + LS) + 148 * Sin(L - LS) - 55 * Sin(2 * F - 2 * D)


Dim S
    S = F + (DL + 412 * Sin(2 * F) + 541 * Sin(LS)) / 206264.8062
    
Dim H
    H = F - 2 * D
    
Dim N
    N = -526 * Sin(H) + 44 * Sin(L + H) - 31 * Sin(-L + H) - 23 * Sin(LS + H) + 11 * Sin(-LS + H) - 25 * Sin(-2 * L + F) + 21 * Sin(-L + F)

If (IsMissing(RadDeg)) Then RadDeg = 1
Select Case LongLat
    Case 1 'longitude
        Dim MoonEclipticLongRad
            MoonEclipticLongRad = 2 * WorksheetFunction.Pi * (meanMoonLong + DL / 1296000 - Int(meanMoonLong + DL / 1296000))
    
        Select Case RadDeg
            Case 1 'radians
                MoonEclipticLongLat = MoonEclipticLongRad
            Case 2 'degrees
                MoonEclipticLongLat = WorksheetFunction.Degrees(MoonEclipticLongRad)
        End Select
    
    
    Case 2 'latitude
        Dim MoonEclipticLatRad
            MoonEclipticLatRad = (18520 * Sin(S) + N) / 206264.8062
                
            Select Case RadDeg
                Case 1 'radians
                    MoonEclipticLongLat = MoonEclipticLatRad
                Case 2 'degrees
                    MoonEclipticLongLat = WorksheetFunction.Degrees(MoonEclipticLatRad)
            End Select
    
    Case 3 'Equatorial Horizontal Parallax
        Select Case RadDeg
            Case 1 'degrees
                MoonEclipticLongLat = 0.950724 + 0.051818 * Cos(L) + 0.009531 * Cos(2 * D - L) + 0.007843 * Cos(2 * D) + 0.002824 * Cos(2 * L) + 0.000857 * Cos(2 * D + L) + 0.000533 * Cos(2 * D - LS) + 0.000401 * Cos(2 * D - LS - L) + 0.00032 * Cos(L - LS) - 0.000271 * Cos(D)
        End Select
End Select
     
End Function

Function MoonDecRA(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, DecRA As Integer) As Double

Dim T
    T = JCsince2000(time, timezone, dte)

Dim MoonEclipticLong
    MoonEclipticLong = MoonEclipticLongLat(time, timezone, dte, 1, 1)
    
Dim MoonEclipticLat
    MoonEclipticLat = MoonEclipticLongLat(time, timezone, dte, 2, 1)

Dim X
    X = Cos(MoonEclipticLong) * Cos(MoonEclipticLat)
        
Dim V
    V = Sin(MoonEclipticLong) * Cos(MoonEclipticLat)
    
Dim W
    W = Sin(MoonEclipticLat)
    
Dim Y
    Y = Cos(WorksheetFunction.Radians(23.4393 - 46.815 * T / 3600)) * V - Sin(WorksheetFunction.Radians(23.4393 - 46.815 * T / 3600)) * W

Dim Z
    Z = Sin(WorksheetFunction.Radians(23.4393 - 46.815 * T / 3600)) * V + Cos(WorksheetFunction.Radians(23.4393 - 46.815 * T / 3600)) * W

Dim RHO
    RHO = Sqr(1 - Z ^ 2)


Select Case DecRA
        
        Case 1 'Declination in degrees
                MoonDecRA = WorksheetFunction.Degrees(Atn(Z / RHO))
            
        Case 2 'RA in hours
                If 24 * Atn(Y / (X + RHO)) / WorksheetFunction.Pi > 0 Then
                        MoonDecRA = 24 * Atn(Y / (X + RHO)) / WorksheetFunction.Pi
                Else:   MoonDecRA = 24 * Atn(Y / (X + RHO)) / WorksheetFunction.Pi + 24
                End If
End Select    
End Function

Function LHA(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, Longitude As Double) As Double 'Local Hour Angle? 'in degrees

Dim MoonRA
    MoonRA = MoonDecRA(time, timezone, dte, 2)
    
Dim LST
    LST = MoonFunctions.LST(time, timezone, dte, Longitude)

If LST - 15 * MoonRA > 0 Then
        LHA = LST - 15 * MoonRA
Else:   LHA = 360 + LST - 15 * MoonRA
End If

End Function

Function EquatorialHorizontalParallax(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double) As Double 'in degrees
         EquatorialHorizontalParallax = MoonEclipticLongLat(time, timezone, dte, 3, 1)

End Function

Function MoonElevation(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double 'in degrees (elev1 - par)


Dim MoonDec
    MoonDec = WorksheetFunction.Radians(MoonDecRA(time, timezone, dte, 1))

    Latitude = WorksheetFunction.Radians(Latitude)
    
Dim elev1
    elev1 = WorksheetFunction.Degrees( _
        WorksheetFunction.Asin( _
        Cos(Latitude) _
        * Cos(MoonDec) _
        * Cos(WorksheetFunction.Radians(LHA(time, timezone, dte, Longitude))) _
        + Sin(Latitude) * Sin(MoonDec) _
        ))
    
Dim par
    par = WorksheetFunction.Degrees( _
    WorksheetFunction.Asin( _
    (0.9983271 + 0.0016764 * Cos(2 * Latitude)) _
    * (Cos(WorksheetFunction.Radians(elev1)) _
    * Sin(WorksheetFunction.Radians(EquatorialHorizontalParallax(time, timezone, dte)))) _
    ))
        
MoonElevation = elev1 - par
            
End Function

Function MoonElevationRefr(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double 'in degrees

Dim MoonElevation
    MoonElevation = MoonFunctions.MoonElevation(time, timezone, dte, Latitude, Longitude)

MoonElevationRefr = MoonElevation + 1.02 / (Tan(WorksheetFunction.Radians(MoonElevation + 10.3 / (MoonElevation + 5.11))) * 60)
    
End Function

Function MoonElongation(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double) As Double 'in degrees

Dim T
    T = JCsince2000(time, timezone, dte)
        
Dim LS
    LS = 2 * WorksheetFunction.Pi _
    * (0.99312619 + 99.99735956 * T - 0.00000044 * T ^ 2 - Int(0.99312619 + 99.99735956 * T - 0.00000044 * T ^ 2))

Dim Lsun
    If 280.4664567 + 360007.6982779 * T / 10 + 0.03032028 * T ^ 2 / 100 + T ^ 3 / 49931000 < 0 Then
            Lsun = XLMod(280.4664567 + 360007.6982779 * T / 10 + 0.03032028 * T ^ 2 / 100 + T ^ 3 / 49931000 + 360, 360)
    Else:   Lsun = XLMod(280.4664567 + 360007.6982779 * T / 10 + 0.03032028 * T ^ 2 / 100 + T ^ 3 / 49931000, 360)
    End If

Dim LST2
    LST2 = Lsun + (1.9146 - 0.004817 * T - 0.000014 * T ^ 2) * Sin(LS) + (0.019993 - 0.000101 * T) * Sin(2 * LS) + 0.00029 * Sin(3 * LS)


MoonElongation = WorksheetFunction.Degrees( _
                    WorksheetFunction.Acos( _
                    Cos(MoonEclipticLongLat(time, timezone, dte, 1, 1) - WorksheetFunction.Radians(LST2)) * Cos(MoonEclipticLongLat(time, timezone, dte, 2, 1)) _
                    ))

End Function

Function MoonIllumination(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double) As Double 'in percentage

Dim T
    T = JCsince2000(time, timezone, dte)

Dim MoonElongation
    MoonElongation = MoonFunctions.MoonElongation(time, timezone, dte)

Dim LS
    LS = 2 * WorksheetFunction.Pi _
    * (0.99312619 + 99.99735956 * T - 0.00000044 * T ^ 2 - Int(0.99312619 + 99.99735956 * T - 0.00000044 * T ^ 2))

Dim Lsun
    If 280.4664567 + 360007.6982779 * T / 10 + 0.03032028 * T ^ 2 / 100 + T ^ 3 / 49931000 < 0 Then
            Lsun = XLMod(280.4664567 + 360007.6982779 * T / 10 + 0.03032028 * T ^ 2 / 100 + T ^ 3 / 49931000 + 360, 360)
    Else:   Lsun = XLMod(280.4664567 + 360007.6982779 * T / 10 + 0.03032028 * T ^ 2 / 100 + T ^ 3 / 49931000, 360)
    End If

Dim LST2
    LST2 = Lsun + (1.9146 - 0.004817 * T - 0.000014 * T ^ 2) * Sin(LS) + (0.019993 - 0.000101 * T) * Sin(2 * LS) + 0.00029 * Sin(3 * LS)

Dim i
    i = 180 - MoonElongation - 0.1468 * (1 - 0.0549 * Sin(LS)) * Sin(WorksheetFunction.Radians(MoonElongation)) / (1 - 0.0167 * Sin(WorksheetFunction.Radians(LST2)))

MoonIllumination = 100 * (1 + Cos(WorksheetFunction.Radians(i))) / 2

End Function

Function MoonAzimuth(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double 'in degrees
    
Dim LHA
    LHA = WorksheetFunction.Radians(MoonFunctions.LHA(time, timezone, dte, Longitude))
    
Dim MoonDec
    MoonDec = MoonDecRA(time, timezone, dte, 1)
    
Dim one
    one = Sin(LHA)
    
Dim two
    two = Cos(LHA) * Sin(WorksheetFunction.Radians(Latitude)) - Tan(WorksheetFunction.Radians(MoonDec)) * Cos(WorksheetFunction.Radians(Latitude))
    
Dim arctan2of1and2
    If (one * two > 0) Or (one < 0 And two > 0) Then
            arctan2of1and2 = XLMod(WorksheetFunction.Degrees(WorksheetFunction.Atan2(two, one)) + 360, 360)
    Else:   arctan2of1and2 = WorksheetFunction.Degrees(WorksheetFunction.Atan2(two, one))
    End If
    
If WorksheetFunction.Degrees(LHA) > 180 Then
        MoonAzimuth = arctan2of1and2 - 180
Else:   MoonAzimuth = arctan2of1and2 + 180
End If

End Function

Function MoonDiameter(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double 'in arc minutes

Dim MoonElev
    MoonElev = WorksheetFunction.Radians(MoonElevation(time, timezone, dte, Latitude, Longitude))
    
Dim EquHorPar
    EquHorPar = WorksheetFunction.Radians(EquatorialHorizontalParallax(time, timezone, dte))

MoonDiameter = WorksheetFunction.Degrees((1 + Sin(MoonElev) * Sin(EquHorPar)) * 120 * WorksheetFunction.Asin(0.272481 * Sin(EquHorPar)))


End Function

Function MoonDistance(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double) As Double 'in KM

Dim T
    T = JCsince2000(time, timezone, dte)
    
Dim L
    L = 2 * WorksheetFunction.Pi _
    * (0.374897 + 1325.55241 * T - Int(0.374897 + 1325.55241 * T))
    
Dim D
    D = 2 * WorksheetFunction.Pi _
    * (0.827361 + 1236.853086 * T - Int(0.827361 + 1236.853086 * T))
    
Dim F
    F = 2 * WorksheetFunction.Pi _
    * (0.259086 + 1342.227825 * T - Int(0.259086 + 1342.227825 * T))

Dim LS
    LS = 2 * WorksheetFunction.Pi _
    * (0.99312619 + 99.99735956 * T - 0.00000044 * T ^ 2 - Int(0.99312619 + 99.99735956 * T - 0.00000044 * T ^ 2))

Dim one
    one = -20905355 * Cos(L) _
        - 3699111 * Cos(2 * D - L) _
        - 2955968 * Cos(2 * D) _
        - 569925 * Cos(2 * L) _
        + (1 - 0.002516 * T) * 48888 * Cos(LS) _
        - 3149 * Cos(2 * F) _
        + 246158 * Cos(2 * D - 2 * L) _
        - (1 - 0.002516 * T) * 152138 * Cos(2 * D - LS - L) _
        - 170733 * Cos(2 * D + L) _
        - (1 - 0.002516 * T) * 204586 * Cos(2 * D - LS) _
        - (1 - 0.002516 * T) * 129620 * Cos(LS - L) _
        + 108743 * Cos(D) _
        + (1 - 0.002516 * T) * 104755 * Cos(LS + L) _
        + 10321 * Cos(2 * D - 2 * F) _
        + 79661 * Cos(L - 2 * F) _
        - 34782 * Cos(4 * D - L) _
        - 23210 * Cos(3 * L) _
        - 21636 * Cos(4 * D - 2 * L) _
        + (1 - 0.002516 * T) * 24208 * Cos(2 * D + LS - L) _
        + (1 - 0.002516 * T) * 30824 * Cos(2 * D + LS) _
        - 8379 * Cos(D - L) _
        - (1 - 0.002516 * T) * 16675 * Cos(D + LS) _
    
Dim two
    two = -(1 - 0.002516 * T) * 12831 * Cos(2 * D - LS + L) _
        - 10445 * Cos(2 * D + 2 * L) - 11650 * Cos(4 * D) _
        + 14403 * Cos(2 * D - 3 * L) _
        - (1 - 0.002516 * T) * 7003 * Cos(LS - 2 * L) _
        + (1 - 0.002516 * T) * 10056 * Cos(2 * D - LS - 2 * L) _
        + 6322 * Cos(D + L) _
        - (1 - 0.002516 * T) * (1 - 0.002516 * T) * 9884 * Cos(2 * D - 2 * LS) _
        + (1 - 0.002516 * T) * 5751 * Cos(LS + 2 * L) _
        - (1 - 0.002516 * T) * (1 - 0.002516 * T) * 4950 * Cos(2 * D - 2 * LS - L) _
        + 4130 * Cos(2 * D + L - 2 * F) _
        - (1 - 0.002516 * T) * 3958 * Cos(4 * D - LS - L) _
        + 3258 * Cos(3 * D - L) _
        + (1 - 0.002516 * T) * 2616 * Cos(2 * D + LS + L) _
        - (1 - 0.002516 * T) * 1897 * Cos(4 * D - LS - 2 * L) _
        - (1 - 0.002516 * T) * (1 - 0.002516 * T) * 2117 * Cos(2 * LS - L) _
        + (1 - 0.002516 * T) * (1 - 0.002516 * T) * 2354 * Cos(2 * D + 2 * LS - L) _
        - 1423 * Cos(4 * D + L) - 1117 * Cos(4 * L) _
        - (1 - 0.002516 * T) * 1571 * Cos(4 * D - LS) _
        - 1739 * Cos(D - 2 * L) _
        - 4421 * Cos(2 * L - 2 * F) _
        + (1 - 0.002516 * T) * (1 - 0.002516 * T) * 1165 * Cos(2 * LS + L) _
        + 8752 * Cos(2 * D - L - 2 * F) _
    
Dim three
    three = 385000.56 + (one + two) / 1000
    
    
MoonDistance = three
    
'R / km (BE2) R stands for "Range" I believe, JPL uses "range" as well instead of the word "distance".

End Function

Module Name: SunFunctions

Function SunDecRA(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, DecRA As Integer)

Dim T
    T = JCsince2000(time, timezone, dte)
    
Dim M
    M = XLMod(357.5291 + 35999.0503 * T - 0.0001559 * T ^ 2 - 0.00000048 * T ^ 3, 360)

Dim L0
    L0 = XLMod(280.46645 + 36000.76983 * T + 0.0003032 * T ^ 2, 360)

Dim DL
    DL = XLMod((1.9146 - 0.004817 * T - 0.000014 * T ^ 2) * Sin(WorksheetFunction.Radians(M)) + (0.019993 - 0.000101 * T) * Sin(WorksheetFunction.Radians(2 * M)) + 0.00029 * Sin(WorksheetFunction.Radians(3 * M)), 360)

Dim L
    L = XLMod(L0 + DL, 360)

Dim X
    X = Cos(WorksheetFunction.Radians(L))

Dim Y
    Y = Cos(WorksheetFunction.Radians(23.4393 - 46.815 * T / 3600)) * Sin(WorksheetFunction.Radians(L))

Dim Z
    Z = Sin(WorksheetFunction.Radians(23.4393 - 46.815 * T / 3600)) * Sin(WorksheetFunction.Radians(L))

Dim R
    R = Sqr(1 - Z ^ 2)

Select Case DecRA

    Case 1 'Declination in degrees
    SunDecRA = WorksheetFunction.Degrees(Atn(Z / R))
    
    Case 2 'RA in degrees
    If 2 * WorksheetFunction.Degrees(Atn(Y / (X + R))) > 0 Then
            SunDecRA = 2 * WorksheetFunction.Degrees(Atn(Y / (X + R)))
    Else:   SunDecRA = 2 * WorksheetFunction.Degrees(Atn(Y / (X + R))) + 360
    End If

End Select

End Function

Function SunAzimuth(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double

Dim T
    T = JCsince2000(time, timezone, dte)

Dim JDay
    JDay = MoonFunctions.JDay(time, timezone, dte)

Dim SunRA 'degrees
    SunRA = SunDecRA(time, timezone, dte, 2)
    
Dim SunDec 'radians
    SunDec = WorksheetFunction.Radians(SunDecRA(time, timezone, dte, 1))

Dim theta 'degrees
    theta = XLMod(280.46061837 + 360.98564736629 * (JDay - 2451545) + 0.000387933 * T ^ 2 - T ^ 3 / 3871000010# + Longitude, 360)

Dim tau 'radians
    If theta - SunRA > 0 Then
            tau = WorksheetFunction.Radians(theta - SunRA)
    Else:   tau = WorksheetFunction.Radians(theta - SunRA + 360)
    End If

Dim sinofh
    sinofh = Sin(WorksheetFunction.Radians(Latitude)) * Sin(SunDec) + Cos(WorksheetFunction.Radians(Latitude)) * Cos(SunDec) * Cos(tau)

Dim one
    one = Sin(tau)

Dim two
    two = Cos(tau) * Sin(WorksheetFunction.Radians(Latitude)) - Tan(SunDec) * Cos(WorksheetFunction.Radians(Latitude))

Dim arctan2of1and2
    If (one * two > 0) Or (one < 0 And two > 0) Then
            arctan2of1and2 = XLMod(WorksheetFunction.Degrees(WorksheetFunction.Atan2(two, one)) + 360, 360)
    Else:   arctan2of1and2 = WorksheetFunction.Degrees(WorksheetFunction.Atan2(two, one))
    End If

If WorksheetFunction.Degrees(tau) > 180 Then
        SunAzimuth = arctan2of1and2 - 180
Else:   SunAzimuth = arctan2of1and2 + 180
End If
End Function

Function SunElevation(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double

Dim T
    T = JCsince2000(time, timezone, dte)

Dim JDay
    JDay = MoonFunctions.JDay(time, timezone, dte)

Dim SunDec 'radians
    SunDec = WorksheetFunction.Radians(SunDecRA(time, timezone, dte, 1))

Dim SunRA 'degrees
    SunRA = SunDecRA(time, timezone, dte, 2)

Dim theta 'degrees
    theta = XLMod(280.46061837 + 360.98564736629 * (JDay - 2451545) + 0.000387933 * T ^ 2 - T ^ 3 / 3871000010# + Longitude, 360)

Dim tau
    If theta - SunRA > 0 Then
            tau = WorksheetFunction.Radians(theta - SunRA)
    Else:   tau = WorksheetFunction.Radians(theta - SunRA + 360)
    End If

Dim sinofh
    sinofh = Sin(WorksheetFunction.Radians(Latitude)) * Sin(SunDec) + Cos(WorksheetFunction.Radians(Latitude)) * Cos(SunDec) * Cos(tau)


    SunElevation = WorksheetFunction.Degrees(WorksheetFunction.Asin(sinofh))

End Function

Function SunElevationRefr(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double

Dim SunElevation 'degrees
    SunElevation = SunFunctions.SunElevation(time, timezone, dte, Latitude, Longitude)

SunElevationRefr = SunElevation + 1.02 / (Tan(WorksheetFunction.Radians(SunElevation + 10.3 / (SunElevation + 5.11))) * 60)

End Function

Module Name: YallopFormulas

Function DAZ(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double 'in degrees

Dim SunAzimuth 'in degrees
    SunAzimuth = SunFunctions.SunAzimuth(time, timezone, dte, Latitude, Longitude)

Dim MoonAzimuth 'in degrees
    MoonAzimuth = MoonFunctions.MoonAzimuth(time, timezone, dte, Latitude, Longitude)

DAZ = SunAzimuth - MoonAzimuth 'in degrees

End Function

Function ArcV(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double 'in degrees

Dim ArcL 'in radians
    ArcL = WorksheetFunction.Radians(MoonFunctions.MoonElongation(time, timezone, dte))
    
Dim DAZ 'in radians
    DAZ = WorksheetFunction.Radians(YallopFormulas.DAZ(time, timezone, dte, Latitude, Longitude))

ArcV = WorksheetFunction.Degrees(WorksheetFunction.Acos(Cos(ArcL) / Cos(DAZ))) 'in degrees

End Function

Function CrescentWidth(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double 'topographic 'in arcminutes

Dim SemiDiameter 'topoghrapic 'in arcmin
    SemiDiameter = MoonFunctions.MoonDiameter(time, timezone, dte, Latitude, Longitude) / 2

Dim ArcL 'in radians
    ArcL = WorksheetFunction.Radians(MoonFunctions.MoonElongation(time, timezone, dte))

CrescentWidth = SemiDiameter * (1 - Cos(ArcL))

End Function

Function YallopQValue(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As Double

Dim ArcV
    ArcV = YallopFormulas.ArcV(time, timezone, dte, Latitude, Longitude)
    
Dim CrescentWidth
    CrescentWidth = YallopFormulas.CrescentWidth(time, timezone, dte, Latitude, Longitude)

Dim MoonElevation
    MoonElevation = MoonFunctions.MoonElevation(time, timezone, dte, Latitude, Longitude)

If MoonElevation < -0.833 Then
        YallopQValue = 999 '999 = moon has set
Else:   YallopQValue = (ArcV - (11.8371 - 6.3226 * CrescentWidth + 0.7319 * CrescentWidth ^ 2 - 0.1081 * CrescentWidth ^ 3)) / 10
End If

'=(L3-(11.8371-6.3226*N3+0.7319*N3^2-0.1081*N3^3))/10

End Function

Function YallopQCode(ByVal time As Double, ByVal timezone As Integer, ByVal dte As Double, ByVal Latitude As Double, ByVal Longitude As Double) As String

Dim YallopQValue
    YallopQValue = YallopFormulas.YallopQValue(time, timezone, dte, Latitude, Longitude)

Select Case YallopQValue
       
       Case Is = 999
            YallopQCode = "MhS" '= moon has set
               
       Case Is > 0.216
            YallopQCode = "A"
       
       Case Is > -0.014
            YallopQCode = "B"
       
       Case Is > -0.16
            YallopQCode = "C"
         
       Case Is > -0.232
            YallopQCode = "D"
            
       Case Is > -0.293
            YallopQCode = "E"
            
       Case Else
            YallopQCode = "F"
            
End Select

End Function

Have a look at this. I suspect you are experiencing the same bug — Cristian Buse, Commented May 19, 2022 at 18:35
I'm getting "statement too complex" error when importing the two files mentioned. — hidp123, Commented May 20, 2022 at 11:54
Are you importing modules or copy-pasting? Try downloading the ZIP and import the files via the VBE interface. Or clone the repo locally and then import from the clone. — Cristian Buse, Commented May 20, 2022 at 11:57
It works fine for me from the Zip download. Not sure what method you are using to download the file. Use the ZIp or clone the repo. It should work — Cristian Buse, Commented May 20, 2022 at 15:56
Thanks, this has improved speed a lot a lot! Can it be made even faster, as I'm dealing with lots of rows and coloumns. There's a speed issue I've noticed when using column filters, the calculating % at the bottom is very slow. I'll post issue on GitHub. — hidp123, Commented May 20, 2022 at 17:44

FreeMan · Accepted Answer · 2022-05-20 14:36:47Z

At a bare minimum declare your variable types when Dimming variables^*. If you don't, they're automatically Variant. The VBA runtime must then evaluate the variable each and every time it's accessed to determine what kind of value is currently held in the variable. This adds unnecessary processing time.

Unless, of course the variable content type is changing. If it's truly changing, then reevaluate whether you should be using the same variable to hold different data types (hint: you shouldn't except in rare cases) and rewrite the code to have different variables for different types.

Why are you recalculating historical data? Do the sun and moon positions for yesterday, last week, last year change every time you open the workbook? Calculate the historical data once and store it, then run calculations only on new data (daily, weekly, monthly) and it will run much faster.

Maybe have a [Historical Data] worksheet where all the inputs and the outputs are stored. Then have a [Current Data] worksheet where you enter, copy/paste, whatever, an amount of new data. When you've collected a sufficient amount, kick off the calculations. The last step of the calculations would be to Copy, then Paste Special | Values the newly calculated data to the appropriate location in [Historical Data].

Maybe you enter data daily, but need to be able to see the last 30 day's worth of data. Each day could recalculate the last 30 day's worth, then you add today's data and part of the processing deletes the (previously copy/pasted) 31st day's data from the [Current Data] worksheet.

No, this doesn't really help speed up the code, but by rethinking the process, you may well be able to avoid worrying about slow code by minimizing code execution in the first place.

^*NOTE: it's always a good idea to declare data types when declaring variables so the next programmer (possibly future you) knows what the intent is. In this case of slow processing, the time to determine the type at runtime was the key factor involved.

Rick James · Accepted Answer · 2022-05-21 01:26:33Z

Does VBA had a "MOD" or "Fraction" function? Or "common subexpression elimination"? I'm looking at the inefficiency of

meanMoonLong = 0.60643382 + 1336.85522467 * T - 0.00000313 * T ^ 2
         - Int(0.60643382 + 1336.85522467 * T - 0.00000313 * T ^ 2)

Also, unless VBA special cases power-2, T ^ 2 should be written T * T

one and two have (1 - 0.002516 * T) a lot. Pull that out into a temp variable. (That may help with the "too complicated" error.)

Stack Exchange Network

VBA functions to calculate sun and moon positions

2 Answers 2

Not the answer you're looking for? Browse other questions tagged
performance
vba
excel
computational-geometry
physics
or ask your own question.

Hot Network Questions

VBA functions to calculate sun and moon positions

2 Answers 2

Not the answer you're looking for? Browse other questions tagged performancevbaexcelcomputational-geometryphysics or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
performance
vba
excel
computational-geometry
physics
or ask your own question.