include <imhdr.h>

define	SZ_KEYW		11	# size of buffer for header keyword

# fradius -- compute a measure of the radius
# This task computes the average radius as a measure of the spread
# of the psf.  
# The output from fextract is the input to this routine.
#
# Phil Hodge,  1-May-1990  Task created.

procedure fradius()

pointer input			# scr for name of input image
pointer output			# scr for name of output ascii file
#--
pointer sp
pointer rownum			# scratch for list of row numbers
pointer im			# pointer to imhdr struct
pointer v			# pointer to input data
char	keyword[SZ_KEYW]	# keyword for getting refoc position
real	refoc			# position of refoc mirror
real	radius			# characteristic radius of the PSF
int	fd			# fd for output file
int	numimages		# number of focus positions
int	naxis1			# width of input image
int	box1			# separation betw different stars in a slice
int	box2			# separation betw slices at diff focus positions
int	y1, y2			# limits of image section on second axis
int	k			# loop index for focus position
pointer immap(), imgs2r()
int	open(), imgeti(), imaccf()
int	findrow()
real	imgetr()

begin
	call smark (sp)
	call salloc (input, SZ_FNAME, TY_CHAR)
	call salloc (output, SZ_FNAME, TY_CHAR)
	call clgstr ("input", Memc[input], SZ_FNAME)
	call clgstr ("output", Memc[output], SZ_FNAME)

	im = immap (Memc[input], READ_ONLY, NULL)
	fd = open (Memc[output], NEW_FILE, TEXT_FILE)

	# Get the list of numbers of rows that have been filled in.
	call salloc (rownum, SZ_LINE, TY_CHAR)
	call imgstr (im, "rownum", Memc[rownum], SZ_LINE)

	# fextract puts different stars along the horizontal direction
	# and different focus positions (images input to fextract) in
	# the vertical direction.  The slice we are going to take from
	# the input image is of size naxis1 x box2.
	naxis1 = IM_LEN(im,1)
	box1 = imgeti (im, "box1")
	box2 = imgeti (im, "box2")
	numimages = IM_LEN(im,2) / box2		# should divide exactly

	# Process each focus position.
	do k = 1, numimages {

	    # Has this row been filled in?  If not, skip it.
	    if (findrow (Memc[rownum], k) == NO)
		next

	    # Get the refoc mirror position from the input header.
	    call sprintf (keyword, SZ_KEYW, "refoc%d")
		call pargi (k)
	    if (imaccf (im, keyword) == YES)
		refoc = imgetr (im, keyword)
	    else
		refoc = 0.

	    # Get the data (a horizontal slice of size naxis1 x box2)
	    # for all stars for a given focus position.
	    y1 = (k - 1) * box2 + 1
	    y2 = k * box2
	    v = imgs2r (im, 1, naxis1, y1, y2)

	    # Compute the radius of the PSF in this slice.
	    call fr_getrad (Memr[v], naxis1, box1, box2, radius)

	    # Write the refoc position and the radius to the file.
	    call fprintf (fd, "%15.7g %15.7g\n")
		call pargr (refoc)
		call pargr (radius)
	    call flush (fd)
	}

	call close (fd)
	call imunmap (im)
	call sfree (sp)
end


# fr_getrad -- get the radius of the PSF
# This computes the average radius weighted by the pixel values.

procedure fr_getrad (v, naxis1, box1, box2, radius)

real	v[naxis1,box2]		# i: input section
int	naxis1			# i: size of first axis of input image
int	box1, box2		# i: size of one box
real	radius			# o: radius of PSF
#--
double	sumv, sumrv		# for accumulating sums for radius
double	sumvx, sumvy, sumvv	# for getting centroid of PSF
double	x0, y0			# location of center of PSF
double	x, y			# distances from [x0,y0] along x & y axes
double	r			# distance from [x0,y0] to a pixel
int	nstars			# number of stars (boxes) at each focus position
int	star			# loop index for star (box) number
int	off			# offset in v for current star
int	j, k			# loop indexes

begin
	nstars = naxis1 / box1		# should divide exactly

	sumv = 0.
	sumrv = 0.

	# Do for each box (star) in the input slice.
	do star = 1, nstars {

	    # Offset along first axis of v to beginning of current box.
	    off = (star - 1) * box1

	    # Accumulate sums of the x & y positions weighted by data values.
	    sumvx = 0.
	    sumvy = 0.
	    sumvv = 0.
	    do k = 1, box2 {
		do j = off+1, off+box1 {
		    sumvx = sumvx + v[j,k] * j
		    sumvy = sumvy + v[j,k] * k
		    sumvv = sumvv + v[j,k]
		}
	    }
	    if (sumvv <= 0.)				# no data?
		next

	    # Find the centroid for the current box.  The x0 position is
	    # relative to the beginning of v, not the beginning of the box.
	    x0 = sumvx / sumvv
	    y0 = sumvy / sumvv
	    if (x0 <= off || x0 > off+box1)		# bad x position?
		next
	    if (y0 < 1 || y0 > box2)			# bad y position?
		next

	    do k = 1, box2 {
		y = k - y0
		do j = off+1, off+box1 {
		    x = j - x0
		    r = sqrt (x*x + y*y)
		    sumv = sumv + v[j,k]
		    sumrv = sumrv + r * v[j,k]
		}
	    }
	}

	if (sumv > 0.)
	    radius = sumrv / sumv
	else
	    radius = 0.
end