_auto_declare=1;   % So we don't have to declare variables

require("xspec");

% Load the Cycle 10 ARFS: HEG-1, HEG+1, MEG-1, MEG+1

hma = load_arf("aciss_heg-1_cy10.garf");
hpa = load_arf("aciss_heg1_cy10.garf");
mma = load_arf("aciss_meg-1_cy10.garf");
mpa = load_arf("aciss_meg-1_cy10.garf");

% Load the Cycle 10 RMFS: HEG-1, HEG+1, MEG-1, MEG+1

hmr = load_rmf("aciss_heg-1_cy10.grmf");
hpr = load_rmf("aciss_heg1_cy10.grmf");
mmr = load_rmf("aciss_meg-1_cy10.grmf");
mpr = load_rmf("aciss_meg1_cy10.grmf");

% Assign ARFs and RMFs to Dataset IDs

assign_arf(hma,1);  % Dataset 1 = HEG-1
assign_arf(hpa,2);  % Dataset 2 = HEG+1
assign_arf(mma,3);  % Dataset 3 = MEG-1
assign_arf(mpa,4);  % Dataset 4 = MEG+1

assign_rmf(hmr,1);  % Dataset 1 = HEG-1
assign_rmf(hpr,2);  % Dataset 2 = HEG+1
assign_rmf(mmr,3);  % Dataset 3 = MEG-1
assign_rmf(mpr,4);  % Dataset 4 = MEG+1

% Set all ARF exposures to 10x desired value

expos = 65.e3;

set_arf_exposure([hma,hpa,mma,mpa],10*expos);
set_data_exposure(1,10*expos);
set_data_exposure(2,10*expos);
set_data_exposure(3,10*expos);
set_data_exposure(4,10*expos);

% Define the fit function

fit_fun("phabs(1)*(bbodyrad(1)+powerlaw(1)+gaussian(1)-gaussian(2)-gaussian(3))");

% Set parameter values.  Start with no absorption, and very weak lines

set_par("phabs(1).nH",0.);
set_par("bbodyrad(1).norm",41.);
set_par("bbodyrad(1).kT",1.18);
set_par("powerlaw(1).norm",0.15);
set_par("powerlaw(1).PhoIndex",2.02);
set_par("gaussian(1).norm",3.e-5);
set_par("gaussian(1).LineE",6.4);
set_par("gaussian(1).Sigma",0.1);
set_par("gaussian(2).norm",3.e-5);
set_par("gaussian(2).LineE",6.72);
set_par("gaussian(2).Sigma",0.02);
set_par("gaussian(3).norm",3.e-5);
set_par("gaussian(3).LineE",7.00);
set_par("gaussian(3).Sigma",0.02);

% Create the fake data

fakeit;
list_data;

% Define a flux and equivalent width function

public define flux (id, e_lo, e_hi)
{
   variable l_lo, l_hi, m, e_avg, pflux, eflux;

   l_hi = _A(e_lo);        % Convert keV to A, or visa versa
   l_lo = _A(e_hi);        % Convert keV to A, or visa versa
   m = get_model_flux(id);
   pflux = rebin (l_lo, l_hi, m.bin_lo, m.bin_hi, m.value)[0];   
   e_avg = _A(1.)*(1./m.bin_hi+1./m.bin_lo)/2.*1.602e-9;
   eflux = rebin (l_lo, l_hi, m.bin_lo, m.bin_hi, m.value*e_avg)[0];   

   () = printf("\n photons/cm2/sec: %11.3e, ergs/cm2/sec: %11.3e \n \n", pflux, eflux);

   return pflux, eflux;  % Photon flux & ergs flux (per cm^2/sec) 
}

