diff --git a/src/InsolationCalc.jl b/src/InsolationCalc.jl
index 48c78aa..d123935 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,
@@ -42,8 +44,14 @@ function solar_flux_and_cos_sza(
     S0::FT = IP.tot_solar_irrad(param_set)
     d0::FT = IP.orbit_semimaj(param_set)
     # θ = solar zenith angle, ζ = solar azimuth angle, d = earth-sun distance
-    θ, ζ, d =
-        instantaneous_zenith_angle(date, od, longitude, latitude, param_set)
+    θ, ζ, d = instantaneous_zenith_angle(
+        date,
+        date0,
+        od,
+        longitude,
+        latitude,
+        param_set,
+    )
     # 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..f22b75f 100644
--- a/src/ZenithAngleCalc.jl
+++ b/src/ZenithAngleCalc.jl
@@ -30,10 +30,9 @@ compute_orbital_parameters(param_set) = (
 function get_Δt_years(
     param_set::IP.InsolationParameters{FT},
     date,
-    epoch_string = "2000-01-01T11:58:56.816",
+    date0,
 ) 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
@@ -41,8 +40,8 @@ end
 # calculate the distance, declination, and hour angle (at lon=0)
 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,7 +50,7 @@ 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")
+    #date0 = DateTime(epoch_string, dateformat"y-m-dTHH:MM:SS.s")
 
     days_per_year = Ya / day_length
     Δt_years =
@@ -86,7 +85,7 @@ end
 
 function instantaneous_zenith_angle(
     date::DateTime,
-    epoch_string,
+    date0::DateTime,
     (ϖ, γ, e)::Tuple{FT, FT, FT},
     longitude::FT,
     latitude::FT,
@@ -96,11 +95,10 @@ function instantaneous_zenith_angle(
     λ = deg2rad(longitude)
     ϕ = deg2rad(latitude)
 
-
     d, δ, η_UTC = distance_declination_hourangle(
         date,
+        date0,
         (ϖ, γ, e),
-        epoch_string,
         param_set,
         eot_correction,
     )
@@ -143,6 +141,7 @@ Orbital parameters are computed using the `Milankovitch` method.
 """
 function instantaneous_zenith_angle(
     date::DateTime,
+    date0::DateTime,
     od::OrbitalData,
     longitude::FT,
     latitude::FT,
@@ -150,12 +149,12 @@ function instantaneous_zenith_angle(
     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)
+    Δt_years = get_Δt_years(param_set, date, date0)
     instantaneous_zenith_angle(
         date,
-        epoch_string,
+        date0,
         compute_orbital_parameters(od, Δt_years),
         longitude,
         latitude,
@@ -185,17 +184,16 @@ Orbital parameters at the J2000 epoch from CLIMAParameters are used.
 """
 function instantaneous_zenith_angle(
     date::DateTime,
+    date0::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,
+        date0,
         compute_orbital_parameters(param_set),
         longitude,
         latitude,
@@ -228,6 +226,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 +235,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 +244,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..b0ddfb9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,9 @@ FT = Float32
 include(joinpath(pkgdir(Insolation), "parameters", "create_parameters.jl"))
 param_set = create_insolation_parameters(FT)
 
+epoch_string = "2000-01-01T11:58:56.816"
+date0 = DateTime(epoch_string, dateformat"y-m-dTHH:MM:SS.s")
+
 @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..0f20f3a 100644
--- a/test/test_insolation.jl
+++ b/test/test_insolation.jl
@@ -7,32 +7,32 @@ 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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 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 +42,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 +58,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 +75,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..cf70b2a 100644
--- a/test/test_types.jl
+++ b/test/test_types.jl
@@ -5,6 +5,7 @@ od = Insolation.OrbitalData()
 
 sza, azi, d = instantaneous_zenith_angle(
     date,
+    date0,
     lon,
     lat,
     param_set,
@@ -20,6 +21,7 @@ F = insolation(sza, d, param_set)
 
 sza, azi, d = instantaneous_zenith_angle(
     date,
+    date0,
     od,
     lon,
     lat,
@@ -33,8 +35,14 @@ 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)
+sza, azi, d = instantaneous_zenith_angle(
+    date,
+    date0,
+    lon,
+    lat,
+    param_set,
+    eot_correction = true,
+)
 
 F = insolation(sza, d, param_set)
 @test typeof(sza) == FT
@@ -44,6 +52,7 @@ F = insolation(sza, d, param_set)
 
 sza, azi, d = instantaneous_zenith_angle(
     date,
+    date0,
     od,
     lon,
     lat,
diff --git a/test/test_zenith_angle.jl b/test/test_zenith_angle.jl
index ed3a89e..bd301fc 100644
--- a/test/test_zenith_angle.jl
+++ b/test/test_zenith_angle.jl
@@ -3,23 +3,24 @@ 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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 @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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 @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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 @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,
+    date0,
     od,
     lon,
     lat,
@@ -32,6 +33,7 @@ sza, azi, d = instantaneous_zenith_angle(
 date = Dates.DateTime(2020, 2, 20, 18, 0, 0)
 sza, azi, d = instantaneous_zenith_angle(
     date,
+    date0,
     od,
     lon,
     lat,
@@ -44,26 +46,26 @@ sza, azi, d = instantaneous_zenith_angle(
 # 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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 @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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 @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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 @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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 @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)
+sza, azi, d = instantaneous_zenith_angle(date, date0, od, lon, lat, param_set)
 @test d ≈ IP.orbit_semimaj(param_set) rtol = rtol