;+
; NAME:
;     MRDFITS
;
; PURPOSE:
;     Read all standard FITS data types into arrays or structures.
;
; EXPLANATION:
;      Further information on MRDFITS is available at
;      http://idlastro.gsfc.nasa.gov/mrdfits.html 
;
; CALLING SEQUENCE:
;      Result = MRDFITS( Filename/FileUnit,[Exten_no/Exten_name, Header],
;                       /FSCALE , /DSCALE , /UNSIGNED,
;                       ALIAS=strarr[2,n], /USE_COLNUM,
;                       /NO_TDIM, ROWS = [a,b,...], $
;                       /POINTER_VAR, /FIXED_VAR, EXTNUM= 
;                       RANGE=[a,b], COLUMNS=[a,b,...]), ERROR_ACTION=x,
;                       COMPRESS=comp_prog, STATUS=status, /VERSION )
;
; INPUTS:
;      Filename = String containing the name of the file to be read or
;                 file number of an open unit.  If a unit is specified
;                 if will be left open positioned to read the next HDU.
;                 If the file name ends in .gz (or .Z on Unix systems)
;                 the file will be dynamically decompressed.
;      FiluUnit = An integer file unit which has already been
;                 opened for input.  Data will be read from this
;                 unit and the unit will be left pointing immediately
;                 after the HDU that is read.  Thus to read a compressed
;                 file with many HDU's a user might do something like:
;                      lun=fxposit(filename, 3)  ; Skip the first three HDU's
;                      repeat begin
;                          thisHDU = mrdfits(lun, 0, hdr, status=status)
;                          ... process the HDU ...
;                      endrep until status lt 0
;
;      Exten_no= Extension number to be read, 0 for primary array.
;                 Assumed 0 if not specified.
;                 If a unit rather than a filename
;                 is specified in the first argument, this is
;                 the number of HDU's to skip from the current position.
;      Exten_name - Name of the extension to read (as stored in the EXTNAME
;                 keyword).   This is a slightly slower method then specifying
;                 the extension number.
; 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).
;               Bit field columns
;               are stored in byte arrays of the minimum necessary
;               length.  Spaces and invalid characters are replaced by 
;               underscores, and other invalid tag names are converted using
;               the IDL_VALIDNAME(/CONVERT_ALL) function.
;               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.  This additional column will 
;               have a tag name of the form L#_"colname" where # is the column
;               number and colname is the column name of the variable length
;               column.   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 extension.
;
; OPTIONAL INPUT KEYWORDS:
;       ALIAS    The keyword allows the user to specify the column names
;                to be created when reading FITS data.  The value of
;                this keyword should be a 2xn string array.  The first
;                value of each pair of strings should be the desired
;                tag name for the IDL column.  The second should be
;                the FITS TTYPE value.  Note that there are restrictions
;                on valid tag names.  The order of the ALIAS keyword
;                is compatible with MWRFITS.
;       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 or some combination of these two.
;                If USE_COLNUM is specified names 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.
;       COMPRESS - This keyword allows the user to specify a
;                decompression program to use to decompress a file that
;                will not be automatically recognized based upon
;                the file name.
;       /DSCALE - As with FSCALE except that the resulting data is
;                stored in doubles.
;       ERROR_ACTION - Set the on_error action to this value (defaults
;                to 2).
;       /FIXED_VAR- Translate variable length columns into fixed length columns
;                and provide a length column for truly varying columns.
;                This was only behavior prior to V2.5 for MRDFITS and remains
;                the default (see /POINTER_VAR)
;       /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.
;       NO_TDIM  - Disable processing of TDIM keywords.  If NO_TDIM
;                is specified MRDFITS will ignore TDIM keywords in
;                binary tables.
;       /POINTER_VAR- Use pointer arrays for variable length columns.
;                In addition to changing the format in which
;                variable length arrays are stored, if the pointer_var
;                keyword is set to any value other than 1 this also disables
;                the deletion of variable length columns. (See /FIXED_VAR)
;                Note that because pointers may be present in the output
;                structure, the user is responsible for memory management
;                when deleting or reassigning the structure (e.g. use HEAP_FREE
;                first).
;       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].
;       ROWS -  A scalar or vector specifying a  specific row or rows to read 
;               (first row is 0).   For example to read rows 0,
;               12 and 23 only, set ROWS=[0,12,23].   Valid for images, ASCII 
;               and binary tables, but not GROUPed data.   For images
;               the row numbers refer to the last non-unity index in the array.   
;               Cannot be used at the same time as the RANGE keyword
;       /SILENT - Suppress informative messages.
;       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.
;       /UNSIGNED - For integer data with appropriate zero points and scales
;                read the data into unsigned integer arrays.
;       /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.
;       /VERSION Print the current version number
;
; OPTIONAL OUTPUT KEYWORDS:
;       EXTNUM - the number of the extension actually read.   Useful if the
;                 user specified the extension by name.
;       STATUS - A integer status indicating success or failure of
;                the request.  A status of >=0 indicates a successful read.
;                Currently
;                    0 -> successful completion
;                   -1 -> error
;                   -2 -> end of file
;
; EXAMPLES:
;       (1) Read a FITS primary array:
;               a = mrdfits('TEST.FITS')    or
;               a = mrdfits('TEST.FITS', 0, header)
;       The second example also retrieves header information.
;
;       (2) Read rows 10-100 of the second extension of a FITS file.
;               a = mrdfits('TEST.FITS', 2, header, range=[10,100])
;
;       (3) 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)
;
;       (4) Read rows 3, 34 and 52 of a binary table and request that 
;           variable length columns be stored as a pointer variable in the 
;           output structure
;              a = mrdfits('TEST.FITS',1,rows=[3,34,52],/POINTER)
  
