Use Datetime

src/SolarRadiation.jl

@@ -9,7 +9,7 @@ export NoScattering, DaveFurukawaScattering, ChandrasekharScattering
 
 export scattered_radiation
 export elevation_correction
-export solar_geometry, hour_angle
+export solar_geometry, hour_angle, orbital_angular_frequency
 export solar_radiation, solar_radiation!
 export allocate_output_arrays, allocate_buffers
 

src/solar_geometry.jl

@@ -22,34 +22,46 @@ end
 
 abstract type AbstractSolarGeometryModel end
 
+"""
+    orbital_angular_frequency(days_in_year::Real=365)
+
+Compute Earth's orbital angular frequency for a given year length.
+
+Handles leap years (366 days) and non-standard calendars (e.g., 360-day).
+"""
+orbital_angular_frequency(days_in_year::Real=365) = 2π / days_in_year
+
 """
     McCulloughPorterSolarGeometry
 
 Solar geometry model based on McCullough & Porter (1971).
 
 # Fields
 - `reference_day`: Vernal equinox day of year (default: 80)
-- `orbital_angular_frequency`: Earth's orbital angular frequency (default: 2π/365)
 - `orbital_eccentricity`: Earth's orbital eccentricity (default: 0.0167238)
 - `declination_amplitude`: Solar declination amplitude (default: 0.39784993)
+
+Note: `orbital_angular_frequency` is now computed dynamically based on the year length
+to handle leap years and non-standard calendars. Use `orbital_angular_frequency(days_in_year)`.
 """
-@kwdef struct McCulloughPorterSolarGeometry{RD,OAF,OE,DA} <: AbstractSolarGeometryModel
+@kwdef struct McCulloughPorterSolarGeometry{RD,OE,DA} <: AbstractSolarGeometryModel
     reference_day::RD = 80
-    orbital_angular_frequency::OAF = 2π / 365
     orbital_eccentricity::OE = 0.0167238
     declination_amplitude::DA = 0.39784993
 end
 
 """
-    solar_geometry(model::McCulloughPorterSolarGeometry, latitude; day_of_year, hour_angle)
+    solar_geometry(model::McCulloughPorterSolarGeometry, latitude; day_of_year, hour_angle, days_in_year=365)
 
 Compute solar geometry parameters based on McCullough & Porter (1971).
 
 # Arguments
 - `model`: Solar geometry model with orbital parameters
 - `latitude`: Observer latitude (with angle units, e.g. `u"°"` or `u"rad"`)
-- `day_of_year`: Day of year (1–365)
+- `day_of_year`: Day of year (1–365, 1–366 for leap years, or 1–360 for 360-day calendars)
 - `hour_angle`: Hour angle (radians)
+- `days_in_year`: Number of days in the year (default: 365). Use 366 for leap years,
+   or other values for non-standard calendars (e.g., 360 for 360-day calendar).
 
 # Returns
 NamedTuple with:
@@ -65,10 +77,13 @@ solar_geometry(::McCulloughPorterSolarGeometry, ::Missing; kwargs...) = missing
 function solar_geometry(sm::McCulloughPorterSolarGeometry, latitude::Quantity;
     day_of_year::Real,
     hour_angle::Quantity,
+    days_in_year::Real=365,
 )
-    (; reference_day, orbital_angular_frequency, orbital_eccentricity, declination_amplitude) = sm
+    (; reference_day, orbital_eccentricity, declination_amplitude) = sm
+    # Compute orbital angular frequency dynamically based on year length
+    ω = orbital_angular_frequency(days_in_year)
     # Use short aliases for equations (standard notation)
-    d0, ω, ϵ, se = reference_day, orbital_angular_frequency, orbital_eccentricity, declination_amplitude
+    d0, ϵ, se = reference_day, orbital_eccentricity, declination_amplitude
     d, h = day_of_year, hour_angle
 
     ζ = (ω * (d - d0)) + 2.0ϵ * (sin(ω * d) - sin(ω * d0))          # eq.5 McCullough & Porter (1971)

src/solar_radiation.jl

@@ -366,6 +366,7 @@ function solar_radiation!(out, buffers, solar_model::AbstractSolarRadiation;
     year::Real=1975,
     hours::AbstractVector{<:Real}=0:1:23,
     longitude_correction::Real=0.0,
+    days_in_year::Real=365,
 )
     # Unpack model parameters with short aliases for equations
     (; solar_geometry_model, precipitable_water, diffuse_model, mixing_ratio_height,
@@ -390,7 +391,7 @@ function solar_radiation!(out, buffers, solar_model::AbstractSolarRadiation;
         for j in 1:ntimes
             d, t = days[i], hours[j]
             h, tsn = hour_angle(t, longitude_correction)
-            solar_geom = solar_geometry(solar_geometry_model, ϕ; day_of_year=d, hour_angle=h)
+            solar_geom = solar_geometry(solar_geometry_model, ϕ; day_of_year=d, hour_angle=h, days_in_year)
             δ, z, ar² = solar_geom.solar_declination, solar_geom.zenith_angle, solar_geom.sun_distance_factor
             zsl = z
 
@@ -464,7 +465,8 @@ function solar_radiation!(out, buffers, solar_model::AbstractSolarRadiation;
 end
 
 """
