include <error.h>
include <gset.h>
include	<cif.h>
include "centerflux.h"

# Memory management.
define	Flux_real		Memr[flux_real+$1-1]
define	Mx			Memr[mx+$1-1]
define	My			Memr[my+$1-1]
define	New_mask		Memd[new_mask+$1-1]
define	Strval			Memc[strval]

# Dreaded gotos'
define	center_label_		10
define	next_label_		20

#---------------------------------------------------------------------------
.help cf_interactive Apr94 source
.ih
NAME
cf_interactive -- Interactively view the mean fitting.
.endhelp
#---------------------------------------------------------------------------
procedure cf_interactive (input, flux, mask, wave, length, start_pixel,
			  end_pixel, center, size,
			  niter, low_reject, high_reject, gp, mean, stddev,
			  n_points, min_wave, max_wave)

pointer	input			# I:  The Coordinated Input File object.
double	flux[length]		# I:  The flux array.
double	mask[length]		# I:  The mask array.
double	wave[length]		# I:  The wavelength array.
int	length			# I:  Length of the input data.
int	start_pixel		# IO: Starting pixel.
int	end_pixel		# IO: Ending pixel.
int	center			# IO: Center.
int	size			# IO: Pixel range size.
int	niter			# IO:  Number of iterations to perform.
double	low_reject		# IO: Low sigma rejection.
double	high_reject		# IO: High sigma rejection.
pointer	gp			# IO:  Graphics descriptor.
double	mean			# O:  The mean.
double	stddev			# O:  The error.
int	n_points		# O:  Number of points in the fit.
double	min_wave, max_wave	# O:  Minimum/maximum wavelength of mean.

# Declarations
real	area_end		# End of the area.
bool	area_first		# True if first hit of 'a' command.
real	area_start		# Start of the area.
int	cf_wave2pix()		# Wavelength to pixel.
int	clgcur()		# Get value of the graphics cursor.
int	command			# Colon command
int	ctod()			# String to double-valued.
int	ctoi()			# String to integer-valued.
pointer	flux_real		# Real-value version of flux.
bool	full			# True to plot all the flux, just not fit.
int	graph_start, graph_len	# Pixel start and length of the graph.
int	i, j, n			# Generic.
int	key			# Key initiating action.
real	mean_real		# Real version of the mean.
pointer	mx, my			# Arrays.
pointer	new_mask		# Calculated mask.
real	px1, px2		# Range of plot.
int	range			# Number of pixels to average.
bool	redraw			# True to redraw the plot.
real	rx, ry, rz		# Generic.
bool	show_wave		# Plot in wavelength coordinates?
pointer	sp			# Stack pointer
int	start			# Starting pixel to average.
pointer	strdic()		# Get index into string dictionary.
int	strlen()		# Length of string.
pointer	strval			# String value of the graphics cursor.
pointer	sx, sy			# Generic strings.
int	wcs			# Current graphics WCS.
real	x, y			# Graphics cursor position

data	full/false/
data	show_wave/true/