; 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.
;
;       MRDFITS support 64 bit integer data types, which are tentatively
;       included in the FITS standard.
;       http://fits.gsfc.nasa.gov/fits_64bit.html
;    
;
; PROCEDURES USED:
;       The following procedures are contained in the main MRDFITS program.
;           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.
;
;        Other ASTRON Library routines used
;           FXPAR(), FXADDPAR, FXPOSIT, FXMOVE(), MRD_STRUCT(), MRD_SKIP, 
;
; 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.  
;       V1.2b 
;          Changed MRD_HREAD to handle null characters in headers.
;       V2.0 April 1, 1996
;          -Handles FITS tables with an arbitrary number of columns.
;          -Substantial changes to MRD_STRUCT to allow the use of
;          substructures when more than 127 columns are desired.
;          -All references to table columns are now made through the
;          functions MRD_GETC and MRD_PUTC.  See description above.
;          -Use of SILENT will now eliminate compilation messages for
;          temporary functions.
;          -Bugs in handling of variable length columns with either
;          a single row in the table or a maximum of a single element
;          in the column fixed.
;          -Added support for DCOMPLEX numbers in binary tables (M formats) for
;          IDL versions above 4.0.  
;          -Created regression test procedure to check in new versions.
;          -Added error_action parameter to allow user to specify
;          on_error action.  This should allow better interaction with
;          new CHECK facility.  ON_ERROR statements deleted from
;          most called routines.
;          - Modified MRDFITS to read in headers containing null characters
;          with a warning message printed.
;       V2.0a April 16, 1996
;          - Added IS_IEEE_BIG() checks (and routine) so that we don't
;          worry about IEEE to host conversions if the machine's native
;          format is IEEE Big-endian.
;       V2.1 August 24, 1996
;          - Use resolve_routine for dynamically defined functions
;          for versions > 4.0
;          - Fix some processing in random groups format.
;          - Handle cases where the data segment is--legally--null.
;          In this case MRDFITS returns a scalar 0.
;          - Fix bugs with the values for BSCALE and BZERO (and PSCAL and
;          PZERO) parameters set by MRDFITS.
;       V2.1a April 24, 1997  Handle binary tables with zero length columns
;       V2.1b May 13,1997 Remove whitespace from replicate structure definition
;       V2.1c May 28,1997 Less strict parsing of XTENSION keyword
;       V2.1d June 16, 1997 Fixed problem for >32767 entries introduced 24-Apr
;       V2.1e Aug 12, 1997 Fixed problem handling double complex arrays
;       V2.1f Oct 22, 1997 IDL reserved words can't be structure tag names
;       V2.1g Nov 24, 1997 Handle XTENSION keywords with extra blanks.
;       V2.1h Jul 26, 1998 More flexible parsing of TFORM characters
;       V2.2 Dec 14, 1998 Allow fields with longer names for
;                        later versions of IDL.
;                        Fix handling of arrays in scaling routines.
;                        Allow >128 fields in structures for IDL >4.0
;                        Use more efficient structure copying for
;                        IDL>5.0
;       V2.2b June 17, 1999 Fix bug in handling case where
;                           all variable length columns are deleted
;                           because they are empty.
;       V2.3 March 7, 2000 Allow user to supply file handle rather
;                          than file name.
;                          Added status field.
;                          Now needs FXMOVE routine
;       V2.3b April 4, 2000
;                          Added compress option (from D. Palmer)
;       V2.4  July 4, 2000 Added STATUS=-1 for "File access error" (Zarro/GSFC)
;       V2.4a May 2, 2001  Trim binary format string   (W. Landsman)
;       V2.5 December 5, 2001 Add unsigned, alias, 64 bit integers. version, $
;                           /pointer_val, /fixed_var.
;       V2.5a Fix problem when both the first and the last character
;            in a TTYPEnn value are invalid structure tag characters
;       V2.6 February 15, 2002 Fix error in handling unsigned numbers, $
;                           and 64 bit unsigneds. (Thanks to Stephane Beland)
;       V2.6a September 2, 2002 Fix possible conflicting data structure for
;                          variable length arrays (W. Landsman)
;       V2.7 July, 2003  Added Rows keyword (W. Landsman)
;       V2.7a September  2003 Convert dimensions to long64 to handle huge files
;       V2.8 October 2003 Use IDL_VALIDNAME() function to ensure valid tag names
;                         Removed OLD_STRUCT and TEMPDIR keywords W. Landsman
;       V2.9 February 2004 Added internal MRD_FXPAR procedure for faster
;                     processing of binary table headers E. Sheldon
;       V2.9a March 2004 Restore ability to read empty binary table W. Landsman
;             Swallow binary tables with more columns than given in TFIELDS
;       V2.9b Fix to ensure order of TFORMn doesn't matter
;       V2.9c Check if extra degenerate NAXISn keyword are present W.L. Oct 2004 
;       V2.9d Propagate /SILENT to MRD_HREAD, more LONG64 casting W. L. Dec 2004
;       V2.9e Add typarr[good] to fix a problem reading zero-length columns
;             A.Csillaghy, csillag@ssl.berkeley.edu (RHESSI)
;       V2.9f Fix problem with string variable binary tables, possible math 
;             overflow on non-IEEE machines  WL Feb. 2005 
;       V2.9g Fix problem when setting /USE_COLNUM   WL Feb. 2005
;       V2.10 Use faster keywords to BYTEORDER  WL May 2006
;       V2.11  Add ON_IOERROR, CATCH, and STATUS keyword to MRD_READ_IMAGE to
;             trap EOF in compressed files DZ  Also fix handling of unsigned 
;             images when BSCALE not present  K Chu/WL   June 2006 
;       V2.12 Allow extension to be specified by name, added EXTNUM keyword
;                     WL    December 2006
;       V2.12a Convert ASCII table column to DOUBLE if single precision is
;                 insufficient  
;       V2.12b Fixed problem when both /fscale and /unsigned are set 
;                  C. Markwardt    Aug 2007
;       V2.13  Use SWAP_ENDIAN_INPLACE instead of IEEE_TO_HOST and IS_IEEE_BIG
;                W. Landsman Nov 2007
;       V2.13a One element vector allowed for file name W.L. Dec 2007
;-
PRO mrd_fxpar, hdr, xten, nfld, nrow, rsize, fnames, fforms, scales, offsets
;
;  Check for valid header.  Check header for proper attributes.
;
  S = SIZE(HDR)
  IF ( S[0] NE 1 ) OR ( S[2] NE 7 ) THEN $
    MESSAGE,'FITS Header (first parameter) must be a string array'

  xten = fxpar(hdr, 'XTENSION')
  nfld = fxpar(hdr, 'TFIELDS')
  nrow = long64(fxpar(hdr, 'NAXIS2'))
  rsize = long64(fxpar(hdr, 'NAXIS1'))

  ;; will extract these for each
  names = ['TTYPE','TFORM', 'TSCAL', 'TZERO']
  nnames = n_elements(names)

;  Start by looking for the required TFORM keywords.    Then try to extract it 
;  along with names (TTYPE), scales (TSCAL), and offsets (TZERO)
   
 keyword = STRMID( hdr, 0, 8)

;
;  Find all instances of 'TFORM' followed by
;  a number.  Store the positions of the located keywords in mforms, and the
;  value of the number field in n_mforms
;

  mforms = WHERE(STRPOS(keyword,'TFORM') GE 0, n_mforms)
  if n_mforms GT nfld then begin
        message,/CON, $
        'WARNING - More columns found in binary table than specified in TFIELDS'
        n_mforms = nfld
        mforms = mforms[0:nfld-1]
  endif


  IF ( n_mforms GT 0 ) THEN BEGIN
      numst= STRMID(hdr[mforms], 5 ,3)      
      number = INTARR(n_mforms)-1

      FOR i = 0, n_mforms-1 DO         $
        IF VALID_NUM( numst[i], num) THEN number[i] = num

      igood = WHERE(number GE 0, n_mforms)
      IF n_mforms GT 0 THEN BEGIN
          mforms = mforms[igood]
          number = number[igood]
          numst = numst[igood]
      ENDIF

  ENDIF ELSE RETURN              ;No fields in binary table

  ;; The others
  fnames = strarr(n_mforms)
  fforms = strarr(n_mforms) 
  scales = dblarr(n_mforms)
  offsets = dblarr(n_mforms)

  ;;comments = strarr(n_mnames)

  fnames_names  = 'TTYPE'+numst
  scales_names  = 'TSCAL'+numst
  offsets_names = 'TZERO'+numst
  number = number -1    ;Make zero-based


  match, keyword, fnames_names, mkey_names, mnames, count = N_mnames

  match, keyword, scales_names, mkey_scales, mscales, count = N_mscales

  match, keyword, offsets_names, mkey_offsets, moffsets,count = N_moffsets

  FOR in=0L, nnames-1 DO BEGIN 

      CASE names[in] OF
          'TTYPE': BEGIN 
              tmatches = mnames 
              matches = mkey_names
              nmatches = n_mnames
              result = fnames
          END 
          'TFORM': BEGIN 
              tmatches = lindgen(n_mforms)
              matches = mforms
              nmatches = n_mforms
              result = fforms
          END 
          'TSCAL': BEGIN 
              tmatches = mscales
              matches = mkey_scales
              nmatches = n_mscales
              result = scales
          END 
          'TZERO': BEGIN 
              tmatches = moffsets
              matches = mkey_offsets
              nmatches = n_moffsets
              result = offsets
          END 
          ELSE: message,'What?'
      ENDCASE 

      ;;help,matches,nmatches

