include	"../lib/cf_obj.h"
include	"../lib/spec_obj.h"
include	"../lib/ratio_obj.h"

# Memory management.
define	Int_stand		Memd[int_stand+$1-1]
define	Int_stand_err		Memd[int_stand_err+$1-1]
define	Rpos1			Memd[rpos1+$1-1]
define	Rpos2			Memd[rpos2+$1-1]

#---------------------------------------------------------------------------
.help abssenz May94 source
.ih
NAME
abssenz -- Compute sensitivity ratio between flux and standard star.
.endhelp
#---------------------------------------------------------------------------
procedure abssenz (flux, stand, ratio)

pointer	flux			# I:  The flux object.
pointer	stand			# I:  The standard star object.
pointer	ratio			# I:  The ratio object.

double	err			# Error.
int	i			# Generic.
pointer	int_stand		# Integrated flux of the standard star.
pointer	int_stand_err		# Integrated flux of the standard star error.
bool	itob()			# Convert integer to boolean.
double	value			# Normalized value.
double	norm			# Normalization factor.
pointer	rpos1			# Lowest pixel included in integration.
pointer	rpos2			# Highest pixel included in integration.
pointer	sp			# Stack pointer.
int	strlen()		# String length.
bool	strne()			# Strings not equal?
double	wave			# Wavelength.

begin
        call smark (sp)
        call salloc (int_stand, CF_N(flux), TY_DOUBLE)
        call salloc (int_stand_err, CF_N(flux), TY_DOUBLE)
        call aclrd (Int_stand_err(1), CF_N(flux))
        call salloc (rpos1, CF_N(flux), TY_DOUBLE)
        call salloc (rpos2, CF_N(flux), TY_DOUBLE)

        # Check that the required flux data is available.
	if (CF_N(flux) <= 0)
	    call error (1, "abssenz: No flux data.  Check that the flux table is not empty")
	if (CF_EXIST(flux,CF_FLUX) != YES)
	    call error (1, "abssenz: FLUX data not available! Check flux table for a FLUX column")
	if (CF_EXIST(flux,CF_WHIGH) != YES)
	    call error (1, "abssenz: WHIGH data not available! Check flux table for a WHIGH column")
	if (CF_EXIST(flux,CF_WLOW) != YES)
	    call error (1, "abssenz: WLOW data not available! Check flux table for a WLOW column")

	# Check that the required standard star data is available.
	if (SP_N(stand) <= 0)
	    call error (1, "abssenz: No standard star.  Check that the standard star table is not empty")
	if (SP_EXIST(stand,SP_FLUX) != YES)
	    call error (1, "abssenz: FLUX data not available! Check standard star table for a FLUX column")
	if (SP_EXIST(stand,SP_WAVE) != YES)
	    call error (1, "abssenz: WAVE data not available! Check standard star table for a WAVE column")

        # Check on units.  This is for warning purposes only.
        call strcpy (CF_WUNITS(flux), RA_WUNITS(ratio), RA_SZ_WUNITS)
        if (strlen (CF_WUNITS(flux)) >0 && strlen (SP_WUNITS(stand)) > 0) {
            call strlwr (CF_WUNITS(flux))
            call strlwr (SP_WUNITS(stand))
            if (strne (CF_WUNITS(flux), SP_WUNITS(stand)))
                call eprintf ("Warning: abssenz: Wavelength units between observation and standard star\n are not the same.\n")
        }
            
        if (strlen (CF_UNITS(flux)) >0 && strlen (SP_UNITS(stand)) > 0) {
            call strlwr (CF_UNITS(flux))
            call strlwr (SP_UNITS(stand))
            if (strne (CF_UNITS(flux),SP_UNITS(stand)))
                call eprintf ("Warning: abssenz: Flux units between observation and standard star\n are not the same\n")
        }
            
	# Setup output ratio for what data is valid.
	RA_EXIST(ratio,RA_VALUE) = YES
	RA_EXIST(ratio,RA_WAVE) = YES
	if (CF_EXIST(flux,CF_ERR) == YES ||
	    SP_EXIST(stand,SP_ERR) == YES)
	    RA_EXIST(ratio,RA_ERR) = YES
	RA_EXIST(ratio,RA_CARPOS) = CF_EXIST(flux,CF_CARPOS)
	RA_EXIST(ratio,RA_ORDER) = CF_EXIST(flux,CF_ORDER)

        # Copy aperture and grating information from the observation to
        # the ratio.
        call strcpy (CF_APERTURE(flux), RA_APERTURE(ratio), RA_SZ_APERTURE)
        call strcpy (CF_GRATING(flux), RA_GRATING(ratio), RA_SZ_GRATING)

        # If the input flux has no error, set it to zero.
        if (CF_EXIST(flux,CF_ERR) != YES && RA_EXIST(ratio,RA_ERR) == YES)
            call aclrd (CF_ERR_DATA(flux,1), CF_N(flux))

        # Match the standard star to the observed flux by integrating
        # the standard star.
        call ab_integrate (SP_WAVE_DATA(stand,1), SP_FLUX_DATA(stand,1),
			   SP_ERR_DATA(stand,1), itob (SP_EXIST(stand,RA_ERR)),
			   SP_N(stand), CF_WLOW_DATA(flux,1),
			   CF_WHIGH_DATA(flux,1), CF_N(flux),
                           Int_stand(1), Int_stand_err(1),
			   Rpos1(1), Rpos2(1))

        # Keep only standard fluxes that fall within the requested
        # integration and where the standard has measurable flux.
        do i = 1, CF_N(flux) {
            if (Rpos1(i) >= 1 &&
                Rpos2(i) <= SP_N(stand) &&
                Int_stand(i) > 0.d0) {

		# Calculate the effective wavelength.
                wave = (CF_WHIGH_DATA(flux,i) + CF_WLOW_DATA(flux,i)) / 2.d0

                # Normalize the integral to flux versus unit wavelength.
		norm = CF_WHIGH_DATA(flux,i)-CF_WLOW_DATA(flux,i)
                Int_stand(i) = Int_stand(i) / norm
		if (RA_EXIST(ratio,RA_ERR) == YES)
		    Int_stand_err(i) = Int_stand_err(i) / norm
		value = CF_FLUX_DATA(flux,i) / Int_stand(i)

		# Calculate errors.
		if (RA_EXIST(ratio,RA_ERR) == YES && Int_stand(i) != 0.d0 &&
                    CF_FLUX_DATA(flux,i) != 0.d0) {
		    err = ((Int_stand_err(i)/Int_stand(i)) +
			  (CF_ERR_DATA(flux,i)/CF_FLUX_DATA(flux,i))) *
			  value
		} else
		    err = INDEFD

		# Add the new point.
		call ra_add_point (ratio, CF_CARPOS_DATA(flux,i),
				   CF_ORDER_DATA(flux,1), value,
				   wave, err)
            }
        }

        # Check to make sure that there is a ratio.
        if (RA_N(ratio) <= 0) {
	    call ra_dump (ratio)
            call error (1, "abssenz: no wavelength overlap between observation and standard star")
	}

        # That's all folks.
        call sfree (sp)
end
#---------------------------------------------------------------------------
# end of abbsenz
#---------------------------------------------------------------------------