Commit 18319596 authored by Giorgio Calderone's avatar Giorgio Calderone
Browse files

Updated

parent 4bdc7965

run.jl

+84 −29
Original line number Diff line number Diff line
@@ -13,17 +13,33 @@ if readlines(`hostname`)[1] != "giorgio-pc" 
end



function read_epochs()

    out = DataFrame(id=Int[], filename=String[], date=String[], scale=Float64[], OIII_fwhm=Float64[])

    out = DataFrame(id=Int[], filename=String[], date=String[], scale=Float64[], scale_rel_unc=Float64[], OIII_fwhm=Float64[], OIII_fwhm_unc=Float64[])

    for (root, dirs, files) in walkdir("AT2018zf")

        for file in files

            if file[end-3:end] == ".dat"

                filename = root * "/" * file

                date = filename[27:end-4]

                push!(out, (nrow(out)+1, filename, date, NaN, NaN))

                push!(out, (nrow(out)+1, filename, date, NaN, NaN, NaN, NaN))

            end

        end

    end

    @assert issorted(out.date)



    dd = Date.(out.date, Ref("yyyymmdd"))

    out[!, :day] .= getproperty.(dd - dd[1], :value)

    @assert issorted(out.day)



    (d, c) = csvread("misc/tabula-2019-Trakhtenbrot-mod.csv", ',', header_exists=true)

    tmp = DataFrame(collect(d), Symbol.(c))

    tmp[!, :Date] .= string.(tmp.Date)



    out[!, :instr] .= ""

    for i in 1:nrow(out)

        j = findfirst(tmp.Date .== out[i, :date])

        isnothing(j)  &&  continue

        out[i, :instr] = tmp[j, :Inst]

    end



    return out

end
@@ -78,7 +94,7 @@ function estimate_SE_scale!() 
            id = epochs[irow, :id]

            epochs[irow, :scale] = 1.

            spec = read_spec(id)

            epochs[irow, :scale] = median(spec.flux) / 25000

            epochs[irow, :scale] = 1.44665e-15 / 25000 # median(spec.flux) / 25000



            source = QSO{q1927p654}("1ES 1927+654 ($id)", opt.z, ebv=opt.ebv);

            source.options[:min_spectral_coverage][:OIII_5007] = 0.5
@@ -92,7 +108,9 @@ function estimate_SE_scale!() 
            # Update all_scale and store FWHM

            @info res.model[:OIII_5007].norm.val

            epochs[irow, :scale] *= res.model[:OIII_5007].norm.val

            epochs[irow, :scale_rel_unc] = res.model[:OIII_5007].norm.unc / res.model[:OIII_5007].norm.val

            epochs[irow, :OIII_fwhm] = res.model[:OIII_5007].fwhm.val

            epochs[irow, :OIII_fwhm_unc] = res.model[:OIII_5007].fwhm.unc

        end

        @assert all(epochs[.!isfinite.(epochs.OIII_fwhm), :scale] .== NaN)

        JLD2.save_object(file, epochs)
@@ -139,12 +157,13 @@ function analyze_job(job; Nloop = 6) 
    end

    out = out[:,2:end]



    @gp "set grid" :-

    @gp "set grid" "set key outside top horiz" xlab="Iteration" ylab="[OIII]λ5007 norm. (arb. units)" :-

    for i in 1:length(epoch_ids[job])

        id = epoch_ids[job][i]

        @gp :- 1:Nloop out[i, :] "w lp t '$(id)'" :-

        @gp :- 1:Nloop out[i, :] "w lp t '$(epochs[id, :date])'" :-

    end

    @gp

    save(term="pngcairo size 1000,700", output="$(opt.path)/calib_$(job).png")



    return out

end
@@ -162,9 +181,17 @@ const epochs = read_epochs() 
analyze_boller03()

estimate_SE_scale!()



@gp :prenorm "set key bottom right" "set grid" xlab="Epoch" ylab="Value" :-

@gp :- :prenorm epochs.scale./1e-18 "w lp t '[OIII] norm.'" :-

@gp :- :prenorm epochs.OIII_fwhm./1000 "w lp t '[OIII] FWHM (x1000 km/s)"

@gp :prenorm "set key top right" "set grid" xlab="Epoch" ylab="Value" :-

@gp :- :prenorm 1:nrow(epochs) epochs.scale./1e-18 epochs.scale./1e-18                         "w lp   t '[OIII] norm. (arbitrary scale)' lc rgb 'red'" :-

@gp :- :prenorm 1:nrow(epochs) epochs.scale./1e-18 epochs.scale./1e-18 .* epochs.scale_rel_unc "w yerr notit pt 0                         lc rgb 'red'" :-