;
;  Extract the parameter field from the specified header lines.  If one of the
;  special cases, then done.
;
      IF nmatches GT 0 THEN BEGIN 
          
          ;; "matches" is a subscript for hdr and keyword.
          ;; get just the matches in line
          
          line = hdr[matches]
          svalue = STRTRIM( STRMID(line,9,71),2)
          
          FOR i = 0, nmatches-1 DO BEGIN
              IF ( STRMID(svalue[i],0,1) EQ "'" ) THEN BEGIN

                  ;; Its a string
                  test = STRMID( svalue[i],1,STRLEN( svalue[i] )-1)
                  next_char = 0
                  off = 0
                  value = ''
;
;  Find the next apostrophe.
;
NEXT_APOST:
                  endap = STRPOS(test, "'", next_char)
                  IF endap LT 0 THEN MESSAGE,         $
                    'WARNING: Value of '+nam+' invalid in '+ " (no trailing ')", /info
                  value = value + STRMID( test, next_char, endap-next_char )
;
;  Test to see if the next character is also an apostrophe.  If so, then the
;  string isn't completed yet.  Apostrophes in the text string are signalled as
;  two apostrophes in a row.
;
                  IF STRMID( test, endap+1, 1) EQ "'" THEN BEGIN    
                      value = value + "'"
                      next_char = endap+2      
                      GOTO, NEXT_APOST
                  ENDIF

;
; CM 19 Sep 1997
; This is a string that could be continued on the next line.  Check this
; possibility with the following four criteria: *1) Ends with '&'
; (2) Next line is CONTINUE  (3) LONGSTRN keyword is present (recursive call to
;  FXPAR) 4. /NOCONTINE is not set
                  
;
;  If not a string, then separate the parameter field from the comment field.
;
              ENDIF ELSE BEGIN
                  ;; not a string
                  test = svalue[I]
                  slash = STRPOS(test, "/")
                  IF slash GT 0 THEN BEGIN
                      test = STRMID(test, 0, slash)
                  END 
;
;  Find the first word in TEST.  Is it a logical value ('T' or 'F')?
;
                  test2 = test
                  value = GETTOK(test2,' ')
                  test2 = STRTRIM(test2,2)
                  IF ( value EQ 'T' ) THEN BEGIN
                      value = 1
                  END ELSE IF ( value EQ 'F' ) THEN BEGIN
                      value = 0
                  END ELSE BEGIN
;
;  Test to see if a complex number.  It's a complex number if the value and the
;  next word, if any, both are valid numbers.
;
                      IF STRLEN(test2) EQ 0 THEN GOTO, NOT_COMPLEX
                      test2 = GETTOK(test2,' ')
                      IF VALID_NUM(value,val1) AND VALID_NUM(value2,val2) $
                        THEN BEGIN
                          value = COMPLEX(val1,val2)
                          GOTO, GOT_VALUE
                      ENDIF
;
;  Not a complex number.  Decide if it is a floating point, double precision,
;  or integer number.  If an error occurs, then a string value is returned.
;  If the integer is not within the range of a valid long value, then it will 
;  be converted to a double.  
;
NOT_COMPLEX:
                      ON_IOERROR, GOT_VALUE
                      value = test
                      IF NOT VALID_NUM(value) THEN GOTO, GOT_VALUE

                      IF (STRPOS(value,'.') GE 0) OR (STRPOS(value,'E') $
                                                      GE 0) OR (STRPOS(value,'D') GE 0) THEN BEGIN
                          IF ( STRPOS(value,'D') GT 0 ) OR $
                            ( STRLEN(value) GE 8 ) THEN BEGIN
                              value = DOUBLE(value)
                          END ELSE value = FLOAT(value)
                      ENDIF ELSE BEGIN
                          lmax = 2.0D^31 - 1.0D
                          lmin = -2.0D31
                          value = DOUBLE(value)
                          if (value GE lmin) and (value LE lmax) THEN $
                            value = LONG(value)
                      ENDELSE
                      
;
GOT_VALUE:
                      ON_IOERROR, NULL
                  ENDELSE
              ENDELSE           ; if string
;
;  Add to vector if required.
;
              
              result[tmatches[i]] = value
             
          ENDFOR

          CASE names[in] OF
              'TTYPE': fnames[number] = strtrim(result, 2)
              'TFORM': fforms[number] = strtrim(result, 2)
              'TSCAL': scales[number] = result
              'TZERO': offsets[number] = result
              ELSE: message,'What?'
          ENDCASE 

;
;  Error point for keyword not found.
;
      ENDIF
;



  ENDFOR 
END
  

; Get a tag name give the column name and index
function  mrd_dofn, name, index, use_colnum, alias=alias

    ; Check if the user has specified an alias.

    if n_elements(name) eq 0 then name = 'C'+strtrim(index, 2)
    name = strtrim(name)
    if keyword_set(alias) then begin
	sz = size(alias)
	
	if (sz[0] eq 1 or sz[0] eq 2) and sz[1] eq 2 and sz[sz[0]+1] eq 7 then begin
	    w=where(name eq alias[1,*])
	    if (w[0] ne -1) then begin
		name = alias[0,w[0]];
	    endif
	endif
    endif
    ; Convert the string name to a valid variable name.  If name 
    ; is not defined generate the string Cnn when nn is the index 
    ; number.

    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'+strtrim(index,2) 
        endelse 
    endif else begin 
        str = 'C'+strtrim(index,2) 
    endelse 
 
    return, IDL_VALIDNAME(str,/CONVERT_ALL) 
 
end 

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



; Parse the TFORM keyword and return the type and dimension of the 
; data. 
pro mrd_doff, form, dim, type 
 
    ; 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 then return              ;Modified from len-1 on 26-Jul-1998
 
    if i gt 0 then begin 
        dim = long(strmid(form, 0, i)) 
        if dim EQ 0l then dim = -1l
    endif else begin 
        dim = 0 
    endelse 
 
    type = strmid(form, i, 1) 
end 



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

;  Check that this name is unique with regard to other column names.

function mrd_chkfn, name, namelist, index 
 
    ; 
    ; 

    maxlen = 127

    if strlen(name) gt maxlen then name = strmid(name, 0, maxlen) 
    w = where(name eq strmid(namelist, 0, maxlen) )
    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

; Find the appropriate offset for a given unsigned type.
; The type may be given as the bitpix value or the IDL
; variable type.

function mrd_unsigned_offset, type
    	    
    if (type eq 12 or type eq 16) then begin
	return, uint(32768)
    endif else if (type eq 13 or type eq 32) then begin
	return, ulong('2147483648')
    endif else if (type eq 15 or type eq 64) then begin
	return, ulong64('9223372036854775808');
    endif
    return, 0
end



; Can we treat this data as unsigned?

function mrd_chkunsigned, bitpix, scale, zero, unsigned=unsigned

    if not keyword_set(unsigned) then return, 0
    
    ; This is correct but we should note that
    ; FXPAR returns a double rather than a long.
    ; Since the offset is a power of two
    ; it is an integer that is exactly representable
    ; as a double.  However, if a user were to use
    ; 64 bit integers and an offset close to but not
    ; equal to 2^63, we would erroneously assume that
    ; the dataset was unsigned...

    
    if scale eq 1 then begin
	if (bitpix eq 16 and zero eq 32768L)  or                    $
	   (bitpix eq 32 and zero eq ulong('2147483648')) or        $
	   (bitpix eq 64 and zero eq ulong64('9223372036854775808')) then begin
	    return, 1
	endif
    endif
    return, 0