public define eqw(indx,par)
{
   variable eva, eve;
   variable fv_hold = Fit_Verbose, pmin=0;

   Fit_Verbose = -1;
   () = eval_counts;

   variable a = get_model_flux(indx);
   variable phold = get_par_info(par);

   set_par(par,pmin,0,0,phold.max);
   () = eval_counts;
   variable b = get_model_flux(indx);

   set_par(par,phold.value,phold.freeze,phold.min,phold.max);
   () = eval_counts;

   variable iw = [0:length(a.bin_lo)-1];

   variable diff = a.value - b.value;
   variable bdena = b.value/(a.bin_hi-a.bin_lo);
   variable bdene = reverse(reverse(b.value)/(_A(a.bin_lo)-_A(a.bin_hi)));
   variable fdena = diff/(a.bin_hi-a.bin_lo);
   variable fdene = reverse(reverse(diff)/(_A(a.bin_lo) - _A(a.bin_hi)));
   variable flux = sum( diff[iw] );

   variable iwa = min(where( abs(fdena[iw]) == max(abs(fdena[iw])) ));
   variable iwe = min(where( abs(fdene[iw]) == max(abs(fdene[iw])) ));
   
   eva = flux/bdena[iw[iwa]]*1000.;
   eve = flux/bdene[iw[iwa]]*1000.;

   Fit_Verbose = fv_hold;
 
   () = printf("\n EqW (mA): %11.3e, EqW (eV): %11.3e \n \n", eva, eve);

   return eva, eve;
}

% Get the current, too small, equivalent widths for the Fe lines

(pflux,eflux) = flux(3,0.1,8);  % MEG -1 flux, from 0.1 to 8 keV
(g1_a,g1_e) = eqw(1,"gaussian(1).norm"); % HEG -1 Equivalent width
(g2_a,g2_e) = eqw(1,"gaussian(2).norm"); % HEG -1 Equivalent width
(g3_a,g3_e) = eqw(1,"gaussian(3).norm"); % HEG -1 Equivalent width

% Redefine the line normalizations to get the correct equivalent widths

set_par("gaussian(1).norm",get_par("gaussian(1).norm")*95/g1_e);
set_par("gaussian(2).norm",get_par("gaussian(2).norm")*(-14)/g2_e);
set_par("gaussian(3).norm",get_par("gaussian(3).norm")*(-23)/g3_e);

% Fake the data again, and check the new equivalent widths

fakeit;
(g1_a,g1_e) = eqw(1,"gaussian(1).norm"); % HEG -1 Equivalent width
(g2_a,g2_e) = eqw(1,"gaussian(2).norm"); % HEG -1 Equivalent width
(g3_a,g3_e) = eqw(1,"gaussian(3).norm"); % HEG -1 Equivalent width

% Set the neutral column to a realistic value, and save the parameters

set_par("phabs(1).nH",7.5);
save_par("lmxb.par");

% Set the exposure to the desired value

set_arf_exposure([hma,hpa,mma,mpa],expos);
set_data_exposure(1,expos);
set_data_exposure(2,expos);
set_data_exposure(3,expos);
set_data_exposure(4,expos);

% Fake the data a final time, and look at it

fakeit;
list_data;

% Start fitting the data - first HEG lines, then MEG continuum

group_data([1,2],2);         % Group HEG data by a factor of 2
group_data([3,4],4);         % Group MEG data by a factor of 4
heg = combine_datasets(1,2); % In memory combination of HEG data
meg = combine_datasets(3,4); % In memory combination of MEG data

xnotice([1,2],1.65,2.3); % Only notice the 1.65-2.3 A data for HEG
exclude(3,4);            % Ignore MEG data for now

freeze([1,3,5]);         % Freeze Nh, blackbody kT, Photon Index

() = fit_counts;
list_par;

save_par("heg_fit.par");  % Save the HEG fit parameters

putenv("PGPLOT_BACKGROUND=white");
putenv("PGPLOT_FOREGROUND=black");  % Make plots black on white

set_line_width(3); set_frame_line_width(3); charsize(1.1); % Nicer plot defaults
xrange(1.65,2.3);  % Plot x-range to only show what we've noticed
ylog;              % Logarithmic y-axis

rplot_counts(1);  % HEG-1 data (although fit includes HEG+1);
rplot_counts(2);  % HEG+1 data (although fit includes HEG-1);

