run.jl

using Pkg
Pkg.activate(".")

Z = 0.019422
EBV = 0.077

using Revise, Statistics, Serialization, Dierckx, DataFrames
using QSFit, GFit, Gnuplot, GFitViewer, MyAstroUtils
using CL_1ES_1927p654
using Dates


epoch_filenames = Vector{String}()
for (root, dirs, files) in walkdir("AT2018zf")
    for file in files
        if file[end-3:end] == ".dat"
            push!(epoch_filenames, root * "/" * file)
        end
    end
end



function read_spec(epoch_id; kw...)
    if epoch_id == "B03"
        file = "boller2003.txt"
    else
        file = epoch_filenames[epoch_id]
    end
    println(file)
	spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$epoch_id: $file", kw...);
    if isa(epoch_id, Number)
        spec.flux ./= all_scale[epoch_id]
    end
    spec.err .= 0.05 .* median(spec.flux);
    return spec
end


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],
       filename="output/results_boller.html")
println(res.bestfit[:OIII_5007].norm.val )
# 0.00029085925218136
println(res.bestfit[:OIII_5007].norm.val / res.bestfit[:galaxy].norm.val )
# 59.406476855719376



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)
        spec = read_spec(id)
        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);
            (model, bestfit) = fit(source);
            viewer(model, source, bestfit, showcomps=[:qso_cont, :galaxy, :balmer],
                   filename="output/results_$(id).html")

            # Update all_scale and store FWHM
            @info bestfit[:OIII_5007].norm.val
            all_scale[id] *= bestfit[:OIII_5007].norm.val
            all_fwhm[ id]  = bestfit[:OIII_5007].fwhm.val
        catch
        end
    end
    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 .-= median(y)
    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")



dict_chosen_epochs = Dict(
    :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
resolution = Dict(
    :lowres  => 350.,
    :highres => 150.)


job = :all
chosen_epochs = dict_chosen_epochs[job]
Nloop = 6
for loop in 1:Nloop
    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)'"
    end
    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)

    models = Dict()
    for id in 1:length(chosen_epochs)
        models[(id, :x)] = domain(model[id])[:]
        models[(id, :y)] = model[id]()
        models[(id, :l5100)] = 5100 .* Spline1D(domain(model[id])[:], model[id]())(5100.)

        for cname in [:qso_cont, :OIII_5007,
                      :br_Ha, :na_Ha, :bb_Ha,
                      :br_Hb, :na_Hb, :bb_Hb]
            models[(id, cname)] = model[id](cname)
        end
    end
    serialize("output/results_$(job)_$(loop).dat", (models, model, bestfit))

    model2 = deepcopy(model)
    for id in 1:length(chosen_epochs)
        @info id
        for cname in collect(keys(model2[id]))
            if cname == :OIII_5007
                model2[id][cname].norm.fixed = false
            else
                freeze(model2[id], cname)
            end
        end
    end
    mzer = GFit.cmpfit()
    mzer.config.ftol = mzer.config.gtol = mzer.config.xtol = 1.e-6
    bestfit2 = fit!(model2, source.data, minimizer=mzer)
    #=
    @gp([bestfit[ id][:OIII_5007].norm.val for id in 1:length(chosen_epochs)],
        [bestfit2[id][:OIII_5007].norm.val for id in 1:length(chosen_epochs)])
    @gp([bestfit2[id][:OIII_5007].norm.val for id in 1:length(chosen_epochs)],
        ([bestfit[ id][:OIII_5007].fwhm.val for id in 1:length(chosen_epochs)] .-
         [bestfit2[id][:OIII_5007].fwhm.val for id in 1:length(chosen_epochs)]) ./
        [bestfit[ id][:OIII_5007].fwhm.val for id in 1:length(chosen_epochs)])
    =#

    for id in 1:length(chosen_epochs)
        all_scale[chosen_epochs[id]] *= bestfit2[id][:OIII_5007].norm.val
    end
    serialize("output/scale_$(job)_$(loop).dat", all_scale)
end


(scale_evol, _) = deserialize("output/scale_fwhm.dat")
for loop in 1:Nloop
    ss = deserialize("output/scale_$(job)_$(loop).dat")
    scale_evol = hcat(scale_evol, ss)
end
scale_evol = 2 .*
    (scale_evol[:, 2:end] .- scale_evol[:, 1:end-1]) ./
    (scale_evol[:, 2:end] .+ scale_evol[:, 1:end-1])
@gp "set grid" :-
for epoch in chosen_epochs
    @gp :- 1:Nloop scale_evol[epoch, :] "w lp t '$epoch'"
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])
end


