#hide

# This is some boiler plate to make the notebook stand alone
# Creates a temporary project and adds only the packages used
# in the notebook to the Project.toml
using Pkg

Pkg.activate(; temp=true)
Pkg.add([PackageSpec(name="Plots", version="1.39.0"),
         PackageSpec(name="LaTeXStrings", version="1.3.1"),
         PackageSpec(name="Measures", version="0.3.2"),
         PackageSpec(name="Contour", version="0.6.2"),
         PackageSpec(name="GasDispersion", version="0.2.1")])
Pkg.instantiate();

#hide

using Plots, LaTeXStrings, Measures
using GasDispersion

colours = palette(:Paired_8);

#hide

sets = [ (name="Default", eqn=DefaultSet, color=:black, linestyle=:solid, lw=3),
         (name="CCPS Rural", eqn=CCPSRural, color=colours[2], linestyle=:solid, lw=3),
         (name="CCPS Urban", eqn=CCPSUrban, color=colours[4], linestyle=:solid, lw=3),
         (name="ISC3 Rural", eqn=ISC3Rural, color=colours[1], linestyle=:dash, lw=2),
         (name="ISC3 Urban", eqn=ISC3Urban, color=colours[3], linestyle=:dash, lw=2),
         (name="TNO", eqn=TNO, color=colours[6], linestyle=:dash, lw=2),
         (name="Turner", eqn=Turner, color=colours[7], linestyle=:dash, lw=2) ];

#hide

u0, z0 = 1.0, 1.0

xs = range(1.0,10.0,100)

function plot_windspeed!(plt,class)
    plot!(plt, xlabel=L"z/z_R", ylabel=L"u/u_R", title="$class", legend=false)

    for s in sets[1:5]
        ys = [ GasDispersion._windspeed(u0,z0,x,class,s.eqn) for x in xs]
        plot!(plt,xs,ys, color=s.color, linestyle=s.linestyle, lw=s.lw, label=s.name)
    end

end

ws_legend = plot(framestyle=:none, legend=:topleft, foreground_color_legend = nothing);

for s in sets[1:5]
    plot!(ws_legend,(1), color=s.color, linestyle=:solid, lw=2, label=s.name);
end

#hide

wsA, wsB, wsC, wsD = plot(), plot(), plot(), plot()

plot_windspeed!(wsA,ClassA)
plot_windspeed!(wsB,ClassB)
plot_windspeed!(wsC,ClassC)
plot_windspeed!(wsD,ClassD)

l1 = @layout [ grid(2,2) b{0.14w} ]

plot(wsA, wsB, wsC, wsD, ws_legend, layout=l1, size=(600,400), margin=2.5mm)

#hide

wsE, wsF = plot(), plot()

plot_windspeed!(wsE,ClassE)
plot_windspeed!(wsF,ClassF)

l2 = @layout [ grid(1,2) b{0.14w} ]

plot(wsE, wsF, ws_legend, layout=l2, size=(600,200), margin=2.5mm)

#hide

xs = range(1,1e5,1000)

function plot_cw!(plt,class)
    plot!(plt, xaxis=:log10, yaxis=:log10, xlabel=L"x", ylabel=L"\sigma_y", title="$class", legend=false)

    for s in sets
        ys = [ GasDispersion.crosswind_dispersion(x, Plume, class, s.eqn) for x in xs]
        plot!(plt, xs, ys, color=s.color, linestyle=s.linestyle, lw=s.lw, label=s.name)
    end
end

cw_legend = plot(framestyle=:none, legend=:topleft, foreground_color_legend = nothing);
for s in sets
    plot!(cw_legend,(1), color=s.color, linestyle=:solid, lw=2, label=s.name);
end

#hide

cwA, cwB, cwC, cwD, cwE, cwF = plot(), plot(), plot(), plot(), plot(), plot()

plot_cw!(cwA,ClassA)
plot_cw!(cwB,ClassB)
plot_cw!(cwC,ClassC)
plot_cw!(cwD,ClassD)
plot_cw!(cwE,ClassE)
plot_cw!(cwF,ClassF)

l1 = @layout [ grid(3,2) b{0.14w} ]

plot(cwA, cwB, cwC, cwD, cwE, cwF, cw_legend, layout=l1, size=(600,600), margin=2.5mm)

plot(cwF, legend=:topleft)

#hide

using Logging

