Merge pull request #60 from CliMA/sk/pass_date0

Refactor interface for `instantaneous_zenith_angle` function

Project.toml

@@ -1,7 +1,7 @@
 name = "Insolation"
 uuid = "e98cc03f-d57e-4e3c-b70c-8d51efe9e0d8"
 authors = ["Climate Modeling Alliance"]
-version = "0.7.0"
+version = "0.8.0"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

docs/src/find_equinox_perihelion_dates.jl

@@ -8,13 +8,27 @@ FT = Float64
 include(joinpath(pkgdir(Insolation), "parameters", "create_parameters.jl"))
 param_set = create_insolation_parameters(FT)
 
+date0 = DateTime("2000-01-01T11:58:56.816")
+
 # Difference in NH and SH zenith angles at time x in given year
 function zdiff(x, year, od)
     date = xtomarchdate(x, year)
-    theta_s, dist =
-        daily_zenith_angle(date, od, -45.0, param_set, milankovitch = true)
-    theta_n, dist =
-        daily_zenith_angle(date, od, 45.0, param_set, milankovitch = true)
+    theta_s, dist = daily_zenith_angle(
+        date,
+        date0,
+        od,
+        -45.0,
+        param_set,
+        milankovitch = true,
+    )
+    theta_n, dist = daily_zenith_angle(
+        date,
+        date0,
+        od,
+        45.0,
+        param_set,
+        milankovitch = true,
+    )
     return theta_n - theta_s
 end
 
@@ -28,7 +42,8 @@ end
 # Earth-Sun distance
 function edist(x, year, od)
     date = xtojandate(x, year)
-    _, dist = daily_zenith_angle(date, od, 0.0, param_set, milankovitch = true)
+    _, dist =
+        daily_zenith_angle(date, date0, od, 0.0, param_set, milankovitch = true)
     return dist / IP.orbit_semimaj(param_set)
 end
 

docs/src/library.md

@@ -16,7 +16,7 @@ Insolation.OrbitalData
 ## Zenith Angle
 ```@autodocs
 Modules = [Insolation]
-Private = false
+Private = true
 Pages   = ["ZenithAngleCalc.jl"]
 ```
 
@@ -25,4 +25,4 @@ Pages   = ["ZenithAngleCalc.jl"]
 Modules = [Insolation]
 Private = false
 Pages   = ["InsolationCalc.jl"]
-```
+```

docs/src/plot_diurnal_cycle.jl

@@ -15,11 +15,13 @@ function diurnal_cycle(lat, lon, date, od, timezone, filename)
     hours = collect(range(0, stop = 24, length = nhours))
     insol = zeros(nhours)
     sza = zeros(nhours)
+    date0 = DateTime("2000-01-01T11:58:56.816")
+
     for (i, hr) in enumerate(hours)
         h = Int(round(hr + timezone))
         m = Int(round((hr + timezone - h) * 60))
         datetime = date + Dates.Hour(h) + Dates.Minute(m)
-        S, mu = solar_flux_and_cos_sza(datetime, od, lon, lat, param_set)
+        S, mu = solar_flux_and_cos_sza(datetime, date0, od, lon, lat, param_set)
         insol[i] = S * mu
         sza[i] = rad2deg(acos(mu))
     end

docs/src/plot_insolation.jl

