diff --git a/src/ZenithAngleCalc.jl b/src/ZenithAngleCalc.jl index f0e60cf..326cf64 100644 --- a/src/ZenithAngleCalc.jl +++ b/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 \ No newline at end of file +end