diff --git a/Project.toml b/Project.toml
index 9a8c027..5d1282d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Insolation"
 uuid = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8"
 authors = ["Climate Modeling Alliance"]
-version = "0.7.0"
+version = "0.8.0"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
diff --git a/docs/src/find_equinox_perihelion_dates.jl b/docs/src/find_equinox_perihelion_dates.jl
index bcef68b..168d3a3 100644
--- a/docs/src/find_equinox_perihelion_dates.jl
+++ b/docs/src/find_equinox_perihelion_dates.jl
@@ -8,13 +8,27 @@ FT = Float64
 include(joinpath(pkgdir(Insolation), "parameters", "create_parameters.jl"))
 param_set = create_insolation_parameters(FT)
 
+date0 = DateTime("2000-01-01T11:58:56.816")
+
 # Difference in NH and SH zenith angles at time x in given year
 function zdiff(x, year, od)
     date = xtomarchdate(x, year)
-    theta_s, dist =
-        daily_zenith_angle(date, od, -45.0, param_set, milankovitch = true)
-    theta_n, dist =
-        daily_zenith_angle(date, od, 45.0, param_set, milankovitch = true)
+    theta_s, dist = daily_zenith_angle(
+        date,
+        date0,
+        od,
+        -45.0,
+        param_set,
+        milankovitch = true,
+    )
+    theta_n, dist = daily_zenith_angle(
+        date,
+        date0,
+        od,
+        45.0,
+        param_set,
+        milankovitch = true,
+    )
     return theta_n - theta_s
 end
 
@@ -28,7 +42,8 @@ end
 # Earth-Sun distance
 function edist(x, year, od)
     date = xtojandate(x, year)
-    _, dist = daily_zenith_angle(date, od, 0.0, param_set, milankovitch = true)
+    _, dist =
+        daily_zenith_angle(date, date0, od, 0.0, param_set, milankovitch = true)
     return dist / IP.orbit_semimaj(param_set)
 end
 
diff --git a/docs/src/library.md b/docs/src/library.md
index a5b6e57..73e7b5c 100644
--- a/docs/src/library.md
+++ b/docs/src/library.md
@@ -16,7 +16,7 @@ Insolation.OrbitalData
 ## Zenith Angle
 ```@autodocs
 Modules = [Insolation]
-Private = false
+Private = true
 Pages   = ["ZenithAngleCalc.jl"]
 ```
 
@@ -25,4 +25,4 @@ Pages   = ["ZenithAngleCalc.jl"]
 Modules = [Insolation]
 Private = false
 Pages   = ["InsolationCalc.jl"]
-```
\ No newline at end of file
+```
diff --git a/docs/src/plot_diurnal_cycle.jl b/docs/src/plot_diurnal_cycle.jl
index 2aa2a62..90f1185 100644
--- a/docs/src/plot_diurnal_cycle.jl
+++ b/docs/src/plot_diurnal_cycle.jl
@@ -15,11 +15,13 @@ function diurnal_cycle(lat, lon, date, od, timezone, filename)
     hours = collect(range(0, stop = 24, length = nhours))
     insol = zeros(nhours)
     sza = zeros(nhours)
+    date0 = DateTime("2000-01-01T11:58:56.816")
+
     for (i, hr) in enumerate(hours)
         h = Int(round(hr + timezone))
         m = Int(round((hr + timezone - h) * 60))
         datetime = date + Dates.Hour(h) + Dates.Minute(m)
-        S, mu = solar_flux_and_cos_sza(datetime, od, lon, lat, param_set)
+        S, mu = solar_flux_and_cos_sza(datetime, date0, od, lon, lat, param_set)
         insol[i] = S * mu
         sza[i] = rad2deg(acos(mu))
     end
