% Here's a list of the data files we'll use,
    % 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.

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

    % 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);
}

    % Compute the contribution from each order
    % and over-plot

define oplot_orders (m)
{
   assign_rsp (m, m, 1);
   () = eval_counts;
   assign_rsp (1, 1, 1);
   oplot_model_counts(1);
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Here's where the action starts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
    % Load the ARFs and RMFs

variable aid, rid;
(aid, rid) = read_rsps (arfs, rmfs);
    
    % Generate a matching, empty data set:

assign_rsp (aid, rid, 1);

    % We'll fake a powerlaw spectrum 

fit_fun ("phabs(1) * powerlaw(1)");
set_par (1, 0.01, 0, 0.002, 0.05);
set_par (3, 1.5, 0, 0, 3);
fakeit;
list_par;

    % Plot the fake data:

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

ylog;
variable pid = open_plot("mo_plaw_demo.ps/cps");
plot_data_counts(1);

    % Now, randomize the parameters and 
    % do a fit -- we should recover the
    % known parameter values.

randomize;
() = fit_counts;
oplot_model_counts(1);
list_par;

    % Compute just the 1st order contribution
    % and over-plot that.

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

close_plot(pid);

    % Compare contributions from orders 1, 2, 3
    % with the sum of higher orders 4-11.

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

pid = open_plot("mo_plaw2_demo.ps/cps");
ylog;
plot_data_counts(1);

oplot_orders (1);
oplot_orders (2);
oplot_orders (3);
oplot_orders ([4:11]);

close_plot(pid);