C+
C NAME:
C	ice_pack
C PURPOSE:
C	Compresses an integer array using a modified Rice algorithm
C CATEGORY:
C	gen/for/lib
C CALLING SEQUENCE:
	function ice_pack(kmax, kshift, nlenbit, nsign, nperm, nOrig, Orig, nPack, Pack, iPack)
C INPUTS:
C	kmax		integer		cutoff for applying compression (see PROCEDURE)
C					(kmax > 0)
C	kshift		integer		minimum # bits copied for each difference
C					(kshift <kmax)
C	nlenbit		integer		# bits encoded  when the full value instead
C					of a difference is used (see PROCEDURE)
C	nsign		integer		-1: all values in Orig are <  zero
C					+1: all values in Orig are >= zero
C					 0: both positive and negative values are present
C	nperm(0:kmax-kshift)
C			integer		permutation table (see PROCEDURE)
C					(if not available use nperm(i) = i [i=0,kmax])
C	nOrig		integer		# elements in Orig
C	Orig(nOrig)	integer		array to be compressed
C	nPack		integer		# of elements in Pack available for storing compressed data
C OUTPUTS:
C	ice_pack	integer		# bits in Pack used to hold the compressed data.
C
C					On failure ibit = 0 is returned (happens only if the array
C					Pack is not big enough to hold the compressed data). Since Pack
C					is an integer*4 array nPack*nbit (nbit=32) bits are available.
C
C					Bits 0..ibit-1 are filled. The first 1+(ibit-1)/nbit elements
C					of Pack are used. Only ib = 1+mod(ibit-1,nbit) bits (0..ib-1)
C					of the last element are  used.
C	Pack(nPack)	integer		compressed data; array elements 1..iPack are used.
C	iPack		integer		# elements in Pack used (= 1+(ibit-1)/nbit)
C SEE ALSO:
C	ice_unpack, ice_analyze, ice_write, ice_read
C PROCEDURE:
C >	The subroutine href=ice_analyze= provides optimum values for kmax, kshift,
C	nlenbit, nsign and nperm.
C >	The compression factor is ibit/(nOrig*nbit). Use nbit=16 instead of nbit=32
C	if the original data were actually integer*2 stored in integer*4 form.
C >	The technique is based on the algorithm described by Michael W. Richmond
C	and Nancy E. Ellman, in March 1996? paper which we have in preprint form. 
C	Reference?
C
C	COMPRESSION ALGORITHM
C	---------------------
C
C	The Rice algorithm stores an integer in the following bit pattern:
C		0......0 1 x..........x
C		|______|   |__________|
C		  ipad       itop+1
C	i.e. a string of ipad 0-bits, followed by a terminating 1-bit, followed by a
C	string of itop+1 bits encoding the integer value.
C
C	The actual number encoded is either the full value from the input array, or the
C	difference with the previous element in the array (For the first
C	element in the array the 'previous element' is assumed to be zero).
C
C	The two cases are distinguished using the kmax value: if the absolute difference
C	is less than 2^(kmax-1) then the difference is encoded; if not then the full value
C	is encoded.
C
C	CASE A: Non-zero absolute difference less than 2^(kmax-1)
C	---------------------------------------------------------
C	2^(kmax-1) has only bit kmax-1 set. I.e. for absolute differences less than
C	2^(kmax-1) this bit, and all higher bits, are NOT set, i.e. only the first
C	kmax-1 bits 0...kmax-2 might be set.
C
C	Let itop be the highest 1-bit of the absolute difference. For a non-zero difference
C	0 <= itop <= kmax-2. Since bit itop is by definition set, the absolute difference
C	is fully described by the first itop bits 0...itop-1. An additional bit is needed
C	to encode the sign of the difference. So itop+1 bits are needed to store the difference:
C	first ipad=itop+1 0-bits, followed by a terminator 1-bit, indicate how many bits are
C	used. The terminator bit is then followed by the first itop bits of the absolute difference
C	followed by a sign bit. The total number of bits used is 2*itop+3.
C
C		0..........0 1 x..........x
C		|__________|   |__________|
C		ipad=itop+1       itop+1
C	
C	CASE B: Zero difference
C	-----------------------
C	No bit at all is needed to store a zero, except for a terminator bit. Note that this 
C	is CASE A with ipad = 0, itop = -1.
C
C	CASE C:  Absolute difference greater/equal 2^(kmax-1)
C	-----------------------------------------------------
C	In this case the full value instead of the difference is encoded.
C	For case A the maximum number of leading 0-bits is kmax-1. A string of kmax 0-bits
C	is used to indicate that a full value is encoded. After the terminator bit follow
C	the first nlenbit bits. The total number of bits used is kmax+1+nlenbit
C
C		0..........0 1 x..........x
C		|__________|   |__________|
C		ipad=kmax         nlenbit
C
C	MINIMUM NUMBER OF BITS
C	----------------------
C
C	The above describes the the compression algorithm for the case kshift=0.
C	Setting kshift to non-zero value modifies the way differences are stored: for each
C	difference the kshift lowest bits are always stored. The differences are now encoded
C	as follows:
C
C	Zero difference (ipad = 0)
C		1 0
C	The 1-bit is the terminating bit (i.e. no leading 0-bits), followed by a 0-bit to
C	indicate the presence of a zero difference.
C
C	Difference needing 1..kshift bits (1 <= ipad <= kshift)
C		1 1 s x.....x
C		      |_____|
C		       kshift
C	The leading 1-bit is the terminating bit (i.e. no leading 0-bits). It is followed by
C	a 1-bit to be able to distinguish it from a zero difference. Then a bit follows
C	storing the sign of the difference (s-bit). Then follow the lowest kshift bits.
C
C	Differences needing kshift+1...kmax-1 bits (kmax-1 <= ipad <= kmax-1 above)
C		0.......0 1 x........x
C	        |_______|   |________|
C	       ipad-kshift     ipad
C
C	Full value (ipad=kmax)
C		0.......0 1 x........x
C	        |_______|   |________|
C	       kmax-kshift    nlenbit
C
C	For kshift > 0 a smaller number of leading zeroes is needed to store a value. A small
C	penalty is paid for zero differences (1 extra bit); for differences needing
C	1 <= ipad <= kshift bits the penalty is: 2+kshift-2*ipad. This is positive for
C	ipad < 1+kshift/2. The penalty is largest for ipad = 1 (kshift bits).
C
C	PERMUTATION TABLE
C	-----------------
C
C	The number of leading 0-bits varies from izero=0 to izero=kmax-kshift. These differences
C	serve merely as identifiers for the type of information stored after the terminating
C	1-bit. Rather than using them directly as described above, it is more efficient to
C	use a permutation table based on the histogram of the number of bits needed for to store
C	the differences. E.g. if for some sky image 3-bit difference occur most, followed by 2-bit,
C	then 1-bit, then 0-bit, then 4-bit, 5-bit, etc. then the permutation table would be
C		nperm = 3,2,1,0,4,5...kmax-kshift
C MODIFICATION HISTORY:
C	AUG-2000, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C		Rewrite of Andy Buffington's ricearc.for.
C		Main modifications:
C		- large intermediate byte array for storing individual bits is removed
C		- # leading zeroes used to store full values was reduced from kmax+1 to kmax
C		- use of kshift allows for slightly better compression
C-
	    integer	kmax
	    integer	kshift
	    integer	nlenbit
	    integer	nsign
	    integer	nperm(0:kmax-kshift)		! Permutation table
	    integer	nOrig
	    integer	Orig(0:nOrig-1)
	    integer	nPack
	    integer	Pack(0:nPack-1)
	    integer	iPack

	parameter	(nbit  = 32)
	parameter	(nbit1 = nbit-1)

	integer		nrank(0:nbit1)

	itopbit = nlenbit-1

	!-------
	! Initialize ...

	do i=0,kmax-kshift
	    nrank(nperm(i)) = i
	end do

	do i=0,nPack-1
	    Pack(i) = 0					! Clear all bits in Pack
	end do

	!-------
	! idelmax has all bits zero, except bit kmax-1

	idelmax = ibset(0,kmax-1)			! Differences larger than idelmax are not compressed

	ibitmax = nPack*nbit-1				! # of last available bit in nPack (1st bit is bit 0)
	ibit = -1					! Pack is filled already up to bit ibit

	iprev = 0					! The '0th element' of Orig

	!-------
	! Loop over all elements of input array

	do i=0,nOrig-1
	    if (ibit .eq. ibitmax) then			! Pack is full: bail with return value 0
		ice_pack = 0
		return
	    end if

	    idel    = Orig(i)-iprev			! Difference with previous element
	    idelabs = abs(idel)				! Absolute difference

	    !-------
	    ! Test how to encode

	    if (idel .eq. 0) then			! Zero difference

		!-------
		! Zero difference is a special case. No zero bits and no sign bit are
		! needed, so itop+1 = 0

		itop = -1
		ipad =  0

	    else if (idelabs .lt. idelmax) then		! YES: small enough difference
		ienc = idelabs

		!--------
		! Find the number of bits needed to describe the difference. There
		! will be itop+1 bits, including a sign bit. Store with itop+1 zeroes,
		! a one, then the itop+1 bits.
		!
		! Here 1 <= idelabs < idelmax, i.e. at least one bit of idelabs is set,
		! and bits [kmax-1,nbit-1] are zero. Only kmax-1 bits [0,kmax-2] can be
		! non-zero. Find the highest non-zero bit, itop <= kmax-2. Since this
		! bit is always set we don't need to store it. Only the preceding itop bits
		! [0,itop-1] need to be stored. (except for an additional sign bit, see below).

		itop = kmax-2
		do while (.not. btest(ienc, itop))
		    itop = itop-1
		end do

		ipad = itop+1				! ipad > 0

		!--------
		! Copy the sign bit (bit nbit-1) of the non-zero difference idel to
		! bit itop of ienc. Since we know that bit itop of ienc is set (it's
		! the highest non-zero bit), we only need to clear it if necessary.
		! If ipad <= kshift don't move the sign (it is stored separately later).

		if (ipad .gt. kshift .and. .not. btest(idel,nbit1)) ienc = ibclr(ienc,itop)

	    else					! NO Rice Encoding: difference too big

		!--------
		! Mark Pack with kmax zeroes followed by a one, and the full word
		! (i.e. NOT the difference) for a total of kmax+1+nlenbit bits.
		! We always encode the absolute value of Orig, and if required (i.e. if nsign = 0)
		! add a bit for the sign.

		ipad = kmax
		itop = itopbit
		ienc = abs(Orig(i))			! Pick up array value

		!-------
		! If nsign = 0 (both positive and negative values in Orig) then we need to
		! encode the sign in bit itopbit.

		if (nsign .eq. 0) then
		    if (btest(Orig(i),nbit1)) ienc = ibset(ienc,itop)
		end if

	    end if

	    ibit = ibit+nrank(max(0, ipad-kshift))	! Skip ipad-kshift zero bits

	    ibit  = ibit+1				! Set bit ipad+1 to 1
	    iPack = ibit/nbit
	    Pack(iPack) = ibset(Pack(iPack), mod(ibit,nbit))

	    !------
	    ! If kshift > 0 an extra bit is needed to identify zero differences from differences
	    ! described by ipad <= kshift bits. For a non-zero differences an additional bit is
	    ! needed to store the sign.

	    if (kshift .gt. 0 .and. ipad .le. kshift) then
		ibit  = ibit+1
		iPack = ibit/nbit
		if (ipad .ge. 1) then
		    Pack(iPack) = ibset(Pack(iPack), mod(ibit,nbit))
		    ibit  = ibit+1
		    iPack = ibit/nbit
		    if (btest(idel,nbit1)) Pack(iPack) = ibset(Pack(iPack), mod(ibit,nbit))
		end if
	    end if

	    !--------
	    ! Now transfer most-significant itop+1 bits (bits [0,itop]) of ienc to Pack
	    ! Remember that itop=-1 for zero differences. Move at least kshift bits unless itop=-1.
	    ! If itop=-1 then even for kshift > 0 it is not necessay to write kshift zeroes, since the
	    ! preceeding 0-bit already indicates this condition. 

	    if (itop .ne. -1) then			! Required for kshift > 0
		do k=0,max(kshift-1,itop)
		    ibit = ibit+1
		    iPack = ibit/nbit
		    if (btest(ienc,k)) Pack(iPack) = ibset(Pack(iPack), mod(ibit,nbit))
		end do
	    end if

	    iprev = Orig(i)

	end do

	iPack = iPack+1					! # elements in Pack used
	ice_pack = ibit+1				! # bits used

	return
	end