diff --git a/docs/src/plot_insolation.jl b/docs/src/plot_insolation.jl
index 62c9c1a..1bbec54 100644
--- a/docs/src/plot_insolation.jl
+++ b/docs/src/plot_insolation.jl
@@ -22,12 +22,14 @@ function calc_day_lat_insolation(
     d_arr = Array{I}(round.(collect(range(0, stop = 365, length = n_days))))
     l_arr = collect(range(-90, stop = 90, length = n_lats))
     F_arr = zeros(n_days, n_lats)
+    date0 = DateTime("2000-01-01T11:58:56.816")
     # loop over days
     for (i, d) in enumerate(d_arr)
         for (j, lat) in enumerate(l_arr)
             date = Dates.DateTime(2000, 1, 1) + Dates.Day(d)
             θ, dist = daily_zenith_angle(
                 date,
+                date0,
                 od,
                 lat,
                 param_set,
diff --git a/src/InsolationCalc.jl b/src/InsolationCalc.jl
index 48c78aa..791e6a1 100644
--- a/src/InsolationCalc.jl
+++ b/src/InsolationCalc.jl
@@ -22,6 +22,7 @@ end
 
 """
     solar_flux_and_cos_sza(date::DateTime,
+                      date0::DateTime,
                       od::OrbitalData,
                       longitude::FT,
                       latitude::FT,
@@ -34,6 +35,7 @@ param_set is an AbstractParameterSet from CLIMAParameters.jl.
 """
 function solar_flux_and_cos_sza(
     date::DateTime,
+    date0::DateTime,
     od::OrbitalData,
     longitude::FT,
     latitude::FT,
@@ -41,9 +43,18 @@ function solar_flux_and_cos_sza(
 ) where {FT <: Real}
     S0::FT = IP.tot_solar_irrad(param_set)
     d0::FT = IP.orbit_semimaj(param_set)
+    args = (
+        Insolation.helper_instantaneous_zenith_angle(
+            date,
+            date0,
+            od,
+            param_set,
+        )...,
+        longitude,
+        latitude,
+    )
     # θ = solar zenith angle, ζ = solar azimuth angle, d = earth-sun distance
-    θ, ζ, d =
-        instantaneous_zenith_angle(date, od, longitude, latitude, param_set)
+    θ, ζ, d = instantaneous_zenith_angle(args...)
     # set max. zenith angle to π/2, insolation should not be negative
     if θ > FT(π) / 2
         θ = FT(π) / 2
diff --git a/src/ZenithAngleCalc.jl b/src/ZenithAngleCalc.jl
index 2113767..1b38847 100644
--- a/src/ZenithAngleCalc.jl
+++ b/src/ZenithAngleCalc.jl
@@ -29,20 +29,33 @@ compute_orbital_parameters(param_set) = (
 
 function get_Δt_years(
     param_set::IP.InsolationParameters{FT},
-    date,
-    epoch_string = "2000-01-01T11:58:56.816",
+    date::DateTime,
+    date0::DateTime,
 ) where {FT}
     (; year_anom, day) = param_set
-    date0 = DateTime(epoch_string, dateformat"y-m-dTHH:MM:SS.s")
     days_per_year = year_anom / day
     return FT(datetime2julian(date) - datetime2julian(date0)) / days_per_year
 end
 
-# calculate the distance, declination, and hour angle (at lon=0)
+"""
+    distance_declination_hourangle(
+        date::DateTime,
+        date0::DateTime,
+        (ϖ, γ, e)::Tuple{FT, FT, FT},
+        param_set::IP.AIP,
+        eot_correction::Bool,
+    ) where {FT}
+
+Returns the earth-sun distance (m), declination angle (radians) and hour angle 
+(radians), at 0ᵒ longitude, given the current datetime, epoch datetime, 
+`longitude of the perihelion at epoch`, `obliquity at epoch`, `eccentricity at epoch`, 
+and `eot_correction`.
+
+"""
 function distance_declination_hourangle(
     date::DateTime,
+    date0::DateTime,
     (ϖ, γ, e)::Tuple{FT, FT, FT},
-    epoch_string,
     param_set::IP.AIP,
     eot_correction::Bool,
 ) where {FT}
@@ -51,8 +64,6 @@ function distance_declination_hourangle(
     d0 = IP.orbit_semimaj(param_set)
     M0 = IP.mean_anom_epoch(param_set)
 
-    date0 = DateTime(epoch_string, dateformat"y-m-dTHH:MM:SS.s")
-
     days_per_year = Ya / day_length
     Δt_years =
         FT(datetime2julian(date) - datetime2julian(date0)) / days_per_year
@@ -84,27 +95,29 @@ function distance_declination_hourangle(
     return d, δ, η_UTC
 end
 
+"""
+    instantaneous_zenith_angle(
+        d::FT,
+        δ::FT,
+        η_UTC::FT,
+        longitude::FT,
+        latitude::FT,
+    ) where {FT}
+
+Returns the zenith angle (radians), azimuthal angle (radians) and Earth-Sun distance (m)
+at a particular longitude (degrees) and latitude (degrees) given Earth-Sun distance (m), 
+declination angle (radians), and hour angle (radians) at 0ᵒ longitude.
+
+"""
 function instantaneous_zenith_angle(
-    date::DateTime,
-    epoch_string,
-    (ϖ, γ, e)::Tuple{FT, FT, FT},
+    d::FT,
+    δ::FT,
+    η_UTC::FT,
     longitude::FT,
     latitude::FT,
-    param_set::IP.AIP,
-    eot_correction::Bool,
 ) where {FT}
     λ = deg2rad(longitude)
     ϕ = deg2rad(latitude)
-
-
-    d, δ, η_UTC = distance_declination_hourangle(
-        date,
-        (ϖ, γ, e),
-        epoch_string,
-        param_set,
-        eot_correction,
-    )
-
     # hour angle
     η = mod(η_UTC + λ, FT(2π))
 
@@ -121,92 +134,40 @@ function instantaneous_zenith_angle(
     return θ, ζ, d
 end
 
-"""
-    instantaneous_zenith_angle(date::DateTime,
-                               od::OrbitalData,
-                               longitude::FT,
-                               latitude::FT,
-                               param_set::IP.AIP;
-                               eot_correction::Bool=true,
-                               ) where {FT}
-
-Returns the zenith angle and earth-sun distance
-at a particular longitude and latitude on the given date (and time UTC)
-given orbital parameters: obliquity, longitude of perihelion, and eccentricity
-param_set is an AbstractParameterSet from CLIMAParameters.jl.
-
-`eot_correction` is an optional Boolean keyword argument that defaults to true
-when set to true the equation of time correction is turned on.
-This switch functionality is implemented for easy comparisons with reanalyses.
-
-Orbital parameters are computed using the `Milankovitch` method.
-"""
-function instantaneous_zenith_angle(
+function helper_instantaneous_zenith_angle(
     date::DateTime,
+    date0::DateTime,
     od::OrbitalData,
-    longitude::FT,
-    latitude::FT,
-    param_set::IP.AIP;
-    eot_correction::Bool = true,
-) where {FT}
-
-    epoch_string = "2000-01-01T11:58:56.816"
-    # epoch_string::String = IP.epoch(param_set)
-    Δt_years = get_Δt_years(param_set, date, epoch_string)
-    instantaneous_zenith_angle(
+    param_set::AIP;
+    eot_correction = true,
+)
+    Δt_years = get_Δt_years(param_set, date, date0)
+    return distance_declination_hourangle(
         date,
-        epoch_string,
-        compute_orbital_parameters(od, Δt_years),
-        longitude,
-        latitude,
+        date0,
+        Insolation.compute_orbital_parameters(od, Δt_years),
         param_set,
         eot_correction,
     )
 end
 
-"""
-    instantaneous_zenith_angle(date::DateTime,
-                               longitude::FT,
-                               latitude::FT,
-                               param_set::IP.AIP;
-                               eot_correction::Bool=true,
-                               ) where {FT}
-
-Returns the zenith angle and earth-sun distance
-at a particular longitude and latitude on the given date (and time UTC)
-given orbital parameters: obliquity, longitude of perihelion, and eccentricity
-param_set is an AbstractParameterSet from CLIMAParameters.jl.
-
-`eot_correction` is an optional Boolean keyword argument that defaults to true
-when set to true the equation of time correction is turned on.
-This switch functionality is implemented for easy comparisons with reanalyses.
-
-Orbital parameters at the J2000 epoch from CLIMAParameters are used.
-"""
-function instantaneous_zenith_angle(
+helper_instantaneous_zenith_angle(
     date::DateTime,
-    longitude::FT,
-    latitude::FT,
-    param_set::IP.AIP;
-    eot_correction::Bool = true,
-) where {FT}
-    epoch_string = "2000-01-01T11:58:56.816"
-    # epoch_string::String = IP.epoch(param_set)
-
-    return instantaneous_zenith_angle(
-        date,
-        epoch_string,
-        compute_orbital_parameters(param_set),
-        longitude,
-        latitude,
-        param_set,
-        eot_correction,
-    )
+    date0::DateTime,
+    param_set::AIP;
+    eot_correction = true,
+) = distance_declination_hourangle(
+    date,
+    date0,
+    Insolation.compute_orbital_parameters(param_set),
+    param_set,
+    eot_correction,
+)
 
-end
 
 """
     daily_zenith_angle(date::DateTime,
+                       data0::DateTime,
                        od::OrbitalData,
                        latitude::FT,
                        param_set::IP.AIP;
@@ -228,6 +189,7 @@ when set to false the orbital parameters at the J2000 epoch from CLIMAParameters
 """
 function daily_zenith_angle(
     date::DateTime,
+    date0::DateTime,
     od::OrbitalData,
     latitude::FT,
     param_set::IP.AIP;
@@ -236,9 +198,7 @@ function daily_zenith_angle(
 ) where {FT}
     ϕ = deg2rad(latitude)
 
-    epoch_string = "2000-01-01T11:58:56.816"
-    # epoch_string::String = IP.epoch(param_set)
-    Δt_years = get_Δt_years(param_set, date, epoch_string)
+    Δt_years = get_Δt_years(param_set, date, date0)
 
     # calculate orbital parameters or take values at J2000
     (ϖ, γ, e) =
@@ -247,8 +207,8 @@ function daily_zenith_angle(
 
     d, δ, _ = distance_declination_hourangle(
         date,
+        date0,
         (ϖ, γ, e),
-        epoch_string,
         param_set,
         eot_correction,
     )
diff --git a/test/runtests.jl b/test/runtests.jl
index 7c55d80..b147838 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,11 +10,15 @@ using Optim
 
 import CLIMAParameters as CP
 import Insolation.Parameters as IP
+const AIP = IP.AbstractInsolationParams
+import Insolation.OrbitalData
 FT = Float32
 
 include(joinpath(pkgdir(Insolation), "parameters", "create_parameters.jl"))
 param_set = create_insolation_parameters(FT)
 
+date0 = DateTime("2000-01-01T11:58:56.816")
+
 @testset "Orbital Params" begin
     include("test_orbit_param.jl")
 end
diff --git a/test/test_equinox.jl b/test/test_equinox.jl
index 9fe4a48..2549f9f 100644
--- a/test/test_equinox.jl
+++ b/test/test_equinox.jl
@@ -1,8 +1,8 @@
 # Difference in NH and SH zenith angles at time x in given year
 function zdiff(x, year, od)
     date = xtomarchdate(x, year)
-    theta_s, dist = daily_zenith_angle(date, od, FT(-45), param_set)
-    theta_n, dist = daily_zenith_angle(date, od, FT(45), param_set)
+    theta_s, dist = daily_zenith_angle(date, date0, od, FT(-45), param_set)
+    theta_n, dist = daily_zenith_angle(date, date0, od, FT(45), param_set)
     return theta_n - theta_s
 end
 
diff --git a/test/test_insolation.jl b/test/test_insolation.jl
index e206283..9ee81e4 100644
--- a/test/test_insolation.jl
+++ b/test/test_insolation.jl
@@ -7,32 +7,47 @@ od = Insolation.OrbitalData()
 # sunrise at equator
 date = Dates.DateTime(2020, 1, 1, 6, 0, 0)
 lon, lat = [FT(0.0), FT(0.0)]
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 F = insolation(sza, d, param_set)
 @test F ≈ 0.0 atol = atol
 
-S, mu = solar_flux_and_cos_sza(date, od, lon, lat, param_set)
+S, mu = solar_flux_and_cos_sza(date, date0, od, lon, lat, param_set)
 @test S ≈ IP.tot_solar_irrad(param_set) rtol = rtol_insol
 @test mu ≈ 0.0 atol = atol
 
 # polar night NH 1
 date = Dates.DateTime(2020, 12, 20, 11, 0, 0)
 lon, lat = [FT(0.0), FT(80.0)]
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 F = insolation(sza, d, param_set)
 @test F ≈ 0.0 atol = atol
 
-S, mu = solar_flux_and_cos_sza(date, od, lon, lat, param_set)
+S, mu = solar_flux_and_cos_sza(date, date0, od, lon, lat, param_set)
 @test S ≈ IP.tot_solar_irrad(param_set) rtol = rtol_insol
 @test mu ≈ 0.0 atol = atol
 
 # polar night NH 2
 date = Dates.DateTime(2020, 12, 20, 23, 0, 0)
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 F = insolation(sza, d, param_set)
 @test F ≈ 0.0 atol = atol
 
-S, mu = solar_flux_and_cos_sza(date, od, lon, lat, param_set)
+S, mu = solar_flux_and_cos_sza(date, date0, od, lon, lat, param_set)
 @test S ≈ IP.tot_solar_irrad(param_set) rtol = rtol_insol
 @test mu ≈ 0.0 atol = atol
 
@@ -42,7 +57,7 @@ date = Dates.DateTime(2021, 3, 20, 9, 37, 0) # vernal equinox 2021
 l_arr = FT.(collect(range(-90, stop = 90, length = nlats)))
 F_arr = zeros(nlats)
 for (i, lat) in enumerate(l_arr)
-    θ, dist = daily_zenith_angle(date, od, lat, param_set)
+    θ, dist = daily_zenith_angle(date, date0, od, lat, param_set)
     F_arr[i] = insolation(θ, dist, param_set)
 end
 F_NH = sort(F_arr[l_arr .>= 0])
@@ -58,7 +73,7 @@ F_arr = zeros(ndays, nlats)
 for (i, d) in enumerate(d_arr)
     for (j, lat) in enumerate(l_arr)
         datei = Dates.DateTime(2000, 1, 1) + Dates.Day(d)
-        θ, dist = daily_zenith_angle(datei, od, lat, param_set)
+        θ, dist = daily_zenith_angle(datei, date0, od, lat, param_set)
         F_arr[i, j] = insolation(θ, dist, param_set)
     end
 end
@@ -75,7 +90,7 @@ param_set = create_insolation_parameters(FT, (; lon_perihelion_epoch = ϖ0 + π)
 for (i, d) in enumerate(d_arr)
     for (j, lat) in enumerate(l_arr)
         datei = Dates.DateTime(2000, 1, 1) + Dates.Day(d)
-        θ, dist = daily_zenith_angle(datei, od, lat, param_set)
+        θ, dist = daily_zenith_angle(datei, date0, od, lat, param_set)
         F_arr[i, j] = insolation(θ, dist, param_set)
     end
 end
diff --git a/test/test_perihelion.jl b/test/test_perihelion.jl
index ffd7b99..74585d0 100644
--- a/test/test_perihelion.jl
+++ b/test/test_perihelion.jl
@@ -9,7 +9,7 @@ end
 # Earth-Sun distance
 function edist(x, year, od)
     date = xtojandate(x, year)
-    _, dist = daily_zenith_angle(date, od, FT(0), param_set)
+    _, dist = daily_zenith_angle(date, date0, od, FT(0), param_set)
     return dist / IP.orbit_semimaj(param_set)
 end
 
diff --git a/test/test_types.jl b/test/test_types.jl
index ef017e0..6642514 100644
--- a/test/test_types.jl
+++ b/test/test_types.jl
@@ -3,14 +3,17 @@ lon, lat = [FT(80.0), FT(20.0)]
 
 od = Insolation.OrbitalData()
 
-sza, azi, d = instantaneous_zenith_angle(
-    date,
+args = (
+    Insolation.helper_instantaneous_zenith_angle(
+        date,
+        date0,
+        param_set,
+        eot_correction = false,
+    )...,
     lon,
     lat,
-    param_set,
-    eot_correction = false,
 )
-
+sza, azi, d = instantaneous_zenith_angle(args...)
 
 F = insolation(sza, d, param_set)
 @test typeof(sza) == FT
@@ -18,14 +21,18 @@ F = insolation(sza, d, param_set)
 @test typeof(d) == FT
 @test typeof(F) == FT
 
-sza, azi, d = instantaneous_zenith_angle(
-    date,
-    od,
+args = (
+    Insolation.helper_instantaneous_zenith_angle(
+        date,
+        date0,
+        od,
+        param_set,
+        eot_correction = false,
+    )...,
     lon,
     lat,
-    param_set,
-    eot_correction = false,
 )
+sza, azi, d = instantaneous_zenith_angle(args...)
 
 F = insolation(sza, d, param_set)
 @test typeof(sza) == FT
@@ -33,8 +40,17 @@ F = insolation(sza, d, param_set)
 @test typeof(d) == FT
 @test typeof(F) == FT
 
-sza, azi, d =
-    instantaneous_zenith_angle(date, lon, lat, param_set, eot_correction = true)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(
+        date,
+        date0,
+        param_set,
+        eot_correction = true,
+    )...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 
 F = insolation(sza, d, param_set)
 @test typeof(sza) == FT
@@ -42,15 +58,18 @@ F = insolation(sza, d, param_set)
 @test typeof(d) == FT
 @test typeof(F) == FT
 
-sza, azi, d = instantaneous_zenith_angle(
-    date,
-    od,
+args = (
+    Insolation.helper_instantaneous_zenith_angle(
+        date,
+        date0,
+        od,
+        param_set,
+        eot_correction = true,
+    )...,
     lon,
     lat,
-    param_set,
-    eot_correction = true,
 )
-
+sza, azi, d = instantaneous_zenith_angle(args...)
 
 F = insolation(sza, d, param_set)
 @test typeof(sza) == FT
diff --git a/test/test_zenith_angle.jl b/test/test_zenith_angle.jl
index ed3a89e..f7cebd4 100644
--- a/test/test_zenith_angle.jl
+++ b/test/test_zenith_angle.jl
@@ -1,69 +1,118 @@
 rtol = 1e-2
 od = Insolation.OrbitalData()
+
 # sunrise at equator
 date = Dates.DateTime(2020, 2, 20, 6, 11, 0)
 lon, lat = [FT(0.0), FT(0.0)]
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test sza ≈ π / 2 rtol = rtol
 
 # solar noon at equator
 date = Dates.DateTime(2020, 2, 20, 12, 14, 0)
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test azi ≈ 3π / 2 rtol = rtol
 
 # sunset at equator
 date = Dates.DateTime(2020, 2, 20, 18, 17, 0)
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test sza ≈ π / 2 rtol = rtol
 
 # solar noon at equator, eot correction = false
 date = Dates.DateTime(2020, 2, 20, 12, 0, 0)
-sza, azi, d = instantaneous_zenith_angle(
-    date,
-    od,
+args = (
+    Insolation.helper_instantaneous_zenith_angle(
+        date,
+        date0,
+        od,
+        param_set,
+        eot_correction = false,
+    )...,
     lon,
     lat,
-    param_set;
-    eot_correction = false,
 )
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test azi ≈ 3π / 2 rtol = rtol
 
 # sunset at equator, eot correction = false
 date = Dates.DateTime(2020, 2, 20, 18, 0, 0)
-sza, azi, d = instantaneous_zenith_angle(
-    date,
-    od,
+args = (
+    Insolation.helper_instantaneous_zenith_angle(
+        date,
+        date0,
+        od,
+        param_set,
+        eot_correction = false,
+    )...,
     lon,
     lat,
-    param_set;
-    eot_correction = false,
 )
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test sza ≈ π / 2 rtol = rtol
 
 ## Test Polar Night
 # polar night NH 1
 date = Dates.DateTime(2020, 12, 20, 11, 0, 0)
 lon, lat = [FT(0.0), FT(80.0)]
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test sza > π / 2
 
 # polar night NH 2
 date = Dates.DateTime(2020, 12, 20, 23, 0, 0)
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test sza > π / 2
 
 # polar night SH 1
 date = Dates.DateTime(2020, 6, 20, 11, 0, 0)
 lon, lat = [FT(0.0), FT(-80.0)]
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test sza > π / 2
 
 # polar night SH 2
 date = Dates.DateTime(2020, 6, 20, 23, 0, 0)
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test sza > π / 2
 
 ## Test Distance
 date = Dates.DateTime(2000, 3, 22, 0, 0, 0)
-sza, azi, d = instantaneous_zenith_angle(date, od, lon, lat, param_set)
+args = (
+    Insolation.helper_instantaneous_zenith_angle(date, date0, od, param_set)...,
+    lon,
+    lat,
+)
+sza, azi, d = instantaneous_zenith_angle(args...)
 @test d ≈ IP.orbit_semimaj(param_set) rtol = rtol