% First, load an appropriate ARF and RMF

     () = load_arf ("acis_meg_garf.fits.gz");
     () = load_rmf ("acis_meg_rmf.fits.gz");

% By assigning the responses to a non-existent
% data-set (#1) we initialize a matching (but empty) 
% dataset, which we will then populate with the
% 'fakeit' command:

     assign_arf (1, 1);
     assign_rmf (1, 1);
     list_data;

% Adjust the exposure time by changing the ARF value:

     set_arf_exposure (1, 50.e3);
     
% Load the APEC spectroscopy database..

     plasma (aped);

% ..and generate a multi-component model spectrum:

     variable t, wt, n, m, i; 
     variable v_turb = 200.0;     % Turbulent broadening in km/s

     t = [2.5e6, 4.6e6, 8.4e6];   % Array of Kelvin temperatures
     wt = [0.3, 0.5, 0.2];        % Relative weights for each Temp.

     n = length(t);
     m = Struct_Type[n];

     for (i = 0; i < n; i++)
       {
          m[i] = default_plasma_state();
          m[i].temperature = t[i];
          m[i].norm = 0.04 * wt[i];
	  m[i].vturb = v_turb;
       }
     
     define_model (m);
     list_model;

% ISIS supports fitting user-defined fit-functions
% implmented in S-Lang.  Here, we define a
% fit-function to evaluate our spectral model:
   
     define mymodel_fit (lo,hi,par)
     {
	return model_spectrum (lo,hi);
     }
     add_slang_function ("mymodel");
     
% Using this S-Lang fit-function to define our
% spectral model, we now fold that model through
% the instrument response to generate a fake
% data-set including Poisson noise.

     fit_fun ("mymodel(1)");

     % also include turbulent broadening in the line profile
     use_thermal_profile;    
     fakeit ();
     
% .. and plot a closeup of the region around Ne IX 
% labeling some of the brightest lines in the plot.

     variable pid = open_plot ("ex_fakeit.ps/cps");

     xrange (12.5,14.5);
     errorbars(1);
     title (latex2pg("Fake Data: 200 km~s^{-1} broadening"));

     plot_data_counts (1);

     % show brightest 10 model lines 
     % between 12.5 - 14.5 Angstrom
     variable b = brightest (10, where (wl(12.5,14.5)));
     plot_group (b);

     oplot_model_counts (1);

     close_plot (pid);