using CloudScat, Random, Statistics, LinearAlgebra

# Runtime observer replacement; the installed package source is not changed.
# The constructor fov argument is interpreted HERE as a circular half-angle.
# Outputs are per square metre of a flat entrance aperture, per emitted photon.
# Time bins are EXCESS path delay, not the stock absolute arrival time.
@eval CloudScat begin
    @inline function observeone!(o::Observer, w::World, s::NonNull,
                                r, μ, t, weight, params::Params, tid)
        sobs = norm(o.r-r)
        μobs = (o.r-r)/sobs
        cθ = -(o.e[3] ⋅ μobs)
        cθ >= cos(atan(sqrt(2.)*o.umax)) || return
        trel = t+sobs/co.c-delay(params.source,o.r)
        trel >= -1e-12 || error("Noncausal path")
        ind = 1+floor(Int,max(0.,trel)/o.δt)
        1 <= ind <= size(o.obs,1) || return
        ν(r::Point) = miecollrate(r,w,params)
        τmie = optpath(w.cloud,ν,r,o.r,w.quadrule,tid)
        τray = params.νray_ground*params.H/μobs[3]*
               (exp(-r[3]/params.H)-exp(-o.r[3]/params.H))
        alb = s isa Mie ? s.ω0 : 1.
        f = weight*alb*p(s,μ⋅μobs)*exp(-τmie-τray)*
            transmittance(w.absorber,r,o.r,μobs)*cθ/(sobs^2*params.N)
        o.obs[ind,tid] += f/o.δt
        nothing
    end
end

const MCROOT = @__DIR__
const DT = 2e-6
const WINDOW = 0.012
const ICE = Homogeneous(337.1e-9,1.,10e-6;refindex="WarrenBrandt.dat",verbose=false)

function makecase(;N=10000,zsrc=16000.,ztop=18000.,zbottom=12000.,
                  beta=0.001,ray=true,g=ICE.g,albedo=ICE.ω0,
                  cloudradius=20000.,domainradius=60000.,domaintop=80000.,
                  domainbottom=zbottom-100.,halfangle=30.,xsrc=0.,seed=1,
                  dt=DT,window=WINDOW)
    r=Point(xsrc,0.,zsrc)
    params=Params(N=N,λ=337.1e-9,zmin=domainbottom,source_a=r,source_b=r,
                  nair=ray ? 2.6867811e25 : 0.,
                  νraymax=ray ? 9.223146377701234e-5*exp(-domainbottom/7200) : 1e-7)
    comp=Homogeneous(beta/(ICE.Qext*π*ICE.radius^2),ICE.radius,g,albedo,ICE.Qext)
    world=World(Cylinder(zbottom,ztop,0.,0.,cloudradius),
                Cylinder(domainbottom,domaintop,0.,0.,domainradius),comp)
    observer=Observer(position=[0.,0.,550000.],fov=halfangle,pixels=1,
                      tsample=dt,nsamples=round(Int,window/dt))
    Random.seed!(seed)
    elapsed=@elapsed result=only(CloudScat.main(params,world,[observer];verbose=false))
    R=norm(observer.r-r)
    τover=beta*clamp(ztop-zsrc,0.,ztop-zbottom)
    μobs=(observer.r-r)/R
    τray=params.νray_ground*params.H*exp(-zsrc/params.H)/μobs[3]
    τlos=CloudScat.optpath(world.cloud,r1->CloudScat.miecollrate(r1,world,params),r,observer.r,world.quadrule,1)
    direct=μobs[3]>=cosd(halfangle) ? μobs[3]*exp(-τlos-τray)/(4π*R^2) : 0.
    return result,elapsed,direct,τover,τray
end

function writecurve(path,result)
    open(path,"w") do io
        println(io,"excess_delay_s,fluence_density_m-2_s-1_per_emitted_photon")
        for (t,y) in zip(result.time,result.photo)
            println(io,t,',',y)
        end
    end
end

if abspath(PROGRAM_FILE)==@__FILE__
    ENV["CLOUDSCAT_DEBUG_REDUCTION"]="1"
    for tau in (0.,1.,10.,30.)
        r,secs,d,_,_=makecase(N=10000,beta=tau/2000,seed=123)
        println("BENCH tau=",tau," seconds=",secs," fluence=",sum(r.photo)*DT," direct=",d)
        flush(stdout)
    end
end