path: root/noao/onedspec/irsiids/t_flatdiv.x

include	<imhdr.h>
include	<error.h>

define	MAX_NR_BEAMS	100	# Max number of instrument apertures

# T_FLATDIV -- Divide by a flat field spectrum. This is basically
#  a simple division of two vectors but with the following
#  additional functions:
#
#	1. Check the processing flag of the input spectra to avoid
#	   double processing, and set the flag if the processing is
#	   performed.
#	2. Trap division by zero errors
#	3. Optionally apply coincidence corrections

procedure t_flatdiv ()

int	root, start_rec
int	nrecs
int	len_flat
int	ccmode, qd_flag
real	dtime
real	power
bool	coincidence
pointer	sp, image, str, ofile, flat, recs, bstat, flatsp, im

int	clpopni(), clgeti(), clgwrd(), imgeti()
int	get_next_image(), decode_ranges()
real	clgetr()
bool	clgetb()
pointer	immap()
errchk	get_flatsp

begin
	call smark (sp)
	call salloc (image, SZ_FNAME, TY_CHAR)
	call salloc (str, SZ_LINE, TY_CHAR)
	call salloc (ofile, SZ_FNAME, TY_CHAR)
	call salloc (flat, SZ_FNAME, TY_CHAR)
	call salloc (recs, 300, TY_INT)
	call salloc (bstat, MAX_NR_BEAMS, TY_INT)

	# Open input file name template
	root   = clpopni ("input")

	# Get range specification if any
	call clgstr ("records", Memc[str], SZ_LINE)
	if (decode_ranges (Memc[str], Memi[recs], 100, nrecs) == ERR)
	    call error (0, "Bad range specification")

	# Get rootname for output files and starting record
	# Subtract 1 from start_rec because 1 will be added later.
	call clgstr ("output", Memc[ofile], SZ_FNAME)
	start_rec = clgeti ("start_rec") - 1

	# Get flat field spectrum root name
	call clgstr ("flat_file", Memc[flat], SZ_FNAME)

	# Apply coincidence corrections?
	coincidence = clgetb ("coincor")
	if (coincidence) {
	    ccmode = clgwrd ("ccmode", Memc[str], SZ_LINE, ",photo,iids,")
	    dtime  = clgetr ("deadtime")
	    power = clgetr ("power")
	}

	# Initialize beam number status
	call init_bs (Memi[bstat])

	# Initialize range decoder
	call reset_next_image ()

	# Loop over all input images - divide and make new image.
	# The output record number is incremented in all cases.
	while (get_next_image (root, Memi[recs], nrecs, Memc[image],
	    SZ_FNAME) != EOF) {
	    start_rec = start_rec + 1

	    # Open image
	    iferr (im = immap (Memc[image], READ_ONLY, 0)) {
		call erract (EA_WARN)
		next
	    }

	    # Get header
	    iferr (qd_flag = imgeti (im, "QD-FLAG"))
		qd_flag = -1
	    
	    # Verify divide flag
	    if (qd_flag != 0) {

		# Get flat field spectrum if needed
		call get_flatsp (im, flatsp, Memc[flat], Memi[bstat], len_flat)

		# Calibrate the current spectrum and make a calibrated version
		call divide (im, flatsp, len_flat, Memc[image], Memc[ofile],
		    start_rec, coincidence, ccmode, dtime, power)

	    } else {
		call eprintf ("[%s] already divided - ignored\n")
		    call pargstr (Memc[image])
	    }
	}

	# Update record number
	call clputi ("next_rec", start_rec)

	# Free space
	call sfree (sp)
	call clpcls (root)
end

# GET_FLATSP -- Load flat field spectrum for the current beam number

procedure get_flatsp (im, sp, flat_file, beam_stat, len_flat)

pointer	im, sp
char	flat_file[SZ_FNAME]
int	beam_stat[ARB], len_flat

int	i
int	beam, len[MAX_NR_BEAMS]
char	sfname[SZ_FNAME]
pointer	flatsp[MAX_NR_BEAMS], imflat

int	strlen(), imgeti()
pointer	imgl1r(), immap()
errchk	immap