end

; Is this one of the IDL unsigned types?
function mrd_unsignedtype, data
       
    type = size(data,/ type) 
    if (type eq 12) or (type eq 13) or (type eq 15) then return, type $
                                                    else return, 0
    
end
    
; Return the currrent version string for MRDFITS
function mrd_version
    return, '2.13a'
end
;=====================================================================
; END OF GENERAL UTILITY FUNCTIONS ===================================
;=====================================================================


; Parse the TFORM keyword and return the type and dimension of the
; data.
pro mrd_atype, form, type, slen


    ; 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
   if (type EQ 'F') or (type EQ 'E') then $     ;Updated April 2007
        if (slen GE 8) then type = 'D'

end


; Read in the table information.
pro mrd_read_ascii, unit, range, nbytes, nrows, nfld, typarr, posarr,   $
     lenarr, nullarr, table, old_struct=old_struct, rows=rows

    ;
    ; Unit          Unit to read data from.
    ; Range         Range of  to be read
    ; 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.
    ; Old_struct    Should recursive structure format be used?

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

    if range[0] gt 0 then mrd_skip, unit, nbytes*range[0]
    readu,unit, bigstr
    if N_elements(rows) GT 0 then bigstr = bigstr[*,rows-range[0]] 

    ; Skip to the end of the data area.

    nSkipRow = nrows - range[1] - 1
    nskipB = 2880 - (nbytes*nrows) mod 2880
    if nskipB eq 2880 then nskipB = 0

    mrd_skip, unit, nskipRow*nbytes+nskipB

    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
	    
		curr_col = table.(i)
	    
            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
                    curr_col[w] = long(flds[w])
                endif else if typarr[i] eq 'E' or typarr[i] eq 'F' then begin
                    curr_col[w] = float(flds[w])
                endif else if typarr[i] eq 'D' then begin
                    curr_col[w] = double(flds[w])
                endif else if typarr[i] eq 'A' then begin
                    curr_col[w] = flds[w]
                endif
            endif

                table.(i) = curr_col
                
        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


; Define a structure to hold a FITS ASCII table.
pro mrd_ascii, header, structyp, use_colnum,   $
    range, table, $
    nbytes, nrows, nfld, typarr, posarr, lenarr, nullarr, $
    fnames, fvalues, scales, offsets, scaling, status, rows = rows, $
    silent=silent, columns=columns, alias=alias

    ;
    ; 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
    ; Nrows                 Number of rows in table
    ; 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
    ; Fname                 Column names
    ; Fvalues               Formats for columns
    ; Scales/offsets        Scaling factors for columns
    ; Scaling               Do we need to scale?
    ; Status                Return status.

    table = 0

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

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

    nfld = fxpar(header, 'TFIELDS')
    nrows = long64( fxpar(header, 'NAXIS2'))
    nbytes = long64( 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

    if N_elements(rows) EQ 0 then nrows = range[1] - range[0] + 1 else begin
          bad = where(rows GT nrows, Nbad)
          if Nbad GT 0  then begin 
             print,'MRDFITS: Row numbers must be between 0 and ' + $
                      strtrim(nrows-1,2)
             status = -1
             return
          endif      
          nrows = N_elements(rows)
     endelse

    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, count=cnt)
	if cnt eq 0 then xx = temporary(fname)
        fform = fxpar(header, 'TFORM' + suffix)
        fpos = fxpar(header, 'TBCOL' + suffix)
        fnull = fxpar(header, 'TNULL' + suffix, count=cnt)
	if cnt eq 0 then fnull = ''
        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, alias=alias)
        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,'+strtrim(flen,2)+'))'
                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, $
           silent=silent)
    endif else begin
        table = mrd_struct(fnames, fvalues, nrows, silent=silent)
    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


; Eliminate columns from the table that do not match the
; user specification.
pro  mrd_columns, table, columns, fnames, fvalues, $
    vcls, vtpes, scales,  offsets, scaling,        $
    structyp=structyp, silent=silent



    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, silent=silent )
        endif else begin
            tabx = mrd_struct(fnames, fvalues, ndim, structyp=structyp, silent=silent )
        endelse
        
        for i=0, n_elements(tcols)-1 do begin
                tabx.(i) = table.(tcols[i]);
        endfor

        table = tabx
    endelse
    
end


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

    error=0
    catch,error
    if error ne 0 then begin
     catch,/cancel
     status=-2
     return
    endif

    ; If necessary skip to beginning of desired data. 

    if range[0] gt 0 then mrd_skip, unit, range[0]*rsize

    status=-2
    if rsize eq 0 then return 

    on_ioerror,done
    readu, unit, table

    if N_elements(rows) GT 0 then begin
    row1 = rows- range[0]
    case size(table,/n_dimen) of 
    1: table = table[row1]
    2: table = table[*,row1]
    3: table = table[*,*,row1]
    4: table = table[*,*,*,row1]
    5: table = table[*,*,*,*,row1]
    6: table = table[*,*,*,*,*,row1]
    7: table = table[*,*,*,*,*,*,row1]
    8: table = table[*,*,*,*,*,*,*,row1]
    else: begin 
          print,'MRDFITS: Subscripted image must be between 1 and 8 dimensions'
          status = -1
          return
          end
    endcase
    endif

    ; Skip to the end of the data

    skipB = 2880 - (maxd*rsize) mod 2880
    if skipB eq 2880 then skipB = 0

    if range[1] lt maxd-1 then begin
        skipB = skipB + (maxd-range[1]-1)*rsize
    endif

    mrd_skip, unit, skipB
    swap_endian_inplace, table,/swap_if_little

    ; Fix offset for unsigned data
    type = mrd_unsignedtype(table)
    if type gt 0 then begin
	table = table - mrd_unsigned_offset(type)
    endif
    
    status=0
    done:

;-- probably an EOF 

    if status ne 0 then free_lun,unit
 
    return
end 

; Truncate superfluous axes.

pro mrd_axes_trunc,naxis, dims, silent

    mysilent = silent
    for i=naxis-1,1,-1 do begin 

        if dims[i] eq 1 then begin
            if not mysilent then begin
                print, 'MRDFITS: Truncating unused dimensions'
                mysilent = 1
            endif
            dims = dims[0:i-1] 
            naxis = naxis - 1 
        
        endif else return
     
    endfor 
 
    return
end

; Define structure/array to hold a FITS image. 
pro mrd_image, header, range, maxd, rsize, table, scales, offsets, scaling, $
  status, silent=silent, unsigned=unsigned, rows = rows
 
    ; 
    ; 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                Return status
    ; Silent=silent         Suppress info messages?
 
    table = 0

    ; type    0         1           2         3         4         5  6  7  8  9 10 11        12         13          14          15
    lens =  [ 0,        1,          2,        4,        4,        8, 0, 0, 0, 0, 0, 0,        2,         4,          8,          8] 
    typstrs=['',   'Byte',    'Int*2',  'Int*4', 'Real*4', 'Real*8','','','','','','', 'UInt*2',  'Uint*4',    'Int*8',    'Uint*8']
    typarr= ['', 'bytarr',   'intarr', 'lonarr', 'fltarr', 'dblarr','','','','','','','uintarr', 'ulonarr', 'lon64arr', 'ulon64arr'] 
 
    status = 0 
 
 
    naxis = fxpar(header, 'NAXIS') 
    bitpix= fxpar(header, 'BITPIX') 
    if naxis gt 0 then begin 
          dims = long64(fxpar(header, 'NAXIS*', Count = N_axis)) 
          if N_axis GT naxis then begin
