function  mrd_dofn, name, index, use_colnum 
 
; Convert the string name to a valid variable name.  If name 
; is not defined generate the string Cnn when nn is the index 
; number. 

on_error, 2

table = 0
sz = size(name) 
nsz = n_elements(sz) 
if not use_colnum and (sz(nsz-2) ne 0) then begin 
        if sz(nsz-2) eq 7 then begin 
                str = name(0) 
        endif else begin 
                str = 'C'+string(index) 
        endelse 
endif else begin 
        str = 'C'+string(index) 
endelse 
 
str = strcompress(str, /remove_all) 
len = strlen(str) 
c = strmid(str, 0, 1) 
 
; First character must be alphabetic. 
 
if  not (('a' le c  and 'z' ge c) or ('A' le c  and 'Z' ge c)) then begin 
        str = 'X' + str 
endif 
 
; 
; Replace invalid characters with underscores. 
; We assume ASCII ordering.
;
for i=1, len-1 do begin 
        c = strmid(str, i, 1) 
        if not (('a' le c  and 'z' ge c)  or  $ 
                ('A' le c  and 'Z' ge c)  or  $ 
                ('0' le c  and '9' ge c)  or  $ 
                (c eq '_') or (c eq '$')      $ 
               ) then strput, str, '_', i 
endfor 
 
return, str 
end 

;***************************************************************
pro mrd_doff, form, dim, type 
 
; 
; Parse the TFORM keyword and return the type and dimension of the 
; data. 

on_error, 2

; Find the first non-numeric character. 
 
len = strlen(form) 
 
if len le 0 then return 
 
for i=0, len-1 do begin 
 
        c = strmid(form, i, 1) 
        if c lt '0'  or c gt '9' then goto, not_number 
         
endfor 
 
not_number: 
 
if i ge len-1 then return 
 
if i gt 0 then begin 
        dim = long(strmid(form, 0, i)) 
endif else begin 
        dim = 0 
endelse 
 
type = strmid(form, i, 1) 
end 

;*********************************************************************

function mrd_chkfn, name, namelist, index 
 
; 
;  Check that this name is unique with regard to other column names.
; 
 
if strlen(name) gt 15 then name = strmid(name, 0, 15) 
w = where(name eq strmid(namelist, 0, 15) )
if w(0) ne -1 then begin
        ; We have found a name conflict. 
        ; 
                name = 'gen$name_'+strcompress(string(index+1),/remove_all) 
endif 
 
return, name 
 
end

;=====================================================================
; END OF UTILITY FUNCTIONS ===========================================
;=====================================================================
pro mrd_atype, form, type, slen

;
; Parse the TFORM keyword and return the type and dimension of the
; data.

on_error, 2
; Find the first non-numeric character.


; Get rid of blanks.
form = strcompress(form,/remove_all)
len = strlen(form)
if len le 0 then return

type = strmid(form, 0,1)
length = strmid(form,1,len-1)
;
; Ignore the number of decimal places.  We assume that there
; is a decimal point.
;
p = strpos(length, '.')
if p gt 0 then length = strmid(length,0,p)

if strlen(length) gt 0 then slen = fix(length) else slen = 1

end

pro mrd_read_ascii, unit, range, nbytes, nrows, nfld, typarr, posarr,   $
     lenarr, nullarr, table

; Read in the table information.
;
; Unit          Unit to read data from.
; Nbytes        Number of bytes per row.
; Nrows         Number of rows.
; Nfld          Number of fields in structure.
; Typarr        Array indicating type of variable.
; Posarr        Starting position of fields (first char at 0)
; Lenarr        Length of fields
; Nullarr       Array of null values
; Table         Table to read information into.
;


bigstr = bytarr(nbytes, range(1)-range(0)+1)

if range(0) gt 0 then mrd_skip, unit, nbytes*range(0)

readu,unit, bigstr
s1 = posarr-1
s2 = s1 + lenarr - 1
for i=0, nfld-1 do begin

        flds = strtrim( bigstr(s1(i):s2(i),* ) )
	if strtrim(nullarr(i)) ne '' then begin
	        w = where(flds ne strtrim(nullarr(i)))
	        if w(0) ne -1 then begin
	                if N_elements(w) EQ 1 then w = w(0)
                        if typarr(i) eq 'I' then begin
				table(w).(i) = long(flds(w))
			endif else if typarr(i) eq 'E' or typarr(i) eq 'F' then begin
				table(w).(i) = float(flds(w))
			endif else if typarr(i) eq 'D' then begin
				table(w).(i) = double(flds(w))
			endif else if typarr(i) eq 'A' then begin
				table(w).(i) = flds(w)
			endif
		endif
	endif else begin
		if typarr(i) eq 'I' then begin
			table.(i) = long(flds)
		endif else if typarr(i) eq 'E' or typarr(i) eq 'F' then begin
			table.(i) = float(flds)
		endif else if typarr(i) eq 'D' then begin
			table.(i) = double(flds)
		endif else if typarr(i) eq 'A' then begin
			table.(i) = flds
		endif
	endelse
endfor

end

pro mrd_ascii, header, structyp, use_colnum,   $
    range, table, $
    nbytes, nrows, nfld, typarr, posarr, lenarr, nullarr, $
    fnames, fvalues, scales, offsets, scaling, status, $
    silent=silent, tempdir=tempdir, columns=columns

; Define a structure to hold a FITS binary table.
;
; Header                FITS header for table.
; Structyp              IDL structure type to be used for
;                       structure.
; Use_colnum            Use column numbers not names.
; Range                 Range of rows of interest
; Table                 Structure to be defined.
; Nbytes                Bytes per row
; Nfld                  Number of fields
; Typarr                Array of field types
; Posarr                Array of field offsets
; Lenarr                Array of field lengths
; Nullarr               Array of field null values
; Status                Return status.

on_error, 2

types = ['I', 'E', 'F', 'D', 'A']
sclstr = ['0l', '0.0', '0.0', '0.0d0', ' ']
status = 0

if fxpar(header, 'XTENSION') ne 'TABLE   ' then begin
	print, 'MRDFITS: Header is not from ASCII table.'
	status = -1;
	return
endif

nfld = fxpar(header, 'TFIELDS')
nrows = fxpar(header, 'NAXIS2')
nbytes = fxpar(header, 'NAXIS1')

if range(0) ge 0 then begin
	range(0) = range(0) < (nrows-1)
	range(1) = range(1) < (nrows-1)
endif else begin
	range(0) = 0
	range(1) = nrows-1
endelse

nrows = range(1) - range(0) + 1

if nrows le 0 then begin
	if not keyword_set(silent) then begin
		print,'MRDFITS: ASCII table.  ',strcompress(string(nfld)),  $
		 ' columns, no rows'
	endif
	return
endif

;
;  Loop over the columns