hd = get_combined(heg, &get_data_counts);  % Get the summed data
hm = get_combined(heg, &get_model_counts); % Get the summed model

label("Angstrom","Counts/Bin","Combined HEG Data");
hplot(hd.bin_lo,hd.bin_hi,hd.value);
ohplot(hm.bin_lo,hm.bin_hi,hm.value);

% Get confidence values.  For the script, the values must be captured in 
% variables.  One might want to use different names for different
% parameters.  (On the command line, the return values are "popped off
% the stack" and printed to the screen.)

(lo,hi) = conf("gaussian(1).norm");
(lo,hi) = conf("gaussian(1).LineE");
(lo,hi) = conf("gaussian(1).Sigma");
(lo,hi) = conf("gaussian(2).norm");
(lo,hi) = conf("gaussian(2).LineE");
(lo,hi) = conf("gaussian(2).Sigma");
set_fit_method("minim");  % Slower fit method, but sometimes more accurate
 (lo,hi) = conf("gaussian(2).Sigma");
set_fit_method("lmdif"); % The ISIS default fit method
(lo,hi) = conf("gaussian(3).norm");
(lo,hi) = conf("gaussian(3).LineE");
(lo,hi) = conf("gaussian(3).Sigma");


% Get equivalent widths

set_fit_method("subplex");  % Slow fit method, but often useful for gratings data

(g1_a,g1_e) = eqw(1,"gaussian(1).norm");
set_par("gaussian(1).norm",0.000814,1);
() = fit_counts();
(g1_a,g1_e) = eqw(1,"gaussian(1).norm");
set_par("gaussian(1).norm",0.00112,1);
() = fit_counts();
(g1_a,g1_e) = eqw(1,"gaussian(1).norm");

load_par("heg_fit.par");  % Restore the HEG fit parameters

(g2_a,g2_e) = eqw(1,"gaussian(2).norm");
set_par("gaussian(2).norm",0.000135569,1);
() = fit_counts();
(g2_a,g2_e) = eqw(1,"gaussian(2).norm");
set_par("gaussian(2).norm",0.000306819,1);
() = fit_counts();
(g2_a,g2_e) = eqw(1,"gaussian(2).norm");

load_par("heg_fit.par");  % Restore the HEG fit parameters

(g3_a,g3_e) = eqw(1,"gaussian(3).norm");
set_par("gaussian(2).norm",0.000187434,1);
() = fit_counts();
(g3_a,g3_e) = eqw(1,"gaussian(3).norm");
set_par("gaussian(3).norm",0.000352399,1);
() = fit_counts();
(g3_a,g3_e) = eqw(1,"gaussian(3).norm");

%  Now start fitting the MEG data

exclude([1,2]);         % Ignore HEG
include(3,4);           % Notice MEG
xnotice([3,4],2.1,7.2); % Only notice 1.7-7.2 Angstrom (1.7-5.9 keV)

load_par("lmxb.par");  % Simulation parameters, including lines
() = eval_counts;

set_fit_method("lmdif");                       % Fast fit method
fit_fun("phabs(1)*(bbodyrad(1)+powerlaw(1))"); % Don't fit lines
() = fit_counts;

list_par;
save_par("meg_fit.par");  % Save the MEG fit parameters

% Get the confidence limits on the neutral column

(lo,hi) = conf("phabs(1).nH");

% Plot the MEG data

title("");       % We set the plot title in "label" above, so remove it
xrange(2.1,7.2); % Only plot 2.1-7.2 Angstroms
rplot_counts(3); % MEG-1 data (although fit combined with MEG+1)
rplot_counts(4); % MEG+1 data (although fit combined with MEG-1)

flux_corr(3);      % Calculated the "unfolded" spectrum for MEG-1
plot_data_flux(3); % Plot the MEG-1 unfolded spectrum 

ylin;                % Linear y-axis
plot_bin_density;    % Create plots as counts/sec/Angstrom
plot_data_counts(4); % MEG+1 is the brightest arm
cursor;              % Measure with the pgplot cursor