begin
	call smark (sp)
	call salloc (flux_real, length, TY_REAL)
	call salloc (mx, length, TY_REAL)
	call salloc (my, length, TY_REAL)
	call salloc (new_mask, length, TY_DOUBLE)
	call salloc (strval, SZ_COMMAND, TY_CHAR)
	call salloc (sx, CF_SZ_TMPSTRING, TY_CHAR)
	call salloc (sy, CF_SZ_TMPSTRING, TY_CHAR)
	
        # Convert the flux to type real.
        call achtdr (flux, Flux_real(1), length)
        
        # The graphics interaction loop.
        call greactivate (gp, 0)
        key = 'r'
	if (IS_INDEFD(wave[1]))
	    show_wave = false
	repeat {
	    switch (key) {
            case ':':
                call sscan (Strval)

		# Get the variable
		call gargwrd (Sx, CF_SZ_TMPSTRING)
		call gargwrd (Sy, CF_SZ_TMPSTRING)
		command = strdic (Sx, Sx, CF_SZ_TMPSTRING, CF_COLON_COMMANDS)
		switch (command) {
		case CF_COLON_CENTER:
		    if (strlen (Sy) <= 0) {
			call printf ("center = %d")
			call pargi (center)
		    } else {
			j = 1
			i = ctoi (Sy, j, center)
			goto center_label_
		    }
		    
		case CF_COLON_END:
		    if (strlen (Sy) <= 0) {
			call printf ("end = %d")
			call pargi (end_pixel)
		    } else {
			j = 1
			i = ctoi (Sy, j, end_pixel)
			redraw = true
			full = false
		    }
		    
		case CF_COLON_LOW_REJECT:
		    if (strlen (Sy) <= 0) {
			call printf ("low reject = %g")
			call pargd (low_reject)
		    } else {
			j = 1
			i = ctod (Sy, j, low_reject)
			redraw = true
		    }

		case CF_COLON_HIGH_REJECT:
		    if (strlen (Sy) <= 0) {
			call printf ("high reject = %g")
			call pargd (high_reject)
		    } else {
			j = 1
			i = ctod (Sy, j, high_reject)
			redraw = true
		    }

		case CF_COLON_ITERATIONS:
		    if (strlen (Sy) <= 0) {
			call printf ("iterations = %g")
			call pargi (niter)
		    } else {
			j = 1
			i = ctoi (Sy, j, niter)
			redraw = true
		    }

		case CF_COLON_SIZE:
		    if (strlen (Sy) <= 0) {
			call printf ("size = %d")
			call pargi (size)
		    } else {
			j = 1
			i = ctoi (Sy, j, size)
			redraw = true
			full = false
		    }
		    
		case CF_COLON_START:
		    if (strlen (Sy) <= 0) {
			call printf ("start = %d")
			call pargi (start_pixel)
		    } else {
			j = 1
			i = ctoi (Sy, j, start_pixel)
			redraw = true
			full = false
		    }
		    
		default:
		    call printf ("No such colon command")
                }
		
            case '?':
                call gpagefile (gp, "develop$doc/centerflux.key", NULL)

	    case 'I':
		call gclose (gp)
		gp = NULL
		break
		
	    case 'c':
		center = x
		if (show_wave)
		    center = cf_wave2pix (x, wave, length)

center_label_
		# Remove the start/end pixels but remember the size.
		if (!IS_INDEFI(start_pixel) && !IS_INDEFI(end_pixel)) {
		    size = end_pixel - start_pixel + 1
		    start_pixel = INDEFI
		    end_pixel = INDEFI
		}

		# Redraw.
		redraw = true
		full =  false

	    case 'f':
		full = !full
		redraw = true
		
	    case 'p':
		if (area_first) {
		    area_start = x
		    area_first = false
		    call printf (" move cursor and hit 'p' again")
		} else {
		    area_end = x
		    area_first = true

		    # If in wavelengths, find the pixels.
		    if (show_wave) {
			area_start = cf_wave2pix (area_start, wave, length)
			area_end = cf_wave2pix (area_end, wave, length)
		    }

		    # Reset the beginning and ending pixels.
		    start_pixel = min (area_start, area_end)
		    end_pixel = max (area_start, area_end)

		    # Redraw the graph.
		    redraw = true
		    full = false
		}
		
            case 'q':
                break

            case 'r':
                redraw = true

	    case 'w':
		if (!IS_INDEFD(wave[1])) {
		    show_wave = !show_wave
		    redraw = true
		}
            }

            if (redraw) {

                # Reset the drawing.
                redraw = false
		area_first = true

		# Get start and range of the pixels to average.
		iferr (call cf_calc_range (start_pixel, end_pixel, center, size,
					   length, start, range)) {
		    call eprintf ("centerflux: specified start/end pixels not within image %s.\nRespecify range using the 'p', 'c', or the start, end, center,\nand size colon commands\n")
		    call pargstr (CIF_p_file(input))
		    goto next_label_
		}
		
                iferr (call cf_mean (flux(start), mask(start), wave(start),
                                     range, niter, low_reject, high_reject,
                                     mean, stddev, n_points, min_wave,
                                     max_wave, New_mask(start)))
                    call eprintf ("centerflux: warning: mean cannot be found.\nModify low_reject, high_reject, and iterations\n")

		# Show full spectrum?
		if (full) {
		    graph_start = 1
		    graph_len = length
		} else {
		    graph_start = start
		    graph_len = range
		}

		# Plot in wavelength coordinates?
		if (show_wave) {
		    px1 = wave[graph_start]
		    px2 = wave[graph_start+graph_len-1]
		} else {
		    px1 = graph_start
		    px2 = graph_start+graph_len-1
		}

		# Plot it.
                call gclear (gp)
		call gploto (gp, Flux_real(graph_start), graph_len, px1, px2,
			     CIF_p_file(input))

		# Plot the mean and deviations
		mean_real = mean
		call gline (gp, px1, mean_real, px2, mean_real)
		rx = 2 * stddev
                call ggscale (gp, (px1+px2)/2., mean_real, ry, rz)
                rx = rx / rz
		call gmark (gp, (px1+px2)/2., mean_real, GM_VEBAR,
			    0.25 * rx, rx)

                # If full plot, show what range is include in the mean.
		if (full) {
		    px1 = start
		    px2 = start+range-1
		    if (show_wave) {
			px1 = wave[px1]
			px2 = wave[px2]
		    }
		    call ggwind (gp, rz, rz, rx, ry)
		    rz = ry - rx
		    call gmark (gp, (px1+px2)/2, rx + (0.1 * rz),
				GM_HEBAR, px1-px2, -0.05 * abs(rz))
		}

                # If not full, show which points where masked out.
                if (!full) {
                    n = 0
                    do i = start, start+range-1 {
                        if (New_mask(i) > 0.d0) {
                            n = n + 1
                            if (show_wave)
                                Mx(n) = wave[i]
                            else
                                Mx(n) = i
                            My(n) = flux[i]
                        }
                    }
                    if (n > 0)
                        call gpmark (gp, Mx(1), My(1), n, GM_CROSS,
                                     0.01, 0.01)
                }
	    }
next_label_
	    x = y
	} until (clgcur ("coords", x, y, wcs, key, Strval, SZ_COMMAND) == EOF)

	# That's all folks.
	call sfree (sp)
end
#---------------------------------------------------------------------------
# end of cf_interactive
#---------------------------------------------------------------------------
int procedure cf_wave2pix (wave, wave_array, length)

real	wave			# I:  Wavelength in question.
double	wave_array[length]	# I:  Wavelength array
int	length			# I:  Length of wavelength array.

int	i, pix			# Pixel positions.

begin
	if (wave < wave_array[1])
	    pix = 1
	else if (wave > wave_array[length])
	    pix = length
	else
	    do i = 1, length-1
		if (wave >= wave_array[i] && wave <= wave_array[i+1]) {
		    pix = i
		    break
		}

	return (pix)
end
#---------------------------------------------------------------------------
# End of cf_wave2pix
#---------------------------------------------------------------------------