using Plots, Plots.PlotMeasures, Distributions, DifferentialMobilityAnalyzers, Random, LinearAlgebra, Printf
plotlyjs();    

t,p = 295.15, 1e5                                # Temperature [K], Pressure [Pa]
qsa,β = 1.66e-5, 1/5                             # Qsample [m3 s-1], Sample-to-sheath ratio,
r₁,r₂,l = 9.37e-3,1.961e-2,0.44369               # DMA geometry [m]
leff = 13.0                                      # DMA effective diffusion length [m]
m = 3                                            # Upper number of charges
Λ = DMAconfig(t,p,qsa,qsa/β,r₁,r₂,l,leff,:-,m,:cylindrical)   # Specify DMA with negative polarity
bins,z₁,z₂ = 128, dtoz(Λ,1000e-9), dtoz(Λ,10e-9) # bins, upper, lower mobility limit
δ = setupDMA(Λ, z₁, z₂, bins);                   # Compute matrices

figure("Nimbus Sans L", 2, 3, 2, 8)
𝕟 = DMALognormalDistribution([[400, 30, 1.2],[500, 110, 1.7]], δ)
p = plot(𝕟.Dp, 𝕟.S, xaxis = :log10, xticks = [10, 100, 1000], lt = :steppre, right_margin = 35px,
    legend = :none, ylabel = "dN/dlnD (cm⁻³)", xlim = (8,1000), color =  :black,
    xlabel = "Mobility Diameter (nm)", fmt = :svg)

Random.seed!(703)  # random number seed; fix for repeatable runs
tscan = 120        # SMPS scan time [s] 
Qcpc = 16.66       # CPC flow rate [cm³ s⁻¹]. 16.66 cm³ s⁻¹ = 1 L min⁻¹
t = tscan./bins    # time in each bin

𝕣 = δ.𝐀*𝕟;      # DMA response function 
c = 𝕣.N*Qcpc*t; # number of counts in each bin
R = Float64[]   # Construct a noisy response function 𝕣+ϵ, where ϵ is counting uncertainty
for i = c
    f = rand(Poisson(i),1)   # draw Poisson random number with mean = counts
    push!(R,f[1]/(Qcpc*t))   # convert back to concentration
end

p = plot(𝕣.Dp, R, xaxis = :log10, xticks = [10, 100, 1000], lt = :steppre, label = "𝕣 = 𝐀*𝕟+ϵ", 
    ylabel = "Concentration (cm⁻³)", xlim = (8,1000), color =  color = RGBA(0,0,0.8,1),
    xlabel = "Apparent +1 Mobility Diam. (nm)", left_margin = 25px)
p = plot!(𝕣.Dp, 𝕣.N, label = "𝕣 = 𝐀*𝕟", color = :black, ls = :dash, fmt = :svg)

λ₁,λ₂ = 1e-3, 1e1  # bounds [λ₁,λ₂] within which the optimal distribution lies

eyeM = Matrix{Float64}(I, bins, bins)
setupRegularization(δ.𝐀,eyeM,R,inv(δ.𝐒)*R)   # setup the system of equations to regularize
@time L1,L2,λs,ii = lcurve(λ₁,λ₂;n=200)      # compute the L-curve for plotting only
@time λopt = lcorner(λ₁,λ₂;n=10,r=3)         # compute the optimal λ using recursive algorithn
N =  clean((reginv(λopt, r = :Nλ))[1])       # find the inverted size distribution 
𝕟ᵢₙᵥ= SizeDistribution([],𝕟.De,𝕟.Dp,𝕟.ΔlnD,N./𝕟.ΔlnD,N,:regularized) # store as AerosolSizeDistribution

@printf("Input number concentration %i\n", sum(𝕟.N))
@printf("Inverted number concentration %i\n", sum(𝕟ᵢₙᵥ.N))

figure("Nimbus Sans L", 2, 6, 4.5, 10)

p1 = plot(L1, L2, xaxis = :log10, yaxis = :log10, xlabel = "𝐀*(𝕟-𝕣)", ylabel = "𝐈*(𝕟-𝕟ᵢ)", color = :black, 
    label = "L-curve"*(@sprintf(" λ ∈ [%.0e %.0e]", λ₁,λ₂)), bottom_margin = 50px)

p1 = plot!([L1[ii], L1[ii]], [L2[ii], L2[ii]], marker = :square, color = RGBA(0.8,0,0,0.5), 
            label = "Optimum λ ="*(@sprintf("%.1e", λopt)),lw = 0, ms = 4)

p2 = plot(𝕟.Dp, 𝕟.S, xaxis = :log10, xticks = [10, 100, 1000],left_margin = 65px, ylabel = "dN/dlnD (cm⁻³)", 
    label = "𝕟 = 𝓛𝓝 (N,Dₘ,σ)", xlabel = "Mobility Diameter (nm)", xlim = (8,1000), ls = :dash, 
    lw = 3, color = :black)

p2 = plot!(𝕟ᵢₙᵥ.Dp, 𝕟ᵢₙᵥ.S, color = RGBA(0.8,0,0,1), lt = :steppre, 
    label = "N<sub>inv</sub>=(𝐀ᵀ𝐀+λ²𝐈)⁻¹(𝐀ᵀR-λ²𝐒⁻¹R)", xlim = (8,1000))

p3 = plot(𝕣.Dp, R, xaxis = :log10, xticks = [10, 100, 1000], lt = :steppre,
    label = "𝕣 = 𝐀*𝕟+ϵ", ylabel = "Concentration (cm⁻³)", xlim = (8,1000),
    color = RGBA(0,0,0.8,1), xlabel = "Apparent +1 Mobility Diam. (nm)", left_margin = 65px)

plot(p1,p2,p3, layout = (l = @layout [a; b c]), right_margin = 0px, legend=:right, 
    top_margin = 15px, fmt = :svg)