subroutine rfits(filename,a,nx,ny,nz,datamin,datamax)
c given FITS filename, returns REAL a(nx,ny,nz) & nx,ny,nz + datamin & max
c -lcfitsio
	integer SIZE,status,unit,readwrite,blocksize,naxes(3),group,firstpix
	parameter(SIZE=1200000)
	real nullval,buffer(SIZE),a(*)
	logical anynull
	character*(*) filename
c=============================================================================
	status = 0
	call ftgiou(unit,status)
	readwrite = 0
	call ftopen(unit,filename,readwrite,blocksize,status)
	call ftgknj(unit,'NAXIS',1,3,naxes,nfound,status)
	nx = 0
	ny = 0
	nz = 0
	datamax = -1e30
	datamin = - datamax
	if ( nfound < 0 ) then
		return
	else if ( nfound == 1 ) then
		nx = naxes(1)
	else if ( nfound == 2 ) then
		nx = naxes(1)
		ny = naxes(2)
	else
		nx = naxes(1)
		ny = naxes(2)
		nz = naxes(3)
	end if
	npixels = nx*max(1,ny)*max(1,nz)
	group = 1
	firstpix = 1
	nullval = -999

	ind = 0
	do while ( npixels > 0 ) 
		nbuffer = min(SIZE,nx,npixels)
		call ftgpve(unit,group,firstpix,nbuffer,nullval,buffer,
     &      anynull,status)
		do i=1,nbuffer
			datamin = min(datamin,buffer(i))
			datamax = max(datamax,buffer(i))
			if ( buffer(i) == 255. ) buffer(i) = 0.
			a(ind+i) = buffer(i)
		end do
		ind = ind+nbuffer
		npixels = npixels - nbuffer
		firstpix = firstpix+nbuffer
	end do
	call ftclos(unit,status)
	return
	end