include <imhdr.h>
include	<cif.h>
include "../lib/cf_obj.h"
include "centerflux.h"

# Memory management.
define	Flux			Memd[flux+$1-1]
define	Mask			Memd[mask+$1-1]
define	Wave			Memd[wave+$1-1]

#---------------------------------------------------------------------------
.help centerflux Apr94 source
.ih
NAME
centerflux -- Find center fluxes for a set of input 1D spectra.
.endhelp
#---------------------------------------------------------------------------
procedure centerflux (input, output, gp, start_pixel, end_pixel,
		      center, size, low_reject, high_reject, niter)

pointer	input			# I:  Coordinated Input File (CIF) object.
pointer	output			# I:  CenterFlux (CF) output table object.
pointer	gp			# I:  Graphics descriptor.
int	start_pixel		# IO: Starting pixel.
int	end_pixel		# IO: End pixel.
int	center			# IO: Center of range.
int	size			# IO: Size of range.
double	low_reject, high_reject	# IO: Sigma cutoffs for mean determination.
int	niter			# IO: Number of iterations for mean det.

# Mask parameters.
pointer	mask			# Mask array.
pointer	mask_im			# Mask image descriptor.
int	mask_length		# Current length of mask array.

# Wavelength parameters.
pointer	wave			# Wavelength array.
pointer	wave_im			# Wavelength image descriptor.
int	wave_length		# Current length of wavelength array.

# Flux parameters.
pointer	flux			# Flux array.
pointer	flux_im			# Flux image descriptor.

# Pixel parameters.
int	length			# Length of the current flux array.
int	range			# Number of pixels to take mean of.
int	start			# Starting pixel in the data arrays.

# Results.
double	mean			# The mean.
double	min_wave, max_wave	# Minimum/maximum wavelength
int	n_points		# Number of pixels included in mean.
double	stddev			# The variance.

# Generic.
bool	cif_next()		# TRUE if another input file exists.
int	errcode()		# Get error code.
bool	first			# TRUE if first input processed.
int	imgeti()		# Get integer-valued header parameter.
pointer	imgl1d()		# Get a 1D line of image data.
pointer	immap()			# Open an image.
int	next_type		# Get next group or file.
pointer	sp			# Stack pointer.
bool	strne()			# Are strings not equal?
pointer	sx			# Generic.

