%; Time-stamp: <2001-10-24 17:40:50 dph> 
%; MIT Directory: ~dph/libisis
%; File: Aped_model_demo.sl
%; Author: D. Huenemoerder
%; Original version: 2001.10.24
%;====================================================================
% version: 1
%
% purpose: Demonstrate use of aped plasma models in isis; test aped
%	   plasma user models
%
%

% Where are the functions?
%
autoload("aped_multi_fit", "aped_fit_models.sl");
autoload("aped_multi_init", "aped_fit_models.sl");
autoload("aped_multi_setup", "aped_fit_models.sl");


% We need the plasma database loaded:
%
plasma(aped);

% A couple window, one w/ 2 panels for residual plots:
%
p2 = open_plot("/xwin"); resize(20,0.6); multiplot([3,1]);
p1 = open_plot("/xwin"); resize(20,0.6);

fpha="./Pha/acisf01451_N004_pha2.fits.gz";   % II Peg; Standard Type II Pha file.

hm=3; hp=4; mm=9; mp=10;              % define HEG, MEG order row variables.

h=load_data(fpha,[hm, hp, mm, mp]);   % load HEG, MEG +-1st order from
	 		              % standard product

list_data;			      % show the histogram info.


% Define ARF files string arrays:
%
fgarfs_heg = [
  "Arf/acisf01451_N005_HEG_-1_garf.fits.gz",
  "Arf/acisf01451_N005_HEG_1_garf.fits.gz"
  ];

fgarfs_meg = [
  "Arf/acisf01451_N005_MEG_-1_garf.fits.gz",
  "Arf/acisf01451_N005_MEG_1_garf.fits.gz"
  ];


% Define and load RMF files:  (also in CALDB)
%
frmf_heg = "Ard/acisheg1D1999-07-22rmfN0002.fits";
frmf_meg = "Ard/acismeg1D1999-07-22rmfN0002.fits";

% Load rmfs, define indices, define number of orders
%

load_rmf(frmf_heg); rmf_heg=1;  nheg = length(fgarfs_heg);
load_rmf(frmf_meg); rmf_meg=2;  nmeg = length(fgarfs_meg);


% Load ARFS, assign ARF, RMF to histograms:
%

for (i=0; i < length(fgarfs_heg); i++)  % loop is overkill here, but
				        % demonstrates how to generalize
				        % to many files. 
{
  ()=load_arf(fgarfs_heg[i]);
  assign_arf(i+1,i+1);
  assign_rmf(rmf_heg, i+1);
}

for (i=0; i < length(fgarfs_meg); i++)
{
  ()=load_arf(fgarfs_meg[i]);
  assign_arf(nheg + i+1, nheg + i+1);
  assign_rmf( rmf_meg, nheg + i+1);
}

list_data;			% list histogram info to show assignments

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set up a multi-temperature model.  3 components for demo:

aped_multi_init(3);                                 % defines the model function.
aped_multi_setup(10.^[6.8,7.2,7.6], [1,2,4]*0.01);  % sets model, init params.

elem_list =  [Mg,   Fe,   Si,   S,    O,    N,   Ar,  Ne];   % tweak abundances
abund_list = [0.5, 0.15, 0.45, 0.65, 1.10, 0.6, 1.0, 2.35];
aped_set_abund(elem_list, abund_list);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

ignore([1,2,4]);         % only MEG -1 for starters.
xnotice(3, 1.7, 25);
xrange(1.7,25);
plot_data_counts(3, 2);  
eval_counts;             % see if in the right ballpark...........
list_model;

oplot_model_counts(3,3); 

% pick a small region and fit...short wave region for hot continuum...
%

% group a bit.........
group_data(1,4); group_data(2,4); group_data(3,2); group_data(4,2);

xmin=1.7; xmax=8;       % fit hot continuum....
xnotice([1,2,3,4], xmin, xmax);

yrange(0);xrange(xmin, xmax);  plot_data_counts(3, 2);  
eval_counts; oplot_model_counts(3,3); 
oplot_data_counts(1, 5);oplot_model_counts(1,6); 
plot_group(brightest(20, where(wl(xmin, xmax))), 5); % brightest() is for current model

fit_counts; list_par;

%  Parameters[Variable] = 6[6]
%             Data bins = 2520
%   Best fit chi-square = 1165
%    Reduced chi-square = 0.4633
% aped_multi(1)
%  idx  param                    tie-to  freeze         value          min          max
%   1  aped_multi[1].norm_1          0     0        0.02830747            0            1
%   2  aped_multi[1].norm_2          0     0       0.002516543            0            2
%   3  aped_multi[1].norm_3          0     0        0.04519434            0            4
%   4  aped_multi[1].Temperature_1   0     0           5409037        10000        3e+08
%   5  aped_multi[1].Temperature_2   0     0      3.478148e+07        10000        3e+08
%   6  aped_multi[1].Temperature_3   0     0      3.298896e+07        10000        3e+08