@gp :- :prenorm 1:nrow(epochs) epochs.OIII_fwhm./1000                            "w lp t '[OIII] FWHM (x1000 km/s)' lc rgb 'blue'" :-

@gp :- :prenorm 1:nrow(epochs) epochs.OIII_fwhm./1000 epochs.OIII_fwhm_unc./1000 "w yerr notit pt 0                 lc rgb 'blue'"

i = findall(epochs.instr .== "FLOYDS")

@gp :- :prenorm (1:nrow(epochs))[i] epochs.scale[i]./1e-18                       "w p t 'FLOYDS spectrograph' pt 6 ps 3 lc rgb 'black'"

@gp :- :prenorm (1:nrow(epochs))[i] epochs.OIII_fwhm[i]./1000                    "w p notit                   pt 6 ps 3 lc rgb 'black'"

i = findall(epochs.date .== "20180727")

@gp :- :prenorm (1:nrow(epochs))[i] epochs.scale[i]./1e-18                       "w p t 'Epoch: 20180727' pt 3 ps 4 lc rgb 'black'"

@gp :- :prenorm (1:nrow(epochs))[i] epochs.OIII_fwhm[i]./1000                    "w p               notit pt 3 ps 4 lc rgb 'black'"

save(:prenorm, term="png size 800,600", output="$(opt.path)/oiii_norm_fwhm.png")



w = 10 .^(2:0.01:log10(400));
@@ -213,7 +240,7 @@ resolution = Dict() 


analyze_job(:test)

analyze_job(:highres)

analyze_job(:all)

# analyze_job(:all)  # This takes a really long time

analyze_job(:floyds)

