Newer
Older
function read_epochs()
out = DataFrame(id=Int[], filename=String[], date=String[], scale=Float64[], OIII_fwhm=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))
end
irow = findfirst(epochs.id .== id)
file = epochs[irow, :filename]
spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$id: $file", kw...);
if isa(id, Number)
spec.flux ./= epochs[irow, :scale]
spec.err ./= epochs[irow, :scale]
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
function analyze_boller03()
file = opt.path * "/boller.dat"
if !isfile(file)
source = QSO{q1927p654}("1ES 1927+654 (Boller03)", opt.z, ebv=opt.ebv);
source.options[:min_spectral_coverage][:OIII_5007] = 0.35
# source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"
#=
Resolution: from Sect 3 it is 6 / 5007 * 3e5 ~ 360 km/s
however with this value the [OIII] FWHM is 100 km/s (limit)
I tried 180, 200 and 250, and the [OIII] luminosity do not change
=#
spec = read_spec("B03", resolution=200.)
add_spec!(source, spec);
res = qsfit(source);
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
filename="$(opt.path)/results_boller.html")
JLD2.save_object(file, res)
else
res = JLD2.load_object(file)
end
@assert res.model[:OIII_5007].norm.val == 0.00029042977121756283
@assert res.model[:OIII_5007].norm.val / res.model[:galaxy].norm.val == 60.07373299362222
return res
end
function estimate_SE_scale!()
file = opt.path * "/" * "SE_scale.dat"
if !isfile(file)
for irow in 1:nrow(epochs)
@info irow
id = epochs[irow, :id]
epochs[irow, :scale] = 1.
spec = read_spec(id)
epochs[irow, :scale] = 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
# source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"
spec = read_spec(id)
add_spec!(source, spec);
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
epochs[irow, :scale] *= res.model[:OIII_5007].norm.val
epochs[irow, :OIII_fwhm] = res.model[:OIII_5007].fwhm.val
end
@assert all(epochs[.!isfinite.(epochs.OIII_fwhm), :scale] .== NaN)
JLD2.save_object(file, epochs)
else
empty!(epochs)
append!(epochs, JLD2.load_object(file))
end
end
function analyze_job(job; Nloop = 6)
out = fill(NaN, length(epoch_ids[job]))
estimate_SE_scale!()
for loop in 1:Nloop
file = "$(opt.path)/results_$(job)_$(loop).dat"
if !isfile(file)
source = QSO{q1927p654}("1ES 1927+654", opt.z, ebv=opt.ebv);
source.options[:min_spectral_coverage][:OIII_5007] = 0.5
source.options[:wavelength_range] = [3500, 6900]
@gp :zoom "set grid" :-
for id in epoch_ids[job]
spec = read_spec(id, resolution=get(resolution, job, NaN))
add_spec!(source, spec);
@gp :- :zoom xr=[4750,5150] spec.λ ./ (1 + source.z) spec.flux "w l t '$(id)'"
end
res = qsfit_multi(source);
JLD2.save_object(file, res)
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], filename="$(opt.path)/results_$(job)_$(loop).html")
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], filename="$(opt.path)/results_$(job)_$(loop)_rebin4.html", rebin=4)
else
res = JLD2.load_object(file);
end
OIII_norm = Float64[]
for i in 1:length(epoch_ids[job])
id = epoch_ids[job][i]
push!(OIII_norm, res.multi[i][:OIII_5007].norm.val)
epochs[epochs.id .== id, :scale] *= res.multi[i][:OIII_5007].norm.val
out = out[:,2:end]
@gp "set grid" :-
for i in 1:length(epoch_ids[job])
id = epoch_ids[job][i]
@gp :- 1:Nloop out[i, :] "w lp t '$(id)'"
end
return out
const opt = (
z = 0.01942,
ebv = 0.077,
path = "output")
try; mkdir(opt.path); catch; end
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)"
save(:prenorm, term="png size 800,600", output="$(opt.path)/oiii_norm_fwhm.png")
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="Actual 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]"
y ./= maximum(y)
@gp :- :allepochs x ./ (1 + opt.z) y "u 1:(\$2+$id):($id) w l notit lc pal" :-
save(:allepochs, term="png size 800,600", output="$(opt.path)/all_epochs.png")
save(:allepochs, term="png size 800,2000", output="$(opt.path)/all_epochs_zoom.png")
:vecchi => [17, 10, 21, 23],
:lowres => [8,9,10,11,12,14,15,16,18,21,22,23,24,25],
:highres => [5, 7, 13, 17, 20],
# deleteat!(epoch_ids[:all], 4) #insufficient coverage in na_Hg
# Can't use resolution smaller than 150 km / s otherwise some line
# will be neglected because of spectral coverage
resolution = Dict()
# :lowres => 350.,
# :highres => 150.)
function cont_at(bestfit, λ)
A1 = bestfit[:qso_cont].norm.val - bestfit[:qso_cont].norm.unc
A = bestfit[:qso_cont].norm.val
A2 = bestfit[:qso_cont].norm.val + bestfit[:qso_cont].norm.unc
if A1 > A2
A1, A2 = A2, A1
end
@assert A1 <= A2
B1 = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val - bestfit[:qso_cont].alpha.unc)
B = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val)
B2 = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val + bestfit[:qso_cont].alpha.unc)
if B1 > B2
B1, B2 = B2, B1
end
@assert B1 <= B2
dd = extrema([A1*B1, A1*B2, A2*B1, A2*B2])
return λ * A * B, λ * (dd[2]-dd[1])
tab = DataFrame(epoch=Int[], date=String[], instr=String[], galaxy=Float64[],
contslope=Float64[],
c3500=Float64[], c5100=Float64[],
bb_Ha_norm=Float64[], bb_Ha_fwhm=Float64[], bb_Ha_voff=Float64[],
br_Ha_norm=Float64[], br_Ha_fwhm=Float64[], br_Ha_voff=Float64[],
na_Ha_norm=Float64[], na_Ha_fwhm=Float64[], na_Ha_voff=Float64[],
bb_Hb_norm=Float64[], bb_Hb_fwhm=Float64[], bb_Hb_voff=Float64[],
br_Hb_norm=Float64[], br_Hb_fwhm=Float64[], br_Hb_voff=Float64[],
na_Hb_norm=Float64[], na_Hb_fwhm=Float64[], na_Hb_voff=Float64[],
oiii_fwhm=Float64[], oiii_voff=Float64[])
res.bestfit[id][:qso_cont].alpha.val,
cont_at(res.bestfit[id], 3500)[1],
cont_at(res.bestfit[id], 5100)[1],
res.bestfit[id][:bb_Ha].norm.patched, res.bestfit[id][:bb_Ha].fwhm.patched, res.bestfit[id][:bb_Ha].voff.patched,
res.bestfit[id][:br_Ha].norm.patched, res.bestfit[id][:br_Ha].fwhm.patched, res.bestfit[id][:br_Ha].voff.patched,
res.bestfit[id][:na_Ha].norm.patched, res.bestfit[id][:na_Ha].fwhm.patched, res.bestfit[id][:na_Ha].voff.patched,
res.bestfit[id][:bb_Hb].norm.patched, res.bestfit[id][:bb_Hb].fwhm.patched, res.bestfit[id][:bb_Hb].voff.patched,
res.bestfit[id][:br_Hb].norm.patched, res.bestfit[id][:br_Hb].fwhm.patched, res.bestfit[id][:br_Hb].voff.patched,
res.bestfit[id][:na_Hb].norm.patched, res.bestfit[id][:na_Hb].fwhm.patched, res.bestfit[id][:na_Hb].voff.patched,
res.bestfit[id][:OIII_5007].fwhm.patched, res.bestfit[id][:OIII_5007].voff.patched))
tab[!, :l3500] .= [3500 .* Spline1D(domain(res.multi[id])[:], res.multi[id]())(3500.) for id in 1:length(chosen_epochs)]
tab[!, :l5100] .= [5100 .* Spline1D(domain(res.multi[id])[:], res.multi[id]())(5100.) for id in 1:length(chosen_epochs)]
tab[!, :Ha] .= tab.br_Ha_norm .+ tab.bb_Ha_norm .+ tab.na_Ha_norm
tab[!, :Hb] .= tab.br_Hb_norm .+ tab.bb_Hb_norm .+ tab.na_Hb_norm
dd = Date.(tab.date, Ref("yyyymmdd"))
tab[!, :day] .= getproperty.(dd - Date("2018-03-06") .+ Day(72), :value)
(d, c) = csvread("misc/tabula-2019-Trakhtenbrot-mod.csv", ',', header_exists=true)
tmp = DataFrame(collect(d), Symbol.(c))
tmp[!, :Date] .= string.(tmp.Date)
tab[!, :instr] .= ""
for i in 1:nrow(tab)
j = findfirst(tmp.Date .== tab[i, :date])
isnothing(j) && continue
tab[i, :instr] = tmp[j, :Inst]
end
f = FITS("$(opt.path)/1ES_1927p654_results_$(job).fits", "w")
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)
@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"
@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 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 tab.day tab.contslope tab.pt "w lp t 'SLope' pt var ps 2"
for i in 1:nrow(tab)
ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"
@gp :- domain(res.model[i])[:] res.model[i](:qso_cont) ss
save(:prenorm, term="png size 800,600", output="$(opt.path)/evolution.png")
for i in 1:length(chosen_epochs)
println(res.bestfit[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
save(:prenorm, term="png size 800,600", output="$(opt.path)/evolution_oiii_norm.png")