; Check if extra NAXISn keywords are present (though this is not legal FITS)
                   nextra = N_axis - naxis
                   dim_extra = dims[naxis:N_axis-1]
                   if total(dim_extra) EQ nextra then $
                        dims = dims[0:naxis-1] else $
                   message,'ERROR - NAXIS = ' + strtrim(naxis,2) +  $
                          ' but NAXIS' + strtrim(N_axis,2) + ' keyword present'
          endif
   endif else dims = 0
    
    gcount = fxpar(header, 'GCOUNT') 
    pcount = fxpar(header, 'PCOUNT')
    isgroup = fxpar(header, 'GROUPS')
    gcount = long(gcount)

    xscale = fxpar(header, 'BSCALE', count=cnt)
    if cnt eq 0 then xscale = 1      ;Corrected 06/29/06
    
    xunsigned = mrd_chkunsigned(bitpix,  xscale, $
				fxpar(header, 'BZERO'), unsigned=unsigned)
    ; Note that type is one less than the type signifier returned in the size call.
    type = -1
    
    if not xunsigned then begin 
 
        if bitpix eq 8        then type = 1     $ 
        else if bitpix eq  16 then type = 2     $ 
        else if bitpix eq  32 then type = 3     $ 
        else if bitpix eq -32 then type = 4     $ 
        else if bitpix eq -64 then type = 5     $
        else if bitpix eq  64 then type = 14

    endif else begin

	if bitpix eq 16       then type = 12     $
	else if bitpix eq  32 then type = 13     $
	else if bitpix eq  64 then type = 15
	
    endelse

    if type eq -1 then begin
	print,'MRDFITS: Error: Invalid BITPIX: '+strtrim(bitpix)
	table = 0
	return
    endif

    ; Note that for random groups data we must ignore the first NAXISn keyword. 
    if isgroup GT 0  then begin 


        range[0] = range[0] > 0
        if (range[1] eq -1) then begin
            range[1] = gcount-1
        endif else begin
            range[1] = range[1] < gcount - 1
        endelse
	
	maxd = gcount
        
        if (n_elements(dims) gt 1) then begin
            dims = dims[1:*]
            naxis = naxis-1
        endif else begin
            print, 'MRDFITS: Warning: No data specified for group data.'
            dims = [0]
            naxis = 0
        endelse
        
        ; The last entry is the scaling for the sample data.
        
        if (pcount gt 0) then begin
            scales  = dblarr(pcount+1)
            offsets = dblarr(pcount+1)
        endif
        
        values = strarr(2)
        
        
        mrd_axes_trunc, naxis, dims, keyword_set(silent)
        
        values[0] = typarr[type] + "("+string(pcount)+")" 
        rsize = dims[0] 
        sarr = "(" + strcompress(string(dims[0]), /remo )
         
        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:', $
          ' Ngroup=',strcompress(string(gcount)),' Npar=',                   $
          strcompress(string(pcount),/remo), ' Group=', 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, $
                           silent=silent)

	if xunsigned then begin
	    fxaddpar,header, 'BZERO', 0, 'Reset by MRDFITS v'+mrd_version()
	endif
           

        for i=0, pcount-1 do begin
	    
            istr = strcompress(string(i+1),/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 
         
        mrd_axes_trunc, naxis, dims, keyword_set(silent)

                
        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 

        Nlast = dims[naxis-1]   
        dims[naxis-1] = range[1]-range[0]+1
        pdims = dims
        if N_elements(rows) GT 0 then begin
             if max(rows) GE Nlast then begin 
               print, 'MRDFITS: Row numbers must be between 0 and ' + $
                      strtrim(Nlast-1,2)
               status = -1 & rsize = 0
               return
             endif
             pdims[naxis-1] = N_elements(rows)
        endif 
 
        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(pdims[i]),/remo)
            endfor
            str = str+')'
            print, 'MRDFITS: Image array ',str, '  Type=', typstrs[type]
        endif
         
        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   
        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)

	if xunsigned then begin
	    fxaddpar,header, 'BZERO', 0, 'Updated by MRDFITS v'+mrd_version()
	endif
	
        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

; Scale an array of pointers
pro mrd_ptrscale, array, scale, offset

    for i=0, n_elements(array)-1 do begin
        if ptr_valid(array[i]) then begin
	    array[i] = ptr_new(*array[i] * scale + offset)
	endif
    endfor
end

