Merge branch 'main' into sk/update_inst_zen_ang_comp

.dev/Project.toml

@@ -0,0 +1,5 @@
+[deps]
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
+
+[compat]
+JuliaFormatter = "1"

.dev/clima_formatter_options.jl

@@ -0,0 +1,8 @@
+clima_formatter_options = (
+    indent = 4,
+    margin = 80,
+    always_for_in = true,
+    whitespace_typedefs = true,
+    whitespace_ops_in_indices = true,
+    remove_extra_newlines = false,
+)

.dev/climaformat.jl

@@ -0,0 +1,85 @@
+#!/usr/bin/env julia
+#
+# This is an adapted version of format.jl from JuliaFormatter with the
+# following license:
+#
+#    MIT License
+#    Copyright (c) 2019 Dominique Luna
+#
+#    Permission is hereby granted, free of charge, to any person obtaining a
+#    copy of this software and associated documentation files (the
+#    "Software"), to deal in the Software without restriction, including
+#    without limitation the rights to use, copy, modify, merge, publish,
+#    distribute, sublicense, and/or sell copies of the Software, and to permit
+#    persons to whom the Software is furnished to do so, subject to the
+#    following conditions:
+#
+#    The above copyright notice and this permission notice shall be included
+#    in all copies or substantial portions of the Software.
+#
+#    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+#    OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+#    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN
+#    NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+#    DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
+#    OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
+#    USE OR OTHER DEALINGS IN THE SOFTWARE.
+#
+using Pkg
+Pkg.activate(@__DIR__)
+Pkg.instantiate()
+
+using JuliaFormatter
+
+include("clima_formatter_options.jl")
+
+help = """
+Usage: climaformat.jl [flags] [FILE/PATH]...
+
+Formats the given julia files using the CLIMA formatting options.  If paths
+are given it will format the julia files in the paths. Otherwise, it will
+format all changed julia files.
+
+    -v, --verbose
+        Print the name of the files being formatted with relevant details.
+
+    -h, --help
+        Print this message.
+"""
+
+function parse_opts!(args::Vector{String})
+    i = 1
+    opts = Dict{Symbol, Union{Int, Bool}}()
+    while i ≤ length(args)
+        arg = args[i]
+        if arg[1] != '-'
+            i += 1
+            continue
+        end
+        if arg == "-v" || arg == "--verbose"
+            opt = :verbose
+        elseif arg == "-h" || arg == "--help"
+            opt = :help
+        else
+            error("invalid option $arg")
+        end
+        if opt in (:verbose, :help)
+            opts[opt] = true
+            deleteat!(args, i)
+        end
+    end
+    return opts
+end
+
+opts = parse_opts!(ARGS)
+if haskey(opts, :help)
+    write(stdout, help)
+    exit(0)
+end
+if isempty(ARGS)
+    filenames = readlines(`git ls-files "*.jl"`)
+else
+    filenames = ARGS
+end
+
+format(filenames; clima_formatter_options..., opts...)

docs/make.jl

@@ -4,11 +4,11 @@ using Insolation, Documenter
 ENV["GKSwstype"] = "100"
 
 pages = Any[
-        "Home" => "index.md"
-        "Zenith Angle Equations" => "ZenithAngleEquations.md"
-        "Insolation Examples" => "InsolationExamples.md"
-        "Milankovitch Cycles" => "Milankovitch.md"
-        "APIs" => "library.md"
+    "Home" => "index.md"
+    "Zenith Angle Equations" => "ZenithAngleEquations.md"
+    "Insolation Examples" => "InsolationExamples.md"
+    "Milankovitch Cycles" => "Milankovitch.md"
+    "APIs" => "library.md"
 ]
 
 format = Documenter.HTML(
@@ -20,7 +20,6 @@ makedocs(
     sitename = "Insolation.jl",
     format = format,
     clean = true,
-    strict = true,
     modules = [Insolation],
     pages = pages,
 )

docs/src/find_equinox_perihelion_dates.jl

@@ -10,30 +10,32 @@ param_set = create_insolation_parameters(FT)
 
 # 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., param_set, milankovitch=true)
-    theta_n, dist = daily_zenith_angle(date, od, 45., param_set, milankovitch=true)
+    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)
     return theta_n - theta_s
 end
 
 # x is date relative to March 1, with 1.00 representing March 1 00:00
 function xtomarchdate(x, year)
     basedate = Dates.DateTime(year, 3, 1)
-    deltat = Dates.Second(round((x-1)*IP.day(param_set)))
+    deltat = Dates.Second(round((x - 1) * IP.day(param_set)))
     return basedate + deltat
 end
 
 # Earth-Sun distance
 function edist(x, year, od)
