include("photometer_model.jl")
ENV["CLOUDSCAT_DEBUG_REDUCTION"]="1"
const OUT=joinpath(MCROOT,"results")
mkpath(OUT)

function validate()
    io=open(joinpath(OUT,"validation.txt"),"w")
    for θ in (0.,20.,40.)
        Z=534000.
        r,secs,_,_,_=makecase(N=5000,beta=0.,ray=false,xsrc=Z*tand(θ))
        measured=sum(r.photo)*DT
        expected=θ<=30 ? cosd(θ)^3/(4π*Z^2) : 0.
        @assert isapprox(measured,expected;rtol=1e-10,atol=1e-28)
        println(io,"Vacuum angle=",θ," expected=",expected," measured=",measured," PASS")
    end
    r,_,direct,_,_=makecase(N=5000,beta=3/2000,ray=false,albedo=0.)
    @assert isapprox(sum(r.photo)*DT,direct;rtol=1e-10)
    @assert all(==(0.),r.photo[2:end])
    println(io,"Zero-albedo cloud tau=3: direct-only test PASS")
    println(io,"Binning: 2 us excess delay; 12 ms record. Debug reduction=1.")
    println(io,"Mie g=",ICE.g," albedo=",ICE.ω0," Qext=",ICE.Qext)
    close(io)
    println("VALIDATION PASS"); flush(stdout)
end

function rungrid()
    cases=NamedTuple[]
    for τ in (0.,.1,.3,1.,3.,10.,30.,100.)
        push!(cases,(id="tau_"*replace(string(τ),'.'=>'p'),group="tau_grid",
                     config=(beta=τ/2000,),N=τ>=100 ? 25000 : 50000))
    end
    for d in (200.,1000.,3000.)
        push!(cases,(id="depth_"*string(Int(d)),group="same_tau10",
                     config=(zsrc=18000-d,beta=10/d),N=50000))
    end
    push!(cases,(id="top",group="top_above",config=(zsrc=18000.,beta=.005),N=50000))
    push!(cases,(id="above_1km",group="top_above",config=(zsrc=18000.,ztop=17000.,zbottom=11000.,beta=.005),N=50000))
    for τ in (0.,3.,10.,30.)
        push!(cases,(id="noray_"*string(Int(τ)),group="ray_test",config=(beta=τ/2000,ray=false),N=50000))
    end
    for g in (.75,.95)
        push!(cases,(id="g_"*replace(string(g),'.'=>'p'),group="phase_test",config=(beta=.005,g=g),N=50000))
    end
    for τ in (0.,10.,30.)
        push!(cases,(id="domain_"*string(Int(τ)),group="domain_test",
              config=(beta=τ/2000,domainradius=200000.,domaintop=120000.),N=50000))
    end
    push!(cases,(id="radius40_tau30",group="radius_test",config=(beta=.015,cloudradius=40000.,domainradius=100000.),N=50000))
    push!(cases,(id="fov90_tau10",group="fov_test",config=(beta=.005,halfangle=45.),N=50000))
    # Fine near-top case tests that small overlying tau is not necessarily a pure direct pulse.
    push!(cases,(id="near_top50m",group="top_above",config=(zsrc=17950.,beta=.005),N=50000))
    open(joinpath(OUT,"batch_summary.csv"),"w") do io
        println(io,"case_id,group,batch,N,seed,zsrc_m,ztop_m,zbottom_m,beta_m-1,g,albedo,ray,cloudradius_m,domainradius_m,domaintop_m,full_fov_deg,tau_over,tau_cloud_total,tau_ray,direct_fluence_m-2,total_fluence_12ms_m-2,elapsed_s")
        for (i,c) in enumerate(cases)
            for batch in 1:4
                seed=20261021+100*i+batch
                kw=merge((zsrc=16000.,ztop=18000.,zbottom=12000.,beta=.001,g=ICE.g,albedo=ICE.ω0,ray=true,cloudradius=20000.,domainradius=60000.,domaintop=80000.,halfangle=30.),c.config)
                r,secs,direct,τover,τray=makecase(;N=c.N,seed=seed,kw...)
                filename=c.id*"_batch"*string(batch)*".csv"
                writecurve(joinpath(OUT,filename),r)
                fields=(c.id,c.group,batch,c.N,seed,kw.zsrc,kw.ztop,kw.zbottom,kw.beta,kw.g,kw.albedo,kw.ray,kw.cloudradius,kw.domainradius,kw.domaintop,2kw.halfangle,τover,kw.beta*(kw.ztop-kw.zbottom),τray,direct,sum(r.photo)*DT,secs)
                println(io,join(string.(fields),',')); flush(io)
                println("DONE ",c.id," batch",batch," N=",c.N," seconds=",round(secs;digits=2)," total/vac=",sum(r.photo)*DT*4π*(550000-kw.zsrc)^2); flush(stdout)
            end
        end
    end
end

validate()
rungrid()