Add tests for guillot, ocean, and extend consts/phys coverage

- test_guillot.jl: new testset covering eval_tau (exact values +
  linearity), stellar temperature geometric relations, T4 positivity,
  T4 monotonicity, and profile helpers (length, positivity, sorting)
- test_ocean.jl: new testset covering empty, basin-fill, planet-flood,
  two-liquid ordering, and depth-additivity scenarios
- test_consts.jl: extend with fundamental constants (σSB, k_B, h_pl,
  c_vac, G_grav, zero_celcius), more mmw/atom-count entries, colour
  lookup completeness, species list consistency, and solar metallicity
- test_phys.jl: add Psat at water boiling point, Lv at boiling point,
  Kc positivity/monotonicity, is_vapour for ideal gas, and
  _pretty_color for known and unknown gases

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>

Author: nichollsh

test/runtests.jl

@@ -57,6 +57,8 @@ end
 LoggingExtras.global_logger(Logging.SimpleLogger(Logging.Warn))
 include("test_consts.jl")
 include("test_phys.jl")
+include("test_guillot.jl")
+include("test_ocean.jl")
 include("test_deep_heating.jl")
 if suite != "fast"
     include("test_integration.jl")

test/test_consts.jl

@@ -2,8 +2,70 @@ using Test
 using AGNI
 
 @testset "consts" begin
-    @test isapprox(AGNI.consts.R_gas, 8.314462618; atol=1e-12)
-    @test isapprox(AGNI.consts._lookup_mmw["H2O"], 1.801530E-02; rtol=1e-8)
+
+    # -------------
+    # Fundamental physical constants (NIST CODATA values)
+    # Update these if the source values are changed.
+    # -------------
+    @test isapprox(AGNI.consts.R_gas,        8.314462618;      atol=1e-12)
+    @test isapprox(AGNI.consts.σSB,          5.670374419e-8;   rtol=1e-9)
+    @test isapprox(AGNI.consts.k_B,          1.380649e-23;     rtol=1e-9)
+    @test isapprox(AGNI.consts.h_pl,         6.62607015e-34;   rtol=1e-9)
+    @test isapprox(AGNI.consts.c_vac,        299792458.0;      atol=0.1)
+    @test isapprox(AGNI.consts.G_grav,       6.67430e-11;      rtol=1e-5)
+    @test isapprox(AGNI.consts.zero_celcius, 273.15;           atol=1e-10)
+
+    # -------------
+    # Mean molecular weights [kg mol-1] for key gases
+    # -------------
+    mmw_cases = [
+        ("H2O", 1.801530e-02),
+        ("CO2", 4.401000e-02),
+        ("N2",  2.801340e-02),
+        ("H2",  2.015880e-03),
+    ]
+    for (gas, expected) in mmw_cases
+        @test isapprox(AGNI.consts._lookup_mmw[gas], expected; rtol=1e-6)
+    end
+
+    # -------------
+    # Atom counts for key molecules
+    # -------------
     @test AGNI.consts._lookup_count_atoms["H2O"]["H"] == 2
-    @test "H2O" in keys(AGNI.consts._lookup_color)
+    @test AGNI.consts._lookup_count_atoms["H2O"]["O"] == 1
+    @test AGNI.consts._lookup_count_atoms["CO2"]["C"] == 1
+    @test AGNI.consts._lookup_count_atoms["CO2"]["O"] == 2
+
+    # -------------
+    # Plotting colours: known gases and key elements must be present
+    # -------------
+    for gas in ["H2O", "CO2", "H2", "N2", "CH4"]
+        @test haskey(AGNI.consts._lookup_color, gas)
+    end
+    for elem in ["H", "C", "O", "N", "S", "Fe", "Si"]
+        @test haskey(AGNI.consts._lookup_color, elem)
+    end
+
+    # -------------
+    # Standard species lists are non-empty and self-consistent
+    # -------------
+    @test length(AGNI.consts.vols_standard) > 0
+    @test length(AGNI.consts.vaps_standard) > 0
+    @test length(AGNI.consts.gases_standard) == length(AGNI.consts.vols_standard) +
+                                                length(AGNI.consts.vaps_standard)
+    @test "H2O" in AGNI.consts.vols_standard
+    @test "H"   in AGNI.consts.elems_standard
+
+    # -------------
+    # Solar metallicities: hydrogen is anchored at 12.00 by definition
+    # -------------
+    @test isapprox(AGNI.consts._solar_metallicity["H"], 12.00; atol=1e-10)
+    @test haskey(AGNI.consts._solar_metallicity, "Fe")
+
+    # -------------
+    # Liquid densities for ocean module
+    # -------------
+    @test isapprox(AGNI.consts._lookup_liquid_rho["H2O"], 958.37; rtol=1e-5)
+    @test AGNI.consts._lookup_liquid_rho["CO2"] > AGNI.consts._lookup_liquid_rho["H2O"]
+
 end

