Initial commit

src/ZenithAngleCalc.jl

@@ -14,6 +14,17 @@ function equation_of_time(e::FT, MA::FT, γ::FT, ϖ::FT) where {FT <: Real}
     return mod(_Δt+FT(π), FT(2π)) - FT(π)
 end
 
+# calculate orbital parameters (Milankovitch)
+compute_orbital_parameters(od, Δt_years::FT) where {FT} = (
+    FT(ϖ_spline(od, Δt_years)),
+    FT(γ_spline(od, Δt_years)),
+    FT(e_spline(od, Δt_years)))
+
+compute_orbital_parameters(param_set) = (
+    IP.lon_perihelion_epoch(param_set),
+    IP.obliq_epoch(param_set),
+    IP.eccentricity_epoch(param_set))
+
 # calculate the distance, declination, and hour angle (at lon=0)
 function distance_declination_hourangle(::Type{FT},
                                         date::DateTime,
@@ -35,17 +46,8 @@ function distance_declination_hourangle(::Type{FT},
 
     # mean anomaly given mean anomaly at epoch
     MA = mod(FT(2π) * (Δt_years) + M0, FT(2π))
-
     # calculate orbital parameters or take values at J2000
-    if milankovitch
-        ϖ = FT(ϖ_spline(od, Δt_years));
-        γ = FT(γ_spline(od, Δt_years));
-        e = FT(e_spline(od, Δt_years));
-    else
-        ϖ = FT(IP.lon_perihelion_epoch(param_set));
-        γ = FT(IP.obliq_epoch(param_set));
-        e = FT(IP.eccentricity_epoch(param_set));
-    end
+    (ϖ, γ, e) = milankovitch ? compute_orbital_parameters(od, Δt_years) : compute_orbital_parameters(param_set)
 
     # true anomaly, radians
     TA = true_anomaly(MA, e)
@@ -163,4 +165,4 @@ function daily_zenith_angle(date::DateTime,
     daily_θ = mod(acos(FT(1/π)*(ηd*sin(ϕ)*sin(δ) + cos(ϕ)*cos(δ)*sin(ηd))), FT(2π))
 
     return daily_θ, d
-end
+end