diff --git a/src/ZenithAngleCalc.jl b/src/ZenithAngleCalc.jl index 32e2562..0aba3f0 100644 --- a/src/ZenithAngleCalc.jl +++ b/src/ZenithAngleCalc.jl @@ -14,22 +14,43 @@ 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), +) + +function get_Δt_years( + param_set::IP.InsolationParameters{FT}, + date, + epoch_string = "2000-01-01T11:58:56.816", +) 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) function distance_declination_hourangle( - ::Type{FT}, date::DateTime, - od::OrbitalData, + (ϖ, γ, e)::Tuple{FT, FT, FT}, + epoch_string, param_set::IP.AIP, eot_correction::Bool, - milankovitch::Bool, -) where {FT <: Real} - Ya::FT = IP.year_anom(param_set) - day_length::FT = IP.day(param_set) - d0::FT = IP.orbit_semimaj(param_set) - M0::FT = IP.mean_anom_epoch(param_set) +) where {FT} + Ya = IP.year_anom(param_set) + day_length = IP.day(param_set) + d0 = IP.orbit_semimaj(param_set) + M0 = IP.mean_anom_epoch(param_set) - epoch_string = "2000-01-01T11:58:56.816" - # epoch_string::String = IP.epoch(param_set) date0 = DateTime(epoch_string, dateformat"y-m-dTHH:MM:SS.s") days_per_year = Ya / day_length @@ -39,17 +60,6 @@ function distance_declination_hourangle( # 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 - # true anomaly, radians TA = true_anomaly(MA, e) @@ -98,23 +108,23 @@ when set to false the orbital parameters at the J2000 epoch from CLIMAParameters """ function instantaneous_zenith_angle( date::DateTime, - od::OrbitalData, + epoch_string, + (ϖ, γ, e)::Tuple{FT, FT, FT}, longitude::FT, latitude::FT, - param_set::IP.AIP; - eot_correction::Bool = true, - milankovitch::Bool = true, -) where {FT <: Real} + param_set::IP.AIP, + eot_correction::Bool, +) where {FT} λ = deg2rad(longitude) ϕ = deg2rad(latitude) + d, δ, η_UTC = distance_declination_hourangle( - FT, date, - od, + (ϖ, γ, e), + epoch_string, param_set, eot_correction, - milankovitch, ) # hour angle @@ -133,6 +143,51 @@ function instantaneous_zenith_angle( return θ, ζ, d end +function instantaneous_zenith_angle( + date::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( + date, + epoch_string, + compute_orbital_parameters(od, Δt_years), + longitude, + latitude, + param_set, + eot_correction, + ) +end + +function 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, + ) + +end + """ daily_zenith_angle(date::DateTime, od::OrbitalData, @@ -161,16 +216,24 @@ function daily_zenith_angle( param_set::IP.AIP; eot_correction::Bool = true, milankovitch::Bool = true, -) where {FT <: Real} +) 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) + + # calculate orbital parameters or take values at J2000 + (ϖ, γ, e) = + milankovitch ? compute_orbital_parameters(od, Δt_years) : + compute_orbital_parameters(param_set) + d, δ, _ = distance_declination_hourangle( - FT, date, - od, + (ϖ, γ, e), + epoch_string, param_set, eot_correction, - milankovitch, ) # sunrise/sunset angle diff --git a/test/test_types.jl b/test/test_types.jl index 10c9ca3..ef017e0 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -5,14 +5,13 @@ od = Insolation.OrbitalData() sza, azi, d = instantaneous_zenith_angle( date, - od, lon, lat, param_set, eot_correction = false, - milankovitch = false, ) + F = insolation(sza, d, param_set) @test typeof(sza) == FT @test typeof(azi) == FT @@ -26,7 +25,6 @@ sza, azi, d = instantaneous_zenith_angle( lat, param_set, eot_correction = false, - milankovitch = true, ) F = insolation(sza, d, param_set) @@ -35,15 +33,8 @@ F = insolation(sza, d, param_set) @test typeof(d) == FT @test typeof(F) == FT -sza, azi, d = instantaneous_zenith_angle( - date, - od, - lon, - lat, - param_set, - eot_correction = true, - milankovitch = false, -) +sza, azi, d = + instantaneous_zenith_angle(date, lon, lat, param_set, eot_correction = true) F = insolation(sza, d, param_set) @test typeof(sza) == FT @@ -58,9 +49,9 @@ sza, azi, d = instantaneous_zenith_angle( lat, param_set, eot_correction = true, - milankovitch = true, ) + F = insolation(sza, d, param_set) @test typeof(sza) == FT @test typeof(azi) == FT