From 482cecf92fa9340d8e4b1dff03a8c4d4c4b3b23e Mon Sep 17 00:00:00 2001 From: sriharshakandala Date: Wed, 1 Nov 2023 10:05:25 -0700 Subject: [PATCH 1/2] Update instantaneous zenith angle computation --- src/ZenithAngleCalc.jl | 161 ++++++++++++++++++++++++++--------------- test/test_types.jl | 24 +++--- 2 files changed, 111 insertions(+), 74 deletions(-) diff --git a/src/ZenithAngleCalc.jl b/src/ZenithAngleCalc.jl index 32e2562..06196d8 100644 --- a/src/ZenithAngleCalc.jl +++ b/src/ZenithAngleCalc.jl @@ -14,23 +14,41 @@ function equation_of_time(e::FT, MA::FT, γ::FT, ϖ::FT) where {FT <: Real} return mod(_Δt + FT(π), FT(2π)) - FT(π) end -# calculate the distance, declination, and hour angle (at lon=0) -function distance_declination_hourangle( - ::Type{FT}, - date::DateTime, - od::OrbitalData, - 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) +# 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 - 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") +# calculate the distance, declination, and hour angle (at lon=0) +function distance_declination_hourangle(date::DateTime, + (ϖ, γ, e)::Tuple{FT,FT,FT}, + epoch_string, + param_set::IP.AIP, + eot_correction::Bool, +) 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) + + date0 = DateTime(epoch_string,dateformat"y-m-dTHH:MM:SS.s") days_per_year = Ya / day_length Δt_years = @@ -39,17 +57,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) @@ -96,26 +103,19 @@ This switch functionality is implemented for easy comparisons with reanalyses. when set to true the orbital parameters are calculated for the given DateTime when set to false the orbital parameters at the J2000 epoch from CLIMAParameters are used. """ -function instantaneous_zenith_angle( - date::DateTime, - od::OrbitalData, - longitude::FT, - latitude::FT, - param_set::IP.AIP; - eot_correction::Bool = true, - milankovitch::Bool = true, -) where {FT <: Real} +function instantaneous_zenith_angle(date::DateTime, + epoch_string, + (ϖ, γ, e)::Tuple{FT,FT,FT}, + longitude::FT, + latitude::FT, + param_set::IP.AIP, + eot_correction::Bool, +) where {FT} λ = deg2rad(longitude) ϕ = deg2rad(latitude) - d, δ, η_UTC = distance_declination_hourangle( - FT, - date, - od, - param_set, - eot_correction, - milankovitch, - ) + + d, δ, η_UTC = distance_declination_hourangle(date, (ϖ, γ, e), epoch_string, param_set, eot_correction) # hour angle η = mod(η_UTC + λ, FT(2π)) @@ -133,6 +133,49 @@ 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, @@ -154,24 +197,22 @@ This switch functionality is implemented for easy comparisons with reanalyses. when set to true the orbital parameters are calculated for the given DateTime, when set to false the orbital parameters at the J2000 epoch from CLIMAParameters are used. """ -function daily_zenith_angle( - date::DateTime, - od::OrbitalData, - latitude::FT, - param_set::IP.AIP; - eot_correction::Bool = true, - milankovitch::Bool = true, -) where {FT <: Real} +function daily_zenith_angle(date::DateTime, + od::OrbitalData, + latitude::FT, + param_set::IP.AIP; + eot_correction::Bool=true, + milankovitch::Bool=true) where {FT} ϕ = deg2rad(latitude) - d, δ, _ = distance_declination_hourangle( - FT, - date, - od, - param_set, - eot_correction, - milankovitch, - ) + 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(date, (ϖ, γ, e), epoch_string, param_set, eot_correction) # sunrise/sunset angle T = tan(ϕ) * tan(δ) diff --git a/test/test_types.jl b/test/test_types.jl index 10c9ca3..786f4a9 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -5,13 +5,12 @@ od = Insolation.OrbitalData() sza, azi, d = instantaneous_zenith_angle( date, - od, lon, lat, param_set, - eot_correction = false, - milankovitch = false, -) + eot_correction=false, + ) + F = insolation(sza, d, param_set) @test typeof(sza) == FT @@ -25,9 +24,8 @@ sza, azi, d = instantaneous_zenith_angle( lon, lat, param_set, - eot_correction = false, - milankovitch = true, -) + eot_correction=false, + ) F = insolation(sza, d, param_set) @test typeof(sza) == FT @@ -37,13 +35,11 @@ F = insolation(sza, d, param_set) sza, azi, d = instantaneous_zenith_angle( date, - od, lon, lat, param_set, - eot_correction = true, - milankovitch = false, -) + eot_correction=true, + ) F = insolation(sza, d, param_set) @test typeof(sza) == FT @@ -57,9 +53,9 @@ sza, azi, d = instantaneous_zenith_angle( lon, lat, param_set, - eot_correction = true, - milankovitch = true, -) + eot_correction=true, + ) + F = insolation(sza, d, param_set) @test typeof(sza) == FT From 4568c03cf39746fa74bf8f33d3add628cc32de5f Mon Sep 17 00:00:00 2001 From: sriharshakandala Date: Wed, 1 Nov 2023 10:06:12 -0700 Subject: [PATCH 2/2] Apply formatter --- src/ZenithAngleCalc.jl | 102 +++++++++++++++++++++++++---------------- test/test_types.jl | 21 ++++----- 2 files changed, 70 insertions(+), 53 deletions(-) diff --git a/src/ZenithAngleCalc.jl b/src/ZenithAngleCalc.jl index 06196d8..0aba3f0 100644 --- a/src/ZenithAngleCalc.jl +++ b/src/ZenithAngleCalc.jl @@ -18,37 +18,40 @@ end 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))) + 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)) + IP.eccentricity_epoch(param_set), +) function get_Δt_years( param_set::IP.InsolationParameters{FT}, date, - epoch_string="2000-01-01T11:58:56.816", + 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") + 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(date::DateTime, - (ϖ, γ, e)::Tuple{FT,FT,FT}, - epoch_string, - param_set::IP.AIP, - eot_correction::Bool, +function distance_declination_hourangle( + date::DateTime, + (ϖ, γ, e)::Tuple{FT, FT, FT}, + epoch_string, + param_set::IP.AIP, + eot_correction::Bool, ) 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) - 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 = @@ -103,19 +106,26 @@ This switch functionality is implemented for easy comparisons with reanalyses. when set to true the orbital parameters are calculated for the given DateTime when set to false the orbital parameters at the J2000 epoch from CLIMAParameters are used. """ -function instantaneous_zenith_angle(date::DateTime, - epoch_string, - (ϖ, γ, e)::Tuple{FT,FT,FT}, - longitude::FT, - latitude::FT, - param_set::IP.AIP, - eot_correction::Bool, +function instantaneous_zenith_angle( + date::DateTime, + epoch_string, + (ϖ, γ, e)::Tuple{FT, FT, 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) + d, δ, η_UTC = distance_declination_hourangle( + date, + (ϖ, γ, e), + epoch_string, + param_set, + eot_correction, + ) # hour angle η = mod(η_UTC + λ, FT(2π)) @@ -133,13 +143,14 @@ function instantaneous_zenith_angle(date::DateTime, 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} +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) @@ -155,12 +166,13 @@ function instantaneous_zenith_angle(date::DateTime, ) end -function instantaneous_zenith_angle(date::DateTime, - longitude::FT, - latitude::FT, - param_set::IP.AIP; - eot_correction::Bool=true, - ) where {FT} +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) @@ -197,22 +209,32 @@ This switch functionality is implemented for easy comparisons with reanalyses. when set to true the orbital parameters are calculated for the given DateTime, when set to false the orbital parameters at the J2000 epoch from CLIMAParameters are used. """ -function daily_zenith_angle(date::DateTime, - od::OrbitalData, - latitude::FT, - param_set::IP.AIP; - eot_correction::Bool=true, - milankovitch::Bool=true) where {FT} +function daily_zenith_angle( + date::DateTime, + od::OrbitalData, + latitude::FT, + param_set::IP.AIP; + eot_correction::Bool = true, + milankovitch::Bool = true, +) 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) + (ϖ, γ, e) = + milankovitch ? compute_orbital_parameters(od, Δt_years) : + compute_orbital_parameters(param_set) - d, δ, _ = distance_declination_hourangle(date, (ϖ, γ, e), epoch_string, param_set, eot_correction) + d, δ, _ = distance_declination_hourangle( + date, + (ϖ, γ, e), + epoch_string, + param_set, + eot_correction, + ) # sunrise/sunset angle T = tan(ϕ) * tan(δ) diff --git a/test/test_types.jl b/test/test_types.jl index 786f4a9..ef017e0 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -8,8 +8,8 @@ sza, azi, d = instantaneous_zenith_angle( lon, lat, param_set, - eot_correction=false, - ) + eot_correction = false, +) F = insolation(sza, d, param_set) @@ -24,8 +24,8 @@ sza, azi, d = instantaneous_zenith_angle( lon, lat, param_set, - eot_correction=false, - ) + eot_correction = false, +) F = insolation(sza, d, param_set) @test typeof(sza) == FT @@ -33,13 +33,8 @@ 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, lon, lat, param_set, eot_correction = true) F = insolation(sza, d, param_set) @test typeof(sza) == FT @@ -53,8 +48,8 @@ sza, azi, d = instantaneous_zenith_angle( lon, lat, param_set, - eot_correction=true, - ) + eot_correction = true, +) F = insolation(sza, d, param_set)