Review index · Download raw source

reference/legacy/photometer_model.jl

1using CloudScat, Random, Statistics, LinearAlgebra23# Runtime observer replacement; the installed package source is not changed.4# The constructor fov argument is interpreted HERE as a circular half-angle.5# Outputs are per square metre of a flat entrance aperture, per emitted photon.6# Time bins are EXCESS path delay, not the stock absolute arrival time.7@eval CloudScat begin8    @inline function observeone!(o::Observer, w::World, s::NonNull,9                                r, μ, t, weight, params::Params, tid)10        sobs = norm(o.r-r)11        μobs = (o.r-r)/sobs12        cθ = -(o.e[3] ⋅ μobs)13        cθ >= cos(atan(sqrt(2.)*o.umax)) || return14        trel = t+sobs/co.c-delay(params.source,o.r)15        trel >= -1e-12 || error("Noncausal path")16        ind = 1+floor(Int,max(0.,trel)/o.δt)17        1 <= ind <= size(o.obs,1) || return18        ν(r::Point) = miecollrate(r,w,params)19        τmie = optpath(w.cloud,ν,r,o.r,w.quadrule,tid)20        τray = params.νray_ground*params.H/μobs[3]*21               (exp(-r[3]/params.H)-exp(-o.r[3]/params.H))22        alb = s isa Mie ? s.ω0 : 1.23        f = weight*alb*p(s,μ⋅μobs)*exp(-τmie-τray)*24            transmittance(w.absorber,r,o.r,μobs)*cθ/(sobs^2*params.N)25        o.obs[ind,tid] += f/o.δt26        nothing27    end28end2930const MCROOT = @__DIR__31const DT = 2e-632const WINDOW = 0.01233const ICE = Homogeneous(337.1e-9,1.,10e-6;refindex="WarrenBrandt.dat",verbose=false)3435function makecase(;N=10000,zsrc=16000.,ztop=18000.,zbottom=12000.,36                  beta=0.001,ray=true,g=ICE.g,albedo=ICE.ω0,37                  cloudradius=20000.,domainradius=60000.,domaintop=80000.,38                  domainbottom=zbottom-100.,halfangle=30.,xsrc=0.,seed=1,39                  dt=DT,window=WINDOW)40    r=Point(xsrc,0.,zsrc)41    params=Params(N=N,λ=337.1e-9,zmin=domainbottom,source_a=r,source_b=r,42                  nair=ray ? 2.6867811e25 : 0.,43                  νraymax=ray ? 9.223146377701234e-5*exp(-domainbottom/7200) : 1e-7)44    comp=Homogeneous(beta/(ICE.Qext*π*ICE.radius^2),ICE.radius,g,albedo,ICE.Qext)45    world=World(Cylinder(zbottom,ztop,0.,0.,cloudradius),46                Cylinder(domainbottom,domaintop,0.,0.,domainradius),comp)47    observer=Observer(position=[0.,0.,550000.],fov=halfangle,pixels=1,48                      tsample=dt,nsamples=round(Int,window/dt))49    Random.seed!(seed)50    elapsed=@elapsed result=only(CloudScat.main(params,world,[observer];verbose=false))51    R=norm(observer.r-r)52    τover=beta*clamp(ztop-zsrc,0.,ztop-zbottom)53    μobs=(observer.r-r)/R54    τray=params.νray_ground*params.H*exp(-zsrc/params.H)/μobs[3]55    τlos=CloudScat.optpath(world.cloud,r1->CloudScat.miecollrate(r1,world,params),r,observer.r,world.quadrule,1)56    direct=μobs[3]>=cosd(halfangle) ? μobs[3]*exp(-τlos-τray)/(4π*R^2) : 0.57    return result,elapsed,direct,τover,τray58end5960function writecurve(path,result)61    open(path,"w") do io62        println(io,"excess_delay_s,fluence_density_m-2_s-1_per_emitted_photon")63        for (t,y) in zip(result.time,result.photo)64            println(io,t,',',y)65        end66    end67end6869if abspath(PROGRAM_FILE)==@__FILE__70    ENV["CLOUDSCAT_DEBUG_REDUCTION"]="1"71    for tau in (0.,1.,10.,30.)72        r,secs,d,_,_=makecase(N=10000,beta=tau/2000,seed=123)73        println("BENCH tau=",tau," seconds=",secs," fluence=",sum(r.photo)*DT," direct=",d)74        flush(stdout)75    end76end77