typarr = strarr(nfld)
lenarr = intarr(nfld)
posarr = intarr(nfld)
nullarr = strarr(nfld)
fnames = strarr(nfld)
fvalues = strarr(nfld)
scales = dblarr(nfld)
offsets = dblarr(nfld)

for i=0, nfld-1 do begin
	suffix = strcompress(string(i+1), /remove_all)
	fname = fxpar(header, 'TTYPE' + suffix)
	fform = fxpar(header, 'TFORM' + suffix)
	fpos = fxpar(header, 'TBCOL' + suffix)
	fnull = fxpar(header, 'TNULL' + suffix)
	scales(i) = fxpar(header, 'TSCAL' + suffix)
	if scales(i) eq 0.0d0 then scales(i) = 1.0d0
	offsets(i) = fxpar(header, 'TZERO'+suffix)
	
	fname = mrd_dofn(fname,i+1, use_colnum)
	fnames(i) = fname

	fname = mrd_chkfn(fname, fnames, i)
	
	mrd_atype, fform, ftype, flen
	typarr(i) = ftype
	lenarr(i) = flen
	posarr(i) = fpos
	nullarr(i) = fnull
	
	for j=0, n_elements(types) - 1 do begin
		if ftype eq types(j) then begin
			if ftype ne 'A' then begin
				val = sclstr(j)
			endif else begin
				val = 'string(replicate(32b,'+string(flen)+'))'
			endelse
			
			fvalues(i) = val
			
			goto, next_col
		endif
	endfor
	
	print, 'MRDFITS: Invalid format code:',ftype, ' for column ', i+1
	status = -1
	return
next_col:
endfor

if scaling then begin
	w = where(scales ne 1.0d0 or offsets ne 0.0d0)
	if w(0) eq -1 then scaling = 0
endif
if not scaling and not keyword_set(columns) then begin
	table = mrd_struct(fnames, fvalues, nrows, structyp=structyp, $
	  tempdir=tempdir)
endif else begin
	table = mrd_struct(fnames, fvalues, nrows, tempdir=tempdir)
endelse

if not keyword_set(silent) then begin
	print,'MRDFITS: ASCII table.  ',strcompress(string(nfld)),  $
	 ' columns by ',strcompress(string(nrows)), ' rows.
endif
status = 0
return

end

pro  mrd_columns, table, columns, fnames, fvalues, $
    vcls, vtpes, scales,  offsets, scaling,        $
    structyp=structyp, tempdir=tempdir

on_error, 2
; Eliminate columns from the table that do not match the
; user specification.


sz = size(columns)

type = sz(sz(0)+1)
nele = sz(sz(0)+2)
if type eq 8 or type eq 6 or type eq 0 then return  ; Can't use structs
						    ; or complex.

if type eq 4 or type eq 5 then tcols = fix(columns)
if type eq 1 or type eq 2 or type eq 3 then tcols = columns

; Convert strings to uppercase and compare with column names.

if type eq 7 then begin
	for i=0, nele-1 do begin
		cname = strupcase(columns(i))
		w = where(cname eq strupcase(fnames))
		if w(0) ne -1 then begin
			if n_elements(tcols) eq 0 then begin
				tcols = w(0)+1
			endif else begin
				tcols = [tcols, w(0)+1]
			endelse
		endif
	endfor
endif

; Subtract one from column indices and check that all indices >= 0.
if n_elements(tcols) gt 0 then begin
	tcols = tcols-1
	w = where(tcols ge 0)
	if w(0) eq -1 then begin
		dummy = temporary(tcols)
	endif
endif

if n_elements(tcols) le 0 then begin
	print, 'MRDFITS:  No columns match'
	
	; Undefine variables.  First ensure they are defined, then
	; use temporary() to undefine them.
	table = 0
	fnames = 0
	fvalues = 0
	vcls = 0
	vtpes = 0
	scales = 0
	offsets = 0
	dummy = temporary(fnames)
	dummy = temporary(fvalues)
	dummy = temporary(vcls)
	dummy = temporary(vtpes)
	dummy = temporary(scales)
	dummy = temporary(offsets)
	scaling = 0
	
endif else begin

	; Replace arrays with only desired columns.
	
	fnames = fnames(tcols)
	fvalues = fvalues(tcols)
	
	; Check if there are still variable length columns.
	if n_elements(vcls) gt 0 then begin
		vcls = vcls(tcols)
		vtpes = vtpes(tcols)
		w = where(vcls eq 1)
		if w(0) eq -1 then begin
			dummy = temporary(vcls)
			dummy = temporary(vtpes)
		endif
	endif
	
	; Check if there are still columns that need scaling.
	if n_elements(scales) gt 0 then begin
		scales = scales(tcols)
		offsets = offsets(tcols)
		w = where(scales ne 1.0d0  or offsets ne 0.0d0)
		if w(0) eq -1 then scaling = 0
	endif
	

	ndim = n_elements(table)
	
	if scaling or n_elements(vcls) gt 0 then begin
		tabx = mrd_struct(fnames, fvalues, ndim, tempdir=tempdir)
	endif else begin
		tabx = mrd_struct(fnames, fvalues, ndim, structyp=structyp, $
		  tempdir=tempdir)
	endelse
	
	for i=0, n_elements(tcols)-1 do begin
		tabx.(i) = table.(tcols(i))
	endfor

	table = tabx
endelse
end

pro mrd_read_image, unit, range, rsize, table 
 
; Read in the table information. 
; 
; Unit          Unit to read data from. 
; Table         Table/array to read information into. 
; 

on_error, 0

; If necessary skip to beginning of desired data. 


if range(0) gt 0 then mrd_skip, unit, range(0)*rsize
 
if rsize eq 0 then return 
 
readu, unit, table 
ieee_to_host, table 
end 
 
pro mrd_image, header, range, rsize, table, scales, offsets, scaling, $
  status, silent=silent, tempdir=tempdir
 
; Define structure/array to hold a FITS image. 
; 
; Header                FITS header for table. 
; Range                 Range of data to be retrieved. 
; Rsize                 Size of a row or group. 
; Table                 Structure to be defined. 
; Status                ; Silent=silent		Suppress info messages?
 
on_error, 2

lens = [1,2,4,4,8] 
typstrs=['Byte', 'Int*2', 'Int*4', 'Real*4', 'Real*8']
 
status = 0 
 
 
naxis = fxpar(header, 'NAXIS') 
bitpix= fxpar(header, 'BITPIX') 
if naxis gt 0 then dims = fxpar(header, 'NAXIS*') else dims = 0 
gcount = fxpar(header, 'GCOUNT') 
pcount = fxpar(header, 'PCOUNT')
gcount = long(gcount)
 
if bitpix eq 8 then type = 0            $ 
else if bitpix eq 16 then type = 1      $ 
else if bitpix eq 32 then type = 2      $ 
else if bitpix eq -32 then type = 3     $ 
else if bitpix eq -64 then type = 4 
 
