%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  First, we'll load data and functions
%  and define some useful functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

require("xspec");  % Load the XSPEC module (for phabs)
plasma(aped);     % Load the APEC database

    % Load the LETG spectrum and its associated ARFs and RMFs:
    %
    % For this example, we'll include the first
    % 11 diffraction orders.
    %
    % The ARF and RMF filenames are numbered sequentially
    % so, just for fun, we generate the names rather than
    % listing them explicitly.

variable orders = [1:11];
variable arfs, rmfs, pha;
arfs = "arfs/leg_m" + array_map(String_Type, &string, orders) + ".arf";
rmfs = "rmfs/leg_m" + array_map(String_Type, &string, orders) + ".rmf";
pha  = "hrc_m1.pha";

    % This function simplifies loading N ARFs and RMFs.
    % It takes two arrays listing the filenames and returns
    % two arrays containing the index-labels of the
    % ARFs and RMFs that got loaded.

define read_rsps (arfs, rmfs)
{
   variable i, n = length(arfs);

   variable aid = Int_Type[n];
   for (i = 0; i < n; i++)
     {
	aid[i] = load_arf (arfs[i]);
     }

   variable rid = Int_Type[n];
   for (i = 0; i < n; i++)
     {
	rid[i] = load_rmf (rmfs[i]);
     }

   return (aid, rid);
}

define init_model ()
{
       % create a single-temperature APED model.
   variable plasma_state = default_plasma_state ();
   create_aped_fun ("xaped", plasma_state);

       % The source model:
   fit_fun ("phabs(1)*xaped(1)");

       % An instrumental background term:
   back_fun (1,"poly(1)");

       % initial parameter values
   set_par("phabs(1).nH", 0,1);           % absorbing column = 0
   set_par("xaped(1).norm", 0.015);
   set_par("xaped(1).temperature", 10.0^6.8);   % Kelvin

   set_par("poly(1).a0", 240.2);         % poly(1) constant term
   set_par("poly(1).a1", 0,1);           %         linear term
   set_par("poly(1).a2", 0,1);           %         quadratic term
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Here's where the action starts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % Load the PHA file, ARFs and RMFs
    % and assign all the responses to the
    % single LETG spectrum...

variable aid, rid, did;
did = load_data (pha);
(aid, rid) = read_rsps (arfs, rmfs);
assign_rsp (aid, rid, did);

    % Initialize the spectrum model.

init_model ();

    % ... and fold the model through all 11 dispersion orders.
    %     (fitting to the 3rd order peak)

variable xmin = 43;
variable xmax = 46;
xnotice(1, xmin, xmax);
() = fit_counts();
notice(1);

    % Zoom in on the 40-50 Angstrom region,
    % plot the data,
    % overplot the model (= sum of orders 1-11)

Label_By_Default =0;
xlabel (latex2pg("m\\lambda [\\A]"));
ylabel ("Counts/bin");

variable pid = open_plot("mo_aped_demo.ps/cps");
xrange (xmin, xmax);
ylog;
plot_data_counts(1);
oplot_model_counts(1);

    % Compute the 1st order contribution
    % and overplot that

assign_rsp (1,1,1);
() = eval_counts();
oplot_model_counts(1);

    % Compute the 3rd order contribution
    % and overplot it _using_the_1st_order_grid_:

assign_rsp (3,3,1);     % change the grid to 3rd order
() = eval_counts();     % compute 3rd order model
assign_rsp (1,1,1);     % change the grid back to 1st order
oplot_model_counts(1);

close_plot (pid);