-    date = xtojandate(x,year)
-    _, dist = daily_zenith_angle(date, od, 0., param_set, milankovitch=true)
-    return dist/IP.orbit_semimaj(param_set)
+    date = xtojandate(x, year)
+    _, dist = daily_zenith_angle(date, od, 0.0, param_set, milankovitch = true)
+    return dist / IP.orbit_semimaj(param_set)
 end
 
 # x is date relative to Jan 1, with 1.00 representing Jan 1 00:00
 function xtojandate(x, year)
     basedate = Dates.DateTime(year, 1, 1)
-    deltat = Dates.Second(round((x-1)*IP.day(param_set)))
+    deltat = Dates.Second(round((x - 1) * IP.day(param_set)))
     date = basedate + deltat
     return date
 end

docs/src/plot_diurnal_cycle.jl

@@ -15,19 +15,36 @@ function diurnal_cycle(lat, lon, date, od, timezone, filename)
     hours = collect(range(0, stop = 24, length = nhours))
     insol = zeros(nhours)
     sza = zeros(nhours)
-    for (i,hr) in enumerate(hours)
+    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)
-        insol[i] = S*mu
+        insol[i] = S * mu
         sza[i] = rad2deg(acos(mu))
     end
-    plot(hours, insol, color=:blue, lw=2, label="", 
-            ylabel="Insolation [W/m2]", yguidefontcolor=:blue, dpi=200, 
-            left_margein=10Plots.mm, right_margin=15Plots.mm)
-    plot!(twinx(), hours, sza, color=:red, lw=2, label="", 
-            ylabel="SZA [deg]", yguidefontcolor=:red)
+    plot(
+        hours,
+        insol,
+        color = :blue,
+        lw = 2,
+        label = "",
+        ylabel = "Insolation [W/m2]",
+        yguidefontcolor = :blue,
+        dpi = 200,
+        left_margein = 10Plots.mm,
+        right_margin = 15Plots.mm,
+    )
+    plot!(
+        twinx(),
+        hours,
+        sza,
+        color = :red,
+        lw = 2,
+        label = "",
+        ylabel = "SZA [deg]",
+        yguidefontcolor = :red,
+    )
     xlabel!("Local Time")
     savefig(filename)
 end

docs/src/plot_insolation.jl

@@ -13,21 +13,30 @@ import Insolation.Parameters as IP
                             n_lats::I,
                             param_set::IP.AIP) where {I<:Int}
 """
-function calc_day_lat_insolation(od, n_days::I,
-                                n_lats::I,
-                                param_set::IP.AIP) where {I<:Int}
-  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)
-  # 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, od, lat, param_set, milankovitch=false)
-        F_arr[i, j] = insolation(θ, dist, param_set)
+function calc_day_lat_insolation(
+    od,
+    n_days::I,
+    n_lats::I,
+    param_set::IP.AIP,
+) where {I <: Int}
+    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)
+    # 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,
+                od,
+                lat,
+                param_set,
+                milankovitch = false,
+            )
+            F_arr[i, j] = insolation(θ, dist, param_set)
+        end
     end
-  end
-  return d_arr, l_arr, F_arr
+    return d_arr, l_arr, F_arr
 end
 
 """
