include "register.h"

#   expixi -- Extracts locations and values of subraster of pixels 
#               from science data.
#
#   Version  1.3#  26-APR-85 16:02:26     BALL           Submitted for system testing
#   ----------
#
#   Description:
#   ------------
#   This subroutine accepts the input image and will, for a specified
#   subimage size, return the values and locations of pixels centered
#   around a specified point. The output arrays are presented in such
#   a way that the NAG routine E02CAE/F can directly use them.
#   This (simple) version of the routine does not handle data masks.
#
#
#   Inputs:
#   -------
#   EXTRACTION_PARAMETERS
#   INPUT_SCIENCE_FILE
#
#   Outputs:
#   --------
#   OUTPUT_SUB_IMAGES
#   STATUS
#
#
#    Version    Date        Author     Description
#      10    17-APR-1984  Pat Murphy   Update PDL
#       7     8-FEB-1984  Pat Murphy   Added code for xmin/xmax (for NAG)
#       3    16-JAN-1984  Pat Murphy   Initial FORTRAN coding
#       2    23-AUG-1983  PAT MURPHY   Add PDL
#       1    23-AUG-1983  PAT MURPHY   Prolog Generation
#            21-JAN-1993    WILLIAMSON   Converted to SPP

procedure expixi (input, dim1,dim2, xcentr, ycentr, extent,
			x_values, y_values, values,
			xmin, xmax, ymin, ymax, status)

int dim2,dim1				# Size of input image
int extent				# Size of output sub-images
int i, ii, j, jj			# Local variables
int firstx, lastx			
int firsty, lasty
int status				# The returned status, usually OK

real values  [extent, extent]		# Pixel values in output
real x_values[extent]			# Output x-positions array
real y_values[extent]			# Output y-positions array
real xcentr, ycentr			# Center of range
real input[dim1, dim2]			# The input array

real xmin, xmax				# Range of subraster
real ymin, ymax


begin

        status = OK
#------------ If the extent asked for exceeds the array size, quit -----------

        if (extent > dim1 || extent > dim2) {
        status = FATAL
          return        
        }

#---------- Get the first and last points of a preliminary range -------------

        firstx = max (1, int (xcentr) - ((extent-1) / 2))
        lastx = firstx + extent - 1

        firsty = max (1, nint (ycentr) - ((extent-1) / 2))
        lasty = firsty + extent - 1

#------------- If we overshot the boundaries, make adjustments ---------------

        if (lastx > dim1) {
          firstx = firstx - (lastx - dim1)
          lastx = dim1
        }

        if (lasty > dim2) {
          firsty = firsty - (lasty - dim2)
          lasty = dim2
        }

#------------------------- Now fill the arrays-------------------------------

        do j = firsty, lasty {

          jj = j + 1 - firsty
          y_values [jj] = j

          do i = firstx, lastx {

            ii = i + 1 - firstx
            values[ii, jj] = input [i, j]

          }
        }

          do i = firstx, lastx {

            ii = i + 1 - firstx
            x_values [ii] = float[i]

          }

        xmin = firstx
        xmax = lastx
        ymin = firsty
        ymax = lasty

end