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

• Cristian Buse: Have a look at this. I suspect you are experiencing the same bug
• hidp123: I'm getting "statement too complex" error when importing the two files mentioned.
• Cristian Buse: 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: 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
• hidp123: 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.

Answer 1

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.

Answer 2

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.)