Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
using Pkg
Pkg.activate(".")
using Revise, Statistics, Serialization
using QSFit, GFit, Gnuplot, GFitViewer
using CL_1ES_1927p654
function fitline(x, y, err=[]; center=NaN, sigma=NaN)
domain = Domain(x)
isfinite(center) || (center = mean(x))
isfinite(sigma) || (sigma = (maximum(x) - minimum(x)) / 10.)
(length(err) == 0) && (err = (y .* 0.1) .+ maximum(y) / 10)
model = Model(domain, Reducer(sum, [:continuum, :line]),
:continuum => GFit.OffsetSlope(0., mean(x), 0.),
:line => GFit.Gaussian(1., center, sigma))
bestfit = fit!(model, Measures(y, err))
return (model, bestfit)
end
epochs = Vector{String}()
for (root, dirs, files) in walkdir("AT2018zf")
for file in files
if file[end-3:end] == ".dat"
push!(epochs, root * "/" * file)
end
end
end
@gp :summary xr=[3e3,1e4] cbr=[1,26] :-
@gp :- :summary cblabel="Epoch" xlab="Rest frame wavelength[A]" ylab="Flux density"
for ii in 1:length(epochs)
file = epochs[ii]
println(file)
s = Spectrum(Val(:ASCII), file, columns=[1,2]);
x = s.λ;
y = s.flux;
y .-= median(y)
y ./= quantile(y, 0.95)
@gp :- :summary x y "u 1:(\$2+$ii):($ii) w l notit lc pal"
end
norm = Vector{Float64}()
fwhm = Vector{Float64}()
for i in 1:length(epochs)
s = Spectrum(Val(:ASCII), epochs[i], columns=[1,2]);
x = s.λ;
j = findall(5060 .< x .< 5150);
x = x[j];
y = s.flux[j];
y .-= median(y);
y ./= quantile(y, 0.95);
j = findall(isfinite.(x));
x = x[j];
y = y[j];
model, bestfit = fitline(x, y)
show(bestfit)
push!(norm, bestfit[:line].norm.val)
push!(fwhm, bestfit[:line].sigma.val / 5007 * 3e5 * 2.355)
@gp tit=string(i) x y "w lp notit" x model() "w l notit"
(readline() == "q") && break
end
@gp :pre "set key bottom right" "set ytics nomirror" "set grid" xlab="Epoch" ylab="[OIII] norm." "set y2label '[OIII] FWHM'" norm "u (\$0+1):1 w lp t '[OIII] norm.'" fwhm "u (\$0+1):1 axes x1y2 w lp t '[OIII] FWHM"
chosen_epoch = 17
file = epochs[chosen_epoch]
source = QSO{q1927p654}("1ES 1927+654 ($chosen_epoch)", 0.019422, ebv=0.077);
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"
add_spec!(source, Spectrum(Val(:ASCII), file, columns=[1,2]));
source.data[1].val .*= 1e17;
source.data[1].unc .= 0.05 .* source.data[1].val;
(model, bestfit) = fit(source);
viewer(model, source.data, bestfit, selcomps=[:qso_cont, :galaxy, :balmer])
chosen_epochs = [17, 10, 21, 23] # vecchi
chosen_epochs = [8,9,10,11,12,14,15,16,18,21,22,23,24,25] # low-res
chosen_epochs = [5, 7, 13, 17, 20] # high-res
source = QSO{q1927p654}("1ES 1927+654", 0.019422, ebv=0.077);
source.options[:min_spectral_coverage][:OIII_5007] = 0.5
for id in 1:length(chosen_epochs)
file = epochs[chosen_epochs[id]]
add_spec!(source, Spectrum(Val(:ASCII), file, columns=[1,2]));
source.data[id].val .*= 1e17;
source.data[id].unc .= 0.05 .* source.data[id].val;
end
(model, bestfit) = multi_fit(source);
viewer(model, source.data, bestfit, selcomps=[:qso_cont, :galaxy, :balmer])
serialize("results_highres.dat", (source, model, bestfit))