; Note that for random groups data we must ignore the first NAXISn keyword. 
if pcount gt 1 then begin 
 
        if range(0) eq -1 then begin 
                range(0) = 0 
                range(1) = gcount 
        endif 
         
        naxis = naxis - 1 
        if naxis le 0 then begin 
                rsize = 0 
                table = 0 
                return 
        endif 
         
        scales = dblarr(pcount+1)
        offsets = dblarr(pcount+1)
        dims = dims(1:naxis) 
        values = strarr(2)         
        typarr=["bytarr", "intarr", "lonarr", "fltarr", "dblarr"] 
        values(0) = typarr(type) + "("+string(pcount)+")" 
        rsize = dims(0) 
        sarr = "(" + strcompress(string(dims(0)))
         
        for i=1, naxis-1 do begin 
                sarr = sarr + "," + strcompress(string(dims(i)),/remo)
                rsize = rsize*dims(i) 
        endfor 
         
        sarr = sarr + ")"

	if not keyword_set(silent) then print,'MRDFITS--Image with groups:', $
	  ' Groups=',strcompress(string(gcount)),' Npar=', 		     $
	  strcompress(string(pcount)), ' Array=', sarr, '  Type=',typstrs(type)

	sarr = typarr(type) + sarr
        values(1) = sarr 
        rsize = (rsize + pcount)*lens(type) 
         
        table = mrd_struct(['params','array'], values, range(1)-range(0)+1, $
	  tempdir=tempdir)
	   
    
	
        for i=0, pcount-1 do begin
                istr = strcompress(string(i),/remo)
                scales(i) = fxpar(header, 'PSCAL'+istr)
                if scales(i) eq 0.0d0 then scales(i) =1.0d0
                offsets(i) = fxpar(header, 'PZERO'+istr)
                scales(pcount) = fxpar(header, 'BSCALE')
                if scales(pcount) eq 0.0d0 then scales(pcount) = 1.0d0
                offsets(pcount) = fxpar(header, 'BZERO')
        endfor
        if scaling then begin
                w = where(scales ne 1.0d0 or offsets ne 0.0d0)
                if w(0) eq -1 then scaling = 0
        endif
endif else begin 
 
        if naxis eq 0 then begin 
                rsize = 0 
                table = 0
		if not keyword_set(silent) then begin
			print, 'MRDFITS: Null image, NAXIS=0'
		endif
                return 
        endif 
         
        if gcount gt 1 then begin 
                dims = [dims, gcount] 
                naxis = naxis + 1 
        endif 
         
         
        for i=naxis-1,1,-1 do begin 
                if dims(i) eq 1 then begin
			if not keyword_set(silent) then begin
				print, 'MRDFITS: Truncating unused dimensions'
			endif
                        dims = dims(0:i-1) 
                        naxis = naxis - 1 
                endif else begin 
                        goto, no_trunc 
                endelse 
        endfor 
 
no_trunc:
	if not keyword_set(silent) then begin
		str = '('
		for i=0, naxis-1 do begin
			if i ne 0 then str = str + ','
			str = str + strcompress(string(dims(i)),/remo)
		endfor
		str = str+')'
		print, 'MRDFITS: Image array ',str, '  Type=', typstrs(type)
	endif
		
        maxd = dims(naxis-1) 
         
        if range(0) ne -1 then begin 
                range(0) = range(0)<(maxd-1) 
                range(1) = range(1)<(maxd-1) 
        endif else begin 
                range(0) = 0 
                range(1) = maxd - 1 
        endelse 
         
        dims(naxis-1) = range(1)-range(0)+1 
         
        rsize = 1 
        if naxis gt 1 then for i=0, naxis - 2 do rsize=rsize*dims(i) 
        rsize = rsize*lens(type) 
         
        sz = lonarr(naxis+3) 
        sz(0) = naxis 
        sz(1:naxis) = dims 
        nele = 1l 
         
        for i=0, naxis-1 do begin 
                nele = nele*dims(i) 
        endfor 
         
        sz(naxis+1) = type + 1   
        sz(naxis+2) = nele 
         
        if nele gt 0 then  begin 
                table = make_array(size=sz) 
        endif else begin 
                table = 0 
        endelse 
        
        scales = dblarr(1)
        offsets = dblarr(1)
        scales(0) = fxpar(header, 'BSCALE')
        offsets(0) = fxpar(header, 'BZERO')
        if scales(0) eq 0.0d0 then scales(0) = 1.0d0
        if scaling and scales(0) eq 1.0d0 and offsets(0) eq 0.0d0 then $
            scaling = 0
         
endelse 
         
status = 0 
return 
 
end 

pro mrd_scale, type, scales, offsets, table, header,  $
               fnames, fvalues, nrec, dscale = dscale, structyp=structyp, $
	       tempdir=tempdir
;
; Scale a FITS array or table.
;
; Type:         FITS file type, 0=image/primary array
;                               1=ASCII table
;                               2=Binary table
;
; scales:       An array of scaling info
; offsets:      An array of offset information
; table:        The FITS data.
; header:       The FITS header.
; dscale:       Should data be scaled to R*8?
; fnames:       Names of table columns.
; fvalues:      Values of table columns.
; nrec:         Number of records used.
; structyp:     Structure name.
;
; Modified: 1-Aug-1995 to fix typo.

on_error, 2

w = where(scales ne 1.d0  or offsets ne 0.d0)
if w(0) eq -1 then return
ww = where(scales eq 1.d0 and offsets eq 0.d0)

; First do ASCII and Binary tables.
if type ne 0 then begin
	
	if type eq 1 then begin
		if keyword_set(dscale) then begin
        	        fvalues(w) = '0.0d0'
        	endif else begin
        	        fvalues(w) = '0.0'
        	endelse
	endif else if type eq 2 then begin

	        if keyword_set(dscale) then begin
        	        sclr = '0.d0'
        	        vc = 'dblarr'
        	endif else begin
        	        sclr = '0.0'
        	        vc = 'fltarr'
		endelse
		
        	for i=0, n_elements(w)-1 do begin
                	col = w(i)
	                sz = size(table(0).(col))
        	        if sz(0) eq 0 then begin
                	        fvalues(col) = sclr
	                endif else begin
        	                str = vc + '('
                	        for j=0, sz(0)-1 do begin
                        	        if j ne 0 then str = str + ','
                                	str = str + string(sz(j+1))
	                        endfor
        	                str = str + ')'
                	        fvalues(col) = sclr
                	endelse
        	endfor
        endif

        tabx = mrd_struct(fnames, fvalues, nrec, structyp=structyp, $
	  tempdir=tempdir)

        if ww(0) ne -1 then begin
		for i=0, n_elements(ww)-1 do begin
	                tabx.(ww(i)) = table.(ww(i))
		endfor
        endif
	
        for i=0, n_elements(w)-1 do begin
	        tabx.(w(i)) = table.(w(i))*scales(w(i)) + offsets(w(i))
                istr = strcompress(string(w(i)+1), /remo)
                fxaddpar, header, 'TSCAL'+istr, 1.0, 'Set by MRD_SCALE'
                fxaddpar, header, 'TZERO'+istr, 0.0, 'Set by MRD_SCALE'
        endfor

	table = temporary(tabx)