tabellone = DataFrame(epoch=Int[], date=String[], galaxy=Float64[],
                      contslope=Float64[], cont5100=Float64[], l5100=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[])

for id in 1:length(bestfit.preds)
    push!(tabellone,
          (id, epoch_filenames[id][27:end-4],
           bestfit[id][:galaxy].norm.val,
           bestfit[id][:qso_cont].alpha.val, cont_at(bestfit[id], 4000)[1], models[(id,:l5100)],
           bestfit[id][:bb_Ha].norm.patched, bestfit[id][:bb_Ha].fwhm.patched, bestfit[id][:bb_Ha].voff.patched,
           bestfit[id][:br_Ha].norm.patched, bestfit[id][:br_Ha].fwhm.patched, bestfit[id][:br_Ha].voff.patched,
           bestfit[id][:na_Ha].norm.patched, bestfit[id][:na_Ha].fwhm.patched, bestfit[id][:na_Ha].voff.patched,
           bestfit[id][:bb_Hb].norm.patched, bestfit[id][:bb_Hb].fwhm.patched, bestfit[id][:bb_Hb].voff.patched,
           bestfit[id][:br_Hb].norm.patched, bestfit[id][:br_Hb].fwhm.patched, bestfit[id][:br_Hb].voff.patched,
           bestfit[id][:na_Hb].norm.patched, bestfit[id][:na_Hb].fwhm.patched, bestfit[id][:na_Hb].voff.patched,
           bestfit[id][:OIII_5007].fwhm.patched, bestfit[id][:OIII_5007].voff.patched))
end

dd = Date.(tabellone.date, Ref("yyyymmdd"))
tabellone[!, :day] .= getproperty.(dd - dd[1], :value)

f = FITS("1ES_1927p654_results_$(job).fits", "w")
write(f, tabellone)
close(f)

@gp( [bestfit[id][:qso_cont].alpha.val for id in 1:length(chosen_epochs)],
     [bestfit[id][:br_Ha].norm.patched for id in 1:length(chosen_epochs)],
     "w l")

@gp(:-,
    [bestfit[id][:qso_cont].alpha.val for id in 1:length(chosen_epochs)],
    [bestfit[id][:br_Hb].norm.patched for id in 1:length(chosen_epochs)],
    "w l")


@gp( [bestfit[id][:qso_cont].alpha.val for id in 1:length(chosen_epochs)],
     [cont_at(bestfit[id], 5100)[1] for id in 1:length(chosen_epochs)],
     "w l")

@gp( [cont_at(bestfit[id], 5100)[1] for id in 1:length(chosen_epochs)],
     [bestfit[id][:br_Ha].norm.patched for id in 1:length(chosen_epochs)],
     "w l")


@gp( [models[(id, :l5100)] for id in 1:length(chosen_epochs)],
     [bestfit[id][:br_Ha].norm.patched for id in 1:length(chosen_epochs)],
     "w l")


@gp    tabellone.day [cont_at(bestfit[id], 5100)[1] for id in 1:length(chosen_epochs)] "w lp" ylog=true
@gp :- tabellone.day [bestfit[id][:br_Ha].norm.patched for id in 1:length(chosen_epochs)] .* 1e2 "w lp"
@gp :- tabellone.day [bestfit[id][:br_Hb].norm.patched for id in 1:length(chosen_epochs)] .* 1e2 "w lp"
@gp :- tabellone.day [models[(id, :l5100)] for id in 1:length(chosen_epochs)] "w lp"

@gp    tabellone.day -3 .+ 1e-3 .*[cont_at(bestfit[id], 5100)[1] for id in 1:length(chosen_epochs)] "w lp"
@gp :- tabellone.day [bestfit[id][:qso_cont].alpha.val for id in 1:length(chosen_epochs)] "w lp"


@gp "set grid"
for id in 1:length(chosen_epochs)
    @gp :- models[(id, :x)] models[(id, :y)] "w l t '$id' lw 3"
    #@gp :- models[(id, :x)] models[(id, :qso_cont)] "w l t '$id' lw 2"
end
save(:prenorm, term="png size 800,600", output="output/evolution.png")


@gp "set grid"
for id in 1:length(chosen_epochs)
    println(bestfit[id][:OIII_5007])
    x = models[(id, :x)]
    y  = models[(id, :y)]
    y0 = y .- models[(id, :OIII_5007)]
    #@gp :- x y0 .- Spline1D(x, y0)(5007.) "w l t '$id'"
    @gp :- x y  .- Spline1D(x, y0)(5007.) "w l t '$id'"
end
save(:prenorm, term="png size 800,600", output="output/evolution_oiii_norm.png")