#! /usr/bin/env isis-script

%; Time-stamp: <2004-12-27 21:21:54 dph> 
%; Directory: ~dph/libisis
%; File: factor_rsp
%; Author: D. Huenemoerder   dph@space.mit.edu
%; Original version: 2004.12.23
%;====================================================================
% version: 1
%
% purpose: factor eefrac from grating rmf and multiply into grating arf.
%
% *********************************************************************
%
%  factor_rsp grmf_in garf_in garf_out [ eefrac_out ]
%
%    REQUIRES isis to be in your path.
%
% Adapted from contamarf 
%

variable rmf_in, arf_in, arf_out, eefrac_out ;

% the command name is the 1st arg.

#ifdef 0
% ***************************************************************
% debug:
() = printf("__argc=%d\n", __argc ) ; 
_for( 0, __argc-1, 1 )
{
    variable i=();
    () = printf("    %d: %S\n", i, __argv[i] ) ;
}
% ***************************************************************
#endif

if ( (__argc < 4 ) or (__argc > 5) )
{
    message("\nUsage:  factor_rsp grmf_in garf_in garf_out [eefrac_out]");
    vmessage(" Factor the eefrac from grmf_in and apply to garf_in");
    vmessage(" Write the result to a new ARF.");
    vmessage(" Optionally write the eefrac in arf format.\n");
    exit ;
}

eefrac_out = NULL ; 

if (__argc >= 4)
{
    rmf_in = __argv[1];
    arf_in = __argv[2];
    arf_out= __argv[3];
}

if (__argc == 5 )  eefrac_out =  __argv[4];

define factorrsp( rmf_in, arf_in, arf_out, eefrac_out)
{
    variable a, r, anew, eefrac ;

    a = load_arf( arf_in ) ;
    if ( a<0 ) verror("Failed to load garf_in: %S", arf_in );

    r = load_rmf( rmf_in ) ;
    if ( r<0 ) verror("Failed to load grmf_in: %S", rmf_in );

% for a test of compatibility, try to assign arf and rmf to a fake dataset:

    assign_rsp( a, r, 1 ) ; % possibly generate an error or warning.
    unassign_rmf( r ) ;     % free it up.

% NOTE: get_arf reverses arrays, so they are ascending  wavelength.
%       The OGIP files are ascending in energy.
%       Reverse on output.
%
    eefrac = get_arf( factor_rsp( r ) ) ;
    anew = get_arf( a ) ; 
    anew.value *=  eefrac.value ; 

    variable col, hs, fp ;
    
    () = system( sprintf("cp %s %s", arf_in, arf_out ) );
    () = system( sprintf("chmod u+w %s", arf_out ) );

    fp = fits_open_file( arf_out + "[SPECRESP]", "w" );
    () = _fits_get_colnum (fp, "SPECRESP", &col);
    () = _fits_write_col (fp, col, 1, 1, reverse(anew.value) );

% Add some history.
% NOTE: there doesn't seem to be an appropriate HDUCLASn qualifier
%       to indicate a "FULL" SPECRESP, as there is for the SPECRESP MATRIX.

    hs = "-------------- factor_rsp --------------------------" ; 
    () = _fits_write_history (fp, hs);

    fits_update_key( fp, "HDUCLAS3", "FULL", "EEFRAC included (non-std specresp key");

    hs = sprintf("factor_rsp %s %s %s %S", 
                    rmf_in, arf_in, arf_out, eefrac_out ) ; 
    () = _fits_write_history (fp, hs);

    hs = " ********* SPECRESP includes EEFRAC ************ " ; 
    () = _fits_write_history (fp, hs);

    hs = "----------------------------------------------------";
    () = _fits_write_history (fp, hs);

    fp = 0 ; 

    if ( eefrac_out != NULL )
    {
	() = system( sprintf("cp %s %s", arf_in, eefrac_out ) );
	() = system( sprintf("chmod u+w %s", eefrac_out ) );

	fp = fits_open_file( eefrac_out + "[SPECRESP]", "w" );

	() = _fits_get_colnum (fp, "SPECRESP", &col);
	() = _fits_write_col (fp, col, 1, 1, reverse(eefrac.value) );


	hs = "-------------- factor_rsp --------------------------";
	() = _fits_write_history (fp, hs);

	fits_update_key( fp, "HDUCLAS3", "EEFRAC", "only EEFRAC response (non-std specresp key)");

	hs = sprintf("factor_rsp %s %s %s %s", 
	                 rmf_in, arf_in, arf_out, eefrac_out ) ; 
	() = _fits_write_history (fp, hs);

	hs = " ********* SPECRESP is only EEFRAC ************ " ; 
	() = _fits_write_history (fp, hs);

	hs = "----------------------------------------------------";
	() = _fits_write_history (fp, hs);

	fp = 0;
    }
}

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

factorrsp( rmf_in, arf_in, arf_out, eefrac_out );