; Scale a FITS array or table.
pro mrd_scale, type, scales, offsets, table, header,  $
               fnames, fvalues, nrec, dscale = dscale, structyp=structyp, silent=silent
    ;
    ; 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.
 
    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))

		; Handle pointer columns
		if sz[sz[0]+1] eq 10 then begin
		    fvalues[col] = 'ptr_new()'

		; Scalar columns
		endif else if sz[0] eq 0 then begin
                    fvalues[col] = sclr

		; Vectors
                endif else begin
                    str = vc + '('
                    for j=0, sz[0]-1 do begin
                        if j ne 0 then str = str + ','
                        str = str + strtrim(sz[j+1],2)
                    endfor
                    str = str + ')'
                    fvalues[col] = str
		    
                endelse
		
            endfor
        endif

        tabx = mrd_struct(fnames, fvalues, nrec, structyp=structyp, silent=silent )

	; Just copy the unscaled columns
        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
	    
 		
		sz = size(tabx.(w[i]))
		if sz[sz[0]+1] eq 10 then begin
		    mrd_ptrscale, table.(w[i]), scales[w[i]], offsets[w[i]]
		endif
                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, 'BSCALE', 1.0, 'Set by MRD_SCALE'
            fxaddpar, header, 'BZERO', 0.0, 'Set by MRD_SCALE'

        endif else begin
            ; Random groups.  Get the number of parameters by looking
            ; at the first element in the table.
            nparam = n_elements(table[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, silent=silent)

            for i=0, nparam-1 do begin
                istr = strcompress(string(i+1),/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

; Read a variable length column into a pointer array.
pro mrd_varcolumn, vtype, array, heap, off, siz

    ; Guaranteed to have at least one non-zero length column
    w   = where(siz gt 0)
    nw  = n_elements(w)
    
    if vtype eq 'X' then siz = 1 + (siz-1)/8
    
    siz = siz[w]
    off = off[w]

    unsigned = 0
    if vtype eq '1' then begin
	unsigned = 12
    endif else if vtype eq '2' then begin
	unsigned = 13
    endif else if vtype eq '3' then begin
	unsigned = 15;
    endif
    unsigned = mrd_unsigned_offset(unsigned)
    

    for j=0, nw-1 do begin

        case vtype of

            'L': array[w[j]] = ptr_new(  byte(heap,off[j],siz[j]) )
            'X': array[w[j]] = ptr_new(  byte(heap,off[j],siz[j]) )
            'B': array[w[j]] = ptr_new(  byte(heap,off[j],siz[j]) )
	
            'I': array[w[j]] = ptr_new(  fix(heap, off[j], siz[j]) )
            'J': array[w[j]] = ptr_new(  long(heap, off[j], siz[j]) )
            'K': array[w[j]] = ptr_new(  long64(heap, off[j], siz[j]) )
			  
            'E': array[w[j]] = ptr_new(  float(heap, off[j], siz[j]) )
            'D': array[w[j]] = ptr_new(  double(heap, off[j], siz[j]) )
			  
            'C': array[w[j]] = ptr_new(  complex(heap, off[j], siz[j]) )
            'M': array[w[j]] = ptr_new(  dcomplex(heap, off[j], siz[j]) )
			   
            '1': array[w[j]] = ptr_new(  uint(heap, off[j], siz[j]) )
            '2': array[w[j]] = ptr_new(  ulong(heap, off[j], siz[j]) )
            '3': array[w[j]] = ptr_new(  ulong64(heap, off[j], siz[j]) )
      
        endcase

	; Fix endianness.
        if vtype ne 'B' and vtype ne 'X' and vtype ne 'L' then begin
	    swap_endian_inplace, *array[w[j]],/swap_if_little
        endif

	; Scale unsigneds.
	if unsigned gt 0 then *array[w[j]] = *array[w[j]] - unsigned
	
    endfor
end

; Read a variable length column into a fixed length array.
pro mrd_fixcolumn, vtype, array, heap, off, siz

    w   = where(siz gt 0)
    if w[0] eq -1 then return
    
    nw  = n_elements(w)

    
    if vtype eq 'X' then siz = 1 + (siz-1)/8
    
    siz = siz[w]
    off = off[w]

    for j=0, nw-1 do begin
        case vtype of
            'L': array[0:siz[j]-1,w[j]] = byte(heap,off[j],siz[j])  
            'X': array[0:siz[j]-1,w[j]] = byte(heap,off[j],siz[j])
            'B': array[0:siz[j]-1,w[j]] = byte(heap,off[j],siz[j]) 
	
            'I': array[0:siz[j]-1,w[j]] = fix(heap, off[j], siz[j]) 
            'J': array[0:siz[j]-1,w[j]] = long(heap, off[j], siz[j]) 
            'K': array[0:siz[j]-1,w[j]] = long64(heap, off[j], siz[j]) 
			  
            'E': begin                  ;Delay conversion until after byteswapping to avoid possible math overflow   Feb 2005
	         temp = heap[off[j]: off[j] + 4*siz[j]-1 ]
		 byteorder, temp, /LSWAP, /SWAP_IF_LITTLE	 
	         array[0:siz[j]-1,w[j]] = float(temp,0,siz[j]) 
                 end			  
           'D': begin 
	         temp = heap[off[j]: off[j] + 8*siz[j]-1 ]
		 byteorder, temp, /L64SWAP, /SWAP_IF_LITTLE		 
	         array[0:siz[j]-1,w[j]] = double(temp,0,siz[j]) 
                 end			  
            'C': array[0:siz[j]-1,w[j]] = complex(heap, off[j], siz[j]) 
            'M': array[0:siz[j]-1,w[j]] = dcomplex(heap, off[j], siz[j]) 
			   
            'A': array[w[j]] = string(byte(heap,off[j],siz[j])) 

            '1': array[0:siz[j]-1,w[j]] = uint(heap, off[j], siz[j]) 
            '2': array[0:siz[j]-1,w[j]] = ulong(heap, off[j], siz[j])
            '3': array[0:siz[j]-1,w[j]] = ulong64(heap, off[j], siz[j])
      
        endcase

    endfor

    ; Fix endianness
    if (vtype ne 'A') and (vtype ne 'B') and (vtype ne 'X') and (vtype ne 'L')  and $
       (vtype NE 'D')  and (vtype NE 'E') then begin
	swap_endian_inplace, array, /swap_if_little
    endif

    ; Scale unsigned data
    unsigned = 0
    if vtype eq '1' then begin
	unsigned = 12
    endif else if vtype eq '2' then begin
	unsigned = 13
    endif else if vtype eq '3' then begin
	unsigned = 15;
    endif
    
    if unsigned gt 0 then begin
        unsigned = mrd_unsigned_offset(unsigned)
    endif
    
    if unsigned gt 0 then begin
        for j=0, nw-1 do begin
            array[0:siz[j]-1,w[j]] = array[0:siz[j]-1,w[j]] - unsigned
	endfor
    endif


end
		
; Read the heap area to get the actual values of variable 
; length arrays. 
pro mrd_read_heap, unit, header, range, fnames, fvalues, vcls, vtpes, table, $ 
   structyp, scaling, scales, offsets, status, silent=silent,                $
   columns=columns, rows = rows, pointer_var=pointer_var, fixed_var=fixed_var

    ; 
    ; 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.
    ;
    typstr = 'LXBIJKAEDCM123' 
    prefix = ['bytarr(', 'bytarr(', 'bytarr(', 'intarr(',     $ 
              'lonarr(', 'lon64arr(', 'string(bytarr(', 'fltarr(',         $ 
              'dblarr(', 'cmplxarr(', 'dblarr(2,',            $
	      'uintarr(', 'ulonarr(', 'ulon64arr(']
    
    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 = long64(fxpar(header, 'THEAP')) 
    sz = fxpar(header, 'NAXIS1')*fxpar(header, 'NAXIS2')
    
    if heapoff ne 0 and heapoff lt sz then begin 
        print, 'MRDFITS: ERROR Heap begins within data area' 
        status = -1 
        return 
    endif

    ; Skip to beginning.
    if (heapoff > sz) then begin
        mrd_skip, unit, heapoff-sz
    endif
 
    ; Get the size of the heap. 
    pc = long64(fxpar(header, 'PCOUNT')) 
    if heapoff eq 0 then heapoff = sz 
    hpsiz = pc - (heapoff-sz) 
 
    if (hpsiz gt 0) then heap = bytarr(hpsiz) 
 
 
    ; Read in the heap 
    readu, unit, heap

    ; Skip to the end of the data area.
    skipB = 2880 - (sz+pc) mod 2880
    if skipB ne 2880 then begin
        mrd_skip, unit, skipB
    endif
 
    ; 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]
        curr_col = table.(col)
        vdims[i] = max(curr_col[0,*])
        w = where(curr_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])
    
        ; Note in the following that we ensure that the array is
        ; at least one element long.

        fvalues[col] = prefix[ind] + string((vdims[i] > 1)) + ')'
        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)

    
    ; Should we get rid of empty columns?
    delete = 1
    if keyword_set(pointer_var) then delete = pointer_var eq 1

    if delete then begin
	
        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

        ; Check if all columns have been deleted...
        wx = where(w gt 0)
        if (wx[0] eq -1) then begin
            if not keyword_set(silent) then begin
                print, 'MRDFITS: All columns have been deleted'
            endif
	    table = 0
	    return
        endif
    

        ; Get rid of unused columns.
        corres = corres[wx]
        fnames = fnames[wx] 
        fvalues = fvalues[wx]
        scales = scales[wx]
        offsets = offsets[wx]

        wx = where(vdims gt 0)
    
        if (wx[0] eq -1) then begin
            vcols=[-9999]
	    x=temporary(vtypes)
	    x=temporary(vdims)
        endif else begin 
            vcols = vcols[wx]
            vtypes = vtypes[wx]
            vdims = vdims[wx]
        endelse
    endif

    if not keyword_set(pointer_var) then begin
        ; 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' ; <-- Originally, '0l'; breaks under the virtual machine!
                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
	
    endif else begin
	
        ; We'll just read data into pointer arrays.
	for i=0,n_elements(lencols)-1 do begin
	    col = lencols[i]
	    if vtpes[col] eq 'A' then begin
	        fvalues[col] = '" "'
	    endif else begin
	        fvalues[col] = 'ptr_new()'
	    endelse
	endfor
	
    endelse
	


    ; 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, $
                            silent=silent)
    endif else begin
        tablex = mrd_struct(fnames, fvalues, n_elements(table), silent=silent)
    endelse


    if N_elements(rows) EQ 0 then nrow = range[1]-range[0]+1 $
                             else nrow = N_elements(rows)
    
    ; 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_elements(fnames)-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
                        
                ; If only one row in table, then
                ; IDL will return curr_col as one-dimensional.
                ; Since this is a variable length pointer column we
                ; know that the dimension of the column is 2.
                    curr_col = table.(col)
		
                if (nrow eq 1) then curr_col = reform(curr_col,2,1)
                siz = curr_col[0,*] 
                off = curr_col[1,*] 
                    
                ; Now process each type.
                    curr_colx = tablex.(i)
                    sz = size(curr_colx)
                    if (sz[0] lt 2) then begin
                        curr_colx = reform(curr_colx, 1, n_elements(curr_colx), /overwrite)
                    endif
                           
                    
                ; As above we have to worry about IDL truncating
                ; dimensions.  This can happen if either
                ; nrow=1 or the max dimension of the column is 1.
                    

                    sz = size(tablex.(i))
 
                nel = sz[sz[0]+2]
                if nrow eq 1 and nel eq 1 then begin
                    curr_colx = make_array(1,1,value=curr_colx)
                endif else if nrow eq 1 then begin
                    curr_colx = reform(curr_colx,[nel, 1], /overwrite)
                endif else if nel eq 1 then begin
                    curr_colx = reform(curr_colx,[1, nrow], /overwrite)
                endif

		vtype = vtypes[vc]
		varying = 0
		if n_elements(lencols) gt 0 then begin
		    varying = where(lencols eq col)
		    if varying[0] eq -1 then varying=0 else varying=1
		endif
		
		if varying and keyword_set(pointer_var) and vtype ne 'A' then begin
		    mrd_varcolumn, vtype, curr_colx, heap, off, siz
		endif else begin
		    mrd_fixcolumn, vtype, curr_colx, heap, off, siz
		endelse


                
                if nel eq 1 and nrow eq 1 then begin
                    curr_colx = curr_colx[0]
                endif else if nrow eq 1 then begin
                    curr_colx = reform(curr_colx, nel, /overwrite)
                endif else if nel eq 1 then begin
                    curr_colx = reform(curr_colx, nrow, /overwrite)
                endif

                    sz = size(curr_colx)
                    if sz[1] eq 1 then begin
                         sz_tablex = size(tablex.(i))
                         sdimen = sz_tablex[1:sz_tablex[0]]
                         tablex.(i) = reform(curr_colx,sdimen)
                    endif else begin
                        tablex.(i) = curr_colx
                    endelse
                 
            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.
                xx = table.(ncol)
                tablex.(i) = reform(xx[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 

; Read in the binary table information. 
pro mrd_read_table, unit, range, rsize, structyp, nrows, nfld, typarr, table, rows = rows
 
    ; 
    ; 
    ; 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.
    ; 

    if range[0] gt 0 then mrd_skip, unit, rsize*range[0]
    readu,unit, table
    if N_elements(rows) GT 0 then table = table[rows- range[0]]

    ; Move to the beginning of the heap -- we may have only read some rows of
    ; the data.
    if range[1] lt nrows-1 then begin
        skip_dist = (nrows-range[1]-1)*rsize
        mrd_skip, unit, skip_dist
    endif

    

    ; If necessary then convert to native format.
    if byte(1L, 0,1) EQ 1 then begin
	
        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 'K' then byteorder, fld, /l64swap
            if typ eq 'E' or typarr[i] eq 'C' then $
	                       byteorder, fld, /LSWAP
	    
            if typ eq 'D' or typarr[i] eq 'M' then byteorder, fld, /L64SWAP
	    
            if n_elements(fld) gt 1 then begin
		
                    table.(i) = fld
           endif else begin
                    table.(i) = fld[0]
            endelse
  nxtfld: 
        endfor 
    endif

    ; Handle unsigned fields.
    for i=0, nfld-1 do begin

 
	    type = mrd_unsignedtype(table.(i))

	    if type gt 0 then begin
	        table.(i) = table.(i) - mrd_unsigned_offset(type)
	    endif
	    
	
    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
 
; Define a structure to hold a FITS binary table. 
pro mrd_table, header, structyp, use_colnum,           $ 
    range, rsize, table, nrows, nfld, typarr, fnames, fvalues,   $ 
    vcls, vtpes, scales, offsets, scaling, status, rows = rows, $
    silent=silent, columns=columns, no_tdim=no_tdim, $
    alias=alias, unsigned=unsigned
 
    ; 
    ; 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.

    table = 0

    types =  ['L', 'X', 'B', 'I', 'J', 'K', 'A', 'E', 'D', 'C', 'M', 'P']
    arrstr = ['bytarr(', 'bytarr(', 'bytarr(', 'intarr(', 'lonarr(', 'lon64arr(',      $ 
              'string(replicate(32b,', 'fltarr(', 'dblarr(', 'complexarr(',            $ 
              'dcomplexarr(', 'lonarr(2*']
    bitpix = [  0,   0,   0,  16,  32,  64,   0,  0,   0,   0,   0,   0]

    sclstr = ["'T'", '0B', '0B', '0', '0L', '0LL', '" "', '0.', '0.d0', 'complex(0.,0.)', $ 
              'dcomplex(0.d0,0.d0)', 'lonarr(2)']

    unsarr = ['', '', '', 'uintarr(', 'ulonarr(', 'ulon64arr('];
    unsscl = ['', '', '', '0U',        '0UL',      '0ULL']
 

    status = 0 

; NEW WAY: E.S.S.

    ;; get info from header. Using vectors is much faster
    ;; when there are many columns

    mrd_fxpar, header, xten, nfld, nrow, rsize, fnames, fforms, scales, offsets
    nnames = n_elements(fnames)

    ;; nrow will change later
    nrows = nrow

    ;; Use scale=1 if not found
    if nnames GT 0 then begin
      wsc=where(scales EQ 0.0d,nwsc)
      IF nwsc NE 0 THEN scales[wsc] = 1.0d
    endif

    xten = strtrim(xten,2)
    if xten ne 'BINTABLE' and xten ne 'A3DTABLE' then begin 
        print, 'MRDFITS: ERROR - Header is not from binary table.' 
        nfld = 0 & status = -1 
        return 
    endif 
 
    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

    if N_elements(rows) EQ 0 then nrowp  = nrow else begin 
          bad = where((rows LT range[0]) or (rows GT range[1]), Nbad)
          if Nbad GT 0 then begin 
             print,'MRDFITS: Row numbers must be between 0 and ' + $
                    strtrim(nrow-1,2)      
             status = -1
             return
           endif
           nrowp = N_elements(rows)
    endelse
;    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)
 
    fnames2 = strarr(nfld)

    for i=0, nfld-1 do begin
	
        istr = strcompress(string(i+1), /remo)

        fname = fnames[i]

        ;; check for a name conflict
        fname = mrd_dofn(fname, i+1, use_colnum, alias=alias)
	
        ;; check for a name conflict
        fname = mrd_chkfn(fname, fnames2, i)

        ;; copy in the valid name
        fnames[i] = fname
        ;; for checking conflicts
        fnames2[i] = fname

;        fname = fxpar(header, 'TTYPE' +  istr) 
;        fform = strtrim( fxpar(header, 'TFORM' +   istr),2)
	
;        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, alias=alias) 
;        fname = mrd_chkfn(fname, fnames, i)
	