Logging.disable_logging(Logging.Warn);

#hide

xs = range(1,1e5,1000)

function plot_vt!(plt,class)
    plot!(plt, xaxis=:log10, yaxis=:log10, xlabel=L"x", ylabel=L"\sigma", title="$class", legend=:topleft)

    for s in sets
        ys = [ GasDispersion.vertical_dispersion(x, Plume, class, s.eqn) for x in xs]
        plot!(plt, xs, ys, color=s.color, linestyle=s.linestyle, lw=s.lw, label=s.name)
    end
end;

vtA = plot()
plot_vt!(vtA,ClassA)
plot(vtA)

vtB = plot()
plot_vt!(vtB,ClassB)
plot(vtB)

vtC = plot()
plot_vt!(vtC,ClassC)
plot(vtC)

vtD = plot()
plot_vt!(vtD,ClassD)
plot(vtD)

vtE = plot()
plot_vt!(vtE,ClassE)
plot(vtE)

vtF = plot()
plot_vt!(vtF,ClassF)
plot(vtF)

# TransAlta Sundance - Stack 2
m = 3200/3600 # mass emission rate: 3200kg/h in kg/s
h = 155.5 # stack height, m
d = 7.3   # stack diameter, m
v = 35.6  # stack exit velocity, m/s
T = 439.7 # stack exit temperature, K

# assumed weather conditions
uᵣ  = 2  # windspeed, m/s
zᵣ = 10 # windspeed elevation, m
stability = ClassD

# standard state
Pₛ = 101325 # Pa
Tₛ = 273.15 # K

r = VerticalJet(m, Inf, d, v, h, Pₛ, T, 0.0)

a = SimpleAtmosphere(pressure=Pₛ, temperature=Tₛ, windspeed=uᵣ, windspeed_height=zᵣ, stability=stability)

# a dummy substance, since I know a gaussian plume doesn't require any material
# properties I have just left them as NaNs
SO2 = Substance(name=:SulfurDioxide,molar_weight=0.064066,liquid_density=1,boiling_temp=1,
                latent_heat=1,gas_heat_capacity=1,liquid_heat_capacity=1)

scn = Scenario(SO2,r,a)

conc = plume(scn, GaussianPlume; plumerise=false);

#hide

xs = range(1,5e4,1000)
kms = xs./1000

plt = plot()    
plot!(plt, xlabel="Downwind distance (km)", ylabel="Concentration SO₂ (ppbv)", title="Downwind concentration at 2m elevation", legend=:topright)
plot!(plt, [0, 50], [172, 172], color=:grey, linestyle=:dot, lw=1, label="1-hr AAQO SO₂")

for s in sets
    pl = plume(scn, GaussianPlume, s.eqn; plumerise=false)
    plot!(plt, kms, pl.(xs,0,2).*1e9, color=s.color, linestyle=s.linestyle, lw=s.lw, label=s.name)
end

plot(plt)

# hide

using Contour

xs = range(1,2e4,1000)
ys = range(-2e3,2e3,1000)
c = 172e-9

plt = plot()    
plot!(plt, xlabel="Downwind distance (km)", ylabel="Crosswind distance (km)", title="172ppbv isopleth at 2m elevation", legend=:topright)

for s in sets
    pl = plume(scn, GaussianPlume, s.eqn; plumerise=false)
    pl_zs = [pl(xi,yi,2.0) for xi in xs, yi in ys]
    ctr = Contour.contour(xs,ys,pl_zs,c)
    pl_xs, pl_ys = coordinates(first(lines(ctr)))
    plot!(plt, pl_xs.*1e-3, pl_ys.*1e-3, color=s.color, linestyle=s.linestyle, lw=s.lw, label=s.name)
end

plot(plt, xlims=(0,20), ylims=(-1,1))

conc = plume(scn, GaussianPlume; plumerise=true);

#hide
xs = range(1,5e5,1000)
kms = xs./1000

plt = plot()    
plot!(plt, xlabel="Downwind distance (km)", ylabel="Concentration (ppbv)", title="Downwind concentration at 2m elevation,\n with plume rise", legend=:topright)

for s in sets
    pl = plume(scn, GaussianPlume, s.eqn; plumerise=true)
    plot!(plt, kms, pl.(xs,0,2).*1e9, color=s.color, linestyle=s.linestyle, lw=s.lw, label=s.name)
end

plot(plt)