endif else begin
; Now process images and random groups.

        sz = size(table(0))
        if sz(sz(0)+1) ne 8 then begin
                ; Not a structure so we just have an array of data.
		if keyword_set(dscale) then begin
	                table = table*scales(0)+offsets(0)
		endif else begin
			table = table*float(scales(0)) + float(offsets(0))
		endelse
                fxaddpar, header, 'BZERO', 1.0, 'Set by MRD_SCALE'
                fxaddpar, header, 'BSCALE', 0.0, 'Set by MRD_SCALE'

        endif else begin
                nparam = n_elements(sz(0).(0))
                if keyword_set(dscale) then typ = 'dbl' else typ='flt'
                s1 = typ+'arr('+string(nparam)+')'
                ngr = n_elements(table)
                sz = size(table(0).(1))
                if sz(0) eq 0 then dims = [1] else dims=sz(1:sz(0))
                s2 = typ + 'arr('
                for i=0, n_elements(dims)-1 do begin 
                        if i ne 0 then s2 = s2+ ','
                        s2 = s2+string(dims(i))
                endfor
                s2 = s2+')'
                tabx = mrd_struct(['params', 'array'],[s1,s2],ngr, $
		  tempdir=tempdir)
                for i=0, nparam-1 do begin
                        istr = strcompress(string(i),/remo)
                        fxaddpar, header, 'PSCAL'+istr, 1.0, 'Added by MRD_SCALE'
                        fxaddpar, header, 'PZERO'+istr, 0.0, 'Added by MRD_SCALE'
                        tabx.(0)(i) = table.(0)(i)*scales(i)+offsets(i)
                endfor
                tabx.(1) = table.(1)*scales(nparam) + offsets(nparam)
                fxaddpar, header, 'BSCALE', 1.0, 'Added by MRD_SCALE'
                fxaddpar, header, 'BZERO', 0.0, 'Added by MRD_SCALE'
                table = temporary(tabx)
	endelse
endelse

end

pro mrd_read_heap, unit, header, range, fnames, fvalues, vcls, vtpes, table, $ 
   structyp, scaling, scales, offsets, status, silent=silent,         $
   tempdir=tempdir, columns=columns

on_error,0
; This program reads the heap area to get the actual values of variable 
; length arrays. 
; 
; Unit:         FITS unit number. 
; header:       FITS header. 
; fnames:       Column names. 
; fvalues:      Column values. 
; vcols:        Column numbers of variable length columns. 
; vtypes:       Actual types of variable length columns 
; table:        Table of data from standard data area, on output 
;               contains the variable length data. 
; structyp:     Structure name. 
; scaling:      Is there going to be scaling of the data?
; status:       Set to -1 if an error occurs. 
;
; Modified 17-Feb-1995 by TAM to fix bug when there were
; multiple actual varying columns.
; Modified 2-June-1992 by TAM to fix bug where there are multiple actual
; varying columns and columns deleted.

typstr = 'LXBIJAEDCM' 
prefix = ['bytarr(', 'bytarr(', 'bytarr(', 'intarr(',     $ 
          'lonarr(', 'string(bytarr(', 'fltarr(',         $ 
          'dblarr(', 'cmplxarr(', 'dblarr(2,'] 
status = 0 

; Convert from a list of indicators of whether a column is variable
; length to pointers to only the variable columns.

vcols = where(vcls eq 1)
vtypes = vtpes(vcols)

nv = n_elements(vcols) 
 
; Find the beginning of the heap area. 
 
heapoff = fxpar(header, 'THEAP') 
sz = fxpar(header, 'NAXIS1')*fxpar(header, 'NAXIS2') 
if heapoff ne 0 and heapoff lt sz then begin 
        print, 'MRDFITS: Heap begins within data area' 
        status = -1 
        return 
endif 
 
; Get the size of the heap. 
pc = fxpar(header, 'PCOUNT') 
if heapoff eq 0 then heapoff = sz 
hpsiz = pc - (heapoff-sz) 
 
if (hpsiz gt 0) then heap = bytarr(hpsiz) 
 
; Move to the beginning of the heap -- we are assuming that we have
; finished reading the data but we may have only read some rows.

skip_dist = fxpar(header, 'NAXIS1')
skip_dist = skip_dist*(fxpar(header, 'NAXIS2') - (range(1)+1))

mrd_skip, unit, skip_dist + heapoff-sz 
 
; Read in the heap 
readu, unit, heap 
 
; Find the maximum dimensions of the arrays. 
; 
; Note that the variable length column currently has fields which 
; are I*4 2-element arrays where the first element is the 
; length of the field on the current row and the second is the 
; offset into the heap. 
 
vdims = lonarr(nv) 
for i=0, nv-1 do begin 
    col = vcols(i) 
    vdims(i) = max(table.(col)(0))
    w = where(table.(col)(0,*) ne vdims(i))
    if w(0) ne -1 then begin
    	if n_elements(lencols) eq 0 then begin
		lencols = [col]
	endif else begin
		lencols=[lencols,col]
	endelse
    endif

    if vtypes(i) eq 'X' then vdims(i)=(vdims(i)+7)/8 
    ind = strpos(typstr, vtypes(i))
    fvalues(col) = prefix(ind) + string(vdims(i)) + ')'     
    if vtypes(i) eq 'A' then fvalues(col) = fvalues(col) + ')' 
    
endfor 
 
nfld = n_elements(fnames) 
 
; Get rid of columns which have no actual data. 
w= intarr(nfld) 
w(*) = 1
corres = indgen(nfld)
ww = where(vdims eq 0) 
if ww(0) ne -1 then  begin
	w(vcols(ww)) = 0
	if not keyword_set(silent) then begin
		print, 'MRDFITS: ', strcompress(string(n_elements(ww))),  $
		  ' unused variable length columns deleted'
	endif
endif

; Get rid of unused columns.
corres = corres(where(w gt 0))
fnames = fnames(where(w gt 0)) 
fvalues = fvalues(where(w gt 0))
scales = scales(where(w gt 0))
offsets = offsets(where(w gt 0))

vcols = vcols(where(vdims gt 0))
vtypes = vtypes(where(vdims gt 0))
vdims = vdims (where(vdims gt 0))

