if (0 == is_defined ("Acis_Contamination_File"))
{
   public variable Acis_Contamination_File 
     = "/nfs/cxc/h1/dph/h3/CXC/ARD/acisD1999-08-13contamN0001.fits";
}

static define shift (x, n)
{
   variable len = length(x);
   variable i = [0:len-1];
   
   % allow n to be negative and large
   n = len + n mod len;
   return x[(i + n)mod len];
}

static define read_contam_file (file)
{
   variable s = fits_read_table (file);
   variable i, j;
   variable num_components = length (s.n_energy);

   variable c = struct
     {
	lambdas, alphas, times, betas, num_components, names
     };

   c.lambdas = Array_Type[num_components];
   c.alphas = Array_Type[num_components];
   c.times = Array_Type[num_components];
   c.betas = Array_Type[num_components];
   c.names = String_Type[num_components];

   variable lambdas = Double_Type[0];

   _for (0, num_components-1, 1)
     {
	i = ();

	j = [0:s.n_energy[i]-1];
	c.lambdas[i] = _A(s.energy[i,j]);
	c.alphas[i] = reverse (s.alpha[i,j]);
	
	j = [0:s.n_time[i]-1];
	c.times[i] = s.time[i,j];
	c.betas[i] = s.beta[i,j];

	lambdas = [lambdas, c.lambdas[i]];
%	c.names[i] = sprintf ("_%d", i);
	c.names[i] = s.feature[i];
     }

   % Use the same energy grid for all components.  Be sure to keep edges.
   % First prune duplicate points
   lambdas = lambdas[array_sort(lambdas)];
   i = where ((shift(lambdas, 1) - lambdas) != 0);
   lambdas = lambdas[i];

   _for (0, num_components-1, 1)
     {
	i = ();
	c.alphas[i] = interpol (lambdas, c.lambdas[i], c.alphas[i]);
     }
   c.lambdas = lambdas;
   c.num_components = num_components;
   return c;
}

static variable Contam_Info = read_contam_file (Acis_Contamination_File);
public define aciscontam_fit (lo, hi, par)
{
   variable t = par[0];
   variable alphas = Contam_Info.alphas;
   variable betas = Contam_Info.betas;
   variable times = Contam_Info.times;
   variable expon = 0.0;
   _for (0, Contam_Info.num_components-1, 1)
     {
	variable i = ();
	expon += (par[i+1]*interpol(t, times[i], betas[i])[0]) * alphas[i];
     }
   expon = interpol (0.5*(lo+hi), Contam_Info.lambdas, expon);
   expon [where (expon < 0.0)] = 0;
   return exp (-expon);
}

public define aciscontam_defaults (i)
{
   if (i == 0)
     return (0, 1, 0, NULL);
   
   return (1, 1, 0, 1);
}

add_slang_function ("aciscontam", ["time", Contam_Info.names]);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% added by dph  2003.01.20
%
% correct an arf for contamination.
%
public define contamarf()
{
  % contamarf(arf_in, arf_out[, tstart]);
  % If tstart is not given, look for TSTART in the arf_file

  variable arf_in, arf_out, tstart;

  if ( (_NARGS < 2) or (_NARGS > 3) )
  {
    message("Usage: contamarf(arf_in, arf_out[, tstart]);");
    return -1;
  }

  if( _NARGS == 3)    tstart = ();
  arf_out = ();
  arf_in = ();

  variable a, wlo, whi;

  a = fits_read_table( arf_in );

  (wlo,whi)=_A(a.energ_lo,a.energ_hi);

  if ( _NARGS == 2 )
    tstart = fits_read_key( arf_in, "TSTART");

  if (tstart == NULL)
  {
    vmessage("ERROR: No TSTART in %s\n", arf_in );
    return -1;
  }

  variable alphas = Contam_Info.alphas;
  variable betas = Contam_Info.betas;
  variable times = Contam_Info.times;
  variable expon = 0.0;

   _for (0, Contam_Info.num_components-1, 1)
     {
	variable i = ();
	expon += interpol(tstart, times[i], betas[i])[0] * alphas[i];
     }
  expon = interpol ( 0.5*(wlo+whi), Contam_Info.lambdas, expon );
  expon [where (expon < 0.0)] = 0;

  variable   arf = reverse(exp (-expon)) * a.specresp;
    
  () = system( sprintf("cp %s %s", arf_in, arf_out ) );

  variable col, hs ;
  variable fp = fits_open_file( arf_out + "[SPECRESP]", "w" );

  () = _fits_get_colnum (fp, "SPECRESP", &col);
  () = _fits_write_col (fp, col, 1, 1, arf );
  hs = sprintf ("contamarf: tstart = %f", tstart);
  () = _fits_write_history (fp, hs);
  hs = sprintf ("contamarf: contam file = %s", Acis_Contamination_File);
  () = _fits_write_history (fp, hs);

  fp = 0;
  return 0;

}