test/test_guillot.jl

@@ -0,0 +1,103 @@
+using Test
+using AGNI
+
+# Access the guillot submodule through setpt
+const guillot = AGNI.setpt.guillot
+
+# Guillot (2010) analytical grey-gas T(p) profile module.
+# Constants inside the module:
+#   κ_th = 7e-2 * 10 = 0.7  m2/kg  (LW opacity)
+#   κ_vs = 4e-3 * 10 = 0.04 m2/kg  (SW opacity)
+#   γ = κ_vs/κ_th = 0.04/0.7       (opacity ratio)
+
+@testset "guillot" begin
+
+    # -------------
+    # Optical depth: τ = p * κ_th / g
+    # -------------
+    @testset "eval_tau" begin
+        p_test = [1e3,   1e5,      1e7     ]  # pressure [Pa]
+        g_test = [10.0,  10.0,     25.0    ]  # gravity [m s-2]
+        v_expt = [70.0,  7000.0,   280000.0]  # τ = p * 0.7 / g [dimensionless]
+        for i in eachindex(p_test)
+            @test isapprox(guillot.eval_tau(p_test[i], g_test[i]), v_expt[i]; rtol=1e-10)
+        end
+        # τ must be linear in pressure at fixed gravity
+        @test isapprox(guillot.eval_tau(2e5, 10.0), 2 * guillot.eval_tau(1e5, 10.0); rtol=1e-10)
+    end
+
+    # -------------
+    # Stellar temperatures: geometric relations
+    #   eval_Tirr(Tstar, Rstar, sep) = Tstar * sqrt(Rstar / sep)
+    #   eval_Teqm(Tstar, Rstar, sep) = Tstar * sqrt(Rstar / (2*sep))
+    # So Tirr / Teqm == sqrt(2) exactly.
+    # -------------
+    @testset "stellar_temps" begin
+        Tstar = 5778.0    # solar effective temperature [K]
+        Rstar = 6.957e8   # solar radius [m]
+        sep   = 1.496e11  # Earth-Sun distance [m]
+
+        Tirr = guillot.eval_Tirr(Tstar, Rstar, sep)
+        Teqm = guillot.eval_Teqm(Tstar, Rstar, sep)
+
+        @test Tirr > 0.0
+        @test Teqm > 0.0
+        @test Teqm < Tirr
+        # exact geometric relationship between the two temperatures
+        @test isapprox(Tirr / Teqm, sqrt(2.0); rtol=1e-10)
+    end
+
+    # -------------
+    # T^4 functions must return positive values for physical inputs
+    # Note: eval_T4_avg uses E2(γτ) which diverges at τ=0 (τ>0 only).
+    # eval_T4_cos does not use E1/E2 so τ=0 is fine.
+    # -------------
+    @testset "T4_positive" begin
+        τ_vals_avg = [0.01, 0.1, 1.0, 10.0, 100.0]  # τ > 0 for avg (E1 singularity at 0)
+        τ_vals_cos = [0.0,  0.1, 1.0, 10.0, 100.0]  # τ = 0 is fine for cos
+        Tint = 100.0    # internal temperature [K]
+        Teqm = 800.0    # equilibrium temperature [K]
+        Tirr = 1000.0   # irradiation temperature [K]
+        θ    = 45.0     # zenith angle [degrees]
+        for τ in τ_vals_avg
+            @test guillot.eval_T4_avg(τ, Tint, Teqm) > 0.0
+        end
+        for τ in τ_vals_cos
+            @test guillot.eval_T4_cos(τ, Tint, Tirr, θ) > 0.0
+        end
+    end
+
+    # -------------
+    # T^4 must increase monotonically with optical depth (deeper = hotter)
+    # for the average profile at sensible parameters
+    # -------------
+    @testset "T4_monotone" begin
+        τ_lo, τ_hi = 0.1, 10.0
+        Tint, Teqm, Tirr, θ = 200.0, 1200.0, 1500.0, 30.0
+        @test guillot.eval_T4_avg(τ_hi, Tint, Teqm) > guillot.eval_T4_avg(τ_lo, Tint, Teqm)
+        @test guillot.eval_T4_cos(τ_hi, Tint, Tirr, θ) > guillot.eval_T4_cos(τ_lo, Tint, Tirr, θ)
+    end
+
+    # -------------
+    # Profile helpers: correct length and all-positive temperatures
+    # -------------
+    @testset "profiles" begin
+        grav  = 10.0   # [m s-2]
+        Tint  = 100.0  # [K]
+        Teqm  = 900.0  # [K]
+        Tirr  = 1200.0 # [K]
+        θ     = 45.0   # [degrees]
+        p_arr = [1e1, 1e2, 1e3, 1e4, 1e5, 1e6]  # pressure grid [Pa]
+
+        T_avg = guillot.calc_profile_avg(p_arr, grav, Tint, Teqm)
+        @test length(T_avg) == length(p_arr)
+        @test all(T_avg .> 0.0)
+        @test issorted(T_avg)  # temperature increases towards higher pressures
+
+        T_cos = guillot.calc_profile_cos(p_arr, grav, Tint, Tirr, θ)
+        @test length(T_cos) == length(p_arr)
+        @test all(T_cos .> 0.0)
+        @test issorted(T_cos)
+    end
+
+end

