#! /usr/bin/env isis-script

%; Time-stamp: <2003-02-05 13:23:55 dph> 
%; Directory: ~dph/libisis
%; File: contamarf
%; Author: D. Huenemoerder   dph@space.mit.edu
%; Original version: 2003.01.21
%;====================================================================
% version: 1
% version: 1.1   switch arg order for file and time.  It is more
%	         likely to use a different time than a different file.
%
% purpose: apply contamination function to an arf.
% Note: this is to be run from the unix prompt.
%
% *********************************************************************
%   ** NOTE: make sure you did NOT make the ARFs to be corrected with
%      the 'N0004' ACIS QE's (approved to LETG/ACIS, but sometimes used
%      for other obs).  The 'N0004' QE contains contamination for one
%      epoch.  Look for this HISTORY line in your ARF to check:
%      "HISTORY ACIS QE File: "
% *********************************************************************
%
%  contamarf arf_in arf_out [tstart [contam_file] ]
%
%    REQUIRES isis to be in your path.
%    ASSUMES the contamination file exists as referred to by
%            Acis_Contamination_File.
%
% most code John Davis' acis_contam.sl
%

%%% -----------------------   set this file name approprately.....   
   variable Acis_Contamination_File = "./acisD1999-08-13contamN0001.fits";
%%% -----------------------

   variable tstart = NULL;
   variable arf_in, arf_out;

   if ( (__argc < 3 ) or (__argc > 5) )
   {
     message("\nUsage:  contamarf arf_in arf_out [tstart [contam_file]]");
     message("\n Default TSTART read from arf_in.");
     vmessage(" Default contam_file=%s\n", Acis_Contamination_File);
     vmessage(" Apply the contamination function to the input ARF file.");
     vmessage(" Write the result to a new ARF.");
     vmessage(" TSTART is in seconds since MJDREF=1998.00\n");
     exit ;
   }

   if (__argc == 5 )  Acis_Contamination_File =  __argv[4];
   if (__argc > 3) () = sscanf(__argv[3], "%f", &tstart );
   arf_in = __argv[1];
   arf_out = __argv[2];


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 = NULL;

    if (NULL == stat_file (file))
      verror("Contamination file, %s, is missing", file);

   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] = 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);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% added by dph  2003.01.20
%
% correct an arf for contamination.
%
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 ) );
  () = system( sprintf("chmod u+w %s", 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;

}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  do the job here........

if (tstart == NULL)
 {
   () = contamarf( arf_in, arf_out );
 }
 else
 {
   () = contamarf( arf_in, arf_out, tstart );
 }