define	LOGPTR	20			# log2(maxpts) (1e6)

# ASRT2 -- Vector Quicksort.  The output vectors may be the same as the
# input vectors.

procedure asrt2$t (a, b, c, d, npix)

PIXEL	a[ARB], b[ARB]		# Primary input, output arrays
PIXEL	c[ARB], d[ARB]		# Secondary input, output arrays
int	npix			# Number of pixels

PIXEL	pivot, temp
int	i, j, k, p, lv[LOGPTR], uv[LOGPTR]
define	swap {temp=$1;$1=$2;$2=temp}

begin
	call amov$t (a, b, npix)		# in place sort
	call amov$t (c, d, npix)		# in place sort

	lv[1] = 1
	uv[1] = npix
	p = 1

	while (p > 0) {
	    if (lv[p] >= uv[p])			# only one elem in this subset
		p = p - 1			# pop stack
	    else {
		# Dummy do loop to trigger the Fortran optimizer.
		do p = p, ARB {
		    i = lv[p] - 1
		    j = uv[p]

		    # Select as the pivot the element at the center of the
		    # array, to avoid quadratic behavior on an already sorted
		    # array.

		    k = (lv[p] + uv[p]) / 2
		    swap (b[j], b[k])
		    swap (d[j], d[k])
		    pivot = b[j]			# pivot line

		    while (i < j) {
		    $if (datatype == x)
			for (i=i+1;  abs(b[i]) < abs(pivot);  i=i+1)
		    $else
			for (i=i+1;  b[i] < pivot;  i=i+1)
		    $endif
			    ;
			for (j=j-1;  j > i;  j=j-1)
		    $if (datatype == x)
			    if (abs(b[j]) <= abs(pivot))
		    $else
			    if (b[j] <= pivot)
		    $endif
				break
			if (i < j) {		# out of order pair
			    swap (b[i], b[j])	# interchange elements
			    swap (d[i], d[j])	# interchange elements
			}
		    }

		    j = uv[p]			# move pivot to position i
		    swap (b[i], b[j])		# interchange elements
		    swap (d[i], d[j])		# interchange elements

		    if (i-lv[p] < uv[p] - i) {	# stack so shorter done first
			lv[p+1] = lv[p]
			uv[p+1] = i - 1
			lv[p] = i + 1
		    } else {
			lv[p+1] = i + 1
			uv[p+1] = uv[p]
			uv[p] = i - 1
		    }

		    break
		}
		p = p + 1			# push onto stack
	    }
	}
end