subroutine hpsort_indx(n,da,indx)
c
c***	Heapsort adapted from Numerical Recipes.  da is double precision
c***	array which is sorted into ascending order.  Indx is an
c***	integer array that defines the corresponding permutation and its
c***	inverse.  That is,
c***
c***		da_sort(i) = da_orig(indx(i))
c***		da_orig(indx(i)) = da_sort(i)
c***
c***	Analogous relations can be applied to associated arrays that 
c***	are to be sorted with da as the key.  For example, for array db,
c***
c***		db_sort(i) = db_orig(indx(i))
c***		db_orig(indx(i)) = db_sort(i)
c***
c***	G. Rybicki 3 Feb 95.
c
	INTEGER n
	DOUBLE PRECISION da(n),datemp
	INTEGER i,ir,j,k,indx(n),indxtemp
c
	if(n.lt.2)then
		return
	endif
c
	do i=1,n
		indx(i)=i
	enddo
c
	k=n/2+1
	ir=n
10	continue
		if(k.gt.1)then
			k=k-1
			datemp=da(k)
			indxtemp=indx(k)
		else
			datemp=da(ir)
			indxtemp=indx(ir)
			da(ir)=da(1)
			indx(ir)=indx(1)
			ir=ir-1
			if(ir.eq.1)then
				da(1)=datemp
				indx(1)=indxtemp
				go to 30
			endif
		endif
		i=k
		j=k+k
20		if(j.le.ir)then
			if(j.lt.ir)then
				if(da(j).lt.da(j+1))then
					j=j+1
				endif
			endif
			if(datemp.lt.da(j))then
				da(i)=da(j)
				indx(i)=indx(j)
				i=j
				j=j+j
			else
				j=ir+1
			endif
		go to 20
		endif
		da(i)=datemp
		indx(i)=indxtemp
	go to 10
30	continue
	return
	end
c
c*****************************************************************************
c
	subroutine trimult(a,b,c,x,y,n)
c
c	Does the multiplication of a vector by a tridiagonal matrix, i.e.,
c
c	   y(i) = a(i)*x(i-1) + b(i)*x(i) + c(i)*x(i+1), i=1,...,n,
c
c	with special cases i=1 and i=n, for which the terms containing
c	a(1) and c(n) are assumed absent.    G. Rybicki
c
	INTEGER i,n,n1
	DOUBLE PRECISION a(*),b(*),c(*),x(*),y(*)
c
	n1=n-1
	y(1)=b(1)*x(1)+c(1)*x(2)
	do i=2,n1
	   y(i)=a(i)*x(i-1)+b(i)*x(i)+c(i)*x(i+1)
	enddo
	y(n)=a(n)*x(n1)+b(n)*x(n)
	return
	end
c
c*****************************************************************************
c
	subroutine tridiag(a,b,c,x,y,n)
c
c***	Solves the tridiagonal system
c***
c***		a(i)*x(i-1)+b(i)*x(i)+c(i)*x(i+1)=y(i),   i=1,n
c***
c***	where the terms containing a(1) and c(n) do not appear.
c***	Arrays d and e are scratch arrays of dimension NDIM, which 
c***	must be at least equal to n-1.   G. Rybicki
c
	INTEGER i,k,n,n1,NDIM
	PARAMETER (NDIM=10000)	!alter to suit problem
	double precision a(*),b(*),bb,c(*),d(NDIM),e(NDIM),x(*),y(*)
	n1=n-1
	if(n1.gt.NDIM)then
		stop 'tridiag: NDIM too small'
	endif
	d(1)=-c(1)/b(1)
	e(1)=y(1)/b(1)
	do i=2,n1
		bb=b(i)+a(i)*d(i-1)
		d(i)=-c(i)/bb
		e(i)=(y(i)-a(i)*e(i-1))/bb
	enddo
	x(n)=(y(n)-a(n)*e(n1))/(b(n)+a(i)*d(n1))
	do k=n1,1,-1
		x(k)=e(k)+d(k)*x(k+1)
	enddo
	return
	end
c
c*****************************************************************************
c
	subroutine tridet(a,b,c,n,logdet,signdet)
c
c***	Calculate the determinant of the tridiagonal matrix with diagonals
c***	a, b, and c.  The sign is provided by signdet and the magnitude is
c***	provided by its logarithm logdet.  The elements a(1) and c(n)
c***	are irrelevant and are not used.    G. Rybicki
c
	INTEGER i,n,n1
	DOUBLE PRECISION a(*),b(*),c(*),logdet,signdet
	DOUBLE PRECISION bb,dd,ONE
	PARAMETER (ONE=1.d0)
c
	n1=n-1
	signdet=sign(ONE,b(1))
	logdet=log(abs(b(1)))
	dd=-c(1)/b(1)
	do i=2,n1
		bb=b(i)+a(i)*dd
		dd=-c(i)/bb
		signdet=signdet*sign(ONE,bb)
		logdet=logdet+log(abs(bb))
	enddo
	bb=b(n)+a(n)*dd
	signdet=signdet*sign(ONE,bb)
	logdet=logdet+log(abs(bb))
	return
	end