begin
	# Determine beam number.

	iferr (beam = imgeti (im, "BEAM-NUM"))
	    beam = 0
	beam = beam + 1

	# Validate beam number
	if (beam < 1 || beam > MAX_NR_BEAMS) {
	    call eprintf (" Beam number out of range: %d - using 0\n")
		call pargi (beam)
	    beam = 1
	}

	# Has this beam already been loaded?
	if (IS_INDEFI (beam_stat[beam])) {

	    # Create file name
	    call strcpy (flat_file, sfname, SZ_FNAME)

	    # Flat field file names have beam number appended
	    call sprintf (sfname[strlen(sfname)+1], SZ_FNAME, ".%04d")
		call pargi (beam-1)

	    # Open spectrum
	    imflat = immap (sfname, READ_ONLY, 0)

	    # Allocate space for this beam's sensitivity spectrum
	    call salloc (flatsp[beam], IM_LEN(imflat,1), TY_REAL)

	    # Copy pixels into space
	    call amovr (Memr[imgl1r(imflat)], Memr[flatsp[beam]],
		IM_LEN(imflat,1))

	    # Must be careful that no division by zero occurs.
	    do i = 1, IM_LEN(imflat,1) {
		if (Memr[flatsp[beam]+i-1] == 0.0)
		    Memr[flatsp[beam]+i-1] = 1.0
	    }

	    # Mark this beam accounted for
	    beam_stat[beam] = 1
	    len[beam] = IM_LEN(imflat,1)

	    call imunmap (imflat)
	}

	# Point to the spectrum
	sp = flatsp[beam]
	len_flat = len[beam]

end

# DIVIDE -- Perform the division and create new spectrum

procedure divide (im, flat, len_flat, ifile, ofile, rec, 
	coincidence, ccmode, dtime, power)

pointer	im, flat
int	len_flat, rec, ccmode
real	dtime, power
char	ifile[ARB], ofile[ARB]
bool	coincidence

real	itm, imgetr()
int	i, co_flag, imgeti()
int	ncols, nlines
char	calfname[SZ_FNAME], original[SZ_FNAME]
pointer	imcal, rawpix, calpix

pointer	immap(), impl2r(), imgl2r()

begin
	# Find smallest length of the two possible spectra
	ncols = min (IM_LEN (im, 1), len_flat)

	# Create new spectrum.  Make up a name
	call sprintf (calfname, SZ_FNAME, "%s.%04d")
	    call pargstr (ofile)
	    call pargi (rec)

	call xt_mkimtemp (ifile, calfname, original, SZ_FNAME)
	imcal = immap (calfname, NEW_COPY, im)

	IM_NDIM(imcal) = IM_NDIM(im)
	IM_LEN (imcal,1) = ncols
	IM_PIXTYPE(imcal) = TY_REAL

	# Check for 2D spectrum
	if (IM_NDIM(im) > 1)
	    nlines = IM_LEN(im,2)
	else
	    nlines = 1

	# Copy across the image title
	call strcpy (IM_TITLE(im), IM_TITLE(imcal), SZ_LINE)

	# Operate on the pixels
	do i = 1, nlines {
	    rawpix = imgl2r (im,i)
	    calpix = impl2r (imcal,i)

	    # Apply coincidence correction if needed
	    co_flag = -1
	    if (coincidence) {
		iferr (co_flag = imgeti (im, "CO-FLAG"))
		    ;
	        if (co_flag < 1) {
		    iferr (itm = imgetr (im, "EXPOSURE"))
			iferr (itm = imgetr (im, "ITIME"))
			    iferr (itm = imgetr (im, "EXPTIME"))
				itm = 1.
		    call coincor (Memr[rawpix], Memr[rawpix], ncols,
			itm, co_flag, dtime, power, ccmode)
		}
	    }

	    call adivr (Memr[rawpix], Memr[flat], Memr[calpix], ncols)
	}

	call imaddi (imcal, "QD-FLAG", 0)
	if (co_flag != -1)
	    call imaddi (imcal, "CO-FLAG", co_flag)

	# Send user report
	call printf ("writing [%s]: %s\n")
	    call pargstr (original)
	    call pargstr (IM_TITLE(imcal))
	call flush (STDOUT)

	call imunmap (im)
	call imunmap (imcal)
	call xt_delimtemp (calfname, original)
end

# INIT_BS -- Initialize beam status flags

procedure init_bs (beam_stat)

int	beam_stat[ARB]

int	i

begin
	do i = 1, MAX_NR_BEAMS
	    beam_stat[i] = INDEFI
end