@@ -38,28 +47,55 @@ end
                             stitle,
                             file_name) where {FT<:AbstractFloat,I<:Int}
 """
-function plot_day_lat_insolation(d_arr::Array{I},
-                                 l_arr::Array{FT},
-                                 F_arr::Array{FT},
-                                 scmap,
-                                 stitle,
-                                 file_name) where {FT<:AbstractFloat,I<:Int}
-  if scmap == "YlOrRd"
-    cmap = :YlOrRd
-    vmin, vmax = 0, ceil(max(F_arr...)/100)*100
-  elseif scmap == "PRGn"
-    cmap = :PRGn
-    vmin, vmax = ceil(max(abs.(F_arr)...)/10)*-10, ceil(max(abs.(F_arr)...)/10)*10
-  end
+function plot_day_lat_insolation(
+    d_arr::Array{I},
+    l_arr::Array{FT},
+    F_arr::Array{FT},
+    scmap,
+    stitle,
+    file_name,
+) where {FT <: AbstractFloat, I <: Int}
+    if scmap == "YlOrRd"
+        cmap = :YlOrRd
+        vmin, vmax = 0, ceil(max(F_arr...) / 100) * 100
+    elseif scmap == "PRGn"
+        cmap = :PRGn
+        vmin, vmax = ceil(max(abs.(F_arr)...) / 10) * -10,
+        ceil(max(abs.(F_arr)...) / 10) * 10
+    end
 
-  p1 = contourf(d_arr, l_arr, F_arr', c=cmap, clims=(vmin,vmax), 
-    title=stitle, xlabel="Days since Jan 1 2000", ylabel="Latitude", colorbar_title="ToA Insolation [W/m2]")
-  meanF = mean(F_arr, dims=1)[1,:]
-  x1,x2 = [floor(minimum(meanF)), ceil(maximum(meanF))]
-  p2 = plot(meanF, l_arr, xlabel="Annual-mean TOA \nInsolation [W/m2]", 
-    ylims=[-90,90], yticks=[-90,-60,-30,0,30,60,90], xlims=[x1,x2], legend = false)
-  plot(p1, p2, size=(800,400), layout = grid(1,2, widths=(0.8,0.2)), dpi=250,
-        left_margin=4Plots.mm, bottom_margin=6Plots.mm, right_margin=6Plots.mm)
+    p1 = contourf(
+        d_arr,
+        l_arr,
+        F_arr',
+        c = cmap,
+        clims = (vmin, vmax),
+        title = stitle,
+        xlabel = "Days since Jan 1 2000",
+        ylabel = "Latitude",
+        colorbar_title = "ToA Insolation [W/m2]",
+    )
+    meanF = mean(F_arr, dims = 1)[1, :]
+    x1, x2 = [floor(minimum(meanF)), ceil(maximum(meanF))]
+    p2 = plot(
+        meanF,
+        l_arr,
+        xlabel = "Annual-mean TOA \nInsolation [W/m2]",
+        ylims = [-90, 90],
+        yticks = [-90, -60, -30, 0, 30, 60, 90],
+        xlims = [x1, x2],
+        legend = false,
+    )
+    plot(
+        p1,
+        p2,
+        size = (800, 400),
+        layout = grid(1, 2, widths = (0.8, 0.2)),
+        dpi = 250,
+        left_margin = 4Plots.mm,
+        bottom_margin = 6Plots.mm,
+        right_margin = 6Plots.mm,
+    )
 
-  savefig(file_name)
+    savefig(file_name)
 end

parameters/create_parameters.jl

@@ -5,7 +5,7 @@ function create_insolation_parameters(FT, overrides::NamedTuple = NamedTuple())
     toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
     aliases = string.(fieldnames(IP.InsolationParameters))
     pairs = CP.get_parameter_values!(toml_dict, aliases, "Insolation")
-    params = merge((;pairs...), overrides) # overrides
+    params = merge((; pairs...), overrides) # overrides
     param_set = IP.InsolationParameters{FT}(; params...)
     return param_set
 end

src/Insolation.jl

@@ -25,17 +25,21 @@ struct OrbitalData{E, G, O}
     γ_spline_etp::G
     ϖ_spline_etp::O
     function OrbitalData()
-        datapath = joinpath(artifact"orb_params_dataset", "INSOL.LA2004.BTL.csv");
-        x, _ = readdlm(datapath, ',', Float64, header=true);
-        t_range = x[1,1]*1e3 : 1e3 : x[end,1]*1e3; # array of every 1 kyr to range of years
-        e_spline_etp = CubicSplineInterpolation(t_range, x[:,2], extrapolation_bc = NaN);
-        γ_spline_etp = CubicSplineInterpolation(t_range, x[:,3], extrapolation_bc = NaN);
-        ϖ_spline_etp = CubicSplineInterpolation(t_range, x[:,4], extrapolation_bc = NaN);
+        datapath =
+            joinpath(artifact"orb_params_dataset", "INSOL.LA2004.BTL.csv")
+        x, _ = readdlm(datapath, ',', Float64, header = true)
+        t_range = (x[1, 1] * 1e3):1e3:(x[end, 1] * 1e3) # array of every 1 kyr to range of years
+        e_spline_etp =
+            CubicSplineInterpolation(t_range, x[:, 2], extrapolation_bc = NaN)
+        γ_spline_etp =
+            CubicSplineInterpolation(t_range, x[:, 3], extrapolation_bc = NaN)
+        ϖ_spline_etp =
+            CubicSplineInterpolation(t_range, x[:, 4], extrapolation_bc = NaN)
 
         E = typeof(e_spline_etp)
         G = typeof(γ_spline_etp)
         O = typeof(ϖ_spline_etp)
-        return new{E,G,O}(e_spline_etp,γ_spline_etp,ϖ_spline_etp)
+        return new{E, G, O}(e_spline_etp, γ_spline_etp, ϖ_spline_etp)
     end
 end
 
@@ -60,4 +64,4 @@ end
 include("ZenithAngleCalc.jl")
 include("InsolationCalc.jl")
 
-end # module
+end # module