begin
	call smark (sp)
        call salloc (sx, CF_SZ_TMPSTRING, TY_CHAR)
	
	# For each input spectrum, calculate the center flux and
	# write results to the table.
	mask_length = 0
	wave_length = 0
	first = true
	flux_im = NULL
	wave_im = NULL
	mask_im = NULL
        next_type = CIF_NEXT_GROUP
	while (cif_next (input, next_type)) {
            next_type = CIF_NEXT_GROUP

            # The current flux file being operated on.
            call printf ("centerflux: %s\n")
            call pargstr (CIF_p_file(input))

	    # Open the primary file as an image.
	    if (CIF_p_status(input) != CIF_SAME) {
                if (flux_im != NULL)
                    iferr (call imunmap (flux_im))
                        ;
		
		iferr (flux_im = immap (CIF_p_file(input), READ_ONLY, NULL)) {
		    call eprintf ("WARNING: centerflux: could not read image '%s'\n    ....skipping\n")
		    call pargstr (CIF_p_file(input))
		    next
		}
	    }
	    
	    # Open the wavelength file as an image.
	    if (CIF_in_status(input,CF_WAVE) != CIF_SAME && wave_im != NULL)
                iferr (call imunmap (wave_im))
                    wave_im = NULL
	    if (CIF_in_status(input,CF_WAVE) == CIF_OK)
		iferr (wave_im = immap (CIF_in_file(input, CF_WAVE),
					READ_ONLY, NULL)) {
		    wave_im = NULL
		    call eprintf ("WARNING: centerflux: could not open wavelength file '%s' for flux '%s'\n    ...ignoring wavelength information\n")
		    call pargstr (CIF_in_file(input,CF_WAVE))
		    call pargstr (CIF_p_file(input))
		}	    
            
	    # Open the mask file as an image.
	    if (CIF_in_status(input,CF_MASK) != CIF_SAME && mask_im != NULL)
		iferr (call imunmap (mask_im))
		    mask_im = NULL
	    if (CIF_in_status(input,CF_MASK) == CIF_OK)
		iferr (mask_im = immap (CIF_in_file(input, CF_MASK),
					READ_ONLY, NULL)) {
		    mask_im = NULL
		    call eprintf ("WARNING: centerflux: could not open mask file '%s' for flux '%s'\n    ...ignoring mask information\n")
		    call pargstr (CIF_in_file(input,CF_MASK))
		    call pargstr (CIF_p_file(input))
		}	    
            
	    # If first pass, set the observation parameters.
	    if (first) {

                # Sorry, not first any more.
                first = false

                # Get the constants.
		iferr (call imgstr (flux_im, "grating",
				    CF_GRATING(output), CF_SZ_GRATING)) {
		    call eprintf ("WARNING: centerflux: no GRATING keyword found in image %s.\nNo check for grating mode will be made.\n")
		    call pargstr (CIF_p_file(input))
		    call strcpy ("", CF_GRATING(output), CF_SZ_GRATING)
		}
		
		iferr (call imgstr (flux_im, "aperture",
                             CF_APERTURE(output), CF_SZ_APERTURE)) {
		    call eprintf ("WARNING: centerflux: no APERTURE keyword found in image %s.\nNo check for aperture will be made.\n")
		    call pargstr (CIF_p_file(input))
		    call strcpy ("", CF_APERTURE(output), CF_SZ_GRATING)
		}

		iferr (call imgstr (flux_im, "bunit",
                                    CF_UNITS(output), CF_SZ_UNITS))
                    ;
                if (wave_im != NULL)
                    iferr (call imgstr (wave_im, "bunit",
                                        CF_WUNITS(output), CF_SZ_WUNITS))
			;	
            }

            # Else check for consistency.
            else {
		ifnoerr (call imgstr (flux_im, "grating", Sx,
			     CF_SZ_TMPSTRING)) {
		    if (strne (Sx, CF_GRATING(output))) {
			call eprintf ("centerflux: warning: Grating %s of image %s\ndoes not match table's....skipping\n")
			call pargstr (Sx)
			call pargstr (CIF_p_file(input))
		    }
		}

		ifnoerr (call imgstr (flux_im, "aperture", Sx,
                             CF_SZ_TMPSTRING) ) {
		    if (strne (Sx, CF_APERTURE(output))) {
			call eprintf ("centerflux: warning: Aperture %s of image %s\ndoes not match table's....skipping\n")
			call pargstr (Sx)
			call pargstr (CIF_p_file(input))
			next_type = CIF_NEXT_FILE
			next
		    }
		}
            }
            
            # Check on dimensionality.
            if (IM_NDIM(flux_im) != 1) {
                call eprintf ("centerflux: warning: image %s is not one dimensional\n")
                call pargstr (CIF_p_file(input))
                next_type = CIF_NEXT_FILE
                next
            }

            # Read the flux data.
            flux = imgl1d (flux_im)
	    length = IM_LEN(flux_im,1)
            
            # Read the mask data if available.
            if (mask_im != NULL) {
                if (IM_LEN(mask_im,1) != length) {
                    call eprintf ("centerflux: warning mask %s is not the same length as flux image %s\n")
                    call pargstr (CIF_in_file(input, CF_MASK))
                    call pargstr (CIF_p_file(input))
                    call eprintf ("centeflux: Will not use the mask file.\n")
                    mask_length = INDEFI
                } else {
                    mask_length = length
                    call realloc (mask, mask_length, TY_DOUBLE)
                    call amovd (Memd[imgl1d (mask_im)],
                                Mask(1), mask_length)
                }
	    }

            # See if a new mask is necessary.  If so, then create an empty one.
            if (mask_length != length) {
                mask_length = length
                call realloc (mask, mask_length, TY_DOUBLE)
                call aclrd (Mask(1), mask_length)
            }

            # Read the wavelength data if available.
            if (wave_im != NULL) {
                if (IM_LEN(wave_im,1) != length) {
                    call eprintf ("centerflux: warning: wavelength image %s is not\nthe same length as flux image %s\n")
                    call pargstr (CIF_in_file(input, CF_WAVE))
                    call pargstr (CIF_p_file(input))
                    call eprintf ("centeflux: Will not use the wave file.\n")
                    wave_length = INDEFI
                } else {
                    wave_length = length
                    call realloc (wave, wave_length, TY_DOUBLE)
                    call amovd (Memd[imgl1d (wave_im)],
                                Wave(1), wave_length)
                }
	    }

            # Use the current wavelength vector if it is the same length.
            # Else, create an undefined vector.
            if (wave_length != length) {
                wave_length = length
                call realloc (wave, wave_length, TY_DOUBLE)
                call amovkd (INDEFD, Wave(1), wave_length)
            }

	    # If interactive, calculate means interactively
            if (gp != NULL)
                call cf_interactive (input, Flux(1),
                                     Mask(1), Wave(1), length, start_pixel,
                                     end_pixel, center, size,
				     niter, low_reject, high_reject, gp,
                                     mean, stddev, n_points, min_wave,
                                     max_wave)

	    # Else, just compute and go.
            else {

		# First calculate the range of pixels to be covered.
		iferr (call cf_calc_range (start_pixel, end_pixel, center,
					   size, length, start, range)) {
		    call eprintf ("WARNING: centerflux: specified pixel range out of bounds of image %s\nCheck start_pixel, end_pixel, center, and size parameters.\nskipping...\n")
		    call pargstr (CIF_p_file(input))
		    next
		}

		# Calculate the mean.
		iferr (call cf_mean (Flux(start), Mask(start), Wave(start),
                                     range, niter, low_reject, high_reject,
                                     mean, stddev, n_points, min_wave,
                                     max_wave, INDEFD)) {
                    if (errcode() == 2) {
                        call eprintf ("centerflux: warning: Mean does not converge for image %s.\nCheck mask, iterations, low_reject and high_reject\n")
                        call pargstr (CIF_p_file(input))
                    }
                }
	    }

            # Write the results to the table.
            call cf_add_point (output,
			       imgeti (flux_im, "carpos"),
			       imgeti (flux_im, "sporder"),
			       min_wave, max_wave, mean, stddev, n_points)

        }

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