public define keV_note(a,b,c)
     {
     variable bi, ci ;
     variable conv=12.398424 ;
     if ( b != 0.) bi = conv/b; else return 0;
     if ( c !=0. ) ci = conv/c; else return 0;
     xnotice (a,ci,bi);
     return 1;
     }
%
% Raised on ftools and XSPEC, the following acts like a
% "group min minv" in grppha, starting the grouping at the
% lowest energy bin, and working towards higher energy (as opposed
% to doing it low to high in wavelength space)
%

public define grppha(a,min)
     {
      () = rebin_all(a,0);
      variable s = get_data_counts(a);
      variable len = length(s.value);
      variable idx = [0:len-1:1];
      idx[*] = 1;
      variable sn = 1;
      variable ssum = 0.;
      variable i;
      for (i=len-1; i >=0; i--)
       {
        idx[i] = sn*idx[i];
        ssum += s.value[i];
        if (ssum >= min)
           {
            sn = sn*(-1);
            ssum = 0.;
           }
        }
       rebin_data(a,idx);
       return 1;
      }

%
% Load the data file in
%

() = load_data("pilesim.pi");

%
% Group to a minimum of 20 counts per bin, and only notice ~0.5-3 keV
%

() = grppha(1,20);
() = keV_note(1,0.52,2.8);

%
% Set the pileup kernel, import XSPEC functions, and define the model
%

set_kernel(1,"pileup");
import("xspec");
fit_fun("phabs(1)*(bbody(1))");

%
% Normally, one can now define the starting values for the fits
% by either reading in a file, or editing them in a (vi, emacs, or
% choose your favorite editor) window.  Here, however, I'll define
% them one parameter at a time. Syntax is: parameter value,
% freeze (1) or thaw (0), min value, max value
%

% N_H
set_par(1,0.2,0,0.05,0.4);

% Bbody Norm
set_par(2,2.e-6,0,1.e-7,1.e-5);

% Bbody Temp
set_par(3,0.2,0,0.05,0.4);

% Pileup parameter alpha, which we freeze
set_par(6,1.,1);

%
% Next we define a function to loop through the values of alpha,
% fit the model, and then write the results to a file
%

public define alphasrch()
   {
   variable info=struct {statistic, num_variable_params, num_bins};
   variable alpha=[1.0:-0.01:-0.05];
   variable xi=@alpha;
   variable nh=@alpha;
   variable bbnorm=@alpha;
   variable kt=@alpha;
   variable i=0;
%
% Pileup fits can be fussy, so we try a couple of methods and
% loop through each a few times.  More clever ways of ensuring
% true minima can be programmed.
%
   foreach(alpha)
   {
% Me just being brain dead with roundoff errors above
      if (alpha[i] < 0.) alpha[i] = 0.;
      set_par(6,alpha[i],1);
      renorm_counts;
      set_fit_method("subplex;maxnfe=1000");
      loop(6)
          {
           fit_counts;
          }
      set_fit_method("marquardt;max_loops=100");
      loop(2)
         {
         fit_counts(&info);
         }
      xi[i]=info.statistic;
      nh[i]=get_par(1);
      bbnorm[i]=get_par(2);
      kt[i]=get_par(3);
       i++;
   }
   variable results=struct{nh, bbnorm, kt, alpha, xi};
   results.nh=@nh;
   results.bbnorm=@bbnorm;
   results.kt=@kt;
   results.alpha=@alpha;
   results.xi=@xi;
%
% Write out the results
%
   variable fp = fopen("pilesim.results","w");
   i=0;
   foreach(alpha)
      {
 ()=fprintf(fp,"%11.2e  %11.2e  %11.2e  %11.2e  %12.3e  \n",
                       alpha[i], nh[i], kt[i], bbnorm[i], xi[i]);
      i++;
      }
   () = fprintf(fp,"Alpha            Nh            kT          BB_norm      Xi \n");
   fclose(fp);
%
% One can return the results as a variable or structure, which can
% then be used in other S-Lang programs within ISIS
%
   return results;
   }

%
% Run the above, and just write the results to a file without doing
% anything further to them.
%

alphasrch();