Commit 1c244280 authored by Giorgio Calderone's avatar Giorgio Calderone
Browse files

Updated

parent 9ca6ed7b

run.jl

+51 −27
Original line number Diff line number Diff line
@@ -8,7 +8,7 @@ using Revise, Statistics, Serialization, Dierckx, DataFrames 
using QSFit, GFit, Gnuplot, GFitViewer, MyAstroUtils

using CL_1ES_1927p654

using Dates



using TextParse



mkdir("output")

epoch_filenames = Vector{String}()
@@ -220,7 +220,7 @@ function cont_at(bestfit, λ) 
end





tab = DataFrame(epoch=Int[], date=String[], galaxy=Float64[],

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[],
@@ -233,7 +233,7 @@ tab = DataFrame(epoch=Int[], date=String[], galaxy=Float64[], 


for id in 1:length(res.bestfit.preds)

    push!(tab,

          (id, epoch_filenames[id][27:end-4],

          (id, epoch_filenames[chosen_epochs[id]][27:end-4], "",

           res.bestfit[id][:galaxy].norm.val,

           res.bestfit[id][:qso_cont].alpha.val,

           cont_at(res.bestfit[id], 3500)[1],
@@ -253,47 +253,71 @@ 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 - dd[1], :value)

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("1ES_1927p654_results_$(job).fits", "w")

write(f, tab)

close(f)



@gp    tab.contslope tab.Ha "w l t 'Ha'"  ylog=true

@gp :- tab.contslope tab.Hb "w l t 'Hb'"

@gp :- tab.contslope tab.c5100 "w l t 'c5100'"

@gp :- tab.contslope tab.l5100 "w l t 'l5100'"

@gp :- tab.contslope tab.c3500 "w lp t 'c3500'"

@gp :- tab.contslope tab.l3500 "w lp t 'l3500'"



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.c5100 tab.Ha "w l t 'c5100'"

@gp :- tab.l5100 tab.Ha "w l t 'l5100'"



@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    tab.day tab.c5100 "w lp t 'c5100'" ylog=true

@gp :- tab.day tab.l5100 "w lp t 'l5100'"

@gp :- tab.day tab.Ha    "w lp t 'Ha'"

@gp :- tab.day tab.Hb    "w lp t 'Hb'"



@gp :- tab.day tab.contslope "w lp"

@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"





@gp "set grid"

for id in 1:length(chosen_epochs)

    #@gp :- domain(res.model[id])[:] res.model[id]() "w l t '$id' lw 3"

    @gp :- domain(res.model[id])[:] res.model[id](:qso_cont) "w l t '$id' lw 3"

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

end

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





@gp "set grid"

for id in 1:length(chosen_epochs)

    println(res.bestfit[id][:OIII_5007])

    x  = domain(res.model[id])[:]

    y  = res.model[id]()

    y0 = y .- res.model[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'"

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

end

save(:prenorm, term="png size 800,600", output="output/evolution_oiii_norm.png")