using Pkg
Pkg.activate(".")
Pkg.instantiate()


using Revise
using QSFit, Gnuplot, GFitViewer
using CL_1ES_1927p654

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




# ====================================================================
# Plot spectrum at all epochs
@gp ylog=true xr=[3e3,1e4] yr=[1e-16, 1e-13] cbr=[1,26] :-
@gp :- cblabel="Epoch" xlab="Rest frame wavelength[A]" ylab="Flux density"
for ii in 1:length(datafiles)
    file = datafiles[ii]
    println(file)
    l = readlines(file)
    @gp :- l "u (\$1/1.017):2:($ii) w l notit lc pal"
    (readline() == "q")  &&  break
end


for ii in 1:length(datafiles)
    file = datafiles[ii]
    source = QSO{q1927p654}("1ES 1927+654 ($ii)", 0.019422, ebv=0.077);
    add_spec!(source, Spectrum(Val(:ASCII), file, columns=[1,2]));
    source.data[1].unc .= 0.05 .* source.data[1].val;
    source.data[1].val .*= 1e17;
    source.data[1].unc .*= 1e17;
    (model, bestfit) = fit(source);
    viewer(model, source.data, bestfit, rebin=2, selcomps=[:qso_cont, :galaxy, :balmer])
    (readline() == "q")  &&  break
end

#=
Need a base below both Ha and Hb
Added BOH line at 4490
Avoid using unknown lines

Blue component clearly visible: 20180309.dat
non e` la galassia, non e` l'AGN (finisce a 4300) e non e` il ferro
La powerlaw e` a -4.5!

Good for multiepoch:
20180423: powerlaw a -3, la blue component non c'e` piu`
20180424: powerlaw a -2
20180507: powerlaw a -3.1
20180528: powerlaw a -3.2
20180603: powerlaw a -3.2

Starting from ~ 201806 the host is no longer used

20180706: Hg very bright

20180717: very bright lines and blue component

20180811: got back host galaxy, Hb asymmetric on the other side

20181113: host galaxy vanished again

20190319: powerlaw at -1.5

20190519: Ha very asymmetric
Try multiepochs
20180811
20180812
=#



source = QSO{q1927p654}("1ES 1927+654", 0.019422, ebv=0.077);
add_spec!(source, Spectrum(Val(:ASCII), "AT2018zf/AT2018zf_optspec_20180811.dat", columns=[1,2]));
add_spec!(source, Spectrum(Val(:ASCII), "AT2018zf/AT2018zf_optspec_20180812.dat", columns=[1,2]));
for id in 1:length(source.domain)
    source.data[id].unc .= 0.05 .* source.data[id].val;
    source.data[id].val .*= 1e17;
    source.data[id].unc .*= 1e17;
end
(model, bestfit) = multiepoch_fit(source);
viewer(model, source.data, bestfit, rebin=2,
       selcomps=[Symbol.(:T1_, [:qso_cont, :galaxy, :balmer])
                 Symbol.(:T2_, [:qso_cont, :galaxy, :balmer])])


@gp    :t1 domain(model, id=1) source.data[1].val "w p"
@gp :- :t1 domain(model, id=1) model(id=1) "w l "

@gp    :t2 domain(model, id=2) source.data[2].val "w p"
@gp :- :t2 domain(model, id=2) model(id=2) "w l"


f = open("data.txt", "w")
for i in 1:length(source.domain[1])
    println(f, source.domain[1][1][i], "  ", source.data[1].val[i], "  ", source.data[1].unc[i])
end
close(f)

#=
cd, '/home/gcalderone/my/work/software/qsfit'
compile
cd, '/home/gcalderone/my/work/2020/ChangingLook/CL_1ES_1927p654'

input = qsfit_input('data.txt', type='ASCII', 0.019422, ebv=0.077)
res = qsfit(input)
qsfit_plot, res
=#