# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

include	<fset.h>
include	<imhdr.h>
include	<ttyset.h>
include <math/curfit.h>
include	<pkg/gtools.h>
include	"curfit.h"

define	VERBOSE_OUTPUT	1
define	LIST_OUTPUT	2
define	DEFAULT_OUTPUT	3

define	CF_UNIFORM	1
define	CF_USER		2
define	CF_STATISTICAL	3
define	CF_INSTRUMENTAL	4

# CF_FIT -- Called once for each curve to be fit.  

$for (rd)
procedure cf_fit$t (ic, gt, x, y, wts, nvalues, nmax, device, interactive, ofmt,
    power)

pointer	ic			# ICFIT pointer
pointer	gt			# Graphics tools pointer
PIXEL	x[nmax]			# X data values
PIXEL	y[nmax]			# Y data values
PIXEL	wts[nmax]		# Weights
int	nvalues			# Number of data points
int	nmax			# Maximum number of data points
char	device[SZ_FNAME]	# Output graphics device
int	interactive		# Fit curve interactively?
int	ofmt			# Type of output listing
bool	power			# Convert coeff to power series?

int	ncoeff, i, fd
PIXEL	xmin, xmax
pointer	sp, fname, gp, cv, coeff
pointer	gopen()
int	open(), $tcvstati()

begin
	call smark (sp)
	call salloc (fname, SZ_FNAME, TY_CHAR)

	# Determine data range and set up curve fitting limits.
	call alim$t (x, nvalues, xmin, xmax)
	call ic_putr (ic, "xmin", real (xmin))
	call ic_putr (ic, "xmax", real (xmax))

	if (interactive == YES) {
	    gp = gopen (device, NEW_FILE, STDGRAPH)
	    call icg_fit$t (ic, gp, "cursor", gt, cv, x, y, wts, nvalues)
	    call gclose (gp)
	} else 
	    # Do fit non-interactively
	    call ic_fit$t (ic, cv, x, y, wts, nvalues, YES, YES, YES, YES)

	# Output answers.
	call clgstr ("output", Memc[fname], SZ_FNAME)
	call ic_vshow$t (ic, Memc[fname], cv, x, y, wts, nvalues, gt)

	# Convert coefficients if requested for legendre or chebyshev
	if (power) {
 	    # Calculate and print coefficients
	    ncoeff = $tcvstati (cv, CVNCOEFF)
	    call salloc (coeff, ncoeff, TY_PIXEL)
	    call $tcvpower (cv, Mem$t[coeff], ncoeff)

	    fd = open (Memc[fname], APPEND, TEXT_FILE)
	    call fprintf (fd, "# Power series coefficients would be:\n")
	    call fprintf (fd, "# \t\tcoefficient\n")
	    do i = 1, ncoeff {
		call fprintf (fd, "# \t%d \t%14.7e\n")
		    call pargi (i)
		    call parg$t (Mem$t[coeff+i-1])
	    }
	    call close (fd)
	}

$if (datatype == r)
	call cvfree (cv)
$else
	call $tcvfree (cv)
$endif
	#call ic_close$t (ic)
	call sfree (sp)
end


# CF_LISTXY -- Print answers to STDOUT as x,y pairs.

procedure cf_listxy$t (cv, xvals, yvals, wts, nvalues)

pointer	cv			# Pointer to curfit structure
int	nvalues			# Number of data values
PIXEL	xvals[nvalues]		# Array of x data values
PIXEL	yvals[nvalues]		# Array of y data values
PIXEL	wts[nvalues]		# Array of weights

int	i
PIXEL	$tcveval()

begin
	do i = 1, nvalues {
	    call printf ("\t%14.7e \t%14.7e \t%14.7e \t%14.7e\n")
		call parg$t (xvals[i])
		call parg$t ($tcveval (cv, xvals[i]))
		call parg$t (yvals[i])
		call parg$t (wts[i])
	}
end

# IM_PROJECTION -- Given an image section of arbitrary dimension, compute
# the projection along a single axis by taking the average over the other
# axes.  We do not know about bad pixels.

procedure im_projection$t (im, x, y, w, npix, weighting, axis)

pointer	im			# Pointer to image header structure
PIXEL	x[npix]			# Index of projection vector
PIXEL	y[npix]			# Receives the projection vector
PIXEL	w[npix]			# Receives the weight vector
int	weighting		# Weighting of the individual points
int	npix			# Length of projection vector
int	axis			# The axis to be projected to (x=1)

int	i, lastv
long	v[IM_MAXDIM], nsum, totpix
pointer	pix
PIXEL	asum$t()
pointer	imgnl$t()
errchk	imgnl$t

begin
	if (im == NULL)
	    call error (1, "Image projection operator called with null im")
	if (axis < 1 || axis > IM_NDIM(im))
	    call error (2, "Attempt to take projection over nonexistent axis")


	# Set the y projection vector
	call aclr$t (y, npix)
	call amovkl (long(1), v, IM_MAXDIM)

	switch (axis) {
	case 1:
	    # Since the image is read line by line, it is easy to compute the
	    # projection along the x-axis (axis 1).  We merely sum all of the
	    # image lines.

	    while (imgnl$t (im, pix, v) != EOF)
		call aadd$t (Mem$t[pix], y, y, npix)

	default:
	    # Projecting along any other axis when reading the image line
	    # by line is a bit difficult to understand.  Basically, the
	    # element 'axis' of the V vector (position of the line in the
	    # image) gives us the index into the appropriate element of
	    # y.  When computing the projection over multiple dimensions,
	    # the same output element will be referenced repeatedly.  All
	    # of the elmenents of the input line are summed and added into
	    # this output element.

	    for (lastv=v[axis];  imgnl$t (im, pix, v) != EOF;  lastv=v[axis]) {
		i = lastv
		if (i <= npix)
		    y[i] = y[i] + asum$t (Mem$t[pix], IM_LEN(im,1))
	    }
	}

	# Now compute the number of pixels contributing to each element
	# of the output vector.  This is the number of pixels in the image
	# divided by the length of the projection.

	totpix = 1
	do i = 1, IM_NDIM(im)
	    if (i == axis)
		totpix = totpix * min (npix, IM_LEN(im,i))
	    else
		totpix = totpix * IM_LEN(im,i)
	nsum = totpix / min (npix, IM_LEN(im,axis))

	# Compute the average by dividing by the number if pixels summed at
	# each point.
	call adivk$t (y, PIXEL (nsum), y, npix)

	# Set the x and weight vectors
	do i = 1, npix {
	    x[i] = i
	    switch (weighting) {
	    case CF_STATISTICAL:
		if (y[i] > 0.0)
		    w[i] = 1.0 / y[i]
		else if (y[i] < 0.0)
		    w[i] = abs (1.0 / y[i])
		else
		    w[i] = 1.0
	    case CF_UNIFORM:
	        w[i] = 1.
	    default:
		w[i] = 1.
	    }
	}
end
$endfor