Merge #51

51: Update instantaneous zenith angle computation r=sriharshakandala a=sriharshakandala Co-authored-by: sriharshakandala <sriharsha.kvs@gmail.com>
CliMA · Nov 1, 2023 · 3020aba · 3020aba
2 parents bda8307 + 4568c03
commit 3020aba
Show file tree

Hide file tree

Showing 2 changed files with 100 additions and 46 deletions.
diff --git 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
@@ -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