;+ ; 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