;+
; 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,[Extension, Header],
; /FSCALE , /DSCALE , /UNSIGNED,
; ALIAS=strarr[2,n], /USE_COLNUM,
; /NO_TDIM, ROWS = [a,b,...], $
; /POINTER_VAR, /FIXED_VAR,
; 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.
; Note that the file name may be of the form
; name.gz or name.Z on UNIX systems. If so
; the file will be dynamically decompressed.
; 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
;
; Extension= 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.
;
; 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:
; 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.
;
; As of version 2.5 MRDFITS support 64 bit integer data types.
; These are not standard FITS. 64 bit data also requires
; IDL version 5.2 or greater.
;
;
; 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, IEEE_TO_HOST, FXPOSIT, FXMOVE(), IS_IEEE_BIG()
; 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
; Converted to IDL V5.0 W. Landsman 2-Nov-1997
; 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
;-
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
sz = size(data)
type = sz[sz[0]+1]
if type eq 12 or type eq 13 or type eq 15 then begin
return, type
endif else begin
return, 0
endelse
end
; Return the currrent version string for MRDFITS
function mrd_version
return, '2.9g'
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
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
;
; Unit Unit to read data from.
; Table Table/array to read information into.
;
; If necessary skip to beginning of desired data.
if range[0] gt 0 then mrd_skip, unit, range[0]*rsize
if rsize eq 0 then return
readu, unit, table
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
if not is_ieee_big() then ieee_to_host, table
; Fix offset for unsigned data
type = mrd_unsignedtype(table)
if type gt 0 then begin
table = table - mrd_unsigned_offset(type)
endif
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 scale = 1
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
ieee_to_host, *array[w[j]]
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, /XDRTOF
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, /XDRTOD
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
ieee_to_host, array
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'
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)
is_ieee = is_ieee_big()
; 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 not is_ieee_big() 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, /xdrtof
if typ eq 'D' or typarr[i] eq 'M' then begin
ieee_to_host, fld
endif
if n_elements(fld) gt 1 then begin
table.(i) = fld
endif else begin
table.(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()
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
; Let user know version if MRDFITS being used.
if keyword_set(version) then begin
print,'MRDFITS: Version '+mrd_version()+' Feb 07, 2005'
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, [extension, 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.
sz = size(file)
if (sz[0] ne 0) then begin
print, 'MRDFITS: Vector input not supported'
return, 0
endif
inputUnit = 0
if sz[1] gt 0 and sz[1] lt 4 then begin
inputUnit = 1
unit = file
if fxmove(unit,extension) lt 0 then begin
return, -1
endif
endif else begin
unit = fxposit(file, extension, compress=compress, /readonly)
if unit lt 0 then begin
print, 'MRDFITS: File access error'
status = -1
return, 0
endif
endelse
if eof(unit) then begin
print,'MRDFITS: Extension past EOF'
if inputUnit eq 0 then begin
free_lun, unit
endif
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 begin
free_lun, unit
endif
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 begin
type = 0
endif else if xten eq 'TABLE' then begin
type = 1
endif else if xten eq 'BINTABLE' or xten eq 'A3DTABLE' then begin
type = 2
endif else begin
print, 'MRDFITS: Unable to process extension type:', xten
if inputUnit eq 0 then begin
free_lun,unit
endif
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 $
mrd_read_image, unit, arange, maxd, rsize, table, rows = rows
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
Return to the shapelets web page or the code help menu.
| Last modified on 2nd Mar 2009 by Richard Massey. |