%get the isisscripts
require("isisscripts");

%we are going to use the model tbabs later, so let's set the proper
%abundances and cross-sections for this model (btw., this also needs
%to be done in xspec and sherpa if you happen to use this family of
%models!)
()=xspec_xsect("vern");
()=xspec_abund("wilm");

%make directory where we are going to save our results
variable path = "results/";
()=mkdir(path);

%%applying the same analysis to multiple datasets
%%this is by far not the only way to do it, but it shows some of the
%%basic tricks. Be creative!

%%the datasets we are going to use here is very simple: RXTE-PCA
%%observations of the rather faint X-ray binary LMC X-1. We are also
%%using just 3 spectra - my student worked on 553 over the summer; and
%%doing all 553 by hand would have been really, really ugly ...

%grab all the individual pha files
variable speclist = glob("lmcx1/9*/*/standard2f_0134off_top.pha");

%I think a loop is nicer for demonstration purposes, but this part of
%the script could also be written as a function that one would call
%onto for every entry in speclist

variable spec;
foreach spec (speclist) {

  variable pid = load_data(spec);

  %group the data between 3 & 25 keV to a signal to noise of five;
  %don't fortget the units - they are intrinsically
  %angstrom! 
  group(pid;bounds=[3,25],min_sn=5,unit="keV");
  %ignore all data below 3 and above 45 keV
  notice_values(pid,3,25;unit="kev");

  %let's define a simple isis function consisting of a diskblackbody
  %and powerlaw, modified by absorption
  fit_fun("tbabs*(diskbb+powerlaw)");

  % another function using a convolution model
  % fit_fun("tbabs*simpl(1,diskbb)");

  %let's renorm counts, i.e., change the norm to get closer to the data
  ()=renorm_counts;

  %let's set the equivalent hydrogen density to the literature value
  %and freeze it there
  set_par("tbabs(1).nH",1.0,1);

  %let's also set the disk temperature to a reasonable value for
  %starting the fit at
  set_par("diskbb(1).Tin",0.5);
  
  %fit the data
  % you can change the verbosity of the fit by setting the parameter
  % Fit_Verbose to a higher value
  % if you want the final output of fit_counts (esp. chi, degrees of
  % freedom etc. to be saved in a variable, define a variable and pass
  % it as a reference to fit_counts:
  %variable cc; ()=fit_counts(&cc);
  ()=fit_counts;
  
  %we can now save the bestfit parameters
  %not the string manipulations to get a nice name for the parameter
  %file; save_par is isis-intrinsic
  save_par(path+strtok(spec,"/")[1]+".par");

  %we can also use a rather nice function from the isisscripts to save
  %the whole fit, including information on the fti-quality, exposure,
  %paths to the data, grouping, noticing, etc. in one fits file
  %you can load the data in some other scripts by
  %fits_read_fit("name.fits")
  fits_save_fit(path+strtok(spec,"/")[1]+".fits");

  %let's do some error calculation - Mike has mentioned a parallelized
  %version for error calculations; here I use a different one that is
  %not parallelized (but should be soon), also from the isisscripts;
  %it returns the results in a rather nicely formatted structure that
  %I am going to save in a variable called cont
  variable cont = fit_pars(;saveoutput=0);

  %let's save this again
  %note that I am explicitly passing the variable containing the
  %confidence intervals to fits_save_fit
  save_par(path+strtok(spec,"/")[1]+".par");
  fits_save_fit(path+strtok(spec,"/")[1]+".fits",cont);

  %%this all also works for multiple datasets, say simultaneous
  %%observations from two different instruments, etc.
  
%  list_data;

  %%at the end of the run we have to delete the data loaded, less on
  %%the next loop we have two datasets loaded
  delete_data(pid);
  delete_rmf(pid);
  %if we had an arf, that would also be the place to delete the arf

}

%the folder results will now contain all the individual par and fits files; it
%is often useful to bundle the fits files into one - for example for
%later plotting correlations between different spectral parameters
%(rather more interesting for 500 spectra than for 3, of course)

fits_add_fit("results/all_bestfits.fits",glob("results/9*fits"));

%get me out of here ...
exit;