; Now add columns for lengths of truly variable length records.
if n_elements(lencols) gt 0 then begin
	if not keyword_set(silent) then begin
		print, 'MRDFITS: ', strcompress(string(n_elements(lencols))), $
		  ' length column[s] added'
	endif
	

	for i=0, n_elements(lencols)-1 do begin
		col = lencols(i)
		w = where(col eq corres)
		ww = where(col eq vcols)
		w = w(0)
		ww = ww(0)
		fvstr = '0l'
		fnstr = 'L'+strcompress(string(col),/remo)+'_'+fnames(w)
		nf = n_elements(fnames)
		
		; Note that lencols and col refer to the index of the
		; column before we started adding in the length
		; columns.
		
		if w eq nf-1 then begin
			; Subtract -1 for the length columns so 0 -> -1 and
			; we can distinguish this column.

			corres = [corres, -col-1]
			fnames = [fnames, fnstr]
			fvalues = [fvalues, fvstr]
			scales = [scales, 1.0d0]
			offsets = [offsets, 0.0d0]
		endif else begin
			corres = [corres(0:w),-col-1,corres(w+1:nf-1)]
			fnames = [fnames(0:w),fnstr,fnames(w+1:nf-1)]
			fvalues = [fvalues(0:w),fvstr,fvalues(w+1:nf-1)]
			scales = [scales(0:w), 1.0d0, scales(w+1:nf-1)]
			offsets = [offsets(0:w),0.0d0, offsets(w+1:nf-1)]
 		endelse
	endfor
endif

; Generate a new table with the appropriate structure definitions 
if not scaling and not keyword_set(columns) then begin
  tablex = mrd_struct(fnames, fvalues, n_elements(table), structyp=structyp, $
     tempdir=tempdir)
endif else begin
  tablex = mrd_struct(fnames, fvalues, n_elements(table), tempdir=tempdir)
endelse
 
; I loops over the new table columns, col loops over the old table.
; When col is negative, it is a length column.
for i=0, n_tags(tablex(0))-1 do begin
	
	col = corres(i)
		
        if col ge 0 then begin
		w = where(vcols eq col)
		
		; First handle the case of a column that is not
		; variable length -- just copy the column.
		
		if w(0) eq -1 then begin
 	                tablex.(i) = table.(col)
			
		 
	        endif else begin 
			vc = w(0)
			; Now handle the variable length columns
		        siz = table.(col)(0) 
		        off = table.(col)(1) 
		        w = where(siz gt 0) 
			
			; Don't process rows where the length is 0.
		        off = off(w) 
		        siz = siz(w) 
		        nw = n_elements(w)-1 
			
		; Now process each type.
	        case vtypes(vc) of 
	          'L':begin 
        	        tablex(w).(i)(0:siz-1) = heap(off:off+siz-1) 
        	      end 
	          'X':begin 
        	        siz = 1+(siz-1)/8 
 	                tablex(w).(i)(0:siz-1) = heap(off:off+siz-1) 
	              end 
          	  'B':begin 
              	  	tablex(w).(i)(0:siz-1) = heap(off:off+siz-1) 
              	      end 
          	  'I':begin 
                	for j=0, nw do begin 
                    	 tablex(w(j)).(i) = fix(heap, off(j), siz(j)) 
                	endfor 
                	fld = tablex.(i) 
	                byteorder, fld, /ntohs 
        	        tablex.(i) = fld 
              	      end 
          	  'J':begin 
                	for j=0, nw do begin 
                   	 tablex(w(j)).(i) = long(heap, off(j), siz(j)) 
                	endfor 
                	fld = tablex.(i) 
	                byteorder, fld, /ntohl 
        	        tablex.(i) = fld 
              	      end 
          	   'E':begin 
                	for j=0, nw do begin 
                    	  tablex(w(j)).(i) = float(heap, off(j), siz(j)) 
                	endfor 
                	fld = tablex.(i) 
	                byteorder, fld, /xdrtof 
        	        tablex.(i) = fld 
	              end 
        	   'D':begin 
                	for j=0, nw do begin 
                  	  tablex(w(j)).(i) = double(heap, off(j), siz(j)) 
                	endfor 
	                fld = tablex.(i) 
        	        ieee_to_host, fld 
                	tablex.(i) = fld 
	               end 
             	   'C':begin 
              	        for j=0, nw do begin 
                    	  tablex(w(j)).(i) = complex(heap, off(j), siz(j)) 
	                endfor 
        	        fld = tablex.(i) 
                	byteorder, fld, /xdrtof 
                	tablex.(i) = fld 
	               end 
          	   'M':begin 
                	for j=0, nw do begin 
                    	  tablex(w(j)).(i) = double(heap, off(j), 2, siz(j)) 
	                endfor 
                	fld = tablex.(i) 
	                ieee_to_host, fld 
        	        tablex.(i) = fld 
	               end 
          	   'A':begin 
              	        tablex(w).(i) = string(heap(off:off+siz)) 
              	       end 
      		endcase
	    	endelse
		
	  endif else begin
	  	; Now handle the added columns which hold the lengths
		; of the variable length columns.
		
		ncol = -col - 1 ; Remember we subtracted an extra one.
		tablex.(i) = table.(ncol)(0)
	  endelse
endfor 
 
; Finally get rid of the initial table and return the table with the 
; variable arrays read in. 
; 
table = temporary(tablex) 
return 
end 

pro mrd_read_table, unit, range, rsize, structyp, nfld, typarr, table 
 
; 
; Read in the binary table information. 
; 
; Unit          Unit to read data from. 
; Range         Desired range 
; Rsize         Size of row. 
; structyp      Structure type. 
; Nfld          Number of fields in structure. 
; Typarr        Field types 
; Table         Table to read information into. 
; 

on_error, 2

if range(0) gt 0 then mrd_skip, unit, rsize*range(0)
 
readu,unit, table 
for i=0, nfld-1 do begin 
 
        typ = typarr(i) 
        if typ eq 'B' or typ eq 'A'  or typ eq 'X' or typ eq 'L' $ 
           then goto,nxtfld 
        fld = table.(i) 
        if typ eq 'I' then byteorder, fld, /htons 
        if typ eq 'J' or typ eq 'P' then byteorder, fld, /htonl 
        if typ eq 'E' or typarr(i) eq 'C' then byteorder, fld, /xdrtof 
        if typ eq 'D' or typarr(i) eq 'M' then begin 
           ieee_to_host, fld 
        endif
	if n_elements(fld) gt 1 then begin
	        table.(i) = fld
	endif else begin
		table(0).(i) = fld(0)
	endelse
         
nxtfld: 
endfor 
 
end 
; Check the values of TDIM keywords to see that they have valid
; dimensionalities.  If the TDIM keyword is not present or valid
; then the a one-dimensional array with a size given in the TFORM
; keyword is used.

pro mrd_tdim, header, index, flen, arrstr, no_tdim=no_tdim

; HEADER	Current header array.
; Index		Index of current parameter
; flen		Len given in TFORM keyword
; arrstr	String returned to be included within paren's in definition.
; no_tdim       Disable TDIM processing

arrstr = strcompress(string(flen),/remo)