test/test_ocean.jl

@@ -0,0 +1,77 @@
+using Test
+using AGNI
+
+# Ocean module: pure geometric functions for distributing condensed liquids
+# across ocean basins and continental shelves.
+# Liquid densities used below come from consts._lookup_liquid_rho:
+#   H2O = 958.37 kg/m^3
+#   CO2 = 1178.4 kg/m^3  (denser → sinks to bottom)
+
+@testset "ocean" begin
+
+    # Shared planetary parameters (easy to update for a different test planet)
+    fOB = 0.3      # ocean basin area fraction
+    hCS = 1000.0   # continental shelf height [m]
+    Rpl = 6.371e6  # planet radius [m]
+
+    # -------------
+    # No liquids → empty layers, no coverage
+    # -------------
+    @testset "empty" begin
+        sigs   = Dict{String, Float64}()
+        layers = AGNI.ocean.dist_surf_liq(sigs, fOB, hCS, Rpl)
+        @test length(layers) == 0
+        @test AGNI.ocean.get_topliq(layers)          == "none"
+        @test AGNI.ocean.get_maxdepth(layers)        == 0.0
+        @test AGNI.ocean.get_areacov(layers, fOB)    == 0.0
+    end
+
+    # -------------
+    # Small amount of H2O fits entirely within the ocean basins
+    # Expected: one layer, area coverage = fOB (continental planet)
+    # -------------
+    @testset "basin_fill" begin
+        sigs   = Dict("H2O" => 100.0)  # 100 kg/m^2 of rainout
+        layers = AGNI.ocean.dist_surf_liq(sigs, fOB, hCS, Rpl)
+        @test length(layers) == 1
+        @test AGNI.ocean.get_topliq(layers)          == "H2O"
+        @test AGNI.ocean.get_maxdepth(layers)        > 0.0
+        @test isapprox(AGNI.ocean.get_areacov(layers, fOB), fOB; atol=1e-10)
+    end
+
+    # -------------
+    # Very large amount of H2O overflows basins onto the continental shelf
+    # Expected: area coverage = 1.0 (aqua planet)
+    # -------------
+    @testset "planet_flood" begin
+        sigs   = Dict("H2O" => 1.0e6)  # huge amount; greatly exceeds basin capacity
+        layers = AGNI.ocean.dist_surf_liq(sigs, fOB, hCS, Rpl)
+        @test AGNI.ocean.get_topliq(layers)       == "H2O"
+        @test AGNI.ocean.get_areacov(layers, fOB) == 1.0
+    end
+
+    # -------------
+    # Two liquids: CO2 is denser (1178 kg/m^3) so it sinks below H2O (958 kg/m^3).
+    # The top of the ocean is the least-dense liquid exposed to the atmosphere.
+    # -------------
+    @testset "two_liquids" begin
+        sigs   = Dict("H2O" => 200.0, "CO2" => 100.0)
+        layers = AGNI.ocean.dist_surf_liq(sigs, fOB, hCS, Rpl)
+        @test length(layers) == 2
+        @test AGNI.ocean.get_topliq(layers) == "H2O"   # H2O floats on CO2
+        @test AGNI.ocean.get_maxdepth(layers) > 0.0
+    end
+
+    # -------------
+    # Depth is additive over layers (basin + shelf components)
+    # -------------
+    @testset "depth_additive" begin
+        # Place enough liquid that both in-basin and shelf portions are non-zero
+        sigs   = Dict("H2O" => 1.0e5)
+        layers = AGNI.ocean.dist_surf_liq(sigs, fOB, hCS, Rpl)
+        depth  = AGNI.ocean.get_maxdepth(layers)
+        manual = sum(la[3] + la[4] for la in layers)
+        @test isapprox(depth, manual; atol=1e-10)
+    end
+
+end