;        fnames[i] = fname
	
        fform = fforms[i]

        mrd_doff, fform, dim, ftype
        
        ; Treat arrays of length 1 as scalars.
        if dim eq 1 then begin
            dim = 0
        endif else if dim EQ -1 then begin 
            dimfld[i] = -1
        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 begin
            dim = (dim+7)/8 
            dimfld[i] = strtrim(string(dim),2)
        endif
         
        ; 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('LXBIJKAEDCM', vf) eq -1 then begin 
                print, 'MRDFITS: Invalid type for variable array column '+string(i+1) 
                status = -1 
                return 
            endif 

            vcls[i] = 1
	    
;	    xscale =fxpar(header,'TSCAL'+istr,count=cnt)
;	    if cnt eq 0 then xscale = 1
	    
;	    xunsigned = mrd_chkunsigned(bitpix[ppos], xscale,       $
;				       fxpar(header, 'TZERO'+istr), $
;				       unsigned=unsigned)

	    xunsigned = mrd_chkunsigned(bitpix[ppos], scales[i],       $
				       offsets[i], $
				       unsigned=unsigned)

	    if (xunsigned) then begin
		
		if      vf eq 'I' then vf = '1' $
		else if vf eq 'J' then vf = '2' $
		else if vf eq 'K' then vf = '3'
		
	    endif
							   
            vtpes[i] = vf
            dim = 0
                         
        endif 
         

        for j=0, n_elements(types) - 1 do begin
	    
            if ftype eq types[j] then begin