if keyword_set(no_tdim) then return

tdstr = fxpar(header, 'TDIM'+strcompress(string(index),/remo))
if tdstr eq '' then return

;
; Parse the string.  It should be of the form '(n1,n2,...nx)' where
; all of the n's are positive integers and the product equals flen.
;
tdstr = strcompress(tdstr,/remo)
len = strlen(tdstr)
if strmid(tdstr,0,1) ne '(' and strmid(tdstr,len-1,1) ne ')' or len lt 3 then begin
	print, 'MRDFITS: Error: invalid TDIM for column', index
	return
endif

; Get rid of parens.
tdstr = strmid(tdstr,1,len-2)
len = len-2

nind = 0
cnum = 0

for nchr=0, len-1 do begin
	c = strmid(tdstr,nchr, 1)
	
	if c ge '0' and c le '9' then begin
		cnum = 10*cnum + long(c)
		
	endif else if c eq ',' then begin
	
		if cnum le 0 then begin
			print,'MRDFITS: Error: invalid TDIM for column', index
			return
		endif
		
		if n_elements(numbs) eq 0 then	$
			numbs = cnum $
		else    numbs = [numbs,cnum]
		
		cnum = 0
		
       endif else begin
       
       		print,'MRDFITS: Error: invalid TDIM for column', index
		return
		
       endelse

endfor

; Handle the last number.
if cnum le 0 then begin
	print,'MRDFITS: Error: invalid TDIM for column', index
        return
endif

if n_elements(numbs) eq 0 then numbs = cnum else numbs = [numbs,cnum]

prod = 1

for i=0, n_elements(numbs)-1 do prod = prod*numbs(i)

if prod ne flen then begin
	print,'MRDFITS: Error: TDIM/TFORM dimension mismatch'
	return
endif

arrstr = tdstr
end
 
pro mrd_table, header, structyp, use_colnum,           $ 
    range, rsize, table, nfld, typarr, fnames, fvalues,   $ 
    vcls, vtpes, scales, offsets, scaling, status, $
    silent=silent, tempdir=tempdir, columns=columns, no_tdim=no_tdim
 
; Define a structure to hold a FITS binary table. 
; 
; Header                FITS header for table. 
; Structyp              IDL structure type to be used for 
;                       structure. 
; N_call                Number of times this routine has been called. 
; Table                 Structure to be defined. 
; Status                Return status.
; No_tdim               Disable TDIM processing.
 
on_error, 2

types = ['L', 'X', 'B', 'I', 'J', 'A', 'E', 'D', 'C', 'M', 'P'] 
arrstr = ['bytarr(', 'bytarr(', 'bytarr(', 'intarr(', 'lonarr(',      $ 
          'string(replicate(32b,', 'fltarr(', 'dblarr(', 'complexarr(',    $ 
          'dblarr(2,', 'lonarr(2*'] 
sclstr = ["'T'", '0B', '0B', '0', '0L', '" "', '0.', '0.d0', 'complex(0.,0.)', $ 
          '[0.d0,0.d0]', 'lonarr(2)'] 
 
status = 0 

xten = fxpar(header, 'XTENSION')
if xten ne 'BINTABLE' and xten ne 'A3DTABLE' then begin 
        print, 'MRDFITS: Header is not from binary table.' 
        status = -1; 
        return 
endif 
 
nfld = fxpar(header, 'TFIELDS') 
nrow = fxpar(header, 'NAXIS2')
 
if range(0) ge 0 then begin 
        range(0) = range(0) < (nrow-1) 
        range(1) = range(1) < (nrow-1) 
endif else begin 
        range(0) = 0 
        range(1) = nrow - 1 
endelse 
nrow = range(1) - range(0) + 1 

if nrow le 0 then begin
	if not keyword_set(silent) then begin
		print, 'MRDFITS: Binary table. ', $
		 strcompress(string(nfld)), ' columns, no rows.'
	endif
	return
endif

rsize = fxpar(header, 'NAXIS1') 
 
; 
;  Loop over the columns           
 
typarr = strarr(nfld) 
fnames = strarr(nfld) 
fvalues = strarr(nfld) 
dimfld = strarr(nfld)
scales = dblarr(nfld)
offsets = dblarr(nfld)
vcls = intarr(nfld)
vtpes = strarr(nfld)
 
for i=0, nfld-1 do begin 
        istr = strcompress(string(i+1), /remo)
        fname = fxpar(header, 'TTYPE' +  istr) 
        fform = fxpar(header, 'TFORM' +   istr)
        scales(i) = fxpar(header, 'TSCAL'+istr)
        if scales(i) eq 0.d0 then scales(i) = 1.d0
        offsets(i) = fxpar(header, 'TZERO'+istr)
        fname = mrd_dofn(fname,i+1, use_colnum) 
        fname = mrd_chkfn(fname, fnames, i) 
        fnames(i) = fname 
        mrd_doff, fform, dim, ftype
	
	; Treat arrays of length 1 as scalars.
	
	if dim eq 1 then begin
		dim = 0
	endif else begin
		mrd_tdim, header, i+1, dim, str, no_tdim=no_tdim
		dimfld(i) = str
	endelse
		
        typarr(i) = ftype 
         
        ; Find the number of bytes in a bit array. 
 
        if ftype eq 'X' and dim gt 0 then dim = (dim+7)/8 
         
        ; Add in the structure label. 
        ; 
         
        ; Handle variable length columns. 
        if ftype eq 'P' then begin 
 
                if dim ne 0  and dim ne 1 then begin 
                    print, 'MRDFITS: Invalid dimension for variable array column '+string(i+1) 
                    status = -1 
                    return 
                endif 
                ppos = strpos(fform, 'P') 
                vf = strmid(fform, ppos+1, 1); 
                if strpos('LXBIJAEDCM', vf) eq -1 then begin 
                    print, 'MRDFITS: Invalid type for variable array column '+string(i+1) 
                    status = -1 
                    return 
                endif 

		vcls(i) = 1
                vtpes(i) = vf
		dim = 0
                         
        endif 
         
         
        for j=0, n_elements(types) - 1 do begin 
                if ftype eq types(j) then begin 
                 
                        if dim eq 0 then begin 
                                fvalues(i) = sclstr(j) 
                        endif else begin 
                                line = arrstr(j) + dimfld(i) + ')'
                                if ftype eq 'A' then line = line + ')' 
                                fvalues(i) = line 
                        endelse 
                        goto, next_col 
                endif 
        endfor 
         
        print, 'MRDFITS: Invalid format code:',ftype, ' for column ', i+1 
        status = -1 
        return 
next_col: 
endfor 

; Check if there are any variable length columns.  If not then
; undefine vcls and vtpes
w = where(vcls eq 1)
if w(0) eq -1 then begin
	dummy = temporary(vcls)
	dummy = temporary(vtpes)
	dummy = 0
endif