test/test_phys.jl

@@ -187,4 +187,77 @@ TEST_DIR        = joinpath(ROOT_DIR,"test/")
         atmosphere.deallocate!(atmos)
         @test test_pass
     end
+
+
+    # -------------
+    # Test saturation pressure lookup
+    # At the boiling point of water (373.15 K) Psat should be ~1 atm = 101325 Pa
+    # -------------
+    @testset "Psat" begin
+        gas_H2O::phys.Gas_t = phys.load_gas("$RES_DIR/thermodynamics/", "H2O", true, false)
+        T_boil  = 373.15  # [K]
+        Psat_e  = 101325.0  # expected ~1 atm [Pa]
+        Psat_o  = phys.get_Psat(gas_H2O, T_boil)
+        @test isfinite(Psat_o)
+        @test Psat_o > 0.0
+        @test isapprox(Psat_o, Psat_e; rtol=0.05)
+    end
+
+
+    # -------------
+    # Test latent heat lookup
+    # At the boiling point of water (373.15 K) Lv ≈ 2.256e6 J/kg
+    # -------------
+    @testset "Lv" begin
+        gas_H2O::phys.Gas_t = phys.load_gas("$RES_DIR/thermodynamics/", "H2O", true, false)
+        T_boil  = 373.15   # [K]
+        Lv_e    = 2.256e6  # expected latent heat [J/kg]
+        Lv_o    = phys.get_Lv(gas_H2O, T_boil)
+        @test isfinite(Lv_o)
+        @test Lv_o > 0.0
+        @test isapprox(Lv_o, Lv_e; rtol=0.05)
+    end
+
+
+    # -------------
+    # Test thermal conductivity (kinetic theory estimate)
+    # Must be positive and in a physically reasonable range for water vapour
+    # -------------
+    @testset "Kc" begin
+        gas_H2O::phys.Gas_t = phys.load_gas("$RES_DIR/thermodynamics/", "H2O", true, false)
+        t_test  = [200.0, 500.0, 1000.0]  # [K]
+        for t in t_test
+            Kc = phys.get_Kc(gas_H2O, t)
+            @test isfinite(Kc) && Kc > 0.0
+        end
+        # conductivity must increase with temperature (∝ sqrt(T))
+        @test phys.get_Kc(gas_H2O, 500.0) > phys.get_Kc(gas_H2O, 200.0)
+    end
+
+
+    # -------------
+    # Test is_vapour
+    # Ideal gas is always in the vapour phase (no real-gas condensation)
+    # -------------
+    @testset "is_vapour" begin
+        gas_N2::phys.Gas_t = phys.load_gas("$RES_DIR/thermodynamics/", "N2", true, false)
+        # ideal / no-sat stub is always vapour
+        @test phys.is_vapour(gas_N2, 300.0, 1e5)
+        @test phys.is_vapour(gas_N2, 100.0, 1e8)
+    end
+
+
+    # -------------
+    # Test _pretty_color
+    # Known gases return their lookup colour; unknown gases return a valid hex string
+    # -------------
+    @testset "pretty_color" begin
+        # known lookup entry
+        @test phys._pretty_color("CO2") == AGNI.consts._lookup_color["CO2"]
+        # unknown molecule: must return a 7-character hex code starting with '#'
+        col_sio = phys._pretty_color("SiO")
+        @test length(col_sio) == 7
+        @test col_sio[1] == '#'
+    end
+
 end