-    solar_radiation(solar_model; solar_terrain, days, year, hours, longitude_correction)
+    solar_radiation(solar_model; solar_terrain, days, hours, ...)
+    solar_radiation(solar_model; solar_terrain, dates, hours, ...)
 
 Compute solar radiation for a given model and terrain configuration.
 
@@ -474,28 +476,85 @@ use `solar_radiation!` with pre-allocated buffers for better performance.
 # Arguments
 - `solar_model::AbstractSolarRadiation`: Solar radiation model parameters
 - `solar_terrain::AbstractTerrain`: Terrain configuration (elevation, slope, etc.)
+
+# Numeric interface (R-compatible)
 - `days`: Vector of days of year (default: mid-month days)
-- `year`: Year for leap year handling (default: 1975)
 - `hours`: Hours of day to compute (default: 0:23)
+- `year`: Year for leap year handling (default: 1975)
 - `longitude_correction`: Longitude correction in hours (default: 0.0)
+- `days_in_year`: Number of days in the year (default: 365). Use 366 for leap years,
+   or other values for non-standard calendars (e.g., 360 for 360-day calendar).
+
+# DateTime interface
+- `dates`: Vector of `AbstractDateTime` instances (one per day to simulate)
+- `hours`: Hour offsets as `Period` values (e.g., `Hour(0):Hour(1):Hour(23)`)
+- `timezone_offset`: Timezone offset from UTC in hours (default: 0.0)
+
+The DateTime interface extracts calendar information (leap years, 360-day calendars)
+from the date type automatically.
+
+# Examples
+```julia
+# Numeric interface
+solar_radiation(model; solar_terrain, days=1:10, hours=0:23, year=2024)
+
+# DateTime interface - mid-month days for a year
+using Dates
+dates = [Date(2024, m, 15) for m in 1:12]
+solar_radiation(model; solar_terrain, dates, hours=Hour(0):Hour(1):Hour(23))
+
+# With CFTime (no dependency required)
+using CFTime
+dates = [DateTimeNoLeap(2024, m, 15) for m in 1:12]
+solar_radiation(model; solar_terrain, dates, hours=Hour(0):Hour(1):Hour(23))
+```
 
 # Returns
 NamedTuple with zenith/azimuth angles, integrated irradiances, and spectral data.
 """
 function solar_radiation(solar_model::AbstractSolarRadiation;
     solar_terrain::AbstractTerrain,
-    days::Vector{<:Real}=[15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349],
+    # Numeric interface
+    days::Union{Nothing, AbstractVector{<:Real}}=nothing,
     year::Real=1975,
-    hours::AbstractVector{<:Real}=0:1:23,
+    hours::AbstractVector=0:1:23,
     longitude_correction::Real=0.0,
+    days_in_year::Real=365,
+    # DateTime interface (accepts Date, DateTime, or CFTime types)
+    dates::Union{Nothing, AbstractVector{<:Dates.TimeType}}=nothing,
+    timezone_offset::Real=0.0,
 )
+    # Determine which interface is being used
+    if !isnothing(dates)
+        # DateTime interface: extract from dates + hour offsets
+        days_in_year = Dates.daysinyear(first(dates))
+        year = Dates.year(first(dates))
+        longitude_correction = timezone_offset
+        # Convert dates to day-of-year, hours to numeric
+        days_numeric = [Dates.dayofyear(d) for d in dates]
+        hours_numeric = [Dates.value(Dates.Millisecond(h)) / 3_600_000 for h in hours]
+    elseif !isnothing(days)
+        # Numeric interface
+        days_numeric = collect(days)
+        hours_numeric = collect(hours)
+    else
+        # Default mid-month days
+        days_numeric = [15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349]
+        hours_numeric = collect(hours)
+    end
+
     nmax = solar_model.wavelength_count
-    ndays, ntimes = length(days), length(hours)
+    ndays, ntimes = length(days_numeric), length(hours_numeric)
     nsteps = ndays * ntimes
 
     out = allocate_output_arrays(nsteps, ndays, nmax)
     buffers = allocate_buffers(nmax, solar_model.diffuse_model)
 
     return solar_radiation!(out, buffers, solar_model;
-        solar_terrain, days, year, hours, longitude_correction)
+        solar_terrain,
+        days=days_numeric,
+        year,
+        hours=hours_numeric,
+        longitude_correction,
+        days_in_year)
 end

test/runtests.jl

@@ -4,4 +4,4 @@ using Aqua, SolarRadiation, DataFrames, CSV, Test, SafeTestsets, Unitful
     Aqua.test_all(SolarRadiation)
 end
 
-@safetestset "Test against NicheMapR outputs" begin include("solar_radiation.jl") end
+@safetestset "Test against NicheMapR outputs" begin include("solar_radiation.jl") end

test/solar_radiation.jl

@@ -74,3 +74,120 @@ day_of_year = repeat(days, inner=length(hours))
     @test diffuse_spectra ≈ diffuse_spectra_nmr_units rtol=1e-3
     @test rayleigh_spectra ≈ rayleigh_spectra_nmr_units rtol=1e-7
 end
+
+@testset "dates interface" begin
+    using Dates
+
+    # Create dates matching the numeric days (mid-month days for 1975)
+    year = 1975
+    # The numeric days are day-of-year values, convert to dates
+    dates_vec = [Date(year, 1, 1) + Day(Int(d) - 1) for d in days]
+
+    # Run with dates interface
+    solar_out_dt = solar_radiation(solar_model;
+        solar_terrain,
+        dates=dates_vec,
+        hours=Hour(0):Hour(1):Hour(23),
+    )
+
+    # Compare zenith angles (apply same capping as line 52)
+    zenith_dt = copy(solar_out_dt.zenith_angle)
+    zenith_dt[zenith_dt .> 90u"°"] .= 90u"°"
+    @test zenith_dt ≈ zenith_angle
+    @test solar_out_dt.global_horizontal ≈ global_horizontal
+    @test solar_out_dt.direct_spectra ≈ direct_spectra
+end
+
+@testset "dates with Date type" begin
+    using Dates
+
+    # Simple test with Date objects
+    dates_vec = [Date(2024, 1, 15), Date(2024, 2, 15)]
+    out = solar_radiation(solar_model;
+        solar_terrain,
+        dates=dates_vec,
+        hours=Hour(0):Hour(6):Hour(18),  # 4 hours per day
+    )
+
+    @test length(out.zenith_angle) == 2 * 4  # 2 days × 4 hours
+    @test out.day_of_year[1] == 15  # Jan 15
+    @test out.day_of_year[5] == 46  # Feb 15
+end
+
+@testset "dates with DateTime type" begin
+    using Dates
+
+    # Test with DateTime objects
+    dates_vec = [DateTime(2024, 6, 1), DateTime(2024, 6, 2), DateTime(2024, 6, 3)]
+    out = solar_radiation(solar_model;
+        solar_terrain,
+        dates=dates_vec,
+        hours=Hour(0):Hour(1):Hour(23),
+    )
+
+    @test length(out.zenith_angle) == 3 * 24  # 3 days × 24 hours
+end
+
+@testset "leap year via dates" begin
+    using Dates
+
+    # Leap year - 2024 has 366 days
+    dates_leap = [Date(2024, 12, 31)]
+    out_leap = solar_radiation(solar_model;
+        solar_terrain,
+        dates=dates_leap,
+        hours=[Hour(12)],  # just noon
+    )
+    @test out_leap.day_of_year[1] == 366
+
+    # Non-leap year - 2023 has 365 days
+    dates_normal = [Date(2023, 12, 31)]
+    out_normal = solar_radiation(solar_model;
+        solar_terrain,
+        dates=dates_normal,
+        hours=[Hour(12)],
+    )
+    @test out_normal.day_of_year[1] == 365
+end
+
+@testset "hour offsets" begin
+    using Dates
+
+    dates_vec = [Date(2024, 6, 15)]
+
+    # Test different hour offset patterns
+    out1 = solar_radiation(solar_model;
+        solar_terrain,
+        dates=dates_vec,
+        hours=Hour(0):Hour(1):Hour(23),  # hourly
+    )
+    @test length(out1.zenith_angle) == 24
+
+    out2 = solar_radiation(solar_model;
+        solar_terrain,
+        dates=dates_vec,
+        hours=Hour(0):Hour(3):Hour(21),  # every 3 hours
+    )
+    @test length(out2.zenith_angle) == 8
+
+    out3 = solar_radiation(solar_model;
+        solar_terrain,
+        dates=dates_vec,
+        hours=[Hour(6), Hour(12), Hour(18)],  # specific hours
+    )
+    @test length(out3.zenith_angle) == 3
+end
+
+@testset "orbital_angular_frequency" begin
+    # Standard year
+    @test orbital_angular_frequency(365) ≈ 2π / 365
+
+    # Leap year
+    @test orbital_angular_frequency(366) ≈ 2π / 366
+
+    # 360-day calendar
+    @test orbital_angular_frequency(360) ≈ 2π / 360
+
+    # Default
+    @test orbital_angular_frequency() ≈ 2π / 365
+end