if scaling then begin 
        w = where(scales ne 1.0d0 or offsets ne 0.0d0)
        if w(0) eq -1 then scaling = 0
endif

if n_elements(vcls) eq 0  and  (not scaling) and not keyword_set(columns) $
  then begin
        table = mrd_struct(fnames, fvalues, nrow, structyp=structyp, $
	  tempdir=tempdir)
endif else begin
        table = mrd_struct(fnames, fvalues, nrow, tempdir=tempdir)
endelse

if not keyword_set(silent) then begin
	print, 'MRDFITS: Binary table. ',strcompress(string(nfld)), ' columns by ',  $
	  strcompress(string(nrow)), ' rows.'
	if n_elements(vcls) gt 0 then begin
		print, 'MRDFITS: Uses variable length arrays'
	endif
endif
	
status = 0 
return 
 
end 


function mrdfits, file, extension, header,      $
	structyp = structyp,                    $
	use_colnum = use_colnum,                $
	range = range,                          $
	dscale = dscale, fscale=fscale,		$
	silent = silent, 			$
	columns = columns,			$
	no_tdim = no_tdim,                      $
	tempdir = tempdir

;+
; NAME:
;       MRDFITS
; PURPOSE:
;       Read all standard FITS data types into arrays or structures.
;
; CALLING SEQUENCE:
;       Result = MRDFITS( Filename,[Extension, Header],
;                       /FSCALE , /DSCALE , /USE_COLNUM, /NO_TDIM,
;                       RANGE=[a,b], Columns=[a,b,...])
;
; INPUTS:
;       Filename = String containing the name of the file to be read.
;		   Note that the file name may be of the form
;                  name.gz or name.Z on UNIX systems.  If so
;                  the file will be dynamically decompressed.
;
;       Extension= Extension number to be read, 0 for primary array.
;                  Assumed 0 if not specified.
;
; OUTPUTS:
;       Result = FITS data array or structure constructed from
;                the designated extension.  The format of result depends
;                upon the type of FITS data read.
;             Non-group primary array or IMAGE extension:
;                A simple multidimensional array is returned with the
;                dimensions given in the NAXISn keywords.
;             Grouped image data with PCOUNT=0.
;                As above but with GCOUNT treated as NAXIS(n+1).
;             Grouped image data with PCOUNT>0.
;               The data is returned as an array of structures.  Each
;               structure has two elements.  The first is a one-dimensional
;               array of the group parameters, the second is a multidimensional
;               array as given by the NAXIS2-n keywords.
;             ASCII and BINARY tables.
;               The data is returned as a structure with one column for
;               each field in the table.  The names of the columns are
;               normally taken from the TTYPE keywords (but see USE_COLNUM).
;               Double precision complex columns are stored as double
;               columns with a prefixed dimension of 2.  Bit field columns
;               are stored in byte arrays of the minimum necessary
;               length.  Column names are truncated to 15 characters
;               if longer, spaces are removed, and invalid characters
;               are replaced by underscores.
;
;               Columns specified as variable length columns are stored
;               with a dimension equal to the largest actual dimension
;               used.  Extra values in rows are filled with 0's or blanks.
;		If the size of the variable length column is not
;               a constant, then an additional column is created
;               giving the size used in the current row.  If the length
;               of each element of a variable length column is 0 then
;               the column is deleted.
;
; OPTIONAL OUTPUT:
;       Header = String array containing the header from the FITS extenion.
;
; OPTIONAL INPUT KEYWORDS:
;       FSCALE - If present and non-zero then scale data to float
;                numbers for arrays and columns which have either
;                non-zero offset or non-unity scale.
;                If scaling parameters are applied, then the corresponding
;                FITS scaling keywords will be modified.
;       DSCALE - As with FSCALE except that the resulting data is
;                stored in doubles.
;       SILENT - Suppress informative messages.
;       RANGE  - A scalar or two element vector giving the start
;                and end rows to be retrieved.  For ASCII and BINARY
;                tables this specifies the row number.  For GROUPed data
;                this will specify the groups.  For array images, this
;                refers to the last non-unity index in the array.  E.g.,
;                for a 3 D image with NAXIS* values = [100,100,1], the
;                range may be specified as 0:99, since the last axis
;                is suppressed.  Note that the range uses IDL indexing
;                So that the first row is row 0.
;                If only a single value, x, is given in the range,
;                the range is assumed to be [0,x-1].
;       USE_COLNUM - When creating column names for binary and ASCII tables
;                MRDFITS attempts to use the appropriate TTYPE keyword
;                values.  If USE_COLNUM is specified and non-zero then
;                column names will be generated as 'C1, C2, ... 'Cn'
;                for the number of columns in the table.
;       STRUCTYP - The structyp keyword specifies the name to be used
;                for the structure defined when reading ASCII or binary
;                tables.  Generally users will not be able to conveniently
;                combine data from multiple files unless the STRUCTYP
;                parameter is specified.  An error will occur if the
;                user specifies the same value for the STRUCTYP keyword
;                in calls to MRDFITS in the same IDL session for extensions
;                which have different structures.
;       NO_TDIM  - Disable processing of TDIM keywords.  If NO_TDIM
;                is specified MRDFITS will ignore TDIM keywords in
;                binary tables.
;	TEMPDIR - The tempdir keyword allows the user to specify
;		 the directory where temporary files may be created.
;		 This directory should be both in the IDL path
;                and writable by the user.
;	COLUMNS - This keyword allows the user to specify that only a
;                subset of columns is to be returned.  The columns
;                may be specified either as number 1,... n or by
;                name. (If USE_COLNUM is specified this should be C1,...Cn)
;		 The use of this keyword will not save time or internal
;                memory since the extraction of specified columns
;                is done after all columns have been retrieved from the
;                FITS file.
;
; EXAMPLE:
;       Read a FITS primary array:
;               a = mrdfits('TEST.FITS')    or
;               a = mrdfits('TEST.FITS', 0, header)
;       The second example also retrieves header information.
;
;       Read rows 10-100 of the second extension of a FITS file.
;               a = mrdfits('TEST.FITS', 2, header, range=[10,100])
;
;       Read a table and ask that any scalings be applied and the
;       scaled data be converted to doubles.  Use simple column names,
;       suppress outputs.
;               a = mrdfits('TEST.FITS', 1, /dscale, /use_colnum, /silent)
;
; RESTRICTIONS:
;       (1)     Cannot handle data in non-standard FITS formats.
;       (2)	Doesn't do anything with BLANK or NULL values or
;               NaN's.  They are just read in.  They may be scaled
;               if scaling is applied.
; NOTES:
;       This multiple format FITS reader is designed to provide a
;       single, simple interface to reading all common types of FITS data.
;       MRDFITS DOES NOT scale data by default.  The FSCALE or DSCALE
;       parameters must be used.
;
; PROCEDURES USED:
;	Primary internal procedures:
;           MRD_IMAGE  		-- Generate array/structure for images.
;           MRD_READ_IMAGE 	-- Read image data.
;           MRD_ASCII  		-- Generate structure for ASCII tables.
;           MRD_READ_ASCII 	-- Read an ASCII table.
;           MRD_TABLE 		-- Generate structure for Binary tables.
;           MRD_READ_TABLE 	-- Read binary table info.
;           MRD_READ_HEAP 	-- Read variable length record info.
;           MRD_SCALE 		-- Apply scaling to data.
;	    MRD_COLUMNS         -- Extract columns.
;           MRD_STRUCT 		-- Define a structure dynamically.
;           MRD_FSTRUCT		-- Define structures with too many columns
;                                  to be definable using EXECUTE.
;           FXPOSIT 		-- Move to the n'th FITS extension -- modified
;                      		   from ASTRON library routines.
;	ASTRON Library routines used:
;	    FXPAR, FXADDPAR, IEEE_TO_HOST
;
; MODIfICATION HISTORY:
;       V1.0 November 9, 1994 ----  Initial release.
;          Creator: Thomas A. McGlynn
; 	V1.1 January 20, 1995 T.A. McGlynn
;	   Fixed bug in variable length records.
;          Added TDIM support -- new routine mrd_tdim in MRD_TABLE.
;       V1.2
;	   Added support for dynamic decompression of files.
;	   Fixed further bugs in variable length record handling.
;       V1.2a
;          Added NO_TDIM keyword to turn off TDIM processing for
;          those who don't want it.
;          Bug fixes: Handle one row tables correctly, use BZERO rather than BOFFSET
;            Fix error in scaling of images.  
;-

   on_error, 2
   
