procedure striptz (image, outFile, rowFile, colFile, sublength, skycalc, skyFile, skyXoff, skyYoff)

string image	{prompt="Mosaiced image to be analysed"}
string outFile	{prompt="Name of output text file"}
string rowFile	{prompt="Text file listing row centres for extraction"}
string colFile	{prompt="Text file listing column centres for extraction"}
int sublength	{5, prompt="Half height and width of sub-array (in pix)"}
bool skycalc	{no, prompt="Do statistics on blank patch of sky too?"}
string skyFile  {prompt="Name of output sky values text file"}
int skyXoff	{10, prompt="X-offset of sky sub-array"}
int skyYoff	{10, prompt="Y-offset of sky sub-array"}
 

struct *flistRow	# list-directed structure to go thru rowFile
struct *flistCol	# list-directed structure to go thru colFile


begin

	# Variable Declarations:
	# struct row		# row value currently being centred upon
	# struct col		# column value currently being centred upon

	real row		# row value currently being centred upon
	real col		# column value currently being centred upon
	int x1, x2		# limits of sub-array in the x-direction
	int y1, y2		# limits of sub-array in the y-direction
	string bounds		# bounds of sub-array
	string subArray		# sub-array of image being sampled
	int j			# counter



	flistRow = rowFile
	flistCol = colFile		
	# putting the contents of the list contained in xxxFile
	# to the list-directed structure *flistxxx (where xxx = row/col)



	# Make sure the IMAGES package is loaded:
	if (! defpac ("images")) {
		print("IMAGES package not loaded")
		bye
	}



	# Getting rid of the ".imh" extension:
	j = strlen(image)
	if (substr (image, j-3, j) == ".imh")
		image = substr (image, 1, j-4)



	# Looping through each row/col centre point:
		while ( fscan(flistRow, row) != EOF) {
			print(row)

			while ( fscan(flistCol, col) != EOF) {
				print(col)

				x1=col-sublength
				x2=col+sublength
				y1=row-sublength
				y2=row+sublength
				# i.e. bounds of subArray

				bounds="["//x1//":"//x2//","//y1//":"//y2//"]"
				subArray=image//bounds
				print(subArray)
				# imstat(images=subArray, fields="image,npix,mean,stddev,min,max", format-, >> outFile)
				imstat(images=subArray, fields="image,npix,mean", format-, >> outFile)


				# ADD SKY-CALC IN HERE: IF/THEN etc
				if (skycalc==yes) {
					x1=x1+skyXoff
					x2=x2+skyXoff
					y1=y1+skyYoff
					y2=y2+skyYoff
					bounds="["//x1//":"//x2//","//y1 //":"//y2//"]"
					subArray=image//bounds
					print(subArray)
					imstat(images=subArray, fields="image,npix,mean", format-, >> skyFile)
				}


				print("-----------")

			}

			flistCol = colFile	
			# re-assigning flistCol for next scan thru col file

		}



end