Newer
Older
include("funcs.jl")
term_hw = "pngcairo size 1000,700 fontscale 1.5"
term_fw = "pngcairo size 1200,600 fontscale 1.3"
z = 0.01942,
ebv = 0.077,
path = "output",
subsets = Dict())
mkpath(opt.path)
boller03_analysis(opt)
estimate_single_epoch_scale!(opt, epochs)
f = open("$(opt.path)/scales_SE.txt", "w")
show(f, epochs)
close(f)
@gp "set key top right" "set grid" xlab="Epoch" ylab="Value" :-
@gp :- 1:nrow(epochs) epochs.scale./1e-18 epochs.scale./1e-18 "w lp t '[OIII]λ5007 norm. (arbitrary scale)' lc rgb 'red'" :-
@gp :- 1:nrow(epochs) epochs.scale./1e-18 epochs.scale./1e-18 .* epochs.scale_rel_unc "w yerr notit pt 0 lc rgb 'red'" :-
@gp :- 1:nrow(epochs) epochs.OIII_fwhm./1000 "w lp t '[OIII]λ5007 FWHM (x1000 km/s)' lc rgb 'blue'" :-
@gp :- 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 :- (1:nrow(epochs))[i] epochs.scale[i]./1e-18 "w p t 'FLOYDS spectrograph' pt 6 ps 3 lc rgb 'black'"
@gp :- (1:nrow(epochs))[i] epochs.OIII_fwhm[i]./1000 "w p notit pt 6 ps 3 lc rgb 'black'"
i = findfirst(epochs.date .== "20180727")
@gp :- yr=[-0.3, 1.7]
@gp :- "set arrow 1 from $i,-0.2 to $i,0 filled size screen 0.03,15,1 lw 3 lc rgb 'black'"
save(term=term_hw, output="$(opt.path)/oiii_norm_fwhm_SE.png")
#=
Check the loci of the true source FWHM compatible with the observed
FWHM of ~400 and ~900 km/s, at different instrument resolutions. The
instrument resolutions (Y axis) when projected to the X axis should
identify a single value for the "true" FWHM".
=#
w = 10 .^(2:0.01:log10(400));
r = sqrt.(400^2 .- w.^2) ./ 2.355;
@gp w r "w l t '400 km/s'"
w = 10 .^(2:0.01:log10(900));
r = sqrt.(900^2 .- w.^2) ./ 2.355;
@gp :- w r "w l t '900 km/s'"
@gp :- "set grid" xlab="True FWHM" ylab="Instr. resolution"
# Ensure no sample is negative
@gp :allepochs "set grid" xr=[3e3,1e4] cbr=[1,29] :-
@gp :- :allepochs cblabel="Epoch" xlab="Rest frame wavelength[A]" ylab="Flux density [arb.units]"
for i in 1:nrow(epochs)
s = read_spec(epochs[i, :])
y ./= maximum(y)
@gp :- :allepochs x ./ (1 + opt.z) y "u 1:(\$2+$i):($i) w l notit lc pal" :-
@gp :- :allepochs yr=[0, 29]
save(:allepochs, term=term_hw, output="$(opt.path)/all_epochs.png")
save(:allepochs, term=term_hw, output="$(opt.path)/all_epochs_zoom.png")
# Define epoch subsets
opt.subsets[:lowres] = [8,9,10,11,12,14,15,16,18,21,22,23,24,25]
opt.subsets[:highres] = [5, 7, 13, 17, 20]
opt.subsets[:test] = [7, 13]
opt.subsets[:floyds] = [6, 8, 9, 10, 11, 12, 14, 15, 16, 18, 19, 21, 22, 23, 24, 25, 26]
multi_epoch_analysis!(opt, epochs, :test)
multi_epoch_analysis!(opt, epochs, :highres)
multi_epoch_analysis!(opt, epochs, :floyds)
# Prepare final table
chosen_job = :floyds
res = JLD2.load_object("output/results_$(chosen_job)_$(chosen_iloop).dat")
tab = DataFrame(epoch=Int[], date=String[], day=Int[], instr=String[], galaxy=Float64[],
Ha_bb_norm=Float64[], Ha_bb_fwhm=Float64[], Ha_bb_voff=Float64[],
Ha_br_norm=Float64[], Ha_br_fwhm=Float64[], Ha_br_voff=Float64[],
Ha_na_norm=Float64[], Ha_na_fwhm=Float64[], Ha_na_voff=Float64[],
Hb_bb_norm=Float64[], Hb_bb_fwhm=Float64[], Hb_bb_voff=Float64[],
Hb_br_norm=Float64[], Hb_br_fwhm=Float64[], Hb_br_voff=Float64[],
Hb_na_norm=Float64[], Hb_na_fwhm=Float64[], Hb_na_voff=Float64[],
j = opt.subsets[chosen_job][id]
(j, epochs[j, :date], epochs[j, :day], epochs[j, :instr],
res.multi[id][:qso_cont].alpha.val,
cont_at(res.multi[id], 3500)[1],
cont_at(res.multi[id], 5100)[1],
res.multi[id][:Ha_bb].norm.patched, res.multi[id][:Ha_bb].fwhm.patched, res.multi[id][:Ha_bb].voff.patched,
res.multi[id][:Ha_br].norm.patched, res.multi[id][:Ha_br].fwhm.patched, res.multi[id][:Ha_br].voff.patched,
res.multi[id][:Ha_na].norm.patched, res.multi[id][:Ha_na].fwhm.patched, res.multi[id][:Ha_na].voff.patched,
res.multi[id][:Hb_bb].norm.patched, res.multi[id][:Hb_bb].fwhm.patched, res.multi[id][:Hb_bb].voff.patched,
res.multi[id][:Hb_br].norm.patched, res.multi[id][:Hb_br].fwhm.patched, res.multi[id][:Hb_br].voff.patched,
res.multi[id][:Hb_na].norm.patched, res.multi[id][:Hb_na].fwhm.patched, res.multi[id][:Hb_na].voff.patched,
res.multi[id][:OIII_5007].fwhm.patched, res.multi[id][:OIII_5007].voff.patched))
tab[!, :l3500] .= [3500 .* Spline1D(domain(res.multi[id])[:], res.multi[id]())(3500.) for id in 1:length(res.multi)]
tab[!, :l5100] .= [5100 .* Spline1D(domain(res.multi[id])[:], res.multi[id]())(5100.) for id in 1:length(res.multi)]
tab[!, :Ha] .= tab.Ha_br_norm .+ tab.Ha_bb_norm .+ tab.Ha_na_norm
tab[!, :Hb] .= tab.Hb_br_norm .+ tab.Hb_bb_norm .+ tab.Hb_na_norm
# Assign different plot symbols to each instrument
tab[!, :pt] .= 0
for i in 1:nrow(tab)
j = findfirst(sort(unique(tab.instr)) .== tab[i, :instr])
isnothing(j) && continue
tab[i, :pt] = j
end
setindex!.(Ref(cm), 1:length(cm), collect(keys(cm))[sortperm(collect(values(cm)), rev=true)])
tab.pt .= getindex.(Ref(cm), tab.instr)
# PLOTS
# Assessment of quality of cross calibration
@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])"
@gp :- domain(res.multi[i])[:] res.multi[i]() ss
end
i = 1
@gp :- domain(res.multi[i])[:] res.multi[i](:galaxy) "w l t 'Host galaxy' lc rgb 'black'"
save(term=term_fw, output="$(opt.path)/evolution.png")
@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])[:]
y = res.multi[i]()
y0 = y .- res.multi[i](:OIII_5007)
ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"
@gp :- x y .- Spline1D(x, y0)(5007.) ss
end
save(term="pngcairo size 1200,600 fontscale 1.2", output="$(opt.path)/evolution_oiii_norm.png")
@gp :- xr=[4900, 5100]
save(term=term_fw, output="$(opt.path)/evolution_oiii_norm_zoom.png")
@gp "set grid" xlab="Cont. slope" ylab="Line 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}' pt 1 ps 3 lc var" ylog=true
@gp :- tab.contslope tab.Hb 1:nrow(tab) "w p t 'H{/Symbol b}' pt 2 ps 3 lc var"
i = findfirst(tab.date .== "20180727")
@gp :- tab.contslope[i] tab.Ha[i] "w p notit pt 6 ps 4 lc rgb 'black'"
@gp :- tab.contslope[i] tab.Hb[i] "w p notit pt 6 ps 4 lc rgb 'black'"
save(term=term_fw, output="$(opt.path)/slope_vs_lines.png")
@gp "set grid" xlab="Cont. slope" ylab="AGN cont. 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 '{/Symbol l}L_{/Symbol l} at 3500A' pt 1 ps 3 lc var"
@gp :- tab.contslope tab.c5100 1:nrow(tab) "w p t '{/Symbol l}L_{/Symbol l} at 5100A' pt 2 ps 3 lc var"
@gp :- yr=[300, 2500]
i = findfirst(tab.date .== "20180727")
@gp :- tab.contslope[i] tab.c3500[i] "w p notit pt 6 ps 4 lc rgb 'black'"
@gp :- tab.contslope[i] tab.c5100[i] "w p notit pt 6 ps 4 lc rgb 'black'"
save(term=term_fw, output="$(opt.path)/slope_vs_cont.png")
@gp "set grid" "set key left" xlab="{/Symbol l}L_{/Symbol l} [arb. units]" ylab="H{/Symbol a} luminosity [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 '{/Symbol l}L_{/Symbol l} at 3500A' pt 1 ps 3 lc var"
@gp :- tab.c5100 tab.Ha 1:nrow(tab) "w p t '{/Symbol l}L_{/Symbol l} at 5100A' pt 2 ps 3 lc var"
@gp :- xr=[300, 2.5e3] yr=[6,40]
i = findfirst(tab.date .== "20180727")
@gp :- tab.c3500[i] tab.Ha[i] "w p notit pt 6 ps 4 lc rgb 'black'"
@gp :- tab.c5100[i] tab.Ha[i] "w p notit pt 6 ps 4 lc rgb 'black'"
save(term=term_fw, output="$(opt.path)/cont_vs_lines.png")
@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"
i = findfirst(tab.date .== "20180727")
@gp :- tab.day[[i,i]] tab.c3500[[i,i]] "w p notit pt 6 ps 4 lc rgb 'black'"
@gp :- tab.day[[i,i]] tab.c5100[[i,i]] "w p notit pt 6 ps 4 lc rgb 'black'"
@gp :- tab.day[[i,i]] tab.Ha[[i,i]] "w p notit pt 6 ps 4 lc rgb 'black'"
@gp :- tab.day[[i,i]] tab.Hb[[i,i]] "w p notit pt 6 ps 4 lc rgb 'black'"
save(term=term_fw, output="$(opt.path)/evolution_day.png")