; Check positional arguments.

   if n_params() le 0  or n_params() gt 3 then begin
   	print, 'MRDFITS: Usage'
	print, '   a=mrdfits(file, [extension, header],         $'
	print, '       /fscale, /dscale, /use_colnum, /silent    $'
	print, '       range=, tempdir=, structyp=, columns= )'
	return, 0
   endif
   
   if n_params() eq 1 then extension = 0
   
; Check optional arguments.
;
;  *** Structure name ***

if keyword_set(structyp) then begin
	sz = size(structyp)
	if sz(0) ne 0 then begin
		; Use first element of array
		structyp = structyp(0)
		sz = size(structyp(0))
	endif
	if sz(1) ne 7 then begin
		print, 'MRDFITS: stucture type must be a string'
		return, 0
	endif
endif

;  *** Use column numbers not names?
if not keyword_set(use_colnum) then use_colnum = 0

;  *** Get only a part of the FITS file.
if keyword_set(range) then begin
	if n_elements(range) eq 2 then arange = range $
	else if n_elements(range) eq 1 then arange = [0,range(0)-1] $
	else if n_elements(range) gt 2 then arange = range(0:1) $
	else if n_elements(range) eq 0 then arange = [-1,-1]
endif else arange = [-1,-1]

arange = long(arange)

; Open the file and position to the appropriate extension then read
; the header.

unit = fxposit(file, extension,/readonly)

if unit lt 0 then begin
  	print, 'MRDFITS: File access error'
	return, 0
endif

if eof(unit) then begin
	print,'MRDFITS: Extension past EOF'
	free_lun, unit
	return, 0
endif

mrd_hread, unit, header, status
if status lt 0 then begin
	print, 'MRDFITS: Unable to read header for extension'
	free_lun, unit
	return, 0
endif

if extension gt 0 then begin
	xten = fxpar(header,'XTENSION')
	if xten eq 'IMAGE   ' then type = 0   $
	else if xten eq 'TABLE   ' then type = 1  $
	else if xten eq 'BINTABLE'  or xten eq 'A3DTABLE' then type = 2 $
	else begin 
		print, 'MRDFITS: Unable to process extension type:', xten
		return, 0
	endelse
endif else type = 0

scaling = keyword_set(fscale) or keyword_set(dscale)

if type eq 0 then begin

	;*** Images/arrays
	
	mrd_image, header, arange, rsize, table, scales, offsets, status, $
	  silent=silent, tempdir=tempdir
	if status ge 0 and rsize gt 0 then mrd_read_image, unit, arange, rsize, table
	size = rsize
	
endif else if type eq 1 then begin

        ;*** ASCII tables.
	
	mrd_ascii, header, structyp, use_colnum,   $
	    arange, table, nbytes, nrows, nfld,  $
	    typarr, posarr, lenarr, nullarr, fnames, fvalues, $
	    scales, offsets, scaling, status, silent=silent, tempdir=tempdir, $
	    columns=columns
	size = nbytes*nrows
	
	if status ge 0   and  size gt 0  then begin
	
	        ;*** Read data.
	        mrd_read_ascii, unit,  arange, nbytes, nrows, $
		  nfld, typarr, posarr, lenarr, nullarr, table
	      
	        ;*** Extract desired columns.
		if status ge 0 and keyword_set(columns) then $
			mrd_columns, table, columns, fnames, fvalues, vcls, vtps, $
			  scales, offsets, $
			  scaling, structyp=structyp, tempdir=tempdir
	endif
	
endif else begin

     ; *** Binary tables.

     mrd_table, header, structyp, use_colnum,   $
	 arange, rsize, table, nfld, typarr,        $ 
	 fnames, fvalues, vcls, vtpes, scales, offsets, scaling, status, $
	 silent=silent, tempdir=tempdir, columns=columns, no_tdim=no_tdim


     size = nfld*(arange(1) - arange(0) + 1)
     if status ge 0  and  size gt 0  then begin
     
        ;*** Read data.
     	mrd_read_table, unit, arange, rsize,  $
	  structyp, nfld, typarr, table

     	if status ge 0 and keyword_set(columns) then $
	
	   ;*** Extract desired columns.
           mrd_columns, table, columns, fnames, fvalues, $
	   	vcls, vtpes, scales, offsets, scaling, structyp=structyp, $
	   	tempdir=tempdir
	 
     
     	if status ge 0 and n_elements(vcls) gt 0 then begin 
	
	        ;*** Get variable length columns
	   	mrd_read_heap, unit, header, arange, fnames, fvalues, $
		 vcls, vtpes, table, structyp, scaling, scales, offsets, status, $
		 silent=silent, tempdir=tempdir
     	endif
      endif

endelse


if unit gt 0 then begin
	free_lun, unit
endif

if  status ge 0  and  scaling  and  size gt 0  then begin
	w = where(scales ne 1.d0  or  offsets ne 0.0d0)
	
	;*** Apply scalings.
	if w(0) ne -1 then mrd_scale, type, scales, offsets, table, header,$
	    fnames, fvalues, 1+arange(1)-arange(0), structyp=structyp, $
	    dscale=dscale, tempdir=tempdir
endif

			      
if status ge 0 then return, table else return, 0

end