ddd
@@ -310,36 +337,64 @@ setindex!.(Ref(cm), 1:length(cm), collect(keys(cm))[sortperm(collect(values(cm)) 
tab.pt .= getindex.(Ref(cm), tab.instr)





@gp    tab.contslope tab.Ha    tab.pt "w lp t 'Ha'    pt var ps 2"  ylog=true

@gp :- tab.contslope tab.Hb    tab.pt "w lp t 'Hb'    pt var ps 2"

@gp :- tab.contslope tab.c5100 tab.pt "w lp t 'c5100' pt var ps 2"

@gp :- tab.contslope tab.l5100 tab.pt "w lp t 'l5100' pt var ps 2"

@gp :- tab.contslope tab.c3500 tab.pt "w lp t 'c3500' pt var ps 2"

@gp :- tab.contslope tab.l3500 tab.pt "w lp t 'l3500' pt var ps 2"

# tab = tab[findall(tab.date .!= "20180727"), :]



# @gp "set grid"

# @gp    tab.contslope tab.Ha    tab.pt "w p t 'Ha'    pt var ps 2" ylog=true

# @gp :- tab.contslope tab.Hb    tab.pt "w p t 'Hb'    pt var ps 2"





@gp "set grid" xlab="Cont. slope" ylab="Luminosity [arb. units]" linetypes(:darkrainbow, lw=3)

@gp :- tab.contslope tab.Ha             "w l notit   lc rgb 'gray'" ylog=true

@gp :- tab.contslope tab.Hb             "w l notit   lc rgb 'gray'"

@gp :- tab.contslope tab.Ha 1:nrow(tab) "w p t 'H{/Symbol a} int. lum.'  pt 1 ps 3 lc var" ylog=true

@gp :- tab.contslope tab.Hb 1:nrow(tab) "w p t 'H{/Symbol b} int. lum.'  pt 2 ps 3 lc var"

save(term="pngcairo size 1000,600 fontscale 1.2", output="$(opt.path)/slope_vs_lines.png")





@gp    xlog=true ylog=true

@gp :- tab.c5100 tab.Ha tab.pt "w lp t 'c5100' pt var ps 2"

@gp :- tab.l5100 tab.Ha tab.pt "w lp t 'l5100' pt var ps 2"

@gp "set grid" xlab="Cont. slope" ylab="Luminosity [arb. units]" linetypes(:darkrainbow, lw=3)

@gp :- tab.contslope tab.c3500             "w l notit       lc rgb 'gray'" ylog=true

@gp :- tab.contslope tab.c5100             "w l notit dt 2  lc rgb 'gray'"

@gp :- tab.contslope tab.c3500 1:nrow(tab) "w p t 'Cont. {/Symbol l}L_{/Symbol l} at 3500A'  pt 1 ps 3 lc var"

@gp :- tab.contslope tab.c5100 1:nrow(tab) "w p t 'Cont. {/Symbol l}L_{/Symbol l} at 5100A'  pt 2 ps 3 lc var"

save(term="pngcairo size 1000,600 fontscale 1.2", output="$(opt.path)/slope_vs_cont.png")





@gp    tab.day tab.c5100 tab.pt "w lp t 'c5100' pt var ps 2" ylog=true

@gp :- tab.day tab.l5100 tab.pt "w lp t 'l5100' pt var ps 2"

@gp :- tab.day tab.Ha    tab.pt "w lp t 'Ha'    pt var ps 2"

@gp :- tab.day tab.Hb    tab.pt "w lp t 'Hb'    pt var ps 2"

@gp "set grid" xlab="{/Symbol l}L_{/Symbol l} [arb. units]" ylab="Line int. lum. [arb. units]"  linetypes(:darkrainbow, lw=3) xlog=true ylog=true

@gp :- tab.c3500 tab.Ha             "w l notit       lc rgb 'gray'"

@gp :- tab.c5100 tab.Ha             "w l notit dt 2  lc rgb 'gray'"

@gp :- tab.c3500 tab.Ha 1:nrow(tab) "w p t 'Cont. {/Symbol l}L_{/Symbol l} at 3500A' pt 1 ps 3 lc var"

@gp :- tab.c5100 tab.Ha 1:nrow(tab) "w p t 'Cont. {/Symbol l}L_{/Symbol l} at 5100A' pt 2 ps 3 lc var"

save(term="pngcairo size 1000,600 fontscale 1.2", output="$(opt.path)/cont_vs_lines.png")



@gp  tab.day tab.contslope tab.pt "w lp t 'SLope' pt var ps 2"



@gp "set grid" xlab="Days since 2018-03-06" ylab="L_{/Symbol l} [arb. units]"  linetypes(:darkrainbow, lw=3) ylog=true

@gp :- tab.day tab.c3500             "w l notit       lc rgb 'gray'"

@gp :- tab.day tab.c5100             "w l notit dt 2  lc rgb 'gray'"

@gp :- tab.day tab.Ha                "w l notit dt 3  lc rgb 'gray'"

@gp :- tab.day tab.Hb                "w l notit dt 4  lc rgb 'gray'"

@gp :- tab.day tab.c3500 1:nrow(tab) "w p t 'Cont. {/Symbol l}L_{/Symbol l} at 3500A' pt 1 ps 3 lc var"

@gp :- tab.day tab.c5100 1:nrow(tab) "w p t 'Cont. {/Symbol l}L_{/Symbol l} at 5100A' pt 2 ps 3 lc var"

@gp :- tab.day tab.Ha    1:nrow(tab) "w p t 'H{/Symbol a} int. lum.' pt 3 ps 3 lc var"

@gp :- tab.day tab.Hb    1:nrow(tab) "w p t 'H{/Symbol b} int. lum.' pt 4 ps 3 lc var"

save(term="pngcairo size 1000,600 fontscale 1.2", output="$(opt.path)/evolution_day.png")



@gp "set grid"







@gp "set grid" xlab="Rest-frame wavelength [A]" ylab="L_{/Symbol l} [arb. units]"

@gp :- "set key out box" linetypes(:darkrainbow, lw=3)

for i in 1:nrow(tab)

    ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"

    @gp :- domain(res.multi[i])[:] res.multi[i](:qso_cont) ss

    ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt])"

    @gp :- domain(res.multi[i])[:] res.multi[i]() ss

end

save(:prenorm, term="png size 800,600", output="$(opt.path)/evolution.png")

i = 1

@gp :- domain(res.multi[i])[:] res.multi[i](:galaxy) "w l t 'Host galaxy' lc rgb 'black'"

save(term="pngcairo size 1200,600 fontscale 1.2", output="$(opt.path)/evolution.png")





@gp "set grid"

@gp "set grid" xlab="Rest-frame wavelength [A]" ylab="L_{/Symbol l} - cont. at 5007A [arb. units]"

@gp :- "set key out box" linetypes(:darkrainbow, lw=3) xr=[4600,5300]

for i in 1:length(res.multi)

    println(res.multi[i][:OIII_5007])

    x  = domain(res.multi[i])[:]
@@ -348,4 +403,4 @@ for i in 1:length(res.multi) 
    ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"

    @gp :- x y .- Spline1D(x, y0)(5007.) ss

end

save(:prenorm, term="png size 800,600", output="$(opt.path)/evolution_oiii_norm.png")
 
save(term="pngcairo size 1200,600 fontscale 1.2", output="$(opt.path)/evolution_oiii_norm.png")