%% Solution to "How to Find and Idenfity Weak Emission Lines in a Stellar HETG Spectrum"
%%  Problem set by D. Huenemoerder
%%
%% ISIS script by L. Corrales
%% 2015.08.12 - lia@space.mit.edu
%%

%%--------------------------------------------------------------------------------------------
%% 1. Use TGcat (http://tgcat.mit.edu) to find observations of Capella

%    Sorting by exposure time, I'll evaluate the longest HETG observation, ObsId 1103 (40 ks)
%    I selected ObsId 1103, then Actions -> Download
%    I selected pha2, arf, and rmf for download
%    After unpacking, I moved to folder obs_1103_tgid_3633/

%%--------------------------------------------------------------------------------------------
%% 2. Load the TGcat data into ISIS

% We can load the MEG and HEG +/- 1 orders by specifying what rows of the PHA file to read in
%   Rows 3 and 4 are HEG -1 and +1; Rows 9 and 10 are MEG -1 and +1

variable HEG1 = load_data("pha2.gz", [3,4]);
variable MEG1 = load_data("pha2.gz", [9,10]);

% Show what we have loaded so far
% Note the "id" column, which shows the integer value used to identify each particular dataset

list_data;

% Notice that the arf and rmf were not automatically loaded.
% They will have to be read in separately and assigned to the proper data set

() = load_arf("heg_-1.arf.gz");
() = load_arf("heg_1.arf.gz");
() = load_arf("meg_-1.arf.gz");
() = load_arf("meg_1.arf.gz");

% Examine the loaded arfs
list_arf;

% Now assign them to the datasets - be sure that the ids match properly!
assign_arf([1:4], [HEG1, MEG1]);
list_data;

% Now do the same for the rmfs

() = load_rmf("heg_-1.rmf.gz");
() = load_rmf("heg_1.rmf.gz");
() = load_rmf("meg_-1.rmf.gz");
() = load_rmf("meg_1.rmf.gz");

list_rmf;
assign_rmf([1:4], [HEG1, MEG1]);
list_data;

%%--------------------------------------------------------------------------------------------
%% 3. Make a plot of the combined MEG first orders in counts/bin vs wavelength from 1.7 to 5.5 Angs, 
%%    uniformly grouped by 4 bins, scaled linearly on both axes.

%% ISIS will plot according to certain environment settings
%%    Read help docs on 'plot_bin_integral', 'plot_bin_density', and 'plot_unit'

% My plot settings

plot_bin_integral;
plot_unit("Angstrom");
xrange(1.7, 5.5);
xlin; ylin;

plot_data_counts(MEG1[0]);     % MEG -1
oplot_data_counts(MEG1[1], 2); % MEG +1, note that "2" refers to the color choice

% See how noisy it is! That's why Dave asked us to group the spectrum
%    To find the correct grouping function, I ran 'apropos group' and read the help docs

group_data(MEG1, 4);

% Look at the data sets and notice how the MEG orders have a different number of bins now
list_data;

% Now replot
plot_data_counts(MEG1[0]);
oplot_data_counts(MEG1[1], 2);

%% If you want to plot the *combined* MEG+/-1 orders, 
%%    we recommend using the Remeis set of isis scripts, 
%%    which can be downloaded separately on the web.
%%    The following commands show how to use these scripts to make the plot.
%%    (Note that these plotting tools will sometimes change the plot environment settings.
%%    For this reason, I recommending setting up the plot environment before each plot command,
%%    if you are not using the contributed plotting tools.)

%% Load the scripts
() = evalfile("isisscripts.sl");

%% To access help docs on the plotting functions, you have to call one with no input
plot_counts();

fancy_plot_unit("Angstrom");  % set up plot units
variable pstyle = struct{dcol=1, decol=15, mcol=0, dsym=0};

plot_counts(MEG1[0], pstyle);                 % MEG-1
plot_counts(MEG1[1], pstyle; dcol=2, oplt=1); % MEG+1

% Separate combined orders plot
plot_counts(MEG1, pstyle; dcol=4);    % combined MEG first order

%%--------------------------------------------------------------------------------------------
%% 4. (Modified from Dave's assignment) - Make a plot of the combined HEG and MEG first orders
%%    in photon flux density (photons/cm^2/s/Angs), linearly scaled from 1.7-5.5 Angs, with a linear y-scale

% Group HEG data similar to MEG, but increase grouping by factor of 2
%    (which is the difference between energy resolution on the two gratings)
group_data(MEG1, 4);
group_data(HEG1, 8);

% Note the change
list_data;

% Now flux correct the data
% Read the help file on 'flux_corr'

flux_corr([MEG1, HEG1]);

% How to plot the flux using default ISIS plot functions:

plot_bin_density; % since we what per unit Angs
plot_unit("Angstrom");
xrange(1.7,5.5);
xlin; ylin;

plot_data_flux(MEG1[0], 2);  % MEG-1 (red)
oplot_data_flux(HEG1[0], 4); % HEG-1 (blue)

%% With the Remeis plotting tools
plot_unfold(MEG1, pstyle; dcol=2, yrange=[0,1.5e-3]); % first order MEG (red)
plot_unfold(HEG1, pstyle; dcol=4, oplt=1);            % first order HEG (blue)


%% If we want to *combine* HEG and MEG orders, the data array need to
%%    be of the same length and they need the same grid.
%%    We can force this by using the 'match_dataset_grids' function

match_dataset_grids([MEG1, HEG1]);

% Now when I bin all datasets by a factor of 4, they will have the same grid
group_data([MEG1, HEG1], 4);

% I verified this by plotting
plot_data_flux(MEG1[0], 2);  % MEG-1 (red)
oplot_data_flux(HEG1[0], 4); % HEG-1 (blue)

% Now combine everything!!
plot_unfold([MEG1, HEG1], pstyle; yrange=[0,1.5e-3]);

%%--------------------------------------------------------------------------------------------
%% 5. Identify some lines

% For this one, I'm going to dramatically expand the range of lines we will look at.
xrange(2, 20); % 2 - 20 Angs
plot_unfold([MEG1, HEG1], pstyle);

% Load AtomDB
plasma(aped);

% Define a single temperature plasma model 
% (see the help entry for 'default_plasma_state')
create_aped_fun( "aped_1T", default_plasma_state );

% Set the plasma model to be your fit model
fit_fun("aped_1T(1)");
set_par("aped_1T(1).temperature", 16e6 );
list_par;

% Note that they are all fixed. Thaw parameters of interest for fitting.
thaw(1, 2); % thawing norm and temperature

% Use xnotice to exclusively notice the wavelength range of interest
% The min and max values are always in angstroms
xnotice([MEG1, HEG1], 2, 20);

% note the change
list_data;

() = renorm_counts;
() = fit_counts;

%% Let's see what it looks like
plot_unfold( [MEG1, HEG1], pstyle; mcol=2 ); % model in red

%% All the major lines are there! We can figure out what they are by browsing the database

% Find the 10 brightest lines betwen 20 and 20 Angs
variable bls = brightest( 30, where(wl(2, 20)) );

page_group(bls);
plot_group(bls, 3);

% Let's say I'm only interested in the Ne and O VIII lines, we can sort on that
variable b_Ne = brightest(10, where( el_ion(Ne) and wl(2,20) ) );
plot_group(b_Ne, 4);

variable b_O = brightest(10, where( el_ion(O, 8) and wl(2,20) ) );
plot_group(b_O, 8);

%%--------------------------------------------------------------------------------------------
%% 6. Use the S XVI and S XV lines to get the plasma temperature
%%    This is the "Advanced" portion of Dave's example

% Define a new clean model
fit_fun("aped_1T(2)");
set_par("aped_1T(2).temperature", 16.e6);

% Set up a wavelength grid just for modelling
variable w1, w2, y;
(w1, w2) = linear_grid( 1.7, 5.5, 8192 );
y = eval_fun(w1, w2);

% Show the theoretical model with the high resolution grid
plot_bin_density;
xlabel("Wavelength (Angstroms)");
ylabel("Flux (phot/s/cm^2/Angs)");
hplot( w1, w2, y, 15 );

% Identify the brightest lines
plot_group( brightest(30, where( wl(1.5,5.5)) ), 4 ); % mark them in blue on plot
page_group( brightest(30, where( wl(1.5,5.5)) ) );    % see the info in table form

% Reading from the table, we can see that the S XVI Ly-alpha is a blended transition with
%   upper level index 3 and 4 dropping to lower level index 1
variable ls16 = where( trans(S, 16, [3:4], 1 ) );

% The S XV Ly-alpha line is also blended, but let's use the strongest one (5.04 Angs)
%    upper level index is 7 and lower level is 1
variable ls15 = where( trans(S, 15, 7, 1) );

%% Calculate the line ratios
variable T = 10^[6.0:8.0:0.1];
variable r = ratio_em( ls16, ls15, T );

% plot the line ratios as a function of temperature
xrange(); xlin; ylog;
plot( log10(T), r );