%
% To run this script, save this file as 'alphasrch_isis.sl' in
% the directory that you plan to run it. (Or choose your own favorite
% filename.)  Start up isis, then type at the isis prompt:
%
%    isis> .load alphasrch_isis
%
% *or*
%
%    isis> () = evalfile("alphasrch_isis.sl");
%
% Followed by (assuming the spectrum is in "pilesim.pi", and you want
% to fit that spectrum between 0.5 and 3 keV)
%
%    isis> results = alphasrch("pilesim.pi",0.5,3.0)
%
% "results" will be a structure (with the results, of course), which
% will also be written to the file: 'alphasrch_pilesim.pi_isis.dat'
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Load the xspec functions (normally done in my .isisrc)
%
import("xspec");
%
% Normally defined in my .isisrc, this does ignoring of data by energy
% range (as opposed to wavelength range)
%
define keV_note(a,b,c)
{
   variable bi, ci , 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;
}
%
% Normally defined in my .isisrc, this acts like ftools grppha,
% and groups the data into minimum counts, starting from 0 keV
% (as opposed to 0 wavelength).  *** ISIS starts off with the data
% ungrouped, and let's you *dynamically bin* without leaving the
% program ***
%
define grppha(a,mincounts)
{
   rebin_data(a,0);
   variable s = get_data_counts(a);
   variable len = length(s.value);
   variable idx = [0:len-1:1];
   idx[*] = 1;
   variable sn = 1, ssum = 0., i;
   for (i=len-1; i >=0; i--)
   {
      idx[i] = sn*idx[i];
      ssum += s.value[i];
      if (ssum >= mincounts)
      {
         sn = sn*(-1);
         ssum = 0.;
      }
   }
   rebin_data(a,idx);
   return 1;
}
%
% This reads the data, sets the model, then fits with several methods
%
define alphasrch(filename,lowe,highe)
{
   variable filenum = load_data(filename);
   () = grppha(filenum,20);
   () = keV_note(filenum,lowe,highe);
%
% Set the fit function, and choose to use the pileup kernel for fitting
%
   fit_fun("phabs(1)*( bbody(1) + powerlaw(1) )");
   set_kernel(1,"pileup");
%
   variable parms = [1,   2,     3,     4,      5,   8];
   variable mins =  [0.0, 0.0,   0.1,   0.0,    0.0, 0.0];
   variable maxs =  [0.4, 1.e-4, 0.3,   1.e-4,  4.0, 1.0];
   variable vals =  [0.1, 1.e-5, 0.192, 1.e-13, 2.,  1.0];
   variable froz =  [0,   0,     0,     0,      1,   1];
%
% Set the parameters
%
   variable i=0;
   loop(length(parms))
   {
      set_par(parms[i],vals[i],froz[i],mins[i],maxs[i]);
      i++;
   }
%
   variable info = struct{ statistic, num_variable_params, num_bins };
   variable r = struct{ nh, bbnorm, kt, plnorm, alpha, xi };
   r.alpha=[100:0:-5];
   r.alpha = r.alpha/100.;
   r.xi = @r.alpha;
   r.nh = @r.alpha;
   r.bbnorm = @r.alpha;
   r.kt = @r.alpha;
   r.plnorm = @r.alpha;
%
% We try several methods a few times each, since pile-up fits
% can be fussier than many models (due to many local minima
% in chi^2 space)
%
   variable method=["subplex;maxnfe=1000","marquardt;max_loops=100"];
%
% Open the file that we are writing to
%
   variable fp = fopen("alphasrch_"+filename+"_isis.dat","w");
%
   variable m;
   i=0;
   loop(length(r.alpha))
   {
      set_par(8,r.alpha[i],1);
%
      foreach(method)
      {
         m=();
     set_fit_method(m);
     loop(3)
         {
        () = renorm_counts;
            () = fit_counts(&info);
         }
      }
      r.xi[i] = info.statistic;
%
      r.nh[i] = get_par(1);
%   
      r.bbnorm[i] = get_par(2);
%
      r.kt[i] = get_par(3);
%
      r.plnorm[i] = get_par(4);
%
      ()=fprintf(fp,"%11.2e  %11.2e  %11.2e  %11.2e  %11.2e  %12.3e \n",
         r.alpha[i], r.nh[i], r.kt[i], r.bbnorm[i], r.plnorm[i], r.xi[i]);
      i++;
   }
   () = fprintf(fp,"Alpha//Nh//kT//BB_norm//PL_norm//Xi\n");
   () = fclose(fp);
%
   return r;
}