using Plots

# Saturated water vapor pressure [Pa]
# See Ham, 2005. Useful Equations and Tables in Micrometeorology
function water_vapor_saturated_pressure(Ta, P)
    Ta  =  Ta - 273.15
    es  =  (1.0007 + (3.46e-5 * (P/1e3))) * 0.61121 *
        exp((17.502 * Ta)/(Ta + 240.97)) * 1e3
    return es

end


function water_vapor_density(e, Tair)
    Rs_v = 461.5 # J/(kg.K) specific gas constant for water vapor
    ρv = e/(Rs_v * Tair)
    return ρv
end


# Calculate dry air density ρa [ kg m-3]
# P_a partial pressure of dry air [Pa]
# Tair air temperature [Kelvin]
function dry_air_density(P_a, Tair)
    Rs_a = 287.058 # J/(kg. K) specific gas constant for dry air
    ρa = P_a / (Rs_a  * Tair)
    return ρa
end


# Wet air density ρ [kg/m3]
# Tair air temperature in [Kelvin]
# P absolute atmospheric pressure [Pa]
function wet_air_density(RH, Tair, P)
    es = water_vapor_saturated_pressure(Tair, P)
    e = es*RH/100
    ρv = water_vapor_density(e, Tair)
    ρa = dry_air_density(P - e, Tair)
    ρ = ρv + ρa
    return ρ
end



# Calculate reference density
ρref = wet_air_density(0, 25 + 273.15, 1e5)

# Define ranges
Trange = (1:1:100) .+ 273.15
RHrange = (0:20:100)


# Create plot
plot(Trange .- 273.15, wet_air_density.(0, Trange, 1e5) ./ ρref, label = "RH = 0 %")
vline!([25], linestyle = :dash, label = missing)
for rh in RHrange[2:end]
    display(plot!(Trange .- 273.15, wet_air_density.(rh, Trange, 1e5) ./ ρref, label = "RH = $(rh) %"))
end
xlabel!("Temperature [deg C]")
ylabel!("Density ratio ρ/ρ_ref")
annotate!(23,0.7, Plots.text("Ref", rotation = 90))


savefig("air_density_changes_on_T_and_RH.svg")