Newer
Older
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
## Pre-requisites
Install Julia (version > 1.5), and prepare the `AT2018zf` directory with spectral data.
## Installation
```
git clone https://www.ict.inaf.it/gitlab/giorgio.calderone/cl_1es_1927p654
cd cl_1es_1927p654
```
Start Julia, then:
```julia
using Pkg
Pkg.activate(".")
Pkg.instantiate()
```
## Usage
### Load packages and read spectrum file names (one for each epoch):
```julia
using Pkg
Pkg.activate(".")
#Pkg.update() # ensure code is up to date
using Revise
using QSFit, Gnuplot, GFitViewer
using CL_1ES_1927p654
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
```
### Plot spectrum at all epochs
([gnuplot](http://www.gnuplot.info/) is required to display the plot)
```julia
@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(epochs)
file = epochs[ii]
println(file)
l = readlines(file)
@gp :- l "u (\$1/1.017):2:($ii) w l notit lc pal"
#(readline() == "q") && break
end
```
### Run analysis on a specific epoch
```julia
chosen_epoch = 10
file = epochs[chosen_epoch]
source = QSO{q1927p654}("1ES 1927+654 ($chosen_epoch)", 0.019422, ebv=0.077);
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, rebin=2, selcomps=[:qso_cont, :galaxy, :balmer])
```
### Multiepoch analysis
```julia
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].val .*= 1e17;
source.data[id].unc .= 0.05 .* source.data[id].val;
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])])
```