@@ -22,12 +22,14 @@ function calc_day_lat_insolation(
     d_arr = Array{I}(round.(collect(range(0, stop = 365, length = n_days))))
     l_arr = collect(range(-90, stop = 90, length = n_lats))
     F_arr = zeros(n_days, n_lats)
+    date0 = DateTime("2000-01-01T11:58:56.816")
     # loop over days
     for (i, d) in enumerate(d_arr)
         for (j, lat) in enumerate(l_arr)
             date = Dates.DateTime(2000, 1, 1) + Dates.Day(d)
             θ, dist = daily_zenith_angle(
                 date,
+                date0,
                 od,
                 lat,
                 param_set,

src/InsolationCalc.jl

@@ -22,6 +22,7 @@ end
 
 """
     solar_flux_and_cos_sza(date::DateTime,
+                      date0::DateTime,
                       od::OrbitalData,
                       longitude::FT,
                       latitude::FT,
@@ -34,16 +35,26 @@ param_set is an AbstractParameterSet from CLIMAParameters.jl.
 """
 function solar_flux_and_cos_sza(
     date::DateTime,
+    date0::DateTime,
     od::OrbitalData,
     longitude::FT,
     latitude::FT,
     param_set::IP.AIP,
 ) where {FT <: Real}
     S0::FT = IP.tot_solar_irrad(param_set)
     d0::FT = IP.orbit_semimaj(param_set)
+    args = (
+        Insolation.helper_instantaneous_zenith_angle(
+            date,
+            date0,
+            od,
+            param_set,
+        )...,
+        longitude,
+        latitude,
+    )
     # θ = solar zenith angle, ζ = solar azimuth angle, d = earth-sun distance
-    θ, ζ, d =
-        instantaneous_zenith_angle(date, od, longitude, latitude, param_set)
+    θ, ζ, d = instantaneous_zenith_angle(args...)
     # set max. zenith angle to π/2, insolation should not be negative
     if θ > FT(π) / 2
         θ = FT(π) / 2

src/ZenithAngleCalc.jl

@@ -29,20 +29,33 @@ compute_orbital_parameters(param_set) = (
 
 function get_Δt_years(
     param_set::IP.InsolationParameters{FT},
-    date,
-    epoch_string = "2000-01-01T11:58:56.816",
+    date::DateTime,
+    date0::DateTime,
 ) 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)
+"""
+    distance_declination_hourangle(
+        date::DateTime,
+        date0::DateTime,
+        (ϖ, γ, e)::Tuple{FT, FT, FT},
+        param_set::IP.AIP,
+        eot_correction::Bool,
+    ) where {FT}
+
+Returns the earth-sun distance (m), declination angle (radians) and hour angle 
+(radians), at 0ᵒ longitude, given the current datetime, epoch datetime, 
+`longitude of the perihelion at epoch`, `obliquity at epoch`, `eccentricity at epoch`, 
+and `eot_correction`.
+
+"""
 function distance_declination_hourangle(
     date::DateTime,
+    date0::DateTime,
     (ϖ, γ, e)::Tuple{FT, FT, FT},
-    epoch_string,
     param_set::IP.AIP,
     eot_correction::Bool,
 ) where {FT}
@@ -51,8 +64,6 @@ function distance_declination_hourangle(
     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 =
         FT(datetime2julian(date) - datetime2julian(date0)) / days_per_year
@@ -84,27 +95,29 @@ function distance_declination_hourangle(
     return d, δ, η_UTC
 end
 
+"""
+    instantaneous_zenith_angle(
+        d::FT,
+        δ::FT,
+        η_UTC::FT,
+        longitude::FT,
+        latitude::FT,
+    ) where {FT}
+
+Returns the zenith angle (radians), azimuthal angle (radians) and Earth-Sun distance (m)
+at a particular longitude (degrees) and latitude (degrees) given Earth-Sun distance (m), 
+declination angle (radians), and hour angle (radians) at 0ᵒ longitude.
+
+"""
 function instantaneous_zenith_angle(
-    date::DateTime,
-    epoch_string,
-    (ϖ, γ, e)::Tuple{FT, FT, FT},
+    d::FT,
+    δ::FT,
+    η_UTC::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,
-    )
-
     # hour angle
     η = mod(η_UTC + λ, FT(2π))
 
@@ -121,92 +134,40 @@ function instantaneous_zenith_angle(
     return θ, ζ, d
 end
 
-"""
-    instantaneous_zenith_angle(date::DateTime,
-                               od::OrbitalData,
-                               longitude::FT,
-                               latitude::FT,
-                               param_set::IP.AIP;
-                               eot_correction::Bool=true,
-                               ) where {FT}
-
-Returns the zenith angle and earth-sun distance
-at a particular longitude and latitude on the given date (and time UTC)
-given orbital parameters: obliquity, longitude of perihelion, and eccentricity
-param_set is an AbstractParameterSet from CLIMAParameters.jl.
-
-`eot_correction` is an optional Boolean keyword argument that defaults to true
-when set to true the equation of time correction is turned on.
-This switch functionality is implemented for easy comparisons with reanalyses.
-
-Orbital parameters are computed using the `Milankovitch` method.
-"""
-function instantaneous_zenith_angle(
+function helper_instantaneous_zenith_angle(
     date::DateTime,
+    date0::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(
+    param_set::AIP;
+    eot_correction = true,
+)
+    Δt_years = get_Δt_years(param_set, date, date0)
+    return distance_declination_hourangle(
         date,
-        epoch_string,
-        compute_orbital_parameters(od, Δt_years),
-        longitude,
-        latitude,
+        date0,
+        Insolation.compute_orbital_parameters(od, Δt_years),
         param_set,
         eot_correction,
     )
 end
 
-"""
-    instantaneous_zenith_angle(date::DateTime,
-                               longitude::FT,
-                               latitude::FT,
-                               param_set::IP.AIP;
-                               eot_correction::Bool=true,
-                               ) where {FT}
-
-Returns the zenith angle and earth-sun distance
-at a particular longitude and latitude on the given date (and time UTC)
-given orbital parameters: obliquity, longitude of perihelion, and eccentricity
-param_set is an AbstractParameterSet from CLIMAParameters.jl.
-
-`eot_correction` is an optional Boolean keyword argument that defaults to true
-when set to true the equation of time correction is turned on.
-This switch functionality is implemented for easy comparisons with reanalyses.
-
-Orbital parameters at the J2000 epoch from CLIMAParameters are used.
-"""
-function instantaneous_zenith_angle(
+helper_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,
-    )
+    date0::DateTime,
+    param_set::AIP;
+    eot_correction = true,
+) = distance_declination_hourangle(
+    date,
+    date0,
+    Insolation.compute_orbital_parameters(param_set),
+    param_set,
+    eot_correction,
+)
 
-end
 
 """
     daily_zenith_angle(date::DateTime,
+                       data0::DateTime,
                        od::OrbitalData,
                        latitude::FT,
                        param_set::IP.AIP;
@@ -228,6 +189,7 @@ when set to false the orbital parameters at the J2000 epoch from CLIMAParameters
 """
 function daily_zenith_angle(
     date::DateTime,
+    date0::DateTime,
     od::OrbitalData,
     latitude::FT,
     param_set::IP.AIP;
@@ -236,9 +198,7 @@ function daily_zenith_angle(
 ) 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)
+    Δt_years = get_Δt_years(param_set, date, date0)
 
     # calculate orbital parameters or take values at J2000
     (ϖ, γ, e) =
@@ -247,8 +207,8 @@ function daily_zenith_angle(
 
     d, δ, _ = distance_declination_hourangle(
         date,
+        date0,
         (ϖ, γ, e),
-        epoch_string,
         param_set,
         eot_correction,
     )

test/runtests.jl

@@ -10,11 +10,15 @@ using Optim
 
 import CLIMAParameters as CP
 import Insolation.Parameters as IP
+const AIP = IP.AbstractInsolationParams
+import Insolation.OrbitalData
 FT = Float32
 
 include(joinpath(pkgdir(Insolation), "parameters", "create_parameters.jl"))
 param_set = create_insolation_parameters(FT)
 
+date0 = DateTime("2000-01-01T11:58:56.816")
+
 @testset "Orbital Params" begin
     include("test_orbit_param.jl")
 end

test/test_equinox.jl

@@ -1,8 +1,8 @@
 # Difference in NH and SH zenith angles at time x in given year
 function zdiff(x, year, od)
     date = xtomarchdate(x, year)
-    theta_s, dist = daily_zenith_angle(date, od, FT(-45), param_set)
-    theta_n, dist = daily_zenith_angle(date, od, FT(45), param_set)
+    theta_s, dist = daily_zenith_angle(date, date0, od, FT(-45), param_set)
+    theta_n, dist = daily_zenith_angle(date, date0, od, FT(45), param_set)
     return theta_n - theta_s
 end