Newer
Older
using Revise, Statistics, Serialization, Dierckx, DataFrames
using QSFit, GFit, Gnuplot, GFitViewer, MyAstroUtils
for (root, dirs, files) in walkdir("AT2018zf")
for file in files
if file[end-3:end] == ".dat"
if epoch_id == "B03"
file = "boller2003.txt"
else
file = epoch_filenames[epoch_id]
end
spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$epoch_id: $file", kw...);
source = QSO{q1927p654}("1ES 1927+654 (Boller03)", Z, ebv=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 = fit(source);
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
println(res.bestfit[:OIII_5007].norm.val / res.bestfit[:galaxy].norm.val )
if !isfile("output/scale_fwhm.dat")
all_scale = fill(1., length(epoch_filenames))
all_fwhm = fill(NaN,length(epoch_filenames))
for id in 1:length(epoch_filenames)
all_scale[id] *= median(spec.flux) / 25000
try
source = QSO{q1927p654}("1ES 1927+654 ($id)", Z, ebv=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);
res = fit(source);
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
filename="output/results_$(id).html")
# Update all_scale and store FWHM
@info res.bestfit[:OIII_5007].norm.val
all_scale[id] *= res.bestfit[:OIII_5007].norm.val
all_fwhm[ id] = res.bestfit[:OIII_5007].fwhm.val
all_scale[.!isfinite.(all_fwhm)] .= NaN
serialize("output/scale_fwhm.dat", (all_scale, all_fwhm))
else
(all_scale, all_fwhm) = deserialize("output/scale_fwhm.dat")
end
@gp :prenorm "set key bottom right" "set grid" xlab="Epoch" ylab="Value" :-
@gp :- :prenorm all_scale./1e-18 "w lp t '[OIII] norm.'" :-
@gp :- :prenorm all_fwhm./1000 "w lp t '[OIII] FWHM (x1000 km/s)"
save(:prenorm, term="png size 800,600", output="output/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]"
for id in 1:length(epoch_filenames)
s = read_spec(id)
x = s.λ;
y = s.flux;
y ./= maximum(y)
@gp :- :allepochs x./(1+Z) y "u 1:(\$2+$id):($id) w l notit lc pal" :-
end
@gp :- :allepochs yr=[1, 29]
save(:allepochs, term="png size 800,600", output="output/all_epochs.png")
@gp :- :allepochs xr=[4900, 5100]
save(:allepochs, term="png size 800,2000", output="output/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],
:limited => [1,2,3,13],
:all => collect(5:26))
# deleteat!(dict_chosen_epochs[: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
(all_scale, all_fwhm) = deserialize("output/scale_fwhm.dat")
file = "output/results_$(job)_$(loop).dat"
if !isfile(file)
source = QSO{q1927p654}("1ES 1927+654", Z, ebv=EBV);
source.options[:min_spectral_coverage][:OIII_5007] = 0.5
@gp :zoom "set grid" :-
for id in chosen_epochs
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)'"
res = multi_fit(source);
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
filename="output/results_$(job)_$(loop).html")
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
filename="output/results_$(job)_$(loop)_rebin4.html", rebin=4)
# Find best normalization for [OIII]
model2 = deepcopy(res.model);
for id in 1:length(chosen_epochs)
for cname in collect(keys(model2[id]))
if cname == :OIII_5007
model2[id][cname].norm.fixed = false
else
freeze(model2[id], cname)
end
mzer = GFit.cmpfit()
mzer.config.ftol = mzer.config.gtol = mzer.config.xtol = 1.e-6
bestfit2 = fit!(model2, source.data, minimizer=mzer);
OIII_norm = fill(0., length(chosen_epochs))
for id in 1:length(chosen_epochs)
OIII_norm[id] = bestfit2[id][:OIII_5007].norm.val
all_scale[chosen_epochs[id]] *= bestfit2[id][:OIII_5007].norm.val
end
file = "output/results_$(job)_$(loop).dat"
(res, all_scale, OIII_norm) = deserialize(file);
OIII_norm_evol = hcat(OIII_norm_evol, OIII_norm)
end
@gp :- 1:Nloop OIII_norm_evol[id, :] "w lp t '$(chosen_epochs[id])'"
end
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.model[id])[:], res.model[id]())(3500.) for id in 1:length(chosen_epochs)]
tab[!, :l5100] .= [5100 .* Spline1D(domain(res.model[id])[:], res.model[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
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
cm = countmap(tab.instr)
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="output/evolution.png")
for i in 1:length(chosen_epochs)
println(res.bestfit[i][:OIII_5007])
x = domain(res.model[i])[:]
y = res.model[i]()
y0 = y .- res.model[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="output/evolution_oiii_norm.png")