% "Hard"copy of fit:

pid = dup_plot("Aped_model_demo_1.gif/gif",p1); resize(6*2.54, 0.6);

Object=" II Peg (ObsID 1451) ";
label("Wavelength", "Flux", Object + " Fit to region");
plot_data_counts(3, 2);  oplot_model_counts(3,3); 
oplot_data_counts(1, 6);oplot_model_counts(1,5); 
plot_group(brightest(30, where(wl(xmin, xmax))), 5); % brightest() is for current model

close_plot(pid); window(p1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% residuals plot:
%

window(p2); yrange; xrange(xmin, xmax); rplot_counts(3);
plot_group(brightest(30, where(wl(xmin, xmax))), 5); % brightest() is for current model
window(p1); xrange(xmin, xmax);


%%% "Refine" fit, but freezing hot component, fitting another region:

xmin=10.;  xmax=12;
xnotice([1:4], xmin, xmax);
freeze(3);  freeze(6);
fit_counts; list_par;

%  Parameters[Variable] = 6[4]
%             Data bins = 800
%   Best fit chi-square = 367.3
%    Reduced chi-square = 0.4615
%
%  aped_multi(1)
%  idx  param                    tie-to  freeze         value          min          max
%   1  aped_multi[1].norm_1          0     0        0.01826794            0            1
%   2  aped_multi[1].norm_2          0     0        0.01065223            0            2
%   3  aped_multi[1].norm_3          0     1        0.04519434            0            4
%   4  aped_multi[1].Temperature_1   0     0           2265638        10000        3e+08
%   5  aped_multi[1].Temperature_2   0     0      1.117411e+07        10000        3e+08
%   6  aped_multi[1].Temperature_3   0     1      3.298896e+07        10000        3e+08


% evaluate current model at full resolution over the whole spectral
% range:
%
for (i=1; i<=4; i++) group_data(i,0);          % regroup
xmin=1.7; xmax=25;                             % notice region
xnotice([1,2], 1.7, 17); xnotice([3,4], 1.7, 25); % notice bins
eval_counts;                                      % compute model

%%%% Plot a region:
%
label("Wavelength", "Flux", Object + " model and counts ");
yrange(-4,120);
xmm=[7,9]+3;xrange(xmm[0], xmm[1]);
plot_data_counts(3, 2);  oplot_model_counts(3,3);
plot_group(brightest(40,where(wl(xmm[0],xmm[1]))),5);
oplot_data_counts(1, 6);  oplot_model_counts(1,5);

%%%%%%%%%%%%%%%%%
% Filecopy:
%
pid = dup_plot("Aped_model_demo_2.gif/gif",p1); resize(6*2.54, 0.6);
label("Wavelength", "Flux", Object + " model and counts ");
xmm=[7,9]+3;xrange(xmm[0], xmm[1]);
yrange(-4, 120);
plot_data_counts(3, 2);  oplot_model_counts(3,3);
plot_group(brightest(40,where(wl(xmm[0],xmm[1]))),5);
oplot_data_counts(1, 6);  oplot_model_counts(1,5);
close_plot(pid); window(p1);
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%% simulate a spectrum:
%
ignore([1:4]);
assign_arf(1, 5, 1); assign_arf(3,6,1);   
assign_rmf(1,5);     assign_rmf(2,6);
xnotice(5, 1.7,17); xnotice(6,1.7,25);
fakeit;

xnotice([1,2,5],1.7,17); xnotice([3,4,6],1.7,25);

label("Wavelength [Angstroms]", "Flux", Object+"SIMULATION and model");
xmm=[7,9]+3;
xrange(xmm[0], xmm[1]); yrange(-4,120);
plot_data_counts(6, 2);  oplot_model_counts(6,3);
plot_group(brightest(40,where(wl(xmm[0],xmm[1]))),5);
oplot_data_counts(5, 6);  oplot_model_counts(5,5);


%%%%%%%%%%% repeat for filecopy:
%
pid = dup_plot("Aped_model_demo_3.gif/gif",p1); resize(6*2.54, 0.6);
label("Wavelength [Angstroms]", "Flux", Object+"SIMULATION and model");
xmm=[7,9]+3;
xrange(xmm[0], xmm[1]);  yrange(-4,120);
plot_data_counts(6, 2);  oplot_model_counts(6,3);
plot_group(brightest(40,where(wl(xmm[0],xmm[1]))),5);
oplot_data_counts(5, 6);  oplot_model_counts(5,5);
close_plot(pid); window(p1);
%