;		xscale = fxpar(header, 'TSCAL'+istr, count=cnt)
;		if cnt eq 0 then xscale = 1
		
;		xunsigned = mrd_chkunsigned(bitpix[j], xscale, $
;					   fxpar(header, 'TZERO'+istr), $
;					   unsigned=unsigned)

                xunsigned = mrd_chkunsigned(bitpix[j], scales[i], $
                                            offsets[i], $
                                            unsigned=unsigned)

		if xunsigned then begin
		    fxaddpar, header, 'TZERO'+istr, 0, 'Modified by MRDFITS V'+mrd_version()
                    offsets[i] = 0 ;; C. Markwardt Aug 2007 - reset to zero so offset is not applied twice'
	        endif

                if dim eq 0 then begin

		    if xunsigned then begin
			fvalues[i] = unsscl[j]
		    endif else begin
                        fvalues[i] = sclstr[j]
		    endelse
		    
                endif else begin

		    if xunsigned then begin
			line = unsarr[j]
		    endif else begin
			line = arrstr[j]
		    endelse
		    
                    line = line + 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

    zero = where(long(dimfld) LT 0L, N_zero)
    if N_zero GT 0 then begin
	
        if N_zero Eq nfld then begin
            print,'MRDFITS: Error - All fields have zero length'
            return
        endif
	
        for i=0, N_zero-1 do begin
	    print,'MRDFITS: Table column ' + fnames[zero[i]] + ' has zero length'
	endfor
	
        nfld    = nfld - N_zero
        good    = where(dimfld GE 0)
        fnames  = fnames[good]
        fvalues = fvalues[good]
        typarr = typarr[good]      ;Added 2005-1-6   (A.Csillaghy)

    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,  silent=silent )
	
    endif else begin
	
        table = mrd_struct(fnames, fvalues, nrow, silent=silent )
	
    endelse

    if not keyword_set(silent) then begin
        print, 'MRDFITS: Binary table. ',strcompress(string(nfld)), ' columns by ',  $
          strcompress(string(nrowp)), ' 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,                      $
        error_action = error_action,            $
 	compress=compress,                      $
	alias=alias,                            $
        rows = rows,                        $
	unsigned=unsigned,                      $
	version=version,                        $
	pointer_var=pointer_var,                $
	fixed_var=fixed_var,                    $
        status=status, extnum = extnum

    compile_opt idl2    
    ;   Let user know version if MRDFITS being used.
    if keyword_set(version) then begin
        print,'MRDFITS: Version '+mrd_version()+' Dec 12, 2007'
    endif
    
    ;
    ;  Can't use keyword_set since default is 2, not 0.

    if n_elements(error_action) eq 0 then begin
        error_action = 2
    endif
    
    on_error, error_action
   

   
    ; Check positional arguments.

    if n_params() le 0  or n_params() gt 3 then begin
	if keyword_set(version) then return, 0
        print, 'MRDFITS: Usage'
        print, '   a=mrdfits(file/unit, [exten_no/exten_name, header], /version $'
        print, '       /fscale, /dscale, /unsigned, /use_colnum, /silent    $'
        print, '       range=, rows= , structyp=, columns=, $'
	print, '       /pointer_var, /fixed_var, error_action=, status= )'
        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 N_elements(rows) GT 0 then begin
        range1 = min(rows,max=range2)
        range = [range1,range2]
    endif
    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 begin
	arange = [-1,-1]
    endelse

    arange = long(arange)

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

    if (N_elements(file) GT 1 ) then begin
        print, 'MRDFITS: Vector input not supported'
        return, 0
    endif

    inputUnit = 0
   
    dtype = size(file,/type)
    if dtype gt 0 and dtype lt 4 then begin    ;File unit number specified
	
        inputUnit = 1
        unit = file
    
        if fxmove(unit,extension) lt 0 then begin
            return, -1
        endif
    
    endif else begin                         ;File name specified
        unit = fxposit(file, extension, compress=compress, $
	               /readonly,extnum=extnum, errmsg= errmsg)

        if unit lt 0 then begin
            message, 'File access error',/CON
	    if errmsg NE '' then message,errmsg,/CON
            status = -1
            return, 0
        endif
    endelse

    if eof(unit) then begin
        print,'MRDFITS: Extension past EOF'
	if inputUnit eq 0 then free_lun,unit 
	status = -2
	return, 0
    endif

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

    ; If this is primary array then XTENSION will have value
    ; 0 which will be converted by strtrim to '0'

    xten = strtrim( fxpar(header,'XTENSION'), 2)
    if xten eq '0' or 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 
        message, 'Unable to process extension type:', xten,/CON
	if inputUnit eq 0 then free_lun,unit
	status = -1
        return, 0
    endelse

    scaling = keyword_set(fscale) or keyword_set(dscale)

    if type eq 0 then begin

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

        ;*** ASCII tables.
        
        mrd_ascii, header, structyp, use_colnum,                              $
            arange, table, nbytes, nrows, nfld, rows=rows,                    $
            typarr, posarr, lenarr, nullarr, fnames, fvalues,                 $
            scales, offsets, scaling, status, silent=silent,                  $
            columns=columns, alias=alias
        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,  rows= rows
              
            ;*** 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, silent=silent
        endif
        
    endif else begin

        ; *** Binary tables.

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

        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,  rows = rows, $
              structyp, nrows, nfld, typarr, table

            if status ge 0 and keyword_set(columns) then begin
        
                ;*** Extract desired columns.
                mrd_columns, table, columns, fnames, fvalues,                  $
                  vcls, vtpes, scales, offsets, scaling, structyp=structyp,    $
                  silent=silent
	    
	    endif
         
     
            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, pointer_var=pointer_var, fixed_var=fixed_var, rows= rows
		
	    endif else begin

	        ; Skip remainder of last data block
	        sz = long64(fxpar(header, 'NAXIS1'))* $
                     long64(fxpar(header,'NAXIS2')) +  $
		       long64(fxpar(header, 'PCOUNT'))
	        skipB = 2880 - sz mod 2880
	        if (skipB ne 2880) then mrd_skip, unit, skipB
            endelse
		     
        endif

    endelse


    ; Don't tie up a unit number that we allocated in this routine.
    if unit gt 0 and inputUnit eq 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, silent=silent
    endif

    ; All done. Check the status to see if we ran into problems on the way.
    
    if status ge 0 then return, table else return,0
 
end