; the following are from "AstroLib" ;+ ; NAME: ; MWRFITS ; PURPOSE: ; Write all standard FITS data types from input arrays or structures. ; ; CALLING SEQUENCE: ; MWRFITS, Input, Filename, [Header], ; /LSCALE , /ISCALE, /BSCALE, ; /USE_COLNUM, /Silent, /Create, /No_comment, /Version, $ ; Alias=, /ASCII, Separator=, Terminator=, Null=, ; /Logical_cols, /Bit_cols, /Nbit_cols, ; Group=, Pscale=, Pzero= ; ; INPUTS: ; Input = Array or structure to be written to FITS file. ; ; -When writing FITS primary data or image extensions ; input should be an array. ; --If data is to be grouped ; the Group keyword should be specified to point to ; a two dimensional array. The first dimension of the ; Group array will be PCOUNT while the second dimension ; should be the same as the last dimension of Input. ; --If Input is undefined, then a dummy primary dataset ; or Image extension is created [This might be done, e.g., ; to put appropriate keywords in a dummy primary ; HDU]. ; ; -When writing an ASCII table extension, Input should ; be a structure array where no element of the structure ; is a structure or array (except see below). ; --A byte array will be written as A field. No checking ; is done to ensure that the values in the byte field ; are valid ASCII. ; --Complex numbers are written to two columns with '_R' and ; '_I' appended to the TTYPE fields (if present). The ; complex number is enclosed in square brackets in the output. ; --Strings are written to fields with the length adjusted ; to accommodate the largest string. Shorter strings are ; blank padded to the right. ; ; -When writing a binary table extension, the input should ; be a structure array with no element of the structure ; being a substructure. ; ; If a structure is specified on input and the output ; file does not exist or the /CREATE keyword is specified ; a dummy primary HDU is created. ; ; Filename = String containing the name of the file to be written. ; By default MWRFITS appends a new extension to existing ; files which are assumed to be valid FITS. The /CREATE ; keyword can be used to ensure that a new FITS file ; is created even if the file already exists. ; ; OUTPUTS: ; ; OPTIONAL INPUTS: ; Header = Header should be a string array. Each element of the ; array is added as a row in the FITS header. No ; parsing is done of this data. MWRFITS will prepend ; required structural (and, if specified, scaling) ; keywords before the rows specified in Header. ; Rows describing columns in the table will be appended ; to the contents of Header. ; Header lines will be extended or truncated to ; 80 characters as necessary. ; If Header is specified then on return Header will have ; the header generated for the specified extension. ; ; OPTIONAL INPUT KEYWORDS: ; ALias= Set up aliases to convert from the IDL structure ; to the FITS column name. The value should be ; a STRARR(2,*) value where the first element of ; each pair of values corresponds to a column ; in the structure and the second is the name ; to be used in the FITS file. ; The order of the alias keyword is compatible with ; use in MRDFITS. ; ASCII - Creates an ASCII table rather than a binary table. ; This keyword may be specified as: ; /ASCII - Use default formats for columns. ; ASCII='format_string' allows the user to specify ; the format of various data types such using the following ; syntax 'column_type:format, column_type:format'. E.g., ; ASCII='A:A1,I:I6,L:I10,B:I4,F:G15.9,D:G23.17,C:G15.9,M:G23.17' ; gives the default formats used for each type. The TFORM ; fields for the real and complex types indicate will use corresponding ; E and D formats when a G format is specified. ; Note that the length of the field for ASCII strings and ; byte arrays is automatically determined for each column. ; BIT_COLS= An array of indices of the bit columns. The data should ; comprise a byte array with the appropriate dimensions. ; If the number of bits per row (see NBIT_COLS) ; is greater than 8, then the first dimension of the array ; should match the number of input bytes per row. ; BSCALE Scale floats, longs, or shorts to unsigned bytes (see LSCALE) ; CREATE If this keyword is non-zero, then a new FITS file will ; be created regardless of whether the file currently ; exists. Otherwise when the file already exists, ; a FITS extension will be appended to the existing file ; which is assumed to be a valid FITS file. ; GROUP= This keyword indicates that GROUPed FITS data is to ; be generated. ; Group should be a 2-D array of the appropriate output type. ; The first dimension will set the number of group parameters. ; The second dimension must agree with the last dimension ; of the Input array. ; ISCALE Scale floats or longs to short integer (see LSCALE) ; LOGICAL_COLS= An array of indices of the logical column numbers. ; These should start with the first column having index 0. ; The structure element should either be an array of characters ; with the values 'T' or 'F', or an array of bytes having the ; values byte('T'), byte('F') or 0b. The use of bytes allows ; the specification of undefined values (0b). ; LSCALE Scale floating point numbers to long integers. ; This keyword may be specified in three ways. ; /LSCALE (or LSCALE=1) asks for scaling to be automatically ; determined. LSCALE=value divides the input by value. ; I.e., BSCALE=value, BZERO=0. Numbers out of range are ; given the value of NULL if specified, otherwise they are given ; the appropriate extremum value. LSCALE=(value,value) ; uses the first value as BSCALE and the second as BZERO ; (or TSCALE and TZERO for tables). ; NBIT_COLS= The number of bits actually used in the bit array. ; This argument must point to an array of the same dimension ; as BIT_COLS. ; NO_TYPES If the NO_TYPES keyword is specified, then no TTYPE ; keywords will be created for ASCII and BINARY tables. ; No_comment Do not write comment keywords in the header ; NULL= Value to be written for integers/strings which are ; undefined or unwritable. ; PSCALE= An array giving scaling parameters for the group keywords. ; It should have the same dimension as the first dimension ; of Group. ; PZERO= An array giving offset parameters for the group keywords. ; It should have the same dimension as the first dimension ; of Group. ; Separator= This keyword can be specified as a string which will ; be used to separate fields in ASCII tables. By default ; fields are separated by a blank. ; SILENT Suppress informative messages. Errors will still ; be reported. ; Terminator= This keyword can be specified to provide a string which ; will be placed at the end of each row of an ASCII table. ; No terminator is used when not specified. ; If a non-string terminator is specified (including ; when the /terminator form is used), a new line terminator ; is appended. ; USE_COLNUM When creating column names for binary and ASCII tables ; MWRFITS attempts to use structure field name ; 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 version number of MWRFITS. ; ; EXAMPLE: ; Write a simple array: ; a=fltarr(20,20) ; mwrfits,a,'test.fits' ; ; Append a 3 column, 2 row, binary table extension to file just created. ; a={name:'M31', coords:(30., 20.), distance:2} ; a=replicate(a, 2); ; mwrfits,a,'test.fits' ; ; Now add on an image extension: ; a=lonarr(10,10,10) ; hdr=("COMMENT This is a comment line to put in the header", $ ; "MYKEY = "Some desired keyword value") ; mwrfits,a,'test.fits',hdr ; ; RESTRICTIONS: ; (1) Variable length columns are not supported for anything ; other than simple types (byte, int, long, float, double). ; (2) Empty strings are converted to 1 element blank strings (because ; IDL refuses to write an emptry string (0b) from a structure) ; NOTES: ; This multiple format FITS writer is designed to provide a ; single, simple interface to writing all common types of FITS data. ; Given the number of options within the program and the ; variety of IDL systems available it is likely that a number ; of bugs are yet to be uncovered. If you find an anomaly ; please send a report to: ; Tom McGlynn ; NASA/GSFC Code 660.2 ; tam@silk.gsfc.nasa.gov (or 301-286-7743) ; ; PROCEDURES USED: ; FXPAR(), FXADDPAR ; MODIfICATION HISTORY: ; Version 0.9: By T. McGlynn 1997-07-23 ; Initial beta release. ; Dec 1, 1997, Lindler, Modified to work under VMS. ; Version 0.91: T. McGlynn 1998-03-09 ; Fixed problem in handling null primary arrays. ; Reconverted to IDL 5.0 format using IDLv4_to_v5 ; Version 0.92: T. McGlynn 1998-09-09 ; Add no_comment flag and keep user comments on fields. ; Fix handling of bit fields. ; Version 0.93: T. McGlynn 1999-03-10 ; Fix table appends on VMS. ; Version 0.93a W. Landsman/D. Schlegel ; Update keyword values in chk_and_upd if data type has changed ; Version 0.94: T. McGlynn 2000-02-02 ; Efficient processing of ASCII tables. ; Use G rather than E formats as defaults for ASCII tables ; and make the default precision long enough that transformations ; binary to/from ASCII are invertible. ; Some loop indices made long. ; Fixed some ends to match block beginnings. ; Version 0.95: T. McGlynn 2000-11-06 ; Several fixes to scaling. Thanks to David Sahnow for ; documenting the problems. ; Added PCOUNT,GCOUNT keywords to Image extensions. ; Version numbers shown in SIMPLE/XTENSION comments ; Version 0.96: T. McGlynn 2001-04-06 ; Changed how files are opened to handle ~ consistently ; Version 1.0: T. McGlynn 2001-12-04 ; Unsigned integers, ; 64 bit integers. ; Aliases ; Variable length arrays ; Some code cleanup ; Version 1.1: T. McGlynn 2002-2-18 ; Fixed major bug in processing of unsigned integers. ; (Thanks to Stephane Beland) ; Version 1.2: Stephane Beland 2003-03-17 ; Fixed problem in creating dummy dataset when passing undefined ; data, caused by an update to FXADDPAR routine. ; Version 1.2.1 Stephane Beland 2003-09-10 ; Exit gracefully if write priveleges unavailable ; Version 1.3 Wayne Landsman 2003-10-24 ; Don't use EXECUTE() statement if on a virtual machine ; Version 1.3a Wayne Landsman 2004-5-21 ; Fix for variable type arrays ; Version 1.4 Wayne Landsman 2004-07-16 ; Use STRUCT_ASSIGN when modifying structure with pointer tags ; Version 1.4a Wayne Landsman 2005-01-03 ; Fix writing of empty strings in binary tables ; Version 1.4b Wayne Landsman 2006-02-23 ; Propagate /SILENT keyword to mwr_tablehdr ; Version 1.5 Wayne Landsman 2006-05-24 ; Open file using /SWAP_IF_LITTLE_ENDIAN keyword ; Convert empty strings to 1 element blank strings before writing ; Version 1.5a Wayne Landsman 2006-06-29 ; Fix problem introduced 2006-05-24 with multidimensional strings ; Version 1.5b K. Tolbert 2006-06-29 ; Make V1.5a fix work pre-V6.0 ; Version 1.5c I.Evans/W.Landsman 2006-08-08 ; Allow logical columns to be specified as bytes ; Version 1,5d K. Tolbert 2006-08-11 ; Make V1.5a fix work for scalar empty string ; Version 1.6 W. Landsman 2006-09-22 ; Assume since V5.5, remove VMS support ;- ; What is the current version of this program? function mwr_version return, '1.6' end ; Find the appropriate offset for a given unsigned type ; or just return 0 if the type is not unsigned. function mwr_unsigned_offset, type if (type eq 12) then begin return, uint(32768) endif else if (type eq 13) then begin return, ulong('2147483648') endif else if (type eq 15) then begin return, ulong64('9223372036854775808') endif return, 0 end ; Add a keyword as non-destructively as possible to a FITS header pro chk_and_upd, header, key, value, comment xcomm = "" if n_elements(comment) gt 0 then xcomm = comment if n_elements(header) eq 0 then begin fxaddpar, header, key, value, xcomm endif else begin oldvalue = fxpar(header, key, count=count, comment=oldcomment) if (count eq 1) then begin qchange = 0 ; Set to 1 if either the type of variable or its ; value changes. size1 = size(oldvalue,/type) & size2 = size(value,/type) if size1 NE size2 then qchange = 1 $ else if (oldvalue ne value) then qchange = 1 if (qchange) then begin if n_elements(oldcomment) gt 0 then xcomm = oldcomment[0] fxaddpar, header, key, value, xcomm endif endif else begin fxaddpar, header, key, value, xcomm endelse endelse end ; Get the column name appropriate for a given tag function mwr_checktype, tag, alias=alias if not keyword_set(alias) then return, tag sz = size(alias) ; 1 or 2 D string array with first dimension of 2 if (sz[0] eq 1 or sz[1] eq 2) and sz[1] eq 2 and sz[sz[0]+1] eq 7 then begin w = where(tag eq alias[0,*]) if (w[0] eq -1) then begin return, tag endif else begin return, alias[1,w[0]] endelse endif else begin print,'MWRFITS: Warning: Alias values not strarr(2) or strarr(2,*)' endelse return, tag end ; Create an ASCII table pro mwr_ascii, input, siz, lun, bof, header, $ ascii=ascii, $ null=null, $ use_colnum = use_colnum, $ lscale=lscale, iscale=iscale, $ bscale=bscale, $ no_types=no_types, $ separator=separator, $ terminator=terminator, $ no_comment=no_comment, $ silent=silent, $ alias=alias ; Write the header and data for a FITS ASCII table extension. types= ['A', 'I', 'L', 'B', 'F', 'D', 'C', 'M', 'K'] formats=['A1', 'I6', 'I10', 'I4', 'G15.9','G23.17', 'G15.9', 'G23.17','I20'] lengths=[1, 6, 10, 4, 15, 23, 15, 23, 20] ; Check if the user is overriding any default formats. sz = size(ascii) if sz[0] eq 0 and sz[1] eq 7 then begin ascii = strupcase(strcompress(ascii,/remo)) for i=0, n_elements(types)-1 do begin p = strpos(ascii,types[i]+':') if p ge 0 then begin q = strpos(ascii, ',', p+1) if q lt p then q = strlen(ascii)+1 formats[i] = strmid(ascii, p+2, (q-p)-2) len = 0 reads, formats[i], len, format='(1X,I)' lengths[i] = len endif endfor endif i0 = input[0] ntag = n_tags(i0) tags = tag_names(i0) ctypes = lonarr(ntag) strmaxs = lonarr(ntag) if not keyword_set(separator) then separator=' ' slen = strlen(separator) offsets = 0 tforms = '' ttypes = '' offset = 0 totalFormat = "" xsep = ""; for i=0, ntag-1 do begin totalFormat = totalFormat + xsep; sz = size(i0.(i)) if sz[0] ne 0 and (sz[sz[0]+1] ne 1) then begin print, 'MWRFITS Error: ASCII table cannot contain arrays' return endif ctypes[i] = sz[1] xtype = mwr_checktype(tags[i], alias=alias) ttypes = [ttypes, xtype+' '] if sz[0] gt 0 then begin ; Byte array to be handled as a string. nelem = sz[sz[0]+2] ctypes[i] = sz[sz[0]+1] tf = 'A'+strcompress(string(nelem)) tforms = [tforms, tf] offsets = [offsets, offset] totalFormat = totalFormat + tf offset = offset + nelem endif else if sz[1] eq 7 then begin ; Use longest string to get appropriate size. strmax = max(strlen(input.(i))) strmaxs[i] = strmax tf = 'A'+strcompress(string(strmax), /remo) tforms = [tforms, tf] offsets = [offsets, offset] totalFormat = totalFormat + tf ctypes[i] = 7 offset = offset + strmax endif else if sz[1] eq 6 or sz[1] eq 9 then begin ; Complexes handled as two floats. offset = offset + 1 if sz[1] eq 6 then indx = where(types eq 'C') if sz[1] eq 9 then indx = where(types eq 'M') indx = indx[0] fx = formats[indx] if (strmid(fx, 0, 1) eq "G" or strmid(fx, 0, 1) eq "g") then begin if (sz[1] eq 6) then begin fx = "E"+strmid(fx,1, 99) endif else begin fx = "D"+strmid(fx,1, 99) endelse endif tforms = [tforms, fx, fx] offsets = [offsets, offset, offset+lengths[indx]+1] nel = n_elements(ttypes) ttypes = [ttypes[0:nel-2], xtype+'_R', xtype+'_I'] offset = offset + 2*lengths[indx] + 1 totalFormat = totalFormat + '"[",'+formats[indx]+',1x,'+formats[indx]+',"]"' offset = offset+1 endif else begin if sz[1] eq 1 then indx = where(types eq 'B') $ else if sz[1] eq 2 or sz[1] eq 12 then indx = where(types eq 'I') $ else if sz[1] eq 3 or sz[1] eq 13 then indx = where(types eq 'L') $ else if sz[1] eq 4 then indx = where(types eq 'F') $ else if sz[1] eq 5 then indx = where(types eq 'D') $ else if sz[1] eq 14 or sz[1] eq 15 then indx = where(types eq 'K') $ else begin print, 'MWRFITS Error: Invalid type in ASCII table' return endelse indx = indx[0] fx = formats[indx] if (strmid(fx, 0, 1) eq 'G' or strmid(fx, 0, 1) eq 'g') then begin if sz[1] eq 4 then begin fx = 'E'+strmid(fx, 1, 99) endif else begin fx = 'D'+strmid(fx, 1, 99) endelse endif tforms = [tforms, fx] offsets = [offsets, offset] totalFormat = totalFormat + formats[indx] offset = offset + lengths[indx] endelse if i ne ntag-1 then begin offset = offset + slen endif xsep = ", '"+separator+"', " endfor if keyword_set(terminator) then begin sz = size(terminator); if sz[0] ne 0 or sz[1] ne 7 then begin terminator= string(10B) endif endif if keyword_set(terminator) then offset = offset+strlen(terminator) ; Write required FITS keywords. chk_and_upd, header, 'XTENSION', 'TABLE', 'ASCII table extension written by MWRFITS '+mwr_version() chk_and_upd, header, 'BITPIX', 8,'Required Value: ASCII characters' chk_and_upd, header, 'NAXIS', 2,'Required Value' chk_and_upd, header, 'NAXIS1', offset, 'Number of characters in a row' chk_and_upd, header, 'NAXIS2', n_elements(input), 'Number of rows' chk_and_upd, header, 'PCOUNT', 0, 'Required value' chk_and_upd, header, 'GCOUNT', 1, 'Required value' chk_and_upd, header, 'TFIELDS', n_elements(ttypes)-1, 'Number of fields' ; Recall that the TTYPES, TFORMS, and OFFSETS arrays have an ; initial dummy element. ; Write the TTYPE keywords. if not keyword_set(no_types) then begin for i=1, n_elements(ttypes)-1 do begin key = 'TTYPE'+ strcompress(string(i),/remo) if keyword_set(use_colnum) then begin value = 'C'+strcompress(string(i),/remo) endif else begin value = ttypes[i]+' ' endelse chk_and_upd, header, key, value endfor if (not keyword_set(no_comment)) then begin fxaddpar, header, 'COMMENT', ' ', before='TTYPE1' fxaddpar, header, 'COMMENT', ' *** Column names ***', before='TTYPE1' fxaddpar, header, 'COMMENT', ' ', before='TTYPE1' endif endif ; Write the TBCOL keywords. for i=1, n_elements(ttypes)-1 do begin key= 'TBCOL'+strcompress(string(i),/remo) chk_and_upd, header, key, offsets[i]+1 endfor if (not keyword_set(no_comment)) then begin fxaddpar, header, 'COMMENT', ' ', before='TBCOL1' fxaddpar, header, 'COMMENT', ' *** Column offsets ***', before='TBCOL1' fxaddpar, header, 'COMMENT', ' ', before='TBCOL1' endif ; Write the TFORM keywords for i=1, n_elements(ttypes)-1 do begin key= 'TFORM'+strcompress(string(i),/remo) chk_and_upd, header, key, tforms[i] endfor if (not keyword_set(no_comment)) then begin fxaddpar, header, 'COMMENT', ' ', before='TFORM1' fxaddpar, header, 'COMMENT', ' *** Column formats ***', before='TFORM1' fxaddpar, header, 'COMMENT', ' ', before='TFORM1' endif ; Write the header. mwr_header, lun, header ; Now loop over the structure and write out the data. totalFormat = "("+totalFormat+")"; start = 0L last = 1023L while (start lt n_elements(input)) do begin if (last ge n_elements(input)) then begin last = n_elements(input) - 1 endif strings = string(input[start:last], format=totalFormat) if keyword_set(terminator) then begin strings = strings+terminator endif writeu, lun, strings start = last + 1 last = last + 1024 endwhile ; Check to see if any padding is required. nbytes = n_elements(input)*offset padding = 2880 - nbytes mod 2880 if padding ne 0 then begin pad = replicate(32b, padding) endif writeu, lun, pad return end ; Write a dummy primary header-data unit. pro mwr_dummy, lun fxaddpar, header, 'SIMPLE', 'T','Dummy Created by MWRFITS v'+mwr_version() fxaddpar, header, 'BITPIX', 8, 'Dummy primary header created by MWRFITS' fxaddpar, header, 'NAXIS', 0, 'No data is associated with this header' fxaddpar, header, 'EXTEND', 'T', 'Extensions may (will!) be present' mwr_header, lun, header end ; Check if this is a valid pointer array for variable length data. function mwr_validptr, vtypes, nfld, index, array type = -1 offset = 0L for i=0, n_elements(array)-1 do begin if ptr_valid(array[i]) then begin sz = size(*array[i]) if sz[0] gt 1 then begin print,'MWRFITS: Error: Multidimensional Pointer array' return, 0 endif if type eq -1 then begin type = sz[sz[0] + 1] endif else begin if sz[sz[0] + 1] ne type then begin print,'MWRFITS: Error: Inconsistent type in pointer array' return, 0 endif endelse xsz = sz[1] if sz[0] eq 0 then xsz = 1 offset = offset + xsz endif endfor if type eq -1 then begin ; If there is no data assume an I*2 type type = 2 endif if (type lt 1 or type gt 5) and (type lt 12 or type gt 15) then begin print,'MWRFITS: Error: Unsupported type for variable length array' endif types = 'BIJED IJKK' sizes = [1,2,4,4,8,0,0,0,0,0,0,2,4,8,8] if n_elements(vtypes) eq 0 then begin vtype = {status:0, data:array, $ type: strmid(types, type-1, 1), $ itype: type, ilen: sizes[type-1], $ offset:offset } vtypes = replicate(vtype, nfld) endif else begin ; This ensures compatible structures without ; having to used named structures. vtype = vtypes[0] vtype.status = 0 vtype.data = array vtype.type = strmid(types, type-1, 1) vtype.itype = type vtype.ilen = sizes[type-1] vtype.offset = offset vtypes[index] = vtype endelse vtypes[index].status = 1; return, 1 end ; Handle the header for a binary table. pro mwr_tablehdr, lun, input, header, vtypes, $ no_types=no_types, $ logical_cols = logical_cols, $ bit_cols = bit_cols, $ nbit_cols= nbit_cols, $ no_comment=no_comment, $ alias=alias, $ silent=silent if not keyword_set(no_types) then no_types = 0 nfld = n_tags(input[0]) if nfld le 0 then begin print, 'MWRFITS Error: Input contains no structure fields.' return endif tags = tag_names(input) ; Get the number of rows in the table. nrow = n_elements(input) dims = lonarr(nfld) tdims = strarr(nfld) types = strarr(nfld) pointers= lonarr(nfld) ; offsets = null... Don't want to define this ; in advance since reference to ulon64 won't word with IDL < 5.2 ; ; Get the type and length of each column. We do this ; by examining the contents of the first row of the structure. ; nbyte = 0 for i=0, nfld-1 do begin a = input[0].(i) sz = size(a) nelem = sz[sz[0]+2] type_ele = sz[sz[0]+1] if type_ele eq 7 then begin maxstr = max(strlen(input.(i)) > 1 ) endif dims[i] = nelem if (sz[0] lt 1) or (sz[0] eq 1 and type_ele ne 7) then begin tdims[i] = '' endif else begin tdims[i] = '(' if type_ele eq 7 then begin tdims[i] = tdims[i] + strcompress(string(maxstr), /remo) + ',' endif for j=1, sz[0] do begin tdims[i] = tdims[i] + strcompress(sz[j]) if j ne sz[0] then tdims[i] = tdims[i] + ',' endfor tdims[i] = tdims[i] + ')' endelse case type_ele of 1: begin types[i] = 'B' nbyte = nbyte + nelem end 2: begin types[i] = 'I' nbyte = nbyte + 2*nelem end 3: begin types[i] = 'J' nbyte = nbyte + 4*nelem end 4: begin types[i] = 'E' nbyte = nbyte + 4*nelem end 5: begin types[i] = 'D' nbyte = nbyte + 8*nelem end 6: begin types[i] = 'C' nbyte = nbyte + 8*nelem end 7: begin types[i] = 'A' nbyte = nbyte + maxstr*nelem dims[i] = maxstr*nelem end 9: begin types[i] = 'M' nbyte = nbyte + 16*nelem end 10: begin if not mwr_validptr(vtypes, nfld, i, input.(i)) then begin return endif types[i] = 'P'+vtypes[i].type nbyte = nbyte + 8 dims[i] = 1 test = mwr_unsigned_offset(vtypes[i].itype) if test gt 0 then begin if (n_elements(offsets) lt 1) then begin offsets = ulon64arr(nfld) endif offsets[i] = test endif end 12: begin types[i] = 'I' if (n_elements(offsets) lt 1) then begin offsets = ulon64arr(nfld) endif offsets[i] = mwr_unsigned_offset(12); nbyte = nbyte + 2*nelem end 13: begin types[i] = 'J' if (n_elements(offsets) lt 1) then begin offsets = ulon64arr(nfld) endif offsets[i] = mwr_unsigned_offset(13); nbyte = nbyte + 4*nelem end ; 8 byte integers became standard FITS in December 2005 14: begin types[i] = 'K' nbyte = nbyte + 8*nelem end 15: begin types[i] = 'K' nbyte = nbyte + 8*nelem if (n_elements(offsets) lt 1) then begin offsets = ulon64arr(nfld) endif offsets[i] = mwr_unsigned_offset(15) end 0: begin print,'MWRFITS Error: Undefined structure element??' return end 8: begin print, 'MWRFITS Error: Nested structures' return end else:begin print, 'MWRFITS Error: Cannot parse structure' return end endcase endfor ; Put in the required FITS keywords. chk_and_upd, header, 'XTENSION', 'BINTABLE', 'Binary table written by MWRFITS v'+mwr_version() chk_and_upd, header, 'BITPIX', 8, 'Required value' chk_and_upd, header, 'NAXIS', 2, 'Required value' chk_and_upd, header, 'NAXIS1', nbyte, 'Number of bytes per row' chk_and_upd, header, 'NAXIS2', n_elements(input), 'Number of rows' chk_and_upd, header, 'PCOUNT', 0, 'Normally 0 (no varying arrays)' chk_and_upd, header, 'GCOUNT', 1, 'Required value' chk_and_upd, header, 'TFIELDS', nfld, 'Number of columns in table' ; ; Handle the special cases. ; if keyword_set(logical_cols) then begin nl = n_elements(logical_cols) for i = 0, nl-1 do begin icol = logical_cols[i] if types[icol-1] ne 'A' and types[icol-1] ne 'B' then begin print,'WARNING: Invalid attempt to create Logical column:',icol goto, next_logical endif types[icol-1] = 'L' next_logical: endfor endif if keyword_set(bit_cols) then begin nb = n_elements(bit_cols) if nb ne n_elements(nbit_cols) then begin print,'WARNING: Bit_cols and Nbit_cols not same size' print,' No bit columns generated.' goto, after_bits endif for i = 0, nb-1 do begin nbyte = (nbit_cols[i]+7)/8 icol = bit_cols[i] if types[icol-1] ne 'B' or (dims[icol-1] ne nbyte) then begin print,'WARNING: Invalid attempt to create bit column:',icol goto, next_bit endif types[icol-1] = 'X' tdims[icol-1] = '' dims[icol-1] = nbit_cols[i] next_bit: endfor after_bits: endif ; Write scaling info as needed. if n_elements(offsets) gt 0 then begin w = where(offsets gt 0) for i=0, n_elements(w) - 1 do begin key = 'TSCAL'+strcompress(string(w[i])+1,/remo) chk_and_upd, header, key, 1 endfor for i=0, n_elements(w) - 1 do begin key = 'TZERO'+strcompress(string(w[i]+1),/remo) chk_and_upd, header, key, offsets[w[i]] endfor if not keyword_set(no_comment) then begin key = 'TSCAL'+strcompress(string(w[0])+1,/remo) fxaddpar, header, 'COMMENT', ' ', before=key fxaddpar, header, 'COMMENT', ' *** Unsigned integer column scalings ***', before=key fxaddpar, header, 'COMMENT', ' ', before=key endif endif ; Now add in the TFORM keywords for i=0, nfld-1 do begin if dims[i] eq 1 then begin form = types[i] endif else begin form=strcompress(string(dims[i]),/remove) + types[i] endelse tfld = 'TFORM'+strcompress(string(i+1),/remove) ; Check to see if there is an existing value for this keyword. ; If it has the proper value we will not modify it. ; This can matter if there is optional information coded ; beyond required TFORM information. oval = fxpar(header, tfld) oval = strcompress(string(oval),/remove_all) if (oval eq '0') or (strmid(oval, 0, strlen(form)) ne form) then begin chk_and_upd, header, tfld, form endif endfor if (not keyword_set(no_comment)) then begin fxaddpar, header, 'COMMENT', ' ', before='TFORM1' fxaddpar, header, 'COMMENT', ' *** Column formats ***', before='TFORM1' fxaddpar, header, 'COMMENT', ' ', before='TFORM1' endif ; Now write TDIM info as needed. for i=nfld-1, 0,-1 do begin if tdims[i] ne '' then begin fxaddpar, header, 'TDIM'+strcompress(string(i+1),/remo), tdims[i],after=tfld endif endfor w=where(tdims ne '') if w[0] ne -1 and not keyword_set(no_comment) then begin fxaddpar, header, 'COMMENT', ' ', after=tfld fxaddpar, header, 'COMMENT', ' *** Column dimensions (2 D or greater) ***', after=tfld fxaddpar, header, 'COMMENT', ' ', after=tfld endif for i=0, nfld-1 do begin if tdims[i] ne '' then begin chk_and_upd, header, 'TDIM'+strcompress(string(i+1),/remo), tdims[i] endif endfor if n_elements(vtypes) gt 0 then begin fxaddpar, header, 'THEAP', nbyte*n_elements(input), 'Offset of start of heap' offset = 0L for i=0,n_elements(vtypes)-1 do begin if vtypes[i].status then offset = offset + vtypes[i].offset*vtypes[i].ilen endfor fxaddpar, header, 'PCOUNT', offset, 'Size of heap' endif ; ; Last add in the TTYPE keywords if desired. ; if not no_types then begin for i=0, nfld - 1 do begin key = 'TTYPE'+strcompress(string(i+1),/remove) if not keyword_set(use_colnums) then begin value= mwr_checktype(tags[i],alias=alias)+' ' endif else begin value = 'C'+strmid(key,5,2) endelse chk_and_upd, header, key, value endfor if (not keyword_set(no_comment)) then begin fxaddpar, header, 'COMMENT', ' ', before='TTYPE1' fxaddpar, header, 'COMMENT', ' *** Column names *** ',before='TTYPE1' fxaddpar, header, 'COMMENT', ' ',before='TTYPE1' endif endif if (not keyword_set(no_comment)) then begin fxaddpar, header, 'COMMENT', ' ', after='TFIELDS' fxaddpar, header, 'COMMENT', ' *** End of mandatory fields ***', after='TFIELDS' fxaddpar, header, 'COMMENT', ' ', after='TFIELDS' endif ; Write to the output device. mwr_header, lun, header end ; Modify the structure to put the pointer column in. function mwr_retable, input, vtypes offset = 0L tags = tag_names(input); ;Create an output structure identical to the input structure but with pointers replaced ; by a 2 word lonarr to point to the heap area if vtypes[0].status then begin output = CREATE_STRUCT(tags[0],lonarr(2)) endif else begin output = CREATE_STRUCT(tags[0],input[0].(0)) endelse for i=1, n_elements(tags) -1 do begin if vtypes[i].status then begin output = CREATE_STRUCT(temporary(output), tags[i], lonarr(2)) endif else begin output = CREATE_STRUCT(temporary(output), tags[i], input[0].(i)) endelse endfor output = replicate(temporary(output), N_elements(input) ) struct_assign, input, output ;Available since V5.1 for i=0, n_elements(tags)-1 do begin if vtypes[i].status then begin for j=0, n_elements(input)-1 do begin ptr = input[j].(i) if ptr_valid(ptr) then begin sz = size(*ptr) if sz[0] eq 0 then xsz = 1 else xsz= sz[1] output[j].(i)[0] = xsz output[j].(i)[1] = offset offset = offset + vtypes[i].ilen*xsz endif endfor endif endfor return,output end ; Write the heap data. function mwr_writeheap, lun, vtypes offset = 0L for i=0, n_elements(vtypes)-1 do begin if vtypes[i].status then begin itype = vtypes[i].itype unsigned = mwr_unsigned_offset(itype) ptrs = vtypes[i].data for j=0,n_elements(ptrs)-1 do begin if ptr_valid(ptrs[j]) then begin if (unsigned gt 0) then begin *ptrs[j] = *ptrs[j] + unsigned endif writeu, lun, *ptrs[j] sz = size(*ptrs[j]) xsz = 1 > sz[1] offset = offset + xsz * vtypes[i].ilen endif endfor endif endfor return, offset end ; Write the brinary table. pro mwr_tabledat, lun, input, header, vtypes ; ; file -- unit to which data is to be written. ; Input -- IDL structure ; Header -- Filled header nfld = n_tags(input) ; Any special processing? typ = intarr(nfld) for i=0, nfld-1 do begin typ[i] = size(input.(i),/type) if (typ[i] eq 7) then begin dim = size(input.(i),/dimen) >1 siz = max(strlen(input.(i))>1) fmt = '(a' +strtrim(siz,2) +')' if total(dim) eq 1 then $ input.(i) = string(input.(i), format=fmt) else $ input.(i) = reform(string(input.(i), format=fmt), dim) endif unsigned = mwr_unsigned_offset(typ[i]) if (unsigned gt 0) then begin input.(i) = input.(i) + unsigned endif endfor if n_elements(vtypes) gt 0 then begin input = mwr_retable(input, vtypes) endif ; Write the data segment. ; writeu, lun, input nbyte = long(fxpar(header, 'NAXIS1')) nrow = n_elements(input) heap = 0 if n_elements(vtypes) gt 0 then begin heap = mwr_writeheap(lun, vtypes) endif siz = nbyte*nrow + heap padding = 2880 - (siz mod 2880) if padding eq 2880 then padding = 0 ; ; If necessary write the padding. ; if padding gt 0 then begin pad = bytarr(padding) ; Should be null-filled by default. writeu, lun, pad endif end ; Scale parameters for GROUPed data. pro mwr_pscale, grp, header, pscale=pscale, pzero=pzero ; This function assumes group is a 2-d array. if not keyword_set(pscale) and not keyword_set(pzero) then return if not keyword_set(pscale) then begin pscale = dblarr(sizg[1]) pscale[*] = 1. endif w = where(pzero eq 0.d0) if w[0] ne 0 then begin print, 'MWRFITS Warning: PSCALE value of 0 found, set to 1.' pscale[w] = 1.d0 endif if keyword_set(pscale) then begin for i=0L, sizg[1]-1 do begin key= 'PSCAL' + strcompress(string(i+1),/remo) chk_and_upd, header, key, pscale[i] endfor endif if not keyword_set(pzero) then begin pzero = dblarr(sizg[1]) pzero[*] = 0. endif else begin for i=0L, sizg[1]-1 do begin key= 'PZERO' + strcompress(string(i+1),/remo) chk_and_upd, header, key, pscale[i] endfor endelse for i=0L, sizg[1]-1 do begin grp[i,*] = grp[i,*]/pscale[i] - pzero[i] endfor end ; Find the appropriate scaling parameters. pro mwr_findscale, flag, array, nbits, scale, offset, error error = 0 if n_elements(flag) eq 2 then begin scale = double(flag[0]) offset = double(flag[1]) endif else if n_elements(flag) eq 1 and flag[0] ne 1 then begin minmum = min(array, max=maxmum) offset = 0.d0 scale = double(flag[0]) endif else if n_elements(flag) ne 1 then begin print, 'MWRFITS Error: Invalid scaling parameters.' error = 1 return endif else begin minmum = min(array, max=maxmum) scale = (maxmum-minmum)/(2.d0^nbits) amin = -(2.d0^(nbits-1)) if (amin gt -130) then amin = 0 ; looking for -128 offset = minmum - scale*amin endelse return end ; Scale and possibly convert array according to information ; in flags. pro mwr_scale, array, scale, offset, lscale=lscale, iscale=iscale, $ bscale=bscale, null=null ; First deallocate scale and offset if n_elements(scale) gt 0 then xx = temporary(scale) if n_elements(offset) gt 0 then xx = temporary(offset) if not keyword_set(lscale) and not keyword_set(iscale) and $ not keyword_set(bscale) then return siz = size(array) if keyword_set(lscale) then begin ; Doesn't make sense to scale data that can be stored exactly. if siz[siz[0]+1] lt 4 then return amin = -2.d0^31 amax = -(amin + 1) mwr_findscale, lscale, array, 32, scale, offset, error endif else if keyword_set(iscale) then begin if siz[siz[0]+1] lt 3 then return amin = -2.d0^15 amax = -(amin + 1) mwr_findscale, iscale, array, 16, scale, offset, error endif else begin if siz[siz[0]+1] lt 2 then return amin = 0 amax = 255 mwr_findscale, bscale, array, 8, scale, offset, error endelse ; Check that there was no error in mwr_findscale if error gt 0 then return if scale le 0.d0 then begin print, 'MWRFITS Error: BSCALE/TSCAL=0' return endif array = round((array-offset)/scale) w=where(array lt 0) w = where(array gt amax) if w[0] ne -1 then begin if keyword_set(null) then array[w] = null else array[w]=amax endif w = where(array lt amin) if w[0] ne -1 then begin if keyword_set(null) then array[w] = null else array[w] = amin endif if keyword_set(lscale) then array = long(array) $ else if keyword_set(iscale) then array = fix(array) $ else array = byte(array) end ; Write a header pro mwr_header, lun, header ; Fill strings to at least 80 characters and then truncate. space = string(replicate(32b, 80)) header = strmid(header+space, 0, 80) w = where(strmid(header,0,8) eq "END ") if w[0] eq -1 then begin header = [header, strmid("END"+space,0,80)] endif else begin if (n_elements(w) gt 1) then begin ; Get rid of extra end keywords; print,"MWRFITS Warning: multiple END keywords found." for irec=0L, n_elements(w)-2 do begin header[w[irec]] = strmid('COMMENT INVALID END REPLACED'+ $ space, 0, 80) endfor endif ; Truncate header array at END keyword. header = header[0:w[n_elements(w)-1]] endelse nrec = n_elements(header) if nrec mod 36 ne 0 then header = [header, replicate(space,36 - nrec mod 36)] writeu, lun, byte(header) end ; Move the group information within the data. pro mwr_groupinfix, data, group, hdr siz = size(data) sizg = size(group) ; Check if group info is same type as data if siz[siz[0]+1] ne sizg[3] then begin case siz[siz[0]+1] of 1: begin mwr_groupscale, 127.d0, group, hdr group = byte(group) end 2: begin mwr_groupscale, 32767.d0, group, hdr group = fix(group) end 3: begin mwr_groupscale, 2147483647.d0, group, hdr group = long(group) end 4: group = float(group) 5: group = double(group) else: begin print,'MWRFITS Internal error: Conversion of group data' return end endcase endif nrow = 1 for i=1, siz[0]-1 do begin nrow = nrow*siz[i] endfor data = reform(data, siz[siz[0]+2]) for i=0L, siz[siz[0]] - 1 do begin if i eq 0 then begin gdata = group[*,0] gdata = reform(gdata) tdata = [ gdata , data[0:nrow-1]] endif else begin start = nrow*i fin = start+nrow-1 gdata = group[*,i] tdata = [tdata, gdata ,data[start:fin]] endelse endfor data = temporary(tdata) end ; If an array is being scaled to integer type, then ; check to see if the group parameters will exceed the maximum ; values allowed. If so scale them and update the header. pro mwr_groupscale, maxval, group, hdr sz = size(group) for i=0L, sz[1]-1 do begin pmax = max(abs(group[i,*])) if (pmax gt maxval) then begin ratio = pmax/maxval psc = 'PSCAL'+strcompress(string(i+1),/remo) currat = fxpar(hdr, psc) if (currat ne 0) then begin fxaddpar, hdr, psc, currat*ratio, 'Scaling overriden by MWRFITS' endif else begin fxaddpar, hdr, psc, ratio, ' Scaling added by MWRFITS' endelse group[i,*] = group[i,*]/ratio endif endfor end ; Write out header and image for IMAGE extensions and primary arrays. pro mwr_image, input, siz, lun, bof, hdr, $ null=null, $ group=group, $ pscale=pscale, pzero=pzero, $ lscale=lscale, iscale=iscale, $ bscale=bscale, $ no_comment=no_comment, $ silent=silent type = siz[siz[0] + 1] bitpixes=[8,8,16,32,-32,-64,-32,0,0,-64,0,0,16,32,64,64] ; Convert complexes to two element real array. if type eq 6 or type eq 9 then begin if not keyword_set(silent) then begin print, "MWRFITS Note: Complex numbers treated as arrays" endif array_dimen=(2) if siz[0] gt 0 then array_dimen=[array_dimen, siz[1:siz[0]]] if siz[siz[0]+1] eq 6 then data = float(input,0,array_dimen) $ else data = double(input,0,array_dimen) ; Convert strings to bytes. endif else if type eq 7 then begin data = input len = max(strlen(input)) if len eq 0 then begin print, 'MWRFITS Error: strings all have zero length' return endif for i=0L, n_elements(input)-1 do begin t = len - strlen(input[i]) if t gt 0 then input[i] = input[i] + string(replicate(32B, len)) endfor ; Note that byte operation works on strings in a special way ; so we don't go through the subterfuge we tried above. data = byte(data) endif else if n_elements(input) gt 0 then data = input ; Convert scalar to 1-d array. if siz[0] eq 0 and siz[1] ne 0 then data=(data) ; Do any scaling of the data. mwr_scale, data, scalval, offsetval, lscale=lscale, $ iscale=iscale, bscale=bscale, null=null ; This may have changed the type. siz = size(data) type = siz[siz[0]+1] ; If grouped data scale the group parameters. if keyword_set(group) then mwr_pscale, group, hdr, pscale=pscale, pzero=pzero if bof then begin chk_and_upd, hdr, 'SIMPLE', 'T','Primary Header created by MWRFITS v'+mwr_version() chk_and_upd, hdr, 'BITPIX', bitpixes[type] chk_and_upd, hdr, 'NAXIS', siz[0] chk_and_upd, hdr, 'EXTEND', 'T', 'Extensions may be present' endif else begin chk_and_upd, hdr, 'XTENSION', 'IMAGE','Image Extension created by MWRFITS v'+mwr_version() chk_and_upd, hdr, 'BITPIX', bitpixes[type] chk_and_upd, hdr, 'NAXIS', siz[0] chk_and_upd, hdr, 'PCOUNT', 0 chk_and_upd, hdr, 'GCOUNT', 1 endelse if keyword_set(group) then begin group_offset = 1 endif else group_offset = 0 if keyword_set(group) then begin chk_and_upd, hdr, 'NAXIS1', 0 endif for i=1L, siz[0]-group_offset do begin chk_and_upd, hdr, 'NAXIS'+strcompress(string(i+group_offset),/remo), siz[i] endfor if keyword_set(group) then begin chk_and_upd, hdr, 'GROUPS', 'T' sizg = size(group) if sizg[0] ne 2 then begin print,'MWRFITS Error: Group data is not 2-d array' return endif if sizg[2] ne siz[siz[0]] then begin print,'MWRFITS Error: Group data has wrong number of rows' return endif chk_and_upd,hdr, 'PCOUNT', sizg[1] chk_and_upd, hdr, 'GCOUNT', siz[siz[0]] endif if n_elements(scalval) gt 0 then begin chk_and_upd, hdr, 'BSCALE', scalval chk_and_upd, hdr, 'BZERO', offsetval endif else begin ; Handle unsigned offsets bzero = mwr_unsigned_offset(type) if bzero gt 0 then begin chk_and_upd,hdr,'BSCALE', 1 chk_and_upd, hdr, 'BZERO', bzero data = data + bzero endif endelse if keyword_set(group) then begin if keyword_set(pscale) then begin if n_elements(pscale) ne sizg[1] then begin print, 'MWRFITS Warning: wrong number of PSCALE values' endif else begin for i=1L, sizg[1] do begin chk_and_upd, hdr, 'PSCALE'+strcompress(string(i),/remo) endfor endelse endif if keyword_set(pzero) then begin if n_elements(pscale) ne sizg[1] then begin print, 'MWRFITS Warning: Wrong number of PSCALE values' endif else begin for i=1L, sizg[1] do begin chk_and_upd, hdr, 'PZERO'+strcompress(string(i),/remo) endfor endelse endif endif bytpix=abs(bitpixes[siz[siz[0]+1]])/8 ; Number of bytes per pixel. npixel = n_elements(data) + n_elements(group) ; Number of pixels. if keyword_set(group) then mwr_groupinfix, data, group, hdr ; Write the FITS header mwr_header, lun, hdr ; This is all we need to do if input is undefined. if (n_elements(input) eq 0) or (siz[0] eq 0) then return ; Write the data. writeu, lun, data nbytes = bytpix*npixel filler = 2880 - nbytes mod 2880 if filler eq 2880 then filler = 0 ; Write any needed filler. if filler gt 0 then writeu, lun, replicate(0B,filler) end ; Main routine -- see documentation at start pro mwrfits, xinput, file, header, $ ascii=ascii, $ separator=separator, $ terminator=terminator, $ create=create, $ null=null, $ group=group, $ pscale=pscale, pzero=pzero, $ alias=alias, $ use_colnum = use_colnum, $ lscale=lscale, iscale=iscale, $ bscale=bscale, $ no_types=no_types, $ silent=silent, $ no_comment=no_comment, $ logical_cols=logical_cols, $ bit_cols=bit_cols, $ nbit_cols=nbit_cols, $ version=version ; Check required keywords. if (keyword_set(Version)) then begin print, "MWRFITS V"+mwr_version()+": January 3, 2005" endif if n_elements(file) eq 0 then begin if (not keyword_set(Version)) then begin print, 'MWRFITS: Usage:' print, ' MWRFITS, struct_name, file, [header,] ' print, ' /CREATE, /SILENT, /NO_TYPES, /NO_COMMENT, ' print, ' GROUP=, PSCALE=, PZERO=,' print, ' LSCALE=, ISCALE=, BSCALE=,' print, ' LOGICAL_COLS=, BIT_COLS=, NBIT_COLS=,' print, ' ASCII=, SEPARATOR=, TERMINATOR=, NULL=' print, ' /USE_COLNUM, ALIAS=' endif return endif ; Save the data into an array/structure that we can modify. if n_elements(xinput) gt 0 then input = xinput on_ioerror, open_error ; Open the input file. ; If the create keyword is not specified we ; try to open the file readonly to see if it ; already exists and if so we append to it. ; An error implies the file does not exist. ; ; We use this rather circuitous route to handle ; the unix ~ construction consistently -- findfile ; doesn't reliably understand that. ; if not keyword_set(create) then begin on_ioerror, not_found openr, lun, file, /get_lun,/swap_if_little free_lun, lun on_ioerror, open_error openu, lun, file, /get_lun, /append,/swap_if_little bof = 0 goto, finished_open endif not_found: on_ioerror, open_error openw, lun, file, /get_lun, /swap_if_little bof = 1 on_ioerror, null finished_open: siz = size(input) if siz[siz[0]+1] ne 8 then begin ; If input is not a structure then call image writing utilities. mwr_image, input, siz, lun, bof, header, $ null=null, $ group=group, $ pscale=pscale, pzero=pzero, $ lscale=lscale, iscale=iscale, $ bscale=bscale, $ no_comment=no_comment, $ silent=silent endif else if keyword_set(ascii) then begin if bof then mwr_dummy, lun ; Create an ASCII table. mwr_ascii, input, siz, lun, bof, header, $ ascii=ascii, $ null=null, $ use_colnum = use_colnum, $ lscale=lscale, iscale=iscale, $ bscale=bscale, $ no_types=no_types, $ separator=separator, $ terminator=terminator, $ no_comment=no_comment, $ alias=alias, $ silent=silent endif else begin if bof then mwr_dummy, lun ; Create a binary table. mwr_tablehdr, lun, input, header, vtypes, $ no_types=no_types, $ logical_cols = logical_cols, $ bit_cols = bit_cols, $ nbit_cols= nbit_cols, $ alias=alias, $ no_comment=no_comment, $ silent=silent mwr_tabledat, lun, input, header, vtypes endelse free_lun, lun return ; Handle error in opening file. open_error: on_ioerror, null print, 'MWRFITS Error: Cannot open output: ', file print,!ERROR_STATE.SYS_MSG if n_elements(lun) gt 0 then free_lun, lun return end ;+ ; 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 ; 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 ;- 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.11' 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,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 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 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 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, /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 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, /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() 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()+' May 01, 2006' 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 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 ; FUNCTION FXPAR, HDR, NAME, ABORT, COUNT=MATCHES, COMMENT=COMMENTS, $ START=START, PRECHECK=PRECHECK, POSTCHECK=POSTCHECK, $ NOCONTINUE = NOCONTINUE, $ DATATYPE=DATATYPE ;+ ; NAME: ; FXPAR() ; PURPOSE: ; Obtain the value of a parameter in a FITS header. ; EXPLANATION: ; The first 8 chacters of each element of HDR are searched for a match to ; NAME. If the keyword is one of those allowed to take multiple values ; ("HISTORY", "COMMENT", or " " (blank)), then the value is taken ; as the next 72 characters. Otherwise, it is assumed that the next ; character is "=", and the value (and optional comment) is then parsed ; from the last 71 characters. An error occurs if there is no parameter ; with the given name. ; ; If the value is too long for one line, it may be continued on to the ; the next input card, using the OGIP CONTINUE convention. For more info, ; http://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/ofwg_recomm/r13.html ; ; Complex numbers are recognized as two numbers separated by one or more ; space characters. ; ; If a numeric value has no decimal point (or E or D) it is returned as ; type LONG. If it contains more than 8 numerals, or contains the ; character 'D', then it is returned as type DOUBLE. Otherwise it is ; returned as type FLOAT. If an integer is too large to be stored as ; type LONG, then it is returned as DOUBLE. ; ; CALLING SEQUENCE: ; Result = FXPAR( HDR, NAME [, ABORT, COUNT=, COMMENT=, /NOCONTINUE ] ) ; ; Result = FXPAR(HEADER,'DATE') ;Finds the value of DATE ; Result = FXPAR(HEADER,'NAXIS*') ;Returns array dimensions as ; ;vector ; REQUIRED INPUTS: ; HDR = FITS header string array (e.g. as returned by FXREAD). Each ; element should have a length of 80 characters ; NAME = String name of the parameter to return. If NAME is of the ; form 'keyword*' then an array is returned containing values ; of keywordN where N is an integer. The value of keywordN ; will be placed in RESULT(N-1). The data type of RESULT will ; be the type of the first valid match of keywordN ; found, unless DATATYPE is given. ; OPTIONAL INPUT: ; ABORT = String specifying that FXPAR should do a RETALL if a ; parameter is not found. ABORT should contain a string to be ; printed if the keyword parameter is not found. If not ; supplied, FXPAR will return with a negative !err if a keyword ; is not found. ; DATATYPE = A scalar value, indicating the type of vector ; data. All keywords will be cast to this type. ; Default: based on first keyword. ; Example: DATATYPE=0.0D (cast data to double precision) ; START = A best-guess starting position of the sought-after ; keyword in the header. If specified, then FXPAR ; first searches for scalar keywords in the header in ; the index range bounded by START-PRECHECK and ; START+POSTCHECK. This can speed up keyword searches ; in large headers. If the keyword is not found, then ; FXPAR searches the entire header. ; ; If not specified then the entire header is searched. ; Searches of the form 'keyword*' also search the ; entire header and ignore START. ; ; Upon return START is changed to be the position of ; the newly found keyword. Thus the best way to ; search for a series of keywords is to search for ; them in the order they appear in the header like ; this: ; ; START = 0L ; P1 = FXPAR('P1', START=START) ; P2 = FXPAR('P2', START=START) ; PRECHECK = If START is specified, then PRECHECK is the number ; of keywords preceding START to be searched. ; Default: 5 ; POSTCHECK = If START is specified, then POSTCHECK is the number ; of keywords after START to be searched. ; Default: 20 ; OUTPUT: ; The returned value of the function is the value(s) associated with the ; requested keyword in the header array. ; ; If the parameter is complex, double precision, floating point, long or ; string, then the result is of that type. Apostrophes are stripped from ; strings. If the parameter is logical, 1 is returned for T, and 0 is ; returned for F. ; ; If NAME was of form 'keyword*' then a vector of values are returned. ; ; OPTIONAL INPUT KEYWORDS: ; /NOCONTINUE = If set, then continuation lines will not be read, even ; if present in the header ; OPTIONAL OUTPUT KEYWORD: ; COUNT = Optional keyword to return a value equal to the number of ; parameters found by FXPAR. ; COMMENTS= Array of comments associated with the returned values. ; ; PROCEDURE CALLS: ; GETTOK(), VALID_NUM ; SIDE EFFECTS: ; ; The system variable !err is set to -1 if parameter not found, 0 for a ; scalar value returned. If a vector is returned it is set to the number ; of keyword matches found. ; ; If a keyword occurs more than once in a header, a warning is given, ; and the first occurence is used. However, if the keyword is "HISTORY", ; "COMMENT", or " " (blank), then multiple values are returned. ; ; NOTES: ; The functions SXPAR() and FXPAR() are nearly identical, although ; FXPAR() has slightly more sophisticated parsing. There is no ; particular reason for having two nearly identical procedures, but ; both are too widely used to drop either one. ; ; REVISION HISTORY: ; Version 1, William Thompson, GSFC, 12 April 1993. ; Adapted from SXPAR ; Version 2, William Thompson, GSFC, 14 October 1994 ; Modified to use VALID_NUM instead of STRNUMBER. Inserted ; additional call to VALID_NUM to trap cases where character ; strings did not contain quotation marks. ; Version 3, William Thompson, GSFC, 22 December 1994 ; Fixed bug with blank keywords, following suggestion by Wayne ; Landsman. ; Version 4, Mons Morrison, LMSAL, 9-Jan-98 ; Made non-trailing ' for string tag just be a warning (not ; a fatal error). It was needed because "sxaddpar" had an ; error which did not write tags properly for long strings ; (over 68 characters) ; Version 5, Wayne Landsman GSFC, 29 May 1998 ; Fixed potential problem with overflow of LONG values ; Version 6, Craig Markwardt, GSFC, 28 Jan 1998, ; Added CONTINUE parsing ; Version 7, Craig Markwardt, GSFC, 18 Nov 1999, ; Added START, PRE/POSTCHECK keywords for better ; performance ; Version 8, Craig Markwardt, GSFC, 08 Oct 2003, ; Added DATATYPE keyword to cast vector keywords type ; Version 9, Paul Hick, 22 Oct 2003, Corrected bug (NHEADER-1) ;- ;------------------------------------------------------------------------------ ; ; Check the number of parameters. ; IF N_PARAMS() LT 2 THEN BEGIN PRINT,'Syntax: result = FXPAR( HDR, NAME [, ABORT ])' RETURN, -1 ENDIF ; ; Determine the abort condition. ; VALUE = 0 IF N_PARAMS() LE 2 THEN BEGIN ABORT_RETURN = 0 ABORT = 'FITS Header' END ELSE ABORT_RETURN = 1 IF ABORT_RETURN THEN ON_ERROR,1 ELSE ON_ERROR,2 ; ; 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' ; ; Convert the selected keyword NAME to uppercase. ; NAM = STRTRIM( STRUPCASE(NAME) ) ; ; Determine if NAME is of form 'keyword*'. If so, then strip off the '*', and ; set the VECTOR flag. One must consider the possibility that NAM is an empty ; string. ; NAMELENGTH1 = (STRLEN(NAM) - 1) > 1 IF STRPOS( NAM, '*' ) EQ NAMELENGTH1 THEN BEGIN NAM = STRMID( NAM, 0, NAMELENGTH1) VECTOR = 1 ;Flag for vector output NAME_LENGTH = STRLEN(NAM) ;Length of name NUM_LENGTH = 8 - NAME_LENGTH ;Max length of number portion IF NUM_LENGTH LE 0 THEN MESSAGE, $ 'Keyword length must be 8 characters or less' ; ; Otherwise, extend NAME with blanks to eight characters. ; ENDIF ELSE BEGIN WHILE STRLEN(NAM) LT 8 DO NAM = NAM + ' ' VECTOR = 0 ENDELSE ; ; If of the form 'keyword*', then find all instances of 'keyword' followed by ; a number. Store the positions of the located keywords in NFOUND, and the ; value of the number field in NUMBER. ; IF N_ELEMENTS(START) EQ 0 THEN START = -1L START = LONG(START[0]) IF NOT VECTOR AND START GE 0 THEN BEGIN IF N_ELEMENTS(PRECHECK) EQ 0 THEN PRECHECK = 5 IF N_ELEMENTS(POSTCHECK) EQ 0 THEN POSTCHECK = 20 NHEADER = N_ELEMENTS(HDR) MN = (START - PRECHECK) > 0 MX = (START + POSTCHECK) < (NHEADER-1) ;Corrected bug KEYWORD = STRMID(HDR[MN:MX], 0, 8) ENDIF ELSE BEGIN RESTART: START = -1L KEYWORD = STRMID( HDR, 0, 8) ENDELSE IF VECTOR THEN BEGIN NFOUND = WHERE(STRPOS(KEYWORD,NAM) GE 0, MATCHES) IF ( MATCHES GT 0 ) THEN BEGIN NUMST= STRMID(HDR[NFOUND], NAME_LENGTH, NUM_LENGTH) NUMBER = INTARR(MATCHES)-1 FOR I = 0, MATCHES-1 DO $ IF VALID_NUM( NUMST[I], NUM) THEN NUMBER[I] = NUM IGOOD = WHERE(NUMBER GE 0, MATCHES) IF MATCHES GT 0 THEN BEGIN NFOUND = NFOUND[IGOOD] NUMBER = NUMBER[IGOOD] ENDIF ENDIF ; ; Otherwise, find all the instances of the requested keyword. If more than ; one is found, and NAME is not one of the special cases, then print an error ; message. ; ENDIF ELSE BEGIN NFOUND = WHERE(KEYWORD EQ NAM, MATCHES) IF MATCHES EQ 0 AND START GE 0 THEN GOTO, RESTART IF START GE 0 THEN NFOUND = NFOUND + MN IF (MATCHES GT 1) AND (NAM NE 'HISTORY ') AND $ (NAM NE 'COMMENT ') AND (NAM NE '') THEN $ MESSAGE,/INFORMATIONAL, 'WARNING- Keyword ' + $ NAM + 'located more than once in ' + ABORT IF (MATCHES GT 0) THEN START = NFOUND[MATCHES-1] ENDELSE ; ; Extract the parameter field from the specified header lines. If one of the ; special cases, then done. ; IF MATCHES GT 0 THEN BEGIN LINE = HDR[NFOUND] SVALUE = STRTRIM( STRMID(LINE,9,71),2) IF (NAM EQ 'HISTORY ') OR (NAM EQ 'COMMENT ') OR $ (NAM EQ ' ') THEN BEGIN VALUE = STRTRIM( STRMID(LINE,8,72),2) COMMENTS = STRARR(N_ELEMENTS(VALUE)) ; ; Otherwise, test to see if the parameter contains a string, signalled by ; beginning with a single quote character (') (apostrophe). ; END ELSE FOR I = 0,MATCHES-1 DO BEGIN IF ( STRMID(SVALUE[I],0,1) EQ "'" ) THEN BEGIN 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 '+NAME+' invalid in '+ABORT+ " (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 ; ; Extract the comment, if any. ; SLASH = STRPOS(TEST, "/", ENDAP) IF SLASH LT 0 THEN COMMENT = '' ELSE $ COMMENT = STRMID(TEST, SLASH+1, STRLEN(TEST)-SLASH-1) ; ; 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 KEYWORD_SET(NOCONTINUE) THEN BEGIN OFF = OFF + 1 VAL = STRTRIM(VALUE,2) IF (STRLEN(VAL) GT 0) AND $ (STRMID(VAL, STRLEN(VAL)-1, 1) EQ '&') AND $ (STRMID(HDR[NFOUND[I]+OFF],0,8) EQ 'CONTINUE') THEN BEGIN IF (SIZE(FXPAR(HDR, 'LONGSTRN',/NOCONTINUE)))[1] EQ 7 THEN BEGIN VALUE = STRMID(VAL, 0, STRLEN(VAL)-1) TEST = HDR[NFOUND[I]+OFF] TEST = STRMID(TEST, 8, STRLEN(TEST)-8) TEST = STRTRIM(TEST, 2) IF STRMID(TEST, 0, 1) NE "'" THEN MESSAGE, $ 'ERROR: Invalidly CONTINUEd string in '+ABORT NEXT_CHAR = 1 GOTO, NEXT_APOST ENDIF ENDIF ENDIF ; ; If not a string, then separate the parameter field from the comment field. ; ENDIF ELSE BEGIN TEST = SVALUE[I] SLASH = STRPOS(TEST, "/") IF SLASH GT 0 THEN BEGIN COMMENT = STRMID(TEST, SLASH+1, STRLEN(TEST)-SLASH-1) TEST = STRMID(TEST, 0, SLASH) END ELSE COMMENT = '' ; ; 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 VALUE2 = 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. ; IF VECTOR THEN BEGIN MAXNUM = MAX(NUMBER) IF ( I EQ 0 ) THEN BEGIN IF N_ELEMENTS(DATATYPE) EQ 0 THEN BEGIN ;; Data type determined from keyword SZ_VALUE = SIZE(VALUE) ENDIF ELSE BEGIN ;; Data type requested by user SZ_VALUE = SIZE(DATATYPE[0]) ENDELSE RESULT = MAKE_ARRAY( MAXNUM, TYPE=SZ_VALUE[1]) COMMENTS = STRARR(MAXNUM) ENDIF RESULT[ NUMBER[I]-1 ] = VALUE COMMENTS[ NUMBER[I]-1 ] = COMMENT ENDIF ELSE BEGIN COMMENTS = COMMENT ENDELSE ENDFOR ; ; Set the value of !ERR for the number of matches for vectors, or simply 0 ; otherwise. ; IF VECTOR THEN BEGIN !ERR = MATCHES RETURN, RESULT ENDIF ELSE !ERR = 0 ; ; Error point for keyword not found. ; ENDIF ELSE BEGIN IF ABORT_RETURN THEN MESSAGE,'Keyword '+NAM+' not found in '+ABORT !ERR = -1 ENDELSE ; RETURN, VALUE END FUNCTION FXMOVE, UNIT, N_EXT, SILENT = Silent ;+ ; NAME: ; FXMOVE ; PURPOSE: ; Skip a specified number of extensions in a FITS file ; ; CALLING SEQUENCE: ; STATUS=FXMOVE(UNIT, N_EXT, /Silent) ; ; INPUT PARAMETERS: ; UNIT = An open unit descriptor for a FITS data stream. ; N_EXT = Number of extensions to skip. ; ; OPTIONAL INPUT PARAMETER: ; /SILENT - If set, then any messages about invalid characters in the ; FITS file are suppressed. ; RETURNS: ; 0 if successful. ; -1 if an error is encountered. ; ; COMMON BLOCKS: ; None. ; SIDE EFFECTS: ; Repositions the file pointer. ; PROCEDURE: ; Each FITS header is read in and parsed, and the file pointer is moved ; to where the next FITS extension header until the desired ; extension is reached. ; PROCEDURE CALLS: ; FXPAR(), MRD_HREAD, MRD_SKIP ; MODIFICATION HISTORY: ; Extracted from FXPOSIT 8-March-2000 by T. McGlynn ; Added /SILENT keyword 14-Dec-2000 by W. Landsman ; Save time by not reading the full header W. Landsman Feb. 2003 ;- FOR EXT = 1, N_EXT DO BEGIN ; ; Read the next header, and get the number of bytes taken up by the data. ; IF EOF(UNIT) THEN RETURN, -1 ; Can't use FXHREAD to read from pipe, since it uses ; POINT_LUN. So we read this in ourselves using mrd_hread MRD_HREAD, UNIT, HEADER, STATUS, SILENT = Silent, /FIRSTBLOCK IF STATUS LT 0 THEN RETURN, -1 ; Get parameters that determine size of data ; region. BITPIX = FXPAR(HEADER,'BITPIX') NAXIS = FXPAR(HEADER,'NAXIS') GCOUNT = FXPAR(HEADER,'GCOUNT') IF GCOUNT EQ 0 THEN GCOUNT = 1 PCOUNT = FXPAR(HEADER,'PCOUNT') IF NAXIS GT 0 THEN BEGIN DIMS = FXPAR(HEADER,'NAXIS*') ;Read dimensions NDATA = DIMS[0] IF NAXIS GT 1 THEN FOR I=2,NAXIS DO NDATA = NDATA*DIMS[I-1] ENDIF ELSE NDATA = 0 NBYTES = (ABS(BITPIX) / 8) * GCOUNT * (PCOUNT + NDATA) ; ; Move to the next extension header in the file. ; NREC = LONG((NBYTES + 2879) / 2880) MRD_SKIP, UNIT, NREC*2880L ENDFOR RETURN, 0 END FUNCTION FXPOSIT, XFILE, EXT_NO, readonly=readonly, COMPRESS=COMPRESS, $ SILENT = Silent ;+ ; NAME: ; FXPOSIT ; PURPOSE: ; Return the unit number of a FITS file positioned at specified extension ; EXPLANATION: ; The FITS file will be ready to be read at the beginning of the ; specified extension. Called by headfits.pro, mrdfits.pro, readfits.pro ; ; CALLING SEQUENCE: ; unit=FXPOSIT(FILE, EXT_NO, /READONLY, COMPRESS=program, /SILENT) ; ; INPUT PARAMETERS: ; FILE = FITS file name, scalar string ; EXT_NO = Extension to be moved to, scalar nonnegative integer ; ; RETURNS: ; Unit number of file or -1 if an error is detected. ; ; OPTIONAL KEYWORD PARAMETER: ; /READONLY - If this keyword is set and non-zero, then OPENR rather ; than OPENU will be used to open the FITS file. ; COMPRESS - If this keyword is set and non-zero, then then treat ; the file as compressed. If 1 assume a gzipped file. ; and use IDLs internal decompression facility. For Unix ; compressed or bzip2 compressed files spawn off a process to ; decompress and use its output as the FITS stream. If the ; keyword is not 1, then use its value as a string giving the ; command needed for decompression. ; /SILENT If set, then suppress any messages about invalid characters ; in the FITS file. ; ; COMMON BLOCKS: ; None. ; SIDE EFFECTS: ; Opens and returns the descriptor of a file. ; PROCEDURE: ; Open the appropriate file, or spawn a command and intercept ; the output. ; Call FXMOVE to get to the appropriate extension. ; PROCEDURE CALLS: ; FXMOVE() ; MODIFICATION HISTORY: ; Derived from William Thompson's FXFINDEND routine. ; Modified by T.McGlynn, 5-October-1994. ; Modified by T.McGlynn, 25-Feb-1995 to handle compressed ; files. Pipes cannot be accessed using FXHREAD so ; MRD_HREAD was written. ; W. Landsman 23-Apr-1997 Force the /bin/sh shell when uncompressing ; W. Landsman 26-May-1997 Non-unix is not just VMS ; T. McGlynn 22-Apr-1999 Add /binary modifier needed for Windows ; T. McGlynn 03-June-1999 Use /noshell option to get rid of processes left by spawn. ; Use findfile to retain ability to use wildcards ; W. Landsman 03-Aug-1999 Use EXPAND_TILDE under Unix to find file ; T. McGlynn 04-Apr-2000 Put reading code into FXMOVE, ; additional support for compression from D.Palmer. ; W. Landsman/D.Zarro 04-Jul-2000 Added test for !VERSION.OS EQ 'Win32' (WinNT) ; W. Landsman 12-Dec-2000 Added /SILENT keyword ; W. Landsman April 2002 Use FILE_SEARCH for V5.5 or later ; W. Landsman Feb 2004 Assume since V5.3 (OPENR,/COMPRESS available) ; W. Landsman,W. Thompson, 2-Mar-2004, Add support for BZIP2 ; W. Landsman Don't leave open file if an error occurs ; W. Landsman Sep 2004 Treat FTZ extension as gzip compressed ; W. Landsman Feb 2006 Removed leading spaces (prior to V5.5) ;- ; ON_ERROR,2 compile_opt idl2 ; ; Check the number of parameters. ; IF N_PARAMS() LT 2 THEN BEGIN PRINT,'SYNTAX: UNIT = FXPOSIT(FILE, EXT_NO, /Readonly, compress=prog)' RETURN,-1 ENDIF FILE = FILE_SEARCH(XFILE, COUNT=COUNT) IF COUNT EQ 0 THEN $ RETURN, -1 ; Don't print anything out, just report an error FILE = FILE[0] ; ; Check if this is a compressed file. ; UNIT = -1 UCMPRS = ' ' IF KEYWORD_SET(compress) THEN BEGIN IF strcompress(string(compress),/remo) eq '1' THEN BEGIN compress = 'gunzip' ENDIF UCMPRS = compress; ENDIF ELSE BEGIN LEN = STRLEN(FILE) IF LEN GT 3 THEN $ TAIL = STRLOWCASE(STRMID(FILE, LEN-3, 3)) $ ELSE TAIL = ' ' IF STRMID(TAIL,1,2) EQ '.z' THEN $ UCMPRS = 'uncompress' $ ELSE IF TAIL EQ '.gz' or tail EQ 'ftz' THEN $ UCMPRS = 'gunzip' $ ELSE IF TAIL EQ 'bz2' THEN $ UCMPRS = 'bunzip2' ENDELSE ; Handle compressed files. IF UCMPRS EQ 'gunzip' THEN BEGIN IF KEYWORD_SET(READONLY) THEN BEGIN OPENR, UNIT, FILE, /COMPRESS, /GET_LUN, ERROR = ERROR ENDIF ELSE BEGIN OPENU, UNIT, FILE, /COMPRESS, /GET_LUN, ERROR = ERROR ENDELSE ENDIF ELSE IF UCMPRS NE ' ' THEN BEGIN IF (!VERSION.OS_FAMILY EQ 'unix') THEN BEGIN SPAWN, [UCMPRS,'-c',FILE], UNIT=UNIT, /NOSHELL ENDIF ELSE BEGIN PRINT, 'MRDFITS: Only Unix IDL supports piped spawns' PRINT, ' File must be uncompressed manually' RETURN, -1 ENDELSE ENDIF ELSE BEGIN ; ; Go to the start of the file. ; IF KEYWORD_SET(READONLY) THEN $ OPENR, UNIT, FILE, /GET_LUN, ERROR = ERROR $ ELSE OPENU, UNIT, FILE, /GET_LUN, ERROR = ERROR IF ERROR NE 0 THEN BEGIN PRINT,!ERROR_STATE.MSG RETURN,-1 ENDIF ENDELSE IF EXT_NO LE 0 THEN RETURN, UNIT STAT = FXMOVE(UNIT, EXT_NO, SILENT = Silent) IF STAT LT 0 THEN BEGIN FREE_LUN, UNIT RETURN, STAT ENDIF ELSE RETURN, UNIT END pro fits_open,filename,fcb,write=write,append=append,update=update, $ no_abort=no_abort,message=message,hprint=hprint ;+ ; NAME: ; FITS_OPEN ; ; PURPOSE: ; Opens a FITS (Flexible Image Transport System) data file. ; ; EXPLANATION: ; FITS_OPEN was modified in Sep 2006 to open with /SWAP_IF_LITTLE_ENDIAN ; It must be used with post-Sep 2006 versions of FITS_READ, FITS_WRITE ; and MODFITS. ;*CATEGORY: ; INPUT/OUTPUT ; ;*CALLING SEQUENCE: ; FITS_OPEN, filename, fcb ; ;*INPUTS: ; filename : name of the FITS file to open, scalar string ; FITS_OPEN can also open gzip compressed (.gz) file *for ; reading only*, although there is a performance penalty ;*OUTPUTS: ; fcb : (FITS Control Block) a IDL structure containing information ; concerning the file. It is an input to FITS_READ, FITS_WRITE ; and FITS_CLOSE ; ; INPUT KEYWORD PARAMETERS: ; /APPEND: Set to append to an existing file. ; /HPRINT - print headers with routine HPRINT as they are read. ; (useful for debugging a strange file) ; /NO_ABORT: Set to return to calling program instead of a RETALL ; when an I/O error is encountered. If set, the routine will ; return a non-null string (containing the error message) in the ; keyword MESSAGE. (For backward compatibility, the obsolete ; system variable !ERR is also set to -1 in case of an error.) ; If /NO_ABORT not set, then FITS_OPEN will print the message and ; issue a RETALL ; /UPDATE Set this keyword to open an existing file for update ; /WRITE: Set this keyword to open a new file for writing. ; ; OUTPUT KEYWORD PARAMETERS: ; MESSAGE = value: Output error message. If the FITS file was opened ; successfully, then message = ''. ; ;*NOTES: ; The output FCB should be passed to the other FITS routines (FITS_OPEN, ; FITS_READ, FITS_HELP, and FITS_WRITE). It has the following structure ; when FITS_OPEN is called without /WRITE or /APPEND keywords set. ; ; FCB.FILENAME - name of the input file ; .UNIT - unit number the file is opened to ; .NEXTEND - number of extensions in the file. ; .XTENSION - string array giving the extension type for each ; extension. ; .EXTNAME - string array giving the extension name for each ; extension. (null string if not defined the extension) ; .EXTVER - vector of extension version numbers (0 if not ; defined) ; .EXTLEVEL - vector of extension levels (0 if not defined) ; .GCOUNT - vector with the number of groups in each extension. ; .PCOUNT - vector with parameter count for each group ; .BITPIX - BITPIX for each extension with values ; 8 byte data ; 16 short word integers ; 32 long word integers ; -32 IEEE floating point ; -64 IEEE double precision floating point ; .NAXIS - number of axes for each extension. (0 for null data ; units) ; .AXIS - 2-D array where axis(*,N) gives the size of each axes ; for extension N ; .START_HEADER - vector giving the starting byte in the file ; where each extension header begins ; .START_DATA - vector giving the starting byte in the file ; where the data for each extension begins ; ; .HMAIN - keyword parameters (less standard required FITS ; keywords) for the primary data unit. ; .OPEN_FOR_WRITE - flag (0= open for read, 1=open for write, ; 2=open for update) ; .LAST_EXTENSION - last extension number read. ; .RANDOM_GROUPS - 1 if the PDU is random groups format, ; 0 otherwise ; .NBYTES - total number of (uncompressed) bytes in the FITS file ; ; When FITS open is called with the /WRITE or /APPEND option, FCB ; contains: ; ; FCB.FILENAME - name of the input file ; .UNIT - unit number the file is opened to ; .NEXTEND - number of extensions in the file. ; .OPEN_FOR_WRITE - flag (1=open for write, 2=open for append ; 3=open for update) ; ; ;*EXAMPLES: ; Open a FITS file for reading: ; FITS_OPEN,'myfile.fits',fcb ; ; Open a new FITS file for output: ; FITS_OPEN,'newfile.fits',fcb,/write ; PROCEDURES USED: ; HPRINT, SXDELPAR, SXPAR() ;*HISTORY: ; Written by: D. Lindler August, 1995 ; July, 1996 NICMOS Modified to allow open for overwrite ; to allow primary header to be modified ; DJL Oct. 15, 1996 corrected to properly extend AXIS when more ; than 100 extensions present ; Converted to IDL V5.0 W. Landsman September 1997 ; Use Message = '' rather than !ERR =1 as preferred signal of normal ; operation W. Landsman November 2000 ; Lindler, Dec, 2001, Modified to use 64 bit words for storing byte ; positions within the file to allow support for very large ; files ; Work with gzip compressed files W. Landsman January 2003 ; Fix gzip compress for V5.4 and earlier W.Landsman/M.Fitzgerald Dec 2003 ; Assume since V5.3 (STRSPLIT, OPENR,/COMPRESS) W. Landsman Feb 2004 ; Treat FTZ extension as gzip compressed W. Landsman Sep 2004 ; Assume since V5.4 fstat.compress available W. Landsman Apr 2006 ; FCB.Filename now expands any wildcards W. Landsman July 2006 ; Make ndata 64bit for very large files B. Garwood/W. Landsman Sep 2006 ; Open with /SWAP_IF_LITTLE_ENDIAN, remove obsolete keywords to OPEN ; W. Landsman Sep 2006 ;- ;-------------------------------------------------------------------- compile_opt idl2 ; if no parameters supplied, print calling sequence ; if N_params() LT 1 then begin print,'Syntax - FITS_OPEN, filename, fcb' print,' Input Keywords: /Append, /Hprint, /No_abort, /Update, /Write' print,' Output Keyword: Message= ' return endif ; ; set default keyword parameters ; message = '' open_for_read = 1 open_for_update = 0 open_for_write = 0 open_for_overwrite = 0 if keyword_set(write) then begin open_for_read = 0 open_for_update = 0 open_for_write = 1 open_for_overwrite = 0 end if keyword_set(append) then begin open_for_read = 0 open_for_write = 0 open_for_update = 1 open_for_overwrite = 0 end if keyword_set(update) then begin open_for_read = 1 open_for_write = 0 open_for_update = 0 open_for_overwrite = 1 end ; ; on I/O errors goto statement ioerror: ; on_ioerror,ioerror ; ; open file ; docompress = 0 if open_for_read and (not open_for_overwrite) then begin len = strlen(filename) ext = strlowcase(strmid(filename, len-3, 3)) docompress = (ext EQ '.gz') or (ext EQ 'ftz') endif ; ; open file ; get_lun,unit if docompress then $ openr,unit,filename, /compress,/swap_if_little else begin case 1 of keyword_set(append): openu,unit,filename,/swap_if_little keyword_set(update): openu,unit,filename,/swap_if_little keyword_set(write) : openw,unit,filename,/swap_if_little else : openr,unit,filename,/swap_if_little endcase endelse file = fstat(unit) fname = file.name ;In case the user input a wildcard docompress = file.compress ; Need to spawn to "gzip -l" to get the number of uncompressed bytes in a gzip ; compressed file. I'm not sure how independent this code is with different ; shells and OS's. if docompress then begin spawn,'gzip -l ' + fname, output output = strtrim(output,2) g = where(strmid(output,0,8) EQ 'compress') nbytes_in_file = long64((strsplit(output[g[0]+1],/extract))[1]) endif else nbytes_in_file = file.size ; ; create vectors needed to store header information for each extension ; n = 100 xtension = strarr(n) extname = strarr(n) extver = lonarr(n) extlevel = lonarr(n) gcount = lonarr(n) pcount = lonarr(n) bitpix = lonarr(n) naxis = lonarr(n) axis = lonarr(20,n) start_header = lon64arr(n) ; starting byte in file for header start_data = lon64arr(n) ; starting byte in file for data position = long64(0) ; current byte position in file skip = long64(0) ; Amount to skip from current position ; ; read and process each header in the file if open for read or update ; extend_number = 0 ; current extension number being ; processed if open_for_read or open_for_update then begin main_header = 1 ; first header in file flag h = bytarr(80,36,/nozero) ; read buffer ; ; loop on headers in the file ; repeat begin if skip GT 0 then point_lun,unit,position start = position ; ; loop on header blocks ; first_block = 1 ; first block in header flag repeat begin if position+2879 ge nbytes_in_file then begin if extend_number eq 0 then begin message = 'EOF encountered while reading header' goto,error_exit endif print,'EOF encountered reading extension header' print,'Only '+strtrim(extend_number-1,2) + $ ' extensions processed' goto,done_headers endif readu,unit,h position = position + 2880 hdr = string(h>32b) endline = where(strmid(hdr,0,8) eq 'END ',nend) if nend gt 0 then hdr = hdr[0:endline[0]] if first_block then begin ; ; check for valid header (SIMPLE keyword must be first for PDU and ; XTENSION keyword for the extensions. ; header = hdr keyword = strmid(header[0],0,8) if (extend_number eq 0) and $ (keyword ne 'SIMPLE ') then begin message = 'Invalid header, no SIMPLE keyword' goto,error_exit endif if (extend_number gt 0) and $ (keyword ne 'XTENSION') then begin print,'Invalid extension header encountered' print,'XTENSION keyword missing' print,'Only '+strtrim(extend_number-1,2) + $ ' extensions processed' goto,done_headers endif end else header = [header,hdr] first_block = 0 end until (nend gt 0) ; ; print header if hprint set ; if keyword_set(hprint) then hprint,header ; ; end of loop on header blocks ; ; Increase size of vectors if needed ; if extend_number ge n then begin xtension = [xtension,strarr(n)] extname = [extname,strarr(n)] extver = [extver,lonarr(n)] extlevel = [extver,lonarr(n)] gcount = [gcount,lonarr(n)] pcount = [pcount,lonarr(n)] bitpix = [bitpix,lonarr(n)] naxis = [naxis,lonarr(n)] old_axis = axis axis = lonarr(20,n*2) axis[0,0] = old_axis start_header = [start_header,lonarr(n)] start_data = [start_data,lonarr(n)] n = n*2 end ; ; extract information from header ; xtension[extend_number] = strtrim(sxpar(header,'xtension')) st = sxpar(header,'extname', Count = N_extname) if N_extname EQ 0 then st = '' extname[extend_number] = strtrim(st,2) extver[extend_number] = sxpar(header,'extver') extlevel[extend_number] = sxpar(header,'extlevel') gcount[extend_number] = sxpar(header,'gcount') pcount[extend_number] = sxpar(header,'pcount') bitpix[extend_number] = sxpar(header,'bitpix') nax = sxpar(header,'naxis') naxis[extend_number] = nax if nax gt 0 then for i=1,nax do $ axis[i-1,extend_number] = sxpar(header,'naxis'+strtrim(i,2)) start_data[extend_number] = position start_header[extend_number] = start ; ; if first header, save without FITS required keywords ; if extend_number eq 0 then begin hmain = header random_groups = sxpar(header,'groups') sxdelpar,hmain,['SIMPLE','BITPIX','NAXIS','NAXIS1', $ 'NAXIS2','NAXIS3','NAXIS4','NAXIS5', $ 'NAXIS6','NAXIS7','NAXIS8','EXTEND', $ 'PCOUNT','GCOUNT','GROUPS','BSCALE', $ 'BZERO','NPIX1','NPIX2','PIXVALUE'] if (pcount[0] gt 0) then for i=1,pcount[0] do $ sxdelpar,hmain,['ptype','pscal','pzero']+strtrim(i,2) endif ; ; skip past data to go to next header ; ndata = long64(axis[0,extend_number]) if naxis[extend_number] gt 1 then $ for i=2,naxis[extend_number] do $ ndata = ndata*axis[i-1,extend_number] nbytes = (abs(bitpix[extend_number])/8) * $ (gcount[extend_number]>1)*(pcount[extend_number] + ndata) skip = (nbytes + 2879)/2880*2880 position = position + skip ; ; end loop on headers ; extend_number = extend_number + 1 end until (position ge nbytes_in_file-2879) end ; ; point at end of file in /extend ; done_headers: if open_for_update then point_lun,unit,nbytes_in_file ; ; number of extensions ; if open_for_write then nextend = -1 $ else nextend = extend_number - 1 ; ; set up blank hmain if open for write ; if open_for_write then begin hmain = strarr(1) hmain[0] = 'END ' end ; ; create output structure for the file control block ; if open_for_write or open_for_update then begin fcb = {filename:fname,unit:unit,nextend:nextend, $ open_for_write:open_for_write + open_for_update*2} end else begin nx = nextend fcb = {filename:fname,unit:unit,nextend:nextend, $ xtension:xtension[0:nx],extname:extname[0:nx], $ extver:extver[0:nx],extlevel:extlevel[0:nx], $ gcount:gcount[0:nx],pcount:pcount[0:nx], $ bitpix:bitpix[0:nx],naxis:naxis[0:nx], $ axis:axis[*,0:nx], $ start_header:start_header[0:nx], $ start_data:start_data[0:nx],hmain:hmain, $ open_for_write:open_for_overwrite*3,$ last_extension:-1, $ random_groups:random_groups, $ nbytes: nbytes_in_file } end !err = 1 ;For obsolete users still using !err return ; ; error exit ; ioerror: message = !ERROR_STATE.msg error_exit: free_lun,unit !err = -1 if keyword_set(no_abort) then return message,' ERROR: '+message,/CON return end ; ; pro fits_write,file_or_fcb,data,header_in,extname=extname,extver=extver, $ xtension=xtension, extlevel=extlevel, $ no_abort=no_abort, message = message, header = header, $ no_data = no_data ;+ ; NAME: ; FITS_WRITE ; ; PURPOSE: ; To write a FITS primary data unit or extension. ; ; EXPLANATION: ; ***NOTE** This version of FITS_READ must be used with a post Sep 2006 ; version of FITS_OPEN. ; ; CALLING SEQUENCE: ; FITS_WRITE, filename_or_fcb, data, [header_in] ; ; INPUTS: ; FILENAME_OR_FCB: name of the output data file or the FITS control ; block returned by FITS_OPEN (called with the /WRITE or ; /APPEND) parameters. ; ; OPTIONAL INPUTS: ; DATA: data array to write. If not supplied or set to a scalar, a ; null image is written. ; HEADER_IN: FITS header keyword. If not supplied, a minimal basic ; header will be created. Required FITS keywords, SIMPLE, ; BITPIX, XTENSION, NAXIS, ... are added by FITS_WRITE and ; do not need to be supplied with the header. If supplied, ; their values will be updated as necessary to reflect DATA. ; ; INPUT KEYWORD PARAMETERS: ; ; XTENSION: type of extension to write (Default="IMAGE"). If not ; supplied, it will be taken from HEADER_IN. If not in either ; place, the default is "IMAGE". This parameter is ignored ; when writing the primary data unit. Note that binary and ; and ASCII table extensions already have a properly formatted ; header (e.g. with TTYPE* keywords) and byte array data. ; EXTNAME: EXTNAME for the extension. If not supplied, it will be taken ; from HEADER_IN. If not supplied and not in HEADER_IN, no ; EXTNAME will be written into the output extension. ; EXTVER: EXTVER for the extension. If not supplied, it will be taken ; from HEADER_IN. If not supplied and not in HEADER_IN, no ; EXTVER will be written into the output extension. ; EXTLEVEL: EXTLEVEL for the extension. If not supplied, it will be taken ; from HEADER_IN. If not supplied and not in HEADER_IN, no ; EXTLEVEL will be written into the output extension. ; /NO_ABORT: Set to return to calling program instead of a RETALL ; when an I/O error is encountered. If set, the routine will ; return a non-null string (containing the error message) in the ; keyword MESSAGE. If /NO_ABORT not set, then FITS_WRITE will ; print the message and issue a RETALL ; /NO_DATA: Set if you only want FITS_WRITE to write a header. The ; header supplied will be written without modification and ; the user is expected to write the data using WRITEU to unit ; FCB.UNIT. When FITS_WRITE is called with /NO_DATA, the user is ; responsible for the validity of the header, and must write ; the correct amount and format of the data. When FITS_WRITE ; is used in this fashion, it will pad the data from a previously ; written extension to 2880 blocks before writting the header. ; ; OUTPUT KEYWORD PARAMETERS: ; MESSAGE: value of the error message for use with /NO_ABORT ; HEADER: actual output header written to the FITS file. ; ; NOTES: ; If the first call to FITS_WRITE is an extension, FITS_WRITE will ; automatically write a null image as the primary data unit. ; ; Keywords and history in the input header will be properly separated ; into the primary data unit and extension portions when constructing ; the output header (See FITS_READ for information on the internal ; Header format which separates the extension and PDU header portions). ; ; EXAMPLES: ; Write an IDL variable to a FITS file with the minimal required header. ; FITS_WRITE,'newfile.fits',ARRAY ; ; Write the same array as an image extension, with a null Primary data ; unit. ; FITS_WRITE,'newfile.fits',ARRAY,xtension='IMAGE' ; ; Write 4 additional image extensions to the same file. ; FITS_OPEN,'newfile.fits',fcb ; FITS_WRITE,fcb,data1,extname='FLUX',extver=1 ; FITS_WRITE,fcb,err1,extname'ERR',extver=1 ; FITS_WRITE,fcb,data2,extname='FLUX',extver=2 ; FITS_WRITE,fcb,err2,extname='ERR',extver=2 ; FITS_CLOSE,FCB ; ; PROCEDURES USED: ; FITS_OPEN, SXADDPAR, SXDELPAR, SXPAR() ; HISTORY: ; Written by: D. Lindler August, 1995 ; Work for variable length extensions W. Landsman August 1997 ; Converted to IDL V5.0 W. Landsman September 1997 ; PCOUNT and GCOUNT added for IMAGE extensions J. Graham October 1999 ; Write unsigned data types W. Landsman December 1999 ; Pad data area with zeros not blanks W. McCann/W. Landsman October 2000 ; Return Message='' to signal normal operation W. Landsman Nov. 2000 ; Ensure that required extension table keywords are in proper order ; W.V. Dixon/W. Landsman March 2001 ; Assume since V5.1, remove NaNValue keyword W. Landsman Nov. 2002 ; Removed obsolete !ERR system variable W. Landsman Feb 2004 ; Check that byte array supplied with table extension W. Landsman Mar 2004 ; Make number of bytes 64bit to avoid possible overflow W.L Apr 2006 ; Asuume FITS_OPEN has opened the file with /SWAP_IF_LITTLE_ENDIAN ; W. Landsman September 2006 ;- ;----------------------------------------------------------------------------- ; ; print calling sequence if no parameters supplied ; if n_params() lt 1 then begin print,'Calling Sequence: FITS_WRITE,file_or_fcb,data,header_in' print,'Input Keywords: extname, extver, xtension, extlevel,' + $ '/no_abort, /no_data' print,'Output Keywords: message, header ' return end ; ; Open file if file name is supplied instead of a FCB ; message = '' s = size(file_or_fcb) & fcbtype = s[s[0]+1] fcbsize = n_elements(file_or_fcb) if (fcbsize ne 1) or ((fcbtype ne 7) and (fcbtype ne 8)) then begin message = 'Invalid Filename or FCB supplied' goto,error_exit end if fcbtype eq 7 then begin if keyword_set(no_data) then begin print,'FITS_WRITE: Must have FCB supplied for NO_DATA' retall endif fits_open,file_or_fcb,fcb,/write, $ no_abort=no_abort,message=message if message NE '' then goto,error_exit end else fcb = file_or_fcb ; ; if user did not pad data to 2880 blocks, pad it now ; point_lun,-fcb.unit,current_position npad = 2880 - (current_position mod 2880) if npad eq 2880 then npad = 0 if npad gt 0 then writeu,fcb.unit,bytarr(npad) ; ; if no_data, just go and write user header as supplied ; if keyword_set(no_data) then begin header = header_in goto,write_header end ; ; if header not supplied then set it to a null header ; if n_elements(header_in) le 1 then begin header = strarr(1) header[0] = 'END ' end else header = header_in ; ; on I/O error go to statement IOERROR ; ; on_ioerror,ioerror ; ; verify file is open for writing ; if fcb.open_for_write eq 0 then begin message,'File is not open for writing' goto,error_exit endif ; ; determine bitpix and axis information ; s = size(data) naxis = s[0] if naxis gt 0 then axis = s[1:naxis] idltype = s[naxis+1] if (idltype gt 5) and (idltype NE 12) and (idltype NE 13) then begin message='Data array is an invalid type' goto,error_exit endif bitpixs = [8,8,16,32,-32,-64,0,0,0,0,0,0,16,32] bitpix = bitpixs[idltype] ; ; determine extname, extver, xtension and extlevel and delete current values ; if n_elements(xtension) gt 0 then begin Axtension = xtension end else begin Axtension = sxpar(header,'xtension', Count = N_Axtension) if N_Axtension EQ 0 then Axtension = '' end if Axtension EQ 'BINTABLE' or (Axtension EQ 'TABLE') then $ if idltype GT 1 then begin message='A Byte array must be supplied with a ' + $ 'BINTABLE or TABLE extension' goto, error_exit endif if n_elements(extname) gt 0 then begin Aextname = extname end else begin Aextname = sxpar(header,'extname', Count = N_Aextname) if N_Aextname EQ 0 then Aextname = '' end if n_elements(extver) gt 0 then $ Aextver = extver $ else Aextver = sxpar(header,'extver') if n_elements(extlevel) gt 0 then $ Aextlevel = extlevel $ else Aextlevel = sxpar(header,'extlevel') sxdelpar,header,['XTENSION','EXTNAME','EXTVER','EXTLEVEL'] ; ; separate header into main and extension header ; keywords = strmid(header,0,8) hpos1 = where(keywords eq 'BEGIN MA') & hpos1 = hpos1[0] ;begin main hpos2 = where(keywords eq 'BEGIN EX') & hpos2 = hpos2[0] ;begin ext. hpos3 = where(keywords eq 'END ') & hpos3 = hpos3[0] ;end of header if (hpos1 gt 0) and (hpos2 lt hpos1) then begin message,'Invalid header BEGIN EXTENSION HEADER ... out of place' goto,error_exit endif if (hpos3 lt 0) then begin print,'FITS_WRITE: END missing from input header and was added' header = [header,'END '] hpos2 = n_elements(header)-1 end ; ; determine if a extension was supplied and no primary data unit (PDU) ; was written ; if (fcb.nextend eq -1) then begin ;no pdu written yet? if (hpos2 gt 0) or (Axtension ne '') or (Aextname ne '') or $ (Aextver ne 0) or (Aextlevel ne 0) then begin ; ; write null image PDU ; if (hpos1 gt 0) and (hpos2 gt (hpos1+1)) then $ hmain = [header[hpos1+1:hpos2-1],'END '] fits_write,fcb,0,hmain,/no_abort,message=message if message NE '' then goto,error_exit end end ; ; For extensions, do not use PDU portion of the header ; if (hpos2 gt 0) then header = header[hpos2+1:hpos3] ; ; create required keywords for the header ; h = strarr(20) h[0] = 'END ' if fcb.nextend eq -1 then begin sxaddpar,h,'SIMPLE','T','image conforms to FITS standard' end else begin if Axtension eq '' then Axtension = 'IMAGE ' sxaddpar,h,'XTENSION',Axtension,'extension type' end sxaddpar,h,'BITPIX',bitpix,'bits per data value' sxaddpar,h,'NAXIS',naxis,'number of axes' if naxis gt 0 then for i=1,naxis do $ sxaddpar,h,'NAXIS'+strtrim(i,2),axis[i-1] if fcb.nextend eq -1 then begin sxaddpar,h,'EXTEND','T','file may contain extensions' end else begin ;PCOUNT, GCOUNT are mandatory for extensions sxaddpar,h,'PCOUNT',0 sxaddpar,h,'GCOUNT',1 if (Axtension eq 'BINTABLE') or $ (Axtension eq 'TABLE ') then begin tfields = sxpar(header,'TFIELDS') > 0 sxaddpar,h,'TFIELDS',tfields endif if Aextname ne '' then sxaddpar,h,'EXTNAME',Aextname if Aextver gt 0 then sxaddpar,h,'EXTVER',Aextver if Aextlevel gt 0 then sxaddpar,h,'EXTLEVEL',Aextlevel endelse if idltype EQ 12 then $ sxaddpar,header,'BZERO',32768,'Data is unsigned integer' if idltype EQ 13 then $ sxaddpar,header,'BZERO',2147483648,'Data is unsigned long' if idltype GE 12 then sxdelpar,header,'BSCALE' ; ; delete special keywords from user supplied header ; pcount = sxpar(header,'pcount') groups = sxpar(header,'groups') sxdelpar,header,['SIMPLE','BITPIX','NAXIS','NAXIS1','NAXIS2','NAXIS3', $ 'NAXIS4','NAXIS5','NAXIS6','NAXIS7','NAXIS8','EXTEND', $ 'PCOUNT','GCOUNT','GROUPS','TFIELDS'] if groups then if (pcount gt 0) then for i=1,pcount do $ sxdelpar,header,['ptype','pscal','pzero']+strtrim(i,2) ; ; combine the two headers ; last = where(strmid(h,0,8) eq 'END ') header = [h[0:last[0]-1],header] ; ; convert header to bytes and write ; write_header: last = where(strmid(header,0,8) eq 'END ') n = last[0] + 1 byte_header = replicate(32b,80,n) for i=0,n-1 do byte_header[0,i] = byte(header[i]) writeu,fcb.unit,byte_header ; ; pad header to 2880 byte records ; npad = 2880 - (80L*n mod 2880) if npad eq 2880 then npad = 0 if (npad gt 0) then writeu,fcb.unit,replicate(32b,npad) nbytes_header = npad + n*80 if keyword_set(no_data) then return ; ; process data ; if naxis gt 0 then begin ; ; convert to IEEE ; unsigned = (idltype EQ 12) or (idltype EQ 13) if idltype EQ 12 then newdata = fix(data - 32768) if idltype EQ 13 then newdata = long(data - 2147483648) ; ; write the data ; nbytes = long64(N_elements(data)) * (abs(bitpix)/8) npad = 2880 - (nbytes mod 2880) if npad eq 2880 then npad = 0 if unsigned then writeu,fcb.unit,newdata else writeu,fcb.unit,data if npad gt 0 then begin if Axtension EQ 'TABLE ' then padnum = 32b else padnum = 0b writeu,fcb.unit,replicate(padnum,npad) endif nbytes_data = nbytes + npad end else begin nbytes_data = 0 end ; ; done, update file control block ; fcb.nextend = fcb.nextend + 1 if fcbtype eq 7 then fits_close,fcb else file_or_fcb = fcb !err = 1 return ; ; error exit ; ioerror: message = !error_state.msg error_exit: if fcbtype eq 7 then free_lun,fcb.unit !err = -1 if keyword_set(no_abort) then return message,' ERROR: '+message,/CON retall end ; function HEADFITS, filename, EXTEN = exten, Compress = compress, $ ERRMSG = errmsg, SILENT = silent ;+ ; NAME: ; HEADFITS ; PURPOSE: ; Read a FITS (primary or extension) header into a string array. ; EXPLANATION: ; HEADFITS() can also read gzip (.gz) or Unix compressed (.Z) FITS files. ; ; CALLING SEQUENCE: ; Result = HEADFITS(Filename/Fileunit ,[ ERRMSG =, EXTEN= , COMPRESS=, ; /SILENT ]) ; ; INPUTS: ; Filename = String containing the name of the FITS file to be read. ; File names ending in '.gz' are assumed to be gzip'ed compressed ; and under Unix file names ending in '.Z' are assumed to be ; Unix compressed, and file names ending in .bz2 are assumed to ; be bzip2 compressed. If this default behaviour is not ; sufficient then use the COMPRESS keyword. ; or ; Fileunit - A scalar integer specifying the unit of an already opened ; FITS file. The unit will remain open after exiting ; HEADFITS(). There are two possible reasons for choosing ; to specify a unit number rather than a file name: ; (1) For a FITS file with many extensions, one can move to the ; desired extensions with FXPOSIT() and then use HEADFITS(). This ; is more efficient that repeatedly starting at the beginning of ; the file. ; (2) For reading a FITS file across a Web http: address after opening ; the unit with the SOCKET procedure (IDL V5.4 or later, ; Unix and Windows only) ; OPTIONAL INPUT KEYWORDS: ; EXTEN = integer scalar, specifying which FITS extension to read. ; For example, to read the header of the first extension set ; EXTEN = 1. Default is to read the primary FITS header ; (EXTEN = 0). The EXTEN keyword cannot be used when a unit ; number is supplied instead of a file name. ; COMPRESS - If this keyword is set and non-zero, then treat the file ; as compressed. If 1 assume a gzipped file. Use IDL's ; internal decompression facilities for gzip files, while for ; Unix or bzip2 compression spawn off a process to decompress and ; use its output as the FITS stream. If the keyword is not 1, ; then use its value as a string giving the command needed for ; decompression. See FXPOSIT for more info. ; /SILENT - If set, then suppress any warning messages about invalid ; characters in the FITS file. ; OPTIONAL KEYWORD OUTPUT: ; ERRMSG = If this keyword is present, then any error messages will be ; returned to the user in this parameter rather than ; depending on the MESSAGE routine in IDL. If no errors are ; encountered, then a null string is returned. ; ; OUTPUTS: ; Result of function = FITS header, string array ; ; EXAMPLE: ; Print the main FITS header of a file 'test.fits' into a string ; variable, h ; ; IDL> print, headfits( 'test.fits') ; ; Print the second extension header of a gzip compressed FITS file ; 'test.fits.gz'. Use HPRINT for pretty format ; ; IDL> hprint, headfits( 'test.fits.gz', ext=2) ; ; PROCEDURES CALLED ; FXPOSIT(), MRD_HREAD ; MODIFICATION HISTORY: ; adapted by Frank Varosi from READFITS by Jim Wofford, January, 24 1989 ; Keyword EXTEN added, K.Venkatakrishna, May 1992 ; Make sure first 8 characters are 'SIMPLE' W. Landsman October 1993 ; Check PCOUNT and GCOUNT W. Landsman December 1994 ; Major rewrite, work for Unix gzip files, W. Landsman April 1996 ; Converted to IDL V5.0 W. Landsman September 1997 ; Added COMPRESS keyword W. Landsman April 2000 ; Added ERRMSG keyword W. Landsman July 2000 ; Added /SILENT keyword W. Landsman December 2000 ; Option to read a unit number rather than file name W.L October 2001 ; Test output status of MRD_HREAD call October 2003 W. Landsman ;- On_error,2 if N_params() LT 1 then begin print,'Syntax - header = headfits( filename,[ EXTEN=, ERRMSG=, ' + $ '/SILENT, COMPRESS= ])' return, -1 endif if arg_present(errmsg) then errmsg = '' if not keyword_set(exten) then exten = 0 unitsupplied = size(filename,/TNAME) NE 'STRING' if unitsupplied then unit = filename else begin unit = FXPOSIT( filename, exten, $ /READONLY,compress = compress, SILENT=silent) if unit EQ -1 then begin message = 'Unable to open file ' + filename if N_elements(errmsg) GT 0 then errmsg = message else $ message,'ERROR - ' + message,/CON return,-1 endif if eof(unit) then begin free_lun,unit message = 'Extension past EOF' if N_elements(errmsg) GT 0 then errmsg = message else $ message,'ERROR - ' + message,/CON return,-1 endif endelse MRD_HREAD, unit, header, status, SILENT = silent if not unitsupplied then free_lun, unit if status LT 0 then begin if N_elements(errmsg) GT 0 then errmsg = !ERROR_STATE.MSG else $ message,'ERROR - ' + !ERROR_STATE.MSG,/CON return, -1 endif else return, header end ;+ ; NAME: ; FSC_COLOR ; ; PURPOSE: ; ; The purpose of this function is to obtain drawing colors ; by name and in a device/decomposition independent way. ; The color names and values may be read in as a file, or 104 color ; names and values are supplied with the program. These colors were ; obtained from the file rgb.txt, found on most X-Window distributions. ; Representative colors were chosen from across the color spectrum. To ; see a list of colors available, type: ; ; Print, FSC_Color(/Names), Format='(6A18)' ; ; AUTHOR: ; ; FANNING SOFTWARE CONSULTING: ; David Fanning, Ph.D. ; 1645 Sheely Drive ; Fort Collins, CO 80526 USA ; Phone: 970-221-0438 ; E-mail: davidf@dfanning.com ; Coyote's Guide to IDL Programming: http://www.dfanning.com ; ; CATEGORY: ; ; Graphics, Color Specification. ; ; CALLING SEQUENCE: ; ; color = FSC_Color(theColor, theColorIndex) ; ; NORMAL CALLING SEQUENCE FOR DEVICE-INDEPENDENT COLOR: ; ; If you write your graphics code *exactly* as it is written below, then ; the same code will work in all graphics devices I have tested. ; These include the PRINTER, PS, and Z devices, as well as X, WIN, and MAC. ; ; In practice, graphics code is seldom written like this. (For a variety of ; reasons, but laziness is high on the list.) So I have made the ; program reasonably tolerant of poor programming practices. I just ; point this out as a place you might return to before you write me ; a nice note saying my program "doesn't work". :-) ; ; axisColor = FSC_Color("Green", !D.Table_Size-2) ; backColor = FSC_Color("Charcoal", !D.Table_Size-3) ; dataColor = FSC_Color("Yellow", !D.Table_Size-4) ; thisDevice = !D.Name ; Set_Plot, 'toWhateverYourDeviceIsGoingToBe', /Copy ; Device, .... ; Whatever you need here to set things up properly. ; IF (!D.Flags AND 256) EQ 0 THEN $ ; POLYFILL, [0,1,1,0,0], [0,0,1,1,0], /Normal, Color=backColor ; Plot, Findgen(11), Color=axisColor, Background=backColor, /NoData, $ ; NoErase= ((!D.Flags AND 256) EQ 0) ; OPlot, Findgen(11), Color=dataColor ; Device, .... ; Whatever you need here to wrap things up properly. ; Set_Plot, thisDevice ; ; OPTIONAL INPUT PARAMETERS: ; ; theColor: A string with the "name" of the color. To see a list ; of the color names available set the NAMES keyword. This may ; also be a vector of color names. Colors available are these: ; ; Active Almond Antique White Aquamarine Beige Bisque ; Black Blue Blue Violet Brown Burlywood Cadet Blue ; Charcoal Chartreuse Chocolate Coral Cornflower Blue Cornsilk ; Crimson Cyan Dark Goldenrod Dark Gray Dark Green Dark Khaki ; Dark Orchid Dark Red Dark Salmon Dark Slate Blue Deep Pink Dodger Blue ; Edge Face Firebrick Forest Green Frame Gold ; Goldenrod Gray Green Green Yellow Highlight Honeydew ; Hot Pink Indian Red Ivory Khaki Lavender Lawn Green ; Light Coral Light Cyan Light Gray Light Salmon Light Sea Green Light Yellow ; Lime Green Linen Magenta Maroon Medium Gray Medium Orchid ; Moccasin Navy Olive Olive Drab Orange Orange Red ; Orchid Pale Goldenrod Pale Green Papaya Peru Pink ; Plum Powder Blue Purple Red Rose Rosy Brown ; Royal Blue Saddle Brown Salmon Sandy Brown Sea Green Seashell ; Selected Shadow Sienna Sky Blue Slate Blue Slate Gray ; Snow Spring Green Steel Blue Tan Teal Text ; Thistle Tomato Turquoise Violet Violet Red Wheat ; White Yellow ; ; In addition, these system colors are available if a connection to the window system is available. ; ; Frame Text Active Shadow Highlight Edge Selected Face ; ; The color WHITE is used if this parameter is absent or a color name is mis-spelled. To see a list ; of the color names available in the program, type this: ; ; IDL> Print, FSC_Color(/Names), Format='(6A18)' ; ; theColorIndex: The color table index (or vector of indices the same length ; as the color name vector) where the specified color is loaded. The color table ; index parameter should always be used if you wish to obtain a color value in a ; color-decomposition-independent way in your code. See the NORMAL CALLING ; SEQUENCE for details. If theColor is a vector, and theColorIndex is a scalar, ; then the colors will be loaded starting at theColorIndex. ; ; When the BREWER keyword is set, you must use more arbitrary and less descriptive color ; names. To see a list of those names, use the command above with the BREWER keyword set, ; or call PICKCOLORNAME with the BREWER keyword set: ; ; IDL> Print, FSC_Color(/Names, /BREWER), Format='(8A10)' ; IDL> color = PickColorName(/BREWER) ; ; Here are the Brewer names: ; ; WT1 WT2 WT3 WT4 WT5 WT6 WT7 WT8 ; TAN1 TAN2 TAN3 TAN4 TAN5 TAN6 TAN7 TAN8 ; BLK1 BLK2 BLK3 BLK4 BLK5 BLK6 BLK7 BLK8 ; GRN1 GRN2 GRN3 GRN4 GRN5 GRN6 GRN7 GRN8 ; BLU1 BLU2 BLU3 BLU4 BLU5 BLU6 BLU7 BLU8 ; ORG1 ORG2 ORG3 ORG4 ORG5 ORG6 ORG7 ORG8 ; RED1 RED2 RED3 RED4 RED5 RED6 RED7 RED8 ; PUR1 PUR2 PUR3 PUR4 PUR5 PUR6 PUR7 PUR8 ; PBG1 PBG2 PBG3 PBG4 PBG5 PBG6 PBG7 PBG8 ; YGB1 YGB2 YGB3 YGB4 YGB5 YGB6 YGB7 YGB8 ; RYB1 RYB2 RYB3 RYB4 RYB5 RYB6 RYB7 RYB8 ; TG1 TG2 TG3 TG4 TG5 TG6 TG7 TG8 ; ; As of 3 July 2008, the Brewer names are also now available to the user without using ; the BREWER keyword. If the BREWER keyword is used, *only* Brewer names are available. ; ; RETURN VALUE: ; ; The value that is returned by FSC_Color depends upon the keywords ; used to call it, on the version of IDL you are using,and on the depth ; of the display device when the program is invoked. In general, ; the return value will be either a color index number where the specified ; color is loaded by the program, or a 24-bit color value that can be ; decomposed into the specified color on true-color systems. (Or a vector ; of such numbers.) ; ; If you are running IDL 5.2 or higher, the program will determine which ; return value to use, based on the color decomposition state at the time ; the program is called. If you are running a version of IDL before IDL 5.2, ; then the program will return the color index number. This behavior can ; be overruled in all versions of IDL by setting the DECOMPOSED keyword. ; If this keyword is 0, the program always returns a color index number. If ; the keyword is 1, the program always returns a 24-bit color value. ; ; If the TRIPLE keyword is set, the program always returns the color triple, ; no matter what the current decomposition state or the value of the DECOMPOSED ; keyword. Normally, the color triple is returned as a 1 by 3 column vector. ; This is appropriate for loading into a color index with TVLCT: ; ; IDL> TVLCT, FSC_Color('Yellow', /Triple), !P.Color ; ; But sometimes (e.g, in object graphics applications) you want the color ; returned as a row vector. In this case, you should set the ROW keyword ; as well as the TRIPLE keyword: ; ; viewobj= Obj_New('IDLgrView', Color=FSC_Color('charcoal', /Triple, /Row)) ; ; If the ALLCOLORS keyword is used, then instead of a single value, modified ; as described above, then all the color values are returned in an array. In ; other words, the return value will be either an NCOLORS-element vector of color ; table index numbers, an NCOLORS-element vector of 24-bit color values, or ; an NCOLORS-by-3 array of color triples. ; ; If the NAMES keyword is set, the program returns a vector of ; color names known to the program. ; ; If the color index parameter is not used, and a 24-bit value is not being ; returned, then colorIndex is typically set to !D.Table_Size-1. However, ; this behavior is changed on 8-bit devices (e.g., the PostScript device, ; or the Z-graphics buffer) and on 24-bit devices that are *not* using ; decomposed color. On these devices, the colors are loaded at an ; offset of !D.Table_Size - ncolors - 2, and the color index parameter reflects ; the actual index of the color where it will be loaded. This makes it possible ; to use a formulation as below: ; ; Plot, data, Color=FSC_Color('Dodger Blue') ; ; on 24-bit displays *and* in PostScript output! Note that if you specify a color ; index (the safest thing to do), then it will always be honored. ; ; INPUT KEYWORD PARAMETERS: ; ; ALLCOLORS: Set this keyword to return indices, or 24-bit values, or color ; triples, for all the known colors, instead of for a single color. ; ; BREWER: Set this keyword if you wish to use the Brewer Colors, as defined ; in this reference: ; ; http://www.personal.psu.edu/cab38/ColorBrewer/ColorBrewer_intro.html ; ; As of 3 July 2008, the BREWER names are always available to the user, with or ; without this keyword. If the keyword is used, only BREWER names are available. ; ; DECOMPOSED: Set this keyword to 0 or 1 to force the return value to be ; a color table index or a 24-bit color value, respectively. ; ; CHECK_CONNECTION: If this keyword is set, the program will check to see if it can obtain ; a window connection before it tries to load system colors (which require one). If you ; think you might be using FSC_COLOR in a cron job, for example, you would want to set this ; keyword. If there is no window connection, the system colors are not available from the program. ; ; FILENAME: The string name of an ASCII file that can be opened to read in ; color values and color names. There should be one color per row ; in the file. Please be sure there are no blank lines in the file. ; The format of each row should be: ; ; redValue greenValue blueValue colorName ; ; Color values should be between 0 and 255. Any kind of white-space ; separation (blank characters, commas, or tabs) are allowed. The color ; name should be a string, but it should NOT be in quotes. A typical ; entry into the file would look like this: ; ; 255 255 0 Yellow ; ; NAMES: If this keyword is set, the return value of the function is ; a ncolors-element string array containing the names of the colors. ; These names would be appropriate, for example, in building ; a list widget with the names of the colors. If the NAMES ; keyword is set, the COLOR and INDEX parameters are ignored. ; ; listID = Widget_List(baseID, Value=GetColor(/Names), YSize=16) ; ; ; NODISPLAY: Normally, FSC_COLOR loads "system" colors as part of its palette of colors. ; In order to do so, it has to create an IDL widget, which in turn has to make ; a connection to the windowing system. If your program is being run without a ; window connection, then this program will fail. If you can live without the system ; colors (and most people don't even know they are there, to tell you the truth), ; then setting this keyword will keep them from being loaded, and you can run ; FSC_COLOR without a display. THIS KEYWORD NOW DEPRECIATED IN FAVOR OF CHECK_CONNECTION. ; ; ROW: If this keyword is set, the return value of the function when the TRIPLE ; keyword is set is returned as a row vector, rather than as the default ; column vector. This is required, for example, when you are trying to ; use the return value to set the color for object graphics objects. This ; keyword is completely ignored, except when used in combination with the ; TRIPLE keyword. ; ; SELECTCOLOR: Set this keyword if you would like to select the color name with ; the PICKCOLORNAME program. Selecting this keyword automaticallys sets ; the INDEX positional parameter. If this keyword is used, any keywords ; appropriate for PICKCOLORNAME can also be used. If this keyword is used, ; the first positional parameter can be either a color name or the color ; table index number. The program will figure out what you want. ; ; TRIPLE: Setting this keyword will force the return value of the function to ; *always* be a color triple, regardless of color decomposition state or ; visual depth of the machine. The value will be a three-element column ; vector unless the ROW keyword is also set. ; ; In addition, any keyword parameter appropriate for PICKCOLORNAME can be used. ; These include BOTTOM, COLUMNS, GROUP_LEADER, INDEX, and TITLE. ; ; OUTPUT KEYWORD PARAMETERS: ; ; CANCEL: This keyword is always set to 0, unless that SELECTCOLOR keyword is used. ; Then it will correspond to the value of the CANCEL output keyword in PICKCOLORNAME. ; ; COLORSTRUCTURE: This output keyword (if set to a named variable) will return a ; structure in which the fields will be the known color names (without spaces) ; and the values of the fields will be either color table index numbers or ; 24-bit color values. If you have specified a vector of color names, then ; this will be a structure containing just those color names as fields. ; ; NCOLORS: The number of colors recognized by the program. It will be 104 by default. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; None. ; ; RESTRICTIONS: ; ; Required programs from the Coyote Library: ; ; http://www.dfanning.com/programs/error_message.pro ; http://www.dfanning.com/programs/pickcolorname.pro ; ; EXAMPLE: ; ; To get drawing colors in a device-decomposed independent way: ; ; axisColor = FSC_Color("Green", !D.Table_Size-2) ; backColor = FSC_Color("Charcoal", !D.Table_Size-3) ; dataColor = FSC_Color("Yellow", !D.Table_Size-4) ; Plot, Findgen(11), Color=axisColor, Background=backColor, /NoData ; OPlot, Findgen(11), Color=dataColor ; ; To set the viewport color in object graphics: ; ; theView = Obj_New('IDLgrView', Color=FSC_Color('Charcoal', /Triple)) ; ; To change the viewport color later: ; ; theView->SetProperty, Color=FSC_Color('Antique White', /Triple) ; ; To load the drawing colors "red", "green", and "yellow" at indices 100-102, type this: ; ; IDL> TVLCT, FSC_Color(["red", "green", and "yellow"], /Triple), 100 ; ; MODIFICATION HISTORY: ; ; Written by: David W. Fanning, 19 October 2000. Based on previous ; GetColor program. ; Fixed a problem with loading colors with TVLCT on a PRINTER device. 13 Mar 2001. DWF. ; Added the ROW keyword. 30 March 2001. DWF. ; Added the PICKCOLORNAME code to the file, since I keep forgetting to ; give it to people. 15 August 2001. DWF. ; Added ability to specify color names and indices as vectors. 5 Nov 2002. DWF. ; Fixed a problem with the TRIPLE keyword when specifying a vector of color names. 14 Feb 2003. DWF. ; Fixed a small problem with the starting index when specifying ALLCOLORS. 24 March 2003. DWF. ; Added system color names. 23 Jan 2004. DWF ; Added work-around for WHERE function "feature" when theColor is a one-element array. 22 July 2004. DWF. ; Added support for 8-bit graphics devices when color index is not specified. 25 August 2004. DWF. ; Fixed a small problem with creating color structure when ALLCOLORS keyword is set. 26 August 2004. DWF. ; Extended the color index fix for 8-bit graphics devices on 25 August 2004 to ; 24-bit devices running with color decomposition OFF. I've concluded most of ; the people using IDL don't have any idea how color works, so I am trying to ; make it VERY simple, and yet still maintain the power of this program. So now, ; in general, for most simple plots, you don't have to use the colorindex parameter ; and you still have a very good chance of getting what you expect in a device-independent ; manner. Of course, it would be *nice* if you could use that 24-bit display you paid ; all that money for, but I understand your reluctance. :-) 11 October 2004. DWF. ; Have renamed the first positional parameter so that this variable doesn't change ; while the program is running. 7 December 2004. DWF. ; Fixed an error I introduced on 7 December 2004. Sigh... 7 January 2005. DWF. ; Added eight new colors. Total now of 104 colors. 11 August 2005. DWF. ; Modified GUI to display system colors and removed PickColorName code. 13 Dec 2005. DWF. ; Fixed a problem with colorIndex when SELECTCOLOR keyword was used. 13 Dec 2005. DWF. ; Fixed a problem with color name synonyms. 19 May 2006. DWF. ; The previous fix broke the ability to specify several colors at once. Fixed. 24 July 2006. DWF. ; Updated program to work with 24-bit Z-buffer in IDL 6.4. 11 June 2007. DWF ; Added the CRONJOB keyword. 07 Feb 2008. DWF. ; Changed the CRONJOB keyword to NODISPLAY to better reflect its purpose. 7 FEB 2008. DWF. ; Added the BREWER keyword to allow selection of Brewer Colors. 15 MAY 2008. DWF. ; Added the CHECK_CONNECTION keyword and depreciated the NODISPLAY keyword for cron jobs. 15 MAY 2008. DWF. ; Added the BREWER names to the program with or without the BREWER keyword set. 3 JULY 2008. DWF. ;- ; ;########################################################################### ; ; LICENSE ; ; This software is OSI Certified Open Source Software. ; OSI Certified is a certification mark of the Open Source Initiative. ; ; Copyright 2000-2008 Fanning Software Consulting. ; ; This software is provided "as-is", without any express or ; implied warranty. In no event will the authors be held liable ; for any damages arising from the use of this software. ; ; Permission is granted to anyone to use this software for any ; purpose, including commercial applications, and to alter it and ; redistribute it freely, subject to the following restrictions: ; ; 1. The origin of this software must not be misrepresented; you must ; not claim you wrote the original software. If you use this software ; in a product, an acknowledgment in the product documentation ; would be appreciated, but is not required. ; ; 2. Altered source versions must be plainly marked as such, and must ; not be misrepresented as being the original software. ; ; 3. This notice may not be removed or altered from any source distribution. ; ; For more information on Open Source Software, visit the Open Source ; web site: http://www.opensource.org. ; ;########################################################################### FUNCTION FSC_Color_Count_Rows, filename, MaxRows = maxrows ; This utility routine is used to count the number of ; rows in an ASCII data file. IF N_Elements(maxrows) EQ 0 THEN maxrows = 500L IF N_Elements(filename) EQ 0 THEN BEGIN filename = Dialog_Pickfile() IF filename EQ "" THEN RETURN, -1 ENDIF OpenR, lun, filename, /Get_Lun Catch, theError IF theError NE 0 THEN BEGIN Catch, /Cancel count = count-1 Free_Lun, lun RETURN, count ENDIF RESTART: count = 1L line = '' FOR j=count, maxrows DO BEGIN ReadF, lun, line count = count + 1 ; Try again if you hit MAXROWS without encountering the ; end of the file. Double the MAXROWS parameter. IF j EQ maxrows THEN BEGIN maxrows = maxrows * 2 Point_Lun, lun, 0 GOTO, RESTART ENDIF ENDFOR RETURN, -1 END ;------------------------------------------------------------------------------- FUNCTION FSC_Color_Color24, color ; This FUNCTION accepts a [red, green, blue] triple that ; describes a particular color and returns a 24-bit long ; integer that is equivalent to (can be decomposed into) ; that color. The triple can be either a row or column ; vector of 3 elements or it can be an N-by-3 array of ; color triples. ON_ERROR, 2 s = Size(color) IF s[0] EQ 1 THEN BEGIN IF s[1] NE 3 THEN Message, 'Input color parameter must be a 3-element vector.' RETURN, color[0] + (color[1] * 2L^8) + (color[2] * 2L^16) ENDIF ELSE BEGIN IF s[2] GT 3 THEN Message, 'Input color parameter must be an N-by-3 array.' RETURN, color[*,0] + (color[*,1] * 2L^8) + (color[*,2] * 2L^16) ENDELSE END ;-------------------------------------------------------------------------------------------- FUNCTION FSC_Color, theColour, colorIndex, $ AllColors=allcolors, $ Brewer=brewer, $ Check_Connection=check_connection, $ ColorStructure=colorStructure, $ Cancel=cancelled, $ Decomposed=decomposedState, $ _Extra=extra, $ Filename=filename, $ Names=names, $ NColors=ncolors, $ NODISPLAY=nodisplay, $ Row=row, $ SelectColor=selectcolor, $ Triple=triple ; Do you have to check the window connection? haveConnection = 1 IF Keyword_Set(check_connection) THEN BEGIN ; In CRON jobs, there is no X connection so system colors cannot be supported. ; Here is a quick test to see if we can connect to a windowing system. Catch, theError IF theError NE 0 THEN BEGIN Catch, /CANCEL haveConnection = 0 GOTO, testConnection ENDIF theWindow = !D.Window Window, /FREE, XSIZE=5, YSIZE=5, /PIXMAP Catch, /CANCEL testConnection: ; Come here if you choke on creating a window. IF !D.Window NE theWindow THEN BEGIN WDelete, !D.Window IF theWindow GE 0 THEN WSet, theWindow ENDIF ENDIF ; Handle depreciated NODISPLAY keyword. IF Keyword_Set(nodisplay) THEN haveConnection = 0 ; Error handling for the rest of the program. Catch, theError IF theError NE 0 THEN BEGIN Catch, /Cancel ok = Error_Message(/Traceback) cancelled = 1 RETURN, !P.Color ENDIF ; Need brewer colors? brewer = Keyword_Set(brewer) ; I don't want to change the original variable. IF N_Elements(theColour) NE 0 THEN theColor = theColour ELSE $ theColor = brewer ? 'WT1' : 'WHITE' ; Load the colors. IF N_Elements(filename) NE 0 THEN BEGIN ; Count the number of rows in the file. ncolors = FSC_Color_Count_Rows(filename) ; Read the data. OpenR, lun, filename, /Get_Lun rvalue = BytArr(NCOLORS) gvalue = BytArr(NCOLORS) bvalue = BytArr(NCOLORS) colors = StrArr(NCOLORS) redvalue = 0B greenvalue = 0B bluevalue = 0B colorvalue = "" FOR j=0L, NCOLORS-1 DO BEGIN ReadF, lun, redvalue, greenvalue, bluevalue, colorvalue rvalue[j] = redvalue gvalue[j] = greenvalue bvalue[j] = bluevalue colors[j] = colorvalue ENDFOR Free_Lun, lun ; Trim the colors array of blank characters. colors = StrTrim(colors, 2) ENDIF ELSE BEGIN ; Set up the color vectors. IF Keyword_Set(Brewer) THEN BEGIN ; Set up the color vectors. colors = [ 'WT1', 'WT2', 'WT3', 'WT4', 'WT5', 'WT6', 'WT7', 'WT8'] rvalue = [ 255, 255, 255, 255, 255, 245, 255, 250 ] gvalue = [ 255, 250, 255, 255, 248, 245, 245, 240 ] bvalue = [ 255, 250, 240, 224, 220, 220, 238, 230 ] colors = [ colors, 'TAN1', 'TAN2', 'TAN3', 'TAN4', 'TAN5', 'TAN6', 'TAN7', 'TAN8'] rvalue = [ rvalue, 250, 255, 255, 255, 255, 245, 222, 210 ] gvalue = [ gvalue, 235, 239, 235, 228, 228, 222, 184, 180 ] bvalue = [ bvalue, 215, 213, 205, 196, 181, 179, 135, 140 ] colors = [ colors, 'BLK1', 'BLK2', 'BLK3', 'BLK4', 'BLK5', 'BLK6', 'BLK7', 'BLK8'] rvalue = [ rvalue, 250, 230, 210, 190, 112, 110, 70, 0 ] gvalue = [ gvalue, 250, 230, 210, 190, 128, 110, 70, 0 ] bvalue = [ bvalue, 250, 230, 210, 190, 128, 110, 70, 0 ] colors = [ colors, 'GRN1', 'GRN2', 'GRN3', 'GRN4', 'GRN5', 'GRN6', 'GRN7', 'GRN8'] rvalue = [ rvalue, 250, 223, 173, 109, 53, 35, 0, 0 ] gvalue = [ gvalue, 253, 242, 221, 193, 156, 132, 97, 69 ] bvalue = [ bvalue, 202, 167, 142, 115, 83, 67, 52, 41 ] colors = [ colors, 'BLU1', 'BLU2', 'BLU3', 'BLU4', 'BLU5', 'BLU6', 'BLU7', 'BLU8'] rvalue = [ rvalue, 232, 202, 158, 99, 53, 33, 8, 8 ] gvalue = [ gvalue, 241, 222, 202, 168, 133, 113, 75, 48 ] bvalue = [ bvalue, 250, 240, 225, 211, 191, 181, 147, 107 ] colors = [ colors, 'ORG1', 'ORG2', 'ORG3', 'ORG4', 'ORG5', 'ORG6', 'ORG7', 'ORG8'] rvalue = [ rvalue, 254, 253, 253, 250, 231, 217, 159, 127 ] gvalue = [ gvalue, 236, 212, 174, 134, 92, 72, 51, 39 ] bvalue = [ bvalue, 217, 171, 107, 52, 12, 1, 3, 4 ] colors = [ colors, 'RED1', 'RED2', 'RED3', 'RED4', 'RED5', 'RED6', 'RED7', 'RED8'] rvalue = [ rvalue, 254, 252, 252, 248, 225, 203, 154, 103 ] gvalue = [ gvalue, 232, 194, 146, 97, 45, 24, 12, 0 ] bvalue = [ bvalue, 222, 171, 114, 68, 38, 29, 19, 13 ] colors = [ colors, 'PUR1', 'PUR2', 'PUR3', 'PUR4', 'PUR5', 'PUR6', 'PUR7', 'PUR8'] rvalue = [ rvalue, 244, 222, 188, 152, 119, 106, 80, 63 ] gvalue = [ gvalue, 242, 221, 189, 148, 108, 82, 32, 0 ] bvalue = [ bvalue, 248, 237, 220, 197, 177, 163, 139, 125 ] colors = [ colors, 'PBG1', 'PBG2', 'PBG3', 'PBG4', 'PBG5', 'PBG6', 'PBG7', 'PBG8'] rvalue = [ rvalue, 243, 213, 166, 94, 34, 3, 1, 1 ] gvalue = [ gvalue, 234, 212, 189, 164, 138, 129, 101, 70 ] bvalue = [ bvalue, 244, 232, 219, 204, 171, 139, 82, 54 ] colors = [ colors, 'YGB1', 'YGB2', 'YGB3', 'YGB4', 'YGB5', 'YGB6', 'YGB7', 'YGB8'] rvalue = [ rvalue, 244, 206, 127, 58, 30, 33, 32, 8 ] gvalue = [ gvalue, 250, 236, 205, 175, 125, 95, 48, 29 ] bvalue = [ bvalue, 193, 179, 186, 195, 182, 168, 137, 88 ] colors = [ colors, 'RYB1', 'RYB2', 'RYB3', 'RYB4', 'RYB5', 'RYB6', 'RYB7', 'RYB8'] rvalue = [ rvalue, 201, 245, 253, 251, 228, 193, 114, 59 ] gvalue = [ gvalue, 35, 121, 206, 253, 244, 228, 171, 85 ] bvalue = [ bvalue, 38, 72, 127, 197, 239, 239, 207, 164 ] colors = [ colors, 'TG1', 'TG2', 'TG3', 'TG4', 'TG5', 'TG6', 'TG7', 'TG8'] rvalue = [ rvalue, 84, 163, 197, 220, 105, 51, 13, 0 ] gvalue = [ gvalue, 48, 103, 141, 188, 188, 149, 113, 81 ] bvalue = [ bvalue, 5, 26, 60, 118, 177, 141, 105, 71 ] ENDIF ELSE BEGIN ; Set up the color vectors. Both original and Brewer colors. colors= ['White'] rvalue = [ 255] gvalue = [ 255] bvalue = [ 255] colors = [ colors, 'Snow', 'Ivory','Light Yellow', 'Cornsilk', 'Beige', 'Seashell' ] rvalue = [ rvalue, 255, 255, 255, 255, 245, 255 ] gvalue = [ gvalue, 250, 255, 255, 248, 245, 245 ] bvalue = [ bvalue, 250, 240, 224, 220, 220, 238 ] colors = [ colors, 'Linen','Antique White','Papaya', 'Almond', 'Bisque', 'Moccasin' ] rvalue = [ rvalue, 250, 250, 255, 255, 255, 255 ] gvalue = [ gvalue, 240, 235, 239, 235, 228, 228 ] bvalue = [ bvalue, 230, 215, 213, 205, 196, 181 ] colors = [ colors, 'Wheat', 'Burlywood', 'Tan', 'Light Gray', 'Lavender','Medium Gray' ] rvalue = [ rvalue, 245, 222, 210, 230, 230, 210 ] gvalue = [ gvalue, 222, 184, 180, 230, 230, 210 ] bvalue = [ bvalue, 179, 135, 140, 230, 250, 210 ] colors = [ colors, 'Gray', 'Slate Gray', 'Dark Gray', 'Charcoal', 'Black', 'Honeydew', 'Light Cyan' ] rvalue = [ rvalue, 190, 112, 110, 70, 0, 240, 224 ] gvalue = [ gvalue, 190, 128, 110, 70, 0, 255, 255 ] bvalue = [ bvalue, 190, 144, 110, 70, 0, 255, 240 ] colors = [ colors,'Powder Blue', 'Sky Blue', 'Cornflower Blue', 'Cadet Blue', 'Steel Blue','Dodger Blue', 'Royal Blue', 'Blue' ] rvalue = [ rvalue, 176, 135, 100, 95, 70, 30, 65, 0 ] gvalue = [ gvalue, 224, 206, 149, 158, 130, 144, 105, 0 ] bvalue = [ bvalue, 230, 235, 237, 160, 180, 255, 225, 255 ] colors = [ colors, 'Navy', 'Pale Green','Aquamarine','Spring Green', 'Cyan' ] rvalue = [ rvalue, 0, 152, 127, 0, 0 ] gvalue = [ gvalue, 0, 251, 255, 250, 255 ] bvalue = [ bvalue, 128, 152, 212, 154, 255 ] colors = [ colors, 'Turquoise', 'Light Sea Green', 'Sea Green','Forest Green', 'Teal','Green Yellow','Chartreuse', 'Lawn Green' ] rvalue = [ rvalue, 64, 143, 46, 34, 0, 173, 127, 124 ] gvalue = [ gvalue, 224, 188, 139, 139, 128, 255, 255, 252 ] bvalue = [ bvalue, 208, 143, 87, 34, 128, 47, 0, 0 ] colors = [ colors, 'Green', 'Lime Green', 'Olive Drab', 'Olive','Dark Green','Pale Goldenrod'] rvalue = [ rvalue, 0, 50, 107, 85, 0, 238 ] gvalue = [ gvalue, 255, 205, 142, 107, 100, 232 ] bvalue = [ bvalue, 0, 50, 35, 47, 0, 170 ] colors = [ colors, 'Khaki', 'Dark Khaki', 'Yellow', 'Gold', 'Goldenrod','Dark Goldenrod'] rvalue = [ rvalue, 240, 189, 255, 255, 218, 184 ] gvalue = [ gvalue, 230, 183, 255, 215, 165, 134 ] bvalue = [ bvalue, 140, 107, 0, 0, 32, 11 ] colors = [ colors,'Saddle Brown', 'Rose', 'Pink', 'Rosy Brown','Sandy Brown', 'Peru'] rvalue = [ rvalue, 139, 255, 255, 188, 244, 205 ] gvalue = [ gvalue, 69, 228, 192, 143, 164, 133 ] bvalue = [ bvalue, 19, 225, 203, 143, 96, 63 ] colors = [ colors,'Indian Red', 'Chocolate', 'Sienna','Dark Salmon', 'Salmon','Light Salmon' ] rvalue = [ rvalue, 205, 210, 160, 233, 250, 255 ] gvalue = [ gvalue, 92, 105, 82, 150, 128, 160 ] bvalue = [ bvalue, 92, 30, 45, 122, 114, 122 ] colors = [ colors, 'Orange', 'Coral', 'Light Coral', 'Firebrick', 'Dark Red', 'Brown', 'Hot Pink' ] rvalue = [ rvalue, 255, 255, 240, 178, 139, 165, 255 ] gvalue = [ gvalue, 165, 127, 128, 34, 0, 42, 105 ] bvalue = [ bvalue, 0, 80, 128, 34, 0, 42, 180 ] colors = [ colors, 'Deep Pink', 'Magenta', 'Tomato', 'Orange Red', 'Red', 'Crimson', 'Violet Red' ] rvalue = [ rvalue, 255, 255, 255, 255, 255, 220, 208 ] gvalue = [ gvalue, 20, 0, 99, 69, 0, 20, 32 ] bvalue = [ bvalue, 147, 255, 71, 0, 0, 60, 144 ] colors = [ colors, 'Maroon', 'Thistle', 'Plum', 'Violet', 'Orchid','Medium Orchid'] rvalue = [ rvalue, 176, 216, 221, 238, 218, 186 ] gvalue = [ gvalue, 48, 191, 160, 130, 112, 85 ] bvalue = [ bvalue, 96, 216, 221, 238, 214, 211 ] colors = [ colors,'Dark Orchid','Blue Violet', 'Purple'] rvalue = [ rvalue, 153, 138, 160] gvalue = [ gvalue, 50, 43, 32] bvalue = [ bvalue, 204, 226, 240] colors = [ colors, 'Slate Blue', 'Dark Slate Blue'] rvalue = [ rvalue, 106, 72] gvalue = [ gvalue, 90, 61] bvalue = [ bvalue, 205, 139] colors = [ colors, 'WT1', 'WT2', 'WT3', 'WT4', 'WT5', 'WT6', 'WT7', 'WT8'] rvalue = [ rvalue, 255, 255, 255, 255, 255, 245, 255, 250 ] gvalue = [ gvalue, 255, 250, 255, 255, 248, 245, 245, 240 ] bvalue = [ bvalue, 255, 250, 240, 224, 220, 220, 238, 230 ] colors = [ colors, 'TAN1', 'TAN2', 'TAN3', 'TAN4', 'TAN5', 'TAN6', 'TAN7', 'TAN8'] rvalue = [ rvalue, 250, 255, 255, 255, 255, 245, 222, 210 ] gvalue = [ gvalue, 235, 239, 235, 228, 228, 222, 184, 180 ] bvalue = [ bvalue, 215, 213, 205, 196, 181, 179, 135, 140 ] colors = [ colors, 'BLK1', 'BLK2', 'BLK3', 'BLK4', 'BLK5', 'BLK6', 'BLK7', 'BLK8'] rvalue = [ rvalue, 250, 230, 210, 190, 112, 110, 70, 0 ] gvalue = [ gvalue, 250, 230, 210, 190, 128, 110, 70, 0 ] bvalue = [ bvalue, 250, 230, 210, 190, 128, 110, 70, 0 ] colors = [ colors, 'GRN1', 'GRN2', 'GRN3', 'GRN4', 'GRN5', 'GRN6', 'GRN7', 'GRN8'] rvalue = [ rvalue, 250, 223, 173, 109, 53, 35, 0, 0 ] gvalue = [ gvalue, 253, 242, 221, 193, 156, 132, 97, 69 ] bvalue = [ bvalue, 202, 167, 142, 115, 83, 67, 52, 41 ] colors = [ colors, 'BLU1', 'BLU2', 'BLU3', 'BLU4', 'BLU5', 'BLU6', 'BLU7', 'BLU8'] rvalue = [ rvalue, 232, 202, 158, 99, 53, 33, 8, 8 ] gvalue = [ gvalue, 241, 222, 202, 168, 133, 113, 75, 48 ] bvalue = [ bvalue, 250, 240, 225, 211, 191, 181, 147, 107 ] colors = [ colors, 'ORG1', 'ORG2', 'ORG3', 'ORG4', 'ORG5', 'ORG6', 'ORG7', 'ORG8'] rvalue = [ rvalue, 254, 253, 253, 250, 231, 217, 159, 127 ] gvalue = [ gvalue, 236, 212, 174, 134, 92, 72, 51, 39 ] bvalue = [ bvalue, 217, 171, 107, 52, 12, 1, 3, 4 ] colors = [ colors, 'RED1', 'RED2', 'RED3', 'RED4', 'RED5', 'RED6', 'RED7', 'RED8'] rvalue = [ rvalue, 254, 252, 252, 248, 225, 203, 154, 103 ] gvalue = [ gvalue, 232, 194, 146, 97, 45, 24, 12, 0 ] bvalue = [ bvalue, 222, 171, 114, 68, 38, 29, 19, 13 ] colors = [ colors, 'PUR1', 'PUR2', 'PUR3', 'PUR4', 'PUR5', 'PUR6', 'PUR7', 'PUR8'] rvalue = [ rvalue, 244, 222, 188, 152, 119, 106, 80, 63 ] gvalue = [ gvalue, 242, 221, 189, 148, 108, 82, 32, 0 ] bvalue = [ bvalue, 248, 237, 220, 197, 177, 163, 139, 125 ] colors = [ colors, 'PBG1', 'PBG2', 'PBG3', 'PBG4', 'PBG5', 'PBG6', 'PBG7', 'PBG8'] rvalue = [ rvalue, 243, 213, 166, 94, 34, 3, 1, 1 ] gvalue = [ gvalue, 234, 212, 189, 164, 138, 129, 101, 70 ] bvalue = [ bvalue, 244, 232, 219, 204, 171, 139, 82, 54 ] colors = [ colors, 'YGB1', 'YGB2', 'YGB3', 'YGB4', 'YGB5', 'YGB6', 'YGB7', 'YGB8'] rvalue = [ rvalue, 244, 206, 127, 58, 30, 33, 32, 8 ] gvalue = [ gvalue, 250, 236, 205, 175, 125, 95, 48, 29 ] bvalue = [ bvalue, 193, 179, 186, 195, 182, 168, 137, 88 ] colors = [ colors, 'RYB1', 'RYB2', 'RYB3', 'RYB4', 'RYB5', 'RYB6', 'RYB7', 'RYB8'] rvalue = [ rvalue, 201, 245, 253, 251, 228, 193, 114, 59 ] gvalue = [ gvalue, 35, 121, 206, 253, 244, 228, 171, 85 ] bvalue = [ bvalue, 38, 72, 127, 197, 239, 239, 207, 164 ] colors = [ colors, 'TG1', 'TG2', 'TG3', 'TG4', 'TG5', 'TG6', 'TG7', 'TG8'] rvalue = [ rvalue, 84, 163, 197, 220, 105, 51, 13, 0 ] gvalue = [ gvalue, 48, 103, 141, 188, 188, 149, 113, 81 ] bvalue = [ bvalue, 5, 26, 60, 118, 177, 141, 105, 71 ] ENDELSE ENDELSE ; Don't load system colors if you are doing this without a window connection. IF haveConnection THEN BEGIN ; Add system color names for IDL version 5.6 and higher. IF Float(!Version.Release) GE 5.6 THEN BEGIN tlb = Widget_Base() sc = Widget_Info(tlb, /System_Colors) Widget_Control, tlb, /Destroy frame = sc.window_frame text = sc.window_text active = sc.active_border shadow = sc.shadow_3d highlight = sc.light_3d edge = sc.light_edge_3d selected = sc.highlight face = sc.face_3d colors = [colors, 'Frame', 'Text', 'Active', 'Shadow'] rvalue = [rvalue, frame[0], text[0], active[0], shadow[0]] gvalue = [gvalue, frame[1], text[1], active[1], shadow[1]] bvalue = [bvalue, frame[2], text[2], active[2], shadow[2]] colors = [colors, 'Highlight', 'Edge', 'Selected', 'Face'] rvalue = [rvalue, highlight[0], edge[0], selected[0], face[0]] gvalue = [gvalue, highlight[1], edge[1], selected[1], face[1]] bvalue = [bvalue, highlight[2], edge[2], selected[2], face[2]] ENDIF ENDIF ; Make sure the color parameter and the colors array are an uppercase strings. varInfo = Size(theColor) IF varInfo[varInfo[0] + 1] NE 7 THEN $ Message, 'The color name parameter must be a string.', /NoName theColor = StrUpCase(StrCompress(StrTrim(theColor,2), /Remove_All)) colors = StrUpCase(StrCompress(StrTrim(colors,2), /Remove_All)) ; Check synonyms of color names. FOR j=0, N_Elements(theColor)-1 DO BEGIN IF StrUpCase(theColor[j]) EQ 'GREY' THEN theColor[j] = 'GRAY' IF StrUpCase(theColor[j]) EQ 'LIGHTGREY' THEN theColor[j] = 'LIGHTGRAY' IF StrUpCase(theColor[j]) EQ 'MEDIUMGREY' THEN theColor[j] = 'MEDIUMGRAY' IF StrUpCase(theColor[j]) EQ 'SLATEGREY' THEN theColor[j] = 'SLATEGRAY' IF StrUpCase(theColor[j]) EQ 'DARKGREY' THEN theColor[j] = 'DARKGRAY' IF StrUpCase(theColor[j]) EQ 'AQUA' THEN theColor[j] = 'AQUAMARINE' IF StrUpCase(theColor[j]) EQ 'SKY' THEN theColor[j] = 'SKYBLUE' IF StrUpCase(theColor[j]) EQ 'NAVYBLUE' THEN theColor[j] = 'NAVY' IF StrUpCase(theColor[j]) EQ 'CORNFLOWER' THEN theColor[j] = 'CORNFLOWERBLUE' IF StrUpCase(theColor[j]) EQ 'BROWN' THEN theColor[j] = 'SIENNA' ENDFOR ; How many colors do we have? ncolors = N_Elements(colors) ; Get the decomposed state of the IDL session right now. IF N_Elements(decomposedState) EQ 0 THEN BEGIN IF Float(!Version.Release) GE 5.2 THEN BEGIN IF (!D.Name EQ 'X' OR !D.Name EQ 'WIN' OR !D.Name EQ 'MAC') THEN BEGIN Device, Get_Decomposed=decomposedState ENDIF ELSE decomposedState = 0 ENDIF ELSE decomposedState = 0 IF (Float(!Version.Release) GE 6.4) AND (!D.NAME EQ 'Z') THEN BEGIN Device, Get_Decomposed=decomposedState ENDIF ENDIF ELSE decomposedState = Keyword_Set(decomposedState) ; Get depth of visual display. IF (!D.Flags AND 256) NE 0 THEN Device, Get_Visual_Depth=theDepth ELSE theDepth = 8 IF (Float(!Version.Release) GE 6.4) AND (!D.NAME EQ 'Z') THEN Device, Get_Pixel_Depth=theDepth IF (theDepth EQ 8) OR (decomposedState EQ 0) THEN offset = !D.Table_Size - ncolors - 2 ELSE offset = 0 ; Did the user want to select a color name? If so, we set ; the color name and color index, unless the user provided ; them. In the case of a single positional parameter, we treat ; this as the color index number as long as it is not a string. cancelled = 0.0 IF Keyword_Set(selectcolor) THEN BEGIN CASE N_Params() OF 0: BEGIN theColor = PickColorName(Filename=filename, _Extra=extra, Cancel=cancelled, BREWER=brewer) IF cancelled THEN RETURN, !P.Color IF theDepth GT 8 AND (decomposedState EQ 1) THEN BEGIN colorIndex = !P.Color < (!D.Table_Size - 1) ENDIF ELSE BEGIN colorIndex = Where(StrUpCase(colors) EQ StrUpCase(StrCompress(theColor, /Remove_All)), count) + offset colorIndex = colorIndex[0] IF count EQ 0 THEN Message, 'Cannot find color: ' + StrUpCase(theColor), /NoName ENDELSE END 1: BEGIN IF Size(theColor, /TName) NE 'STRING' THEN BEGIN colorIndex = theColor theColor = brewer ? 'WT1' : 'White' ENDIF ELSE colorIndex = !P.Color < 255 theColor = PickColorName(theColor, Filename=filename, _Extra=extra, Cancel=cancelled, BREWER=brewer) IF cancelled THEN RETURN, !P.Color END 2: BEGIN theColor = PickColorName(theColor, Filename=filename, _Extra=extra, Cancel=cancelled, BREWER=brewer) IF cancelled THEN RETURN, !P.Color END ENDCASE ENDIF ; Make sure you have a color name and color index. CASE N_Elements(theColor) OF 0: BEGIN theColor = brewer ? 'WT1' : 'White' IF N_Elements(colorIndex) EQ 0 THEN BEGIN IF theDepth GT 8 THEN BEGIN colorIndex = !P.Color < (!D.Table_Size - 1) ENDIF ELSE BEGIN colorIndex = Where(colors EQ theColor, count) + offset colorIndex = colorIndex[0] IF count EQ 0 THEN Message, 'Cannot find color: ' + theColor, /NoName ENDELSE ENDIF ELSE colorIndex = 0 > colorIndex < (!D.Table_Size - 1) ENDCASE 1: BEGIN type = Size(theColor, /TNAME) IF type NE 'STRING' THEN Message, 'The color must be expressed as a color name.' theColor = theColor[0] ; Make it a scalar or you run into a WHERE function "feature". :-( IF N_Elements(colorIndex) EQ 0 THEN BEGIN IF (theDepth GT 8) AND (decomposedState EQ 1) THEN BEGIN colorIndex = !P.Color < (!D.Table_Size - 1) ENDIF ELSE BEGIN colorIndex = Where(colors EQ theColor, count) + offset colorIndex = colorIndex[0] IF count EQ 0 THEN Message, 'Cannot find color: ' + theColor, /NoName ENDELSE ENDIF ELSE colorIndex = 0 > colorIndex < (!D.Table_Size - 1) ENDCASE ELSE: BEGIN type = Size(theColor, /TNAME) IF type NE 'STRING' THEN Message, 'The colors must be expressed as color names.' ncolors = N_Elements(theColor) CASE N_Elements(colorIndex) OF 0: colorIndex = Indgen(ncolors) + (!D.Table_Size - (ncolors + 1)) 1: colorIndex = Indgen(ncolors) + colorIndex ELSE: IF N_Elements(colorIndex) NE ncolors THEN $ Message, 'Index vector must be the same length as color name vector.' ENDCASE ; Did the user want color triples? IF Keyword_Set(triple) THEN BEGIN colors = LonArr(ncolors, 3) FOR j=0,ncolors-1 DO colors[j,*] = FSC_Color(theColor[j], colorIndex[j], Filename=filename, $ Decomposed=decomposedState, /Triple, BREWER=brewer) RETURN, colors ENDIF ELSE BEGIN colors = LonArr(ncolors) FOR j=0,ncolors-1 DO colors[j] = FSC_Color(theColor[j], colorIndex[j], Filename=filename, $ Decomposed=decomposedState, BREWER=brewer) RETURN, colors ENDELSE END ENDCASE ; Did the user ask for the color names? If so, return them now. IF Keyword_Set(names) THEN RETURN, Reform(colors, 1, ncolors) ; Process the color names. theNames = StrUpCase( StrCompress(colors, /Remove_All ) ) ; Find the asked-for color in the color names array. theIndex = Where(theNames EQ StrUpCase(StrCompress(theColor, /Remove_All)), foundIt) theIndex = theIndex[0] ; If the color can't be found, report it and continue with ; the first color in the color names array. IF foundIt EQ 0 THEN BEGIN Message, "Can't find color " + theColor + ". Substituting " + StrUpCase(colors[0]) + ".", /Informational theColor = theNames[0] theIndex = 0 ENDIF ; Get the color triple for this color. r = rvalue[theIndex] g = gvalue[theIndex] b = bvalue[theIndex] ; Did the user want a color triple? If so, return it now. IF Keyword_Set(triple) THEN BEGIN IF Keyword_Set(allcolors) THEN BEGIN IF Keyword_Set(row) THEN RETURN, Transpose([[rvalue], [gvalue], [bvalue]]) ELSE RETURN, [[rvalue], [gvalue], [bvalue]] ENDIF ELSE BEGIN IF Keyword_Set(row) THEN RETURN, [r, g, b] ELSE RETURN, [[r], [g], [b]] ENDELSE ENDIF ; Otherwise, we are going to return either an index ; number where the color has been loaded, or a 24-bit ; value that can be decomposed into the proper color. IF decomposedState THEN BEGIN ; Need a color structure? IF Arg_Present(colorStructure) THEN BEGIN theColors = FSC_Color_Color24([[rvalue], [gvalue], [bvalue]]) colorStructure = Create_Struct(theNames[0], theColors[0]) FOR j=1, ncolors-1 DO colorStructure = Create_Struct(colorStructure, theNames[j], theColors[j]) ENDIF IF Keyword_Set(allcolors) THEN BEGIN RETURN, FSC_Color_Color24([[rvalue], [gvalue], [bvalue]]) ENDIF ELSE BEGIN RETURN, FSC_Color_Color24([r, g, b]) ENDELSE ENDIF ELSE BEGIN IF Keyword_Set(allcolors) THEN BEGIN ; Need a color structure? IF Arg_Present(colorStructue) THEN BEGIN allcolorIndex = !D.Table_Size - ncolors - 2 IF allcolorIndex LT 0 THEN $ Message, 'Number of colors exceeds available color table values. Returning.', /NoName IF (allcolorIndex + ncolors) GT 255 THEN $ Message, 'Number of colors exceeds available color table indices. Returning.', /NoName theColors = IndGen(ncolors) + allcolorIndex colorStructure = Create_Struct(theNames[0], theColors[0]) FOR j=1, ncolors-1 DO colorStructure = Create_Struct(colorStructure, theNames[j], theColors[j]) ENDIF IF N_Elements(colorIndex) EQ 0 THEN colorIndex = !D.Table_Size - ncolors - 2 IF colorIndex LT 0 THEN $ Message, 'Number of colors exceeds available color table values. Returning.', /NoName IF (colorIndex + ncolors) GT 255 THEN BEGIN colorIndex = !D.Table_Size - ncolors - 2 ENDIF IF !D.Name NE 'PRINTER' THEN TVLCT, rvalue, gvalue, bvalue, colorIndex RETURN, IndGen(ncolors) + colorIndex ENDIF ELSE BEGIN ; Need a color structure? IF Arg_Present(colorStructure) THEN BEGIN colorStructure = Create_Struct(theColor, colorIndex) ENDIF IF !D.Name NE 'PRINTER' THEN TVLCT, rvalue[theIndex], gvalue[theIndex], bvalue[theIndex], colorIndex RETURN, colorIndex ENDELSE ENDELSE END ;------------------------------------------------------------------------------------------------------- PRO redraw restore, 'fftstruct.xdr' WSet, fftstruct.imageIndex m=239-fftstruct.nshadow pspec=fftstruct.fft_image workarray=pspec[0:319,0:m] workarray=shift(workarray,160,m/2) workarray=alog(workarray) ;log pspec for display widget_control, fftstruct.fieldIndex, set_value=' ' WSet, fftstruct.imageIndex Widget_Control, fftstruct.draw1, Draw_Motion_Events=0 fftstruct.imagetype=2 pspec[0:319,0:m]=workarray pspec=rebin(pspec,640,480) ;change color filter loadct, 13 tvscl, pspec, 0 if fftstruct.setshadow EQ 2 then begin sx = fftstruct.inbox[0] sy = fftstruct.inbox[1] dx = fftstruct.inbox[2] dy = fftstruct.inbox[3] PlotS, [sx, sx, dx, dx, sx], [sy, dy, dy, sy, sy], /Device, $ Color=FSC_COlor('Orange'), thick=2 endif if fftstruct.nbox EQ 0 then return FOR i=0, fftstruct.nbox-1 DO BEGIN sx = fftstruct.boxes[0,i] sy = fftstruct.boxes[1,i] dx = fftstruct.boxes[2,i] dy = fftstruct.boxes[3,i] PlotS, [sx, sx, dx, dx, sx], [sy, dy, dy, sy, sy], /Device, $ Color=FSC_COlor('Red'), thick=2 ENDFOR END PRO makemask restore, 'fftstruct.xdr' workarray=fltarr(320,240) workarray=shift(workarray, -160, -120) workarray=workarray+1. if fftstruct.nshadow GT 0 then workarray[0:319,239-fftstruct.nshadow+1:239]=0. if fftstruct.nbox GT 0 then begin FOR i=0, fftstruct.nbox-1 DO BEGIN sx = fftstruct.boxes[0,i]/2 sy = fftstruct.boxes[1,i]/2 dx = fftstruct.boxes[2,i]/2 dy = fftstruct.boxes[3,i]/2 workarray[sx:dx,sy:dy]=0. ENDFOR endif fftstruct.datamask=shift(workarray,160,120) save, fftstruct, filename = 'fftstruct.xdr' END ; handle cursor events in the image window here PRO imageevents, event ; parts of this are copied from David Fanning's "zimage.pro" routine... ; What type of an event is this? possibleEventTypes = [ 'DOWN', 'UP', 'MOTION', 'SCROLL', 'EXPOSE', 'KEY' ] thisEvent = possibleEventTypes(event.type) restore, 'fftstruct.xdr' if fftstruct.imagetype EQ 2 then begin ;-------------------------------------------------------------------------- ; the following is an adaptation of "drawbox_widget.pro" by David ; Fanning (see http://www.dfanning.com/widget_tips/rubberband_widget.html) ; with the small changes for my fftstruct being used rather than his ; info structure that is passed via the callback.... ; ; This is the event handler for the draw widget graphics window. ; Normally deal only with DOWN, UP, and MOTION events. ; Deal with keyboard events if trying ; to generate a level cut in the power spectrum window IF event.type GT 2 and event.type NE 5 THEN RETURN CASE thisEvent OF 'DOWN': BEGIN ; Turn motion events on for the image display widget. Widget_Control, event.id, Draw_Motion_Events=1 ; Copy window contents t9o pixmap. WSet, fftstruct.pixIndex Device, Copy=[0, 0, 640, 480, 0, 0, fftstruct.imageIndex] ; Get and store the static corner of the box. fftstruct.inbox[0] = event.x fftstruct.inbox[1] = event.y save, fftstruct, filename = 'fftstruct.xdr' ENDCASE 'UP': BEGIN ; Erase the last box drawn. Destroy the pixmap. WSet, fftstruct.imageIndex Device, Copy=[0, 0, 640, 480, 0, 0, fftstruct.pixIndex] ; WDelete, fftstruct.pixIndex ; Turn draw motion events off. Clear any events queued for widget. Widget_Control, event.id, Draw_Motion_Events=0, Clear_Events=1 ; Order the box coordinates. sx = Min([fftstruct.inbox[0], event.x], Max=dx) sy = Min([fftstruct.inbox[1], event.y], Max=dy) ; Here is where you do something useful with the box. ; For example purposes, I'll compute the average pixel ; value of the pixels enclosed by the box and print it. fftstruct.inbox[0]=sx fftstruct.inbox[1]=sy fftstruct.inbox[2]=dx fftstruct.inbox[3]=dy if fftstruct.setshadow LT 2 then begin ; now put the box information in the stack if there is room ; ; make sure we have room in the array (in practical terms, the limit ; is too high to be exceeded, but if needed drop the last "redo" box ; in the list to make room for the new box ; nbox=fftstruct.nbox mbox=fftstruct.mbox if nbox LT 998 then begin if mbox GT nbox then begin if mbox LT 1000 then begin fftstruct.boxes[0:3,nbox+2:mbox+2]=fftstruct.boxes[0:3,nbox:mbox] endif else begin fftstruct.boxes[0:3,nbox+2:999]=fftstruct.boxes[0:3,nbox:997] endelse endif fftstruct.boxes[0:3,nbox:nbox]=fftstruct.inbox nbox=nbox+1 newbox=fftstruct.inbox newbox[2]=639-fftstruct.inbox[0] newbox[0]=639-fftstruct.inbox[2] newbox[3]=479-fftstruct.inbox[1]-2*fftstruct.nshadow newbox[1]=479-fftstruct.inbox[3]-2*fftstruct.nshadow fftstruct.boxes[0:3,nbox:nbox]=newbox nbox=nbox+1 fftstruct.nbox=nbox fftstruct.mbox=mbox+2 endif else begin print, 'Out of room for boxes; save mask and restart to use more....' endelse save, fftstruct, filename = 'fftstruct.xdr' redraw makemask endif else begin ; turn on keyboard events so the user can mark a level to cut at Widget_Control, event.id, Draw_Keyboard_Events=1, Clear_Events=1 save, fftstruct, filename = 'fftstruct.xdr' redraw endelse ENDCASE 'MOTION': BEGIN ; Here is where the actual box is drawn and erased. ; First, erase the last box. WSet, fftstruct.imageIndex Device, Copy=[0, 0, 640, 480, 0, 0, fftstruct.pixIndex] ; Get the coodinates of the new box and draw it. sx = fftstruct.inbox[0] sy = fftstruct.inbox[1] dx = event.x dy = event.y PlotS, [sx, sx, dx, dx, sx], [sy, dy, dy, sy, sy], /Device, $ Color=FSC_Color('White'), thick=2 ENDCASE 'KEY': BEGIN inbox=fftstruct.inbox/2 ; check that the position is inside the defined box... if event.x LT fftstruct.inbox[0] OR event.x GT fftstruct.inbox[2] OR $ event.y LT fftstruct.inbox[1] OR event.y GT fftstruct.inbox[3] then return workarray=shift(fftstruct.fft_image,160,120) subarray=workarray[inbox[0]:inbox[2],inbox[1]:inbox[3]] value=workarray[event.x/2,event.y/2] vector=where(subarray GT value, count) maskarray=shift(fftstruct.datamask,160,120) if count NE 0 then begin subarray=subarray*0.+1. subarray[vector]=0. maskarray[inbox[0]:inbox[2],inbox[1]:inbox[3]]=subarray endif ; do the reflected box newbox=fftstruct.inbox newbox[2]=639-fftstruct.inbox[0] newbox[0]=639-fftstruct.inbox[2] newbox[3]=479-fftstruct.inbox[1]-2*fftstruct.nshadow newbox[1]=479-fftstruct.inbox[3]-2*fftstruct.nshadow newbox=newbox/2 workarray=shift(fftstruct.fft_image,160,120) subarray=workarray[newbox[0]:newbox[2],newbox[1]:newbox[3]] value=workarray[(639-event.x)/2,(479-event.y)/2] dims=size(subarray) vector=where(subarray GT value, count) if count NE 0 then begin subarray=subarray*0.+1. subarray[vector]=0. maskarray[newbox[0]:newbox[2],newbox[1]:newbox[3]]=subarray endif maskarray=shift(maskarray,-160,-120) fftstruct.datamask=maskarray save, fftstruct, filename = 'fftstruct.xdr' print, 'You must save the mask to retain these values; if you "Undo" this will be lost.' fftstruct.setshadow=0 Widget_Control, event.id, Draw_Keyboard_Events=0, Clear_Events=1 ENDCASE ENDCASE ;-------------------------------------------------------------------------- endif else begin str=string(event.x/2,event.y/2,format='("x =",I4," y =",I4)') widget_control, fftstruct.fieldIndex, set_value=str WSet, fftstruct.drawIndex if fftstruct.lineplot EQ 1 then begin if fftstruct.cutplot EQ 1 then begin if fftstruct.imagetype EQ 1 then begin plot,fftstruct.data[0:319,event.y/2], /xst, title = 'Raw Spectrum', ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') endif if fftstruct.imagetype EQ 3 then begin plot,fftstruct.filtered_data[0:319,event.y/2], /xst, title = 'Filtered Spectrum', ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') endif if fftstruct.imagetype EQ 4 then begin diffim=fftstruct.data[0:319,event.y/2]-fftstruct.filtered_data[0:319,event.y/2] plot,diffim, title = 'Difference Image', /xst, ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') endif endif else begin if fftstruct.imagetype EQ 1 then begin plot,fftstruct.data[event.x/2,0:239], /xst, title = 'Raw Spectrum', ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') endif if fftstruct.imagetype EQ 3 then begin plot,fftstruct.filtered_data[event.x/2,0:239], /xst, title = 'Filtered Spectrum', ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') endif if fftstruct.imagetype EQ 4 then begin diffim=fftstruct.data[event.x/2,0:239]-fftstruct.filtered_image[event.x/2,0:239] plot,fftstruct.diffim, title = 'Difference Image', /xst, ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') endif endelse endif if ((event.type EQ 0) AND (fftstruct.setshadow EQ 1)) then begin y=239-event.y/2 fftstruct.setshadow=0 fftstruct.nshadow=y ; the following is to avoid doing anything when nshadow is zero if y GT 0 then begin values=shift(fftstruct.datamask, -160, -120) values[0:319,239-y+1:239]=0. fftstruct.datamask=shift(values, 160, 120) endif save, fftstruct, filename = 'fftstruct.xdr' endif endelse END ; handle cursor events in the plot window here PRO drawevents, event ; no events here yet.... value=0. END PRO Options_Events, event Widget_Control, event.top, Get_UValue=info, /No_Copy Widget_Control, event.id, Get_Value=thisButtonValue CASE thisButtonValue OF '"Select Level" Function': BEGIN restore, 'fftstruct.xdr' fftstruct.setshadow=2 save, fftstruct, filename = 'fftstruct.xdr' ENDCASE 'Set Shadowed Region': BEGIN restore, 'fftstruct.xdr' fftstruct.setshadow=1 save, fftstruct, filename = 'fftstruct.xdr' ENDCASE 'Clear Shadowed Region': BEGIN restore, 'fftstruct.xdr' values=shift(fftstruct.datamask, -160, -120) values[0:319,0:239-fftstruct.nshadow]=1. fftstruct.datamask=shift(values, 160, 120) fftstruct.nshadow=0 save, fftstruct, filename = 'fftstruct.xdr' ENDCASE 'Line Cuts On/OFf': BEGIN restore, 'fftstruct.xdr' fftstruct.lineplot=1-fftstruct.lineplot save, fftstruct, filename = 'fftstruct.xdr' ENDCASE 'Toggle Line/Column Cuts': BEGIN restore, 'fftstruct.xdr' fftstruct.cutplot=1-fftstruct.cutplot save, fftstruct, filename = 'fftstruct.xdr' ENDCASE ENDCASE END ;Menu Bar events handled here PRO ffits2_Events, event Widget_Control, event.top, Get_UValue=info, /No_Copy Widget_Control, event.id, Get_Value=thisButtonValue CASE thisButtonValue OF 'Load Data': BEGIN ;Select file restore, 'fftstruct.xdr' path=fftstruct.lastpath file =Dialog_Pickfile(FILTER='*.fits',PATH=path, GET_PATH=pathname) if (file ne '') then begin fftstruct.lastpath=pathname ;Check to see if it is a MICHELLE file header = headfits(file) inst = SXPAR(header, 'INSTRUME', count=count) IF count EQ 1 THEN data =miche_read(file, /silent) $ ELSE data =readfits(file) ; have to check the image dimensions here and make sure they are OK ; (320 by 240) values=Size(data) if values[0] NE 2 OR values[1] NE 320 OR values[2] NE 240 then begin print, 'The image is the wrong size. It will not be laoded' return endif fftstruct.data=data fftstruct.header=header fftstruct.filtered_data=data fftstruct.fft_image=data*0. fftstruct.datamask=fftstruct.fft_image+1. ;Display array data = rebin(data, 640, 480) WSet, fftstruct.imageIndex fftstruct.imagetype=1 Widget_Control, fftstruct.draw1, Draw_Motion_Events=1 loadct, 13 tvscl, data, 0 ;Save temp file save, fftstruct, filename = 'fftstruct.xdr' endif ENDCASE 'Load Mask': BEGIN ;Select file restore, 'fftstruct.xdr' file =Dialog_Pickfile(FILTER='*.fits') if (file NE '') then begin ;Check to see if it is a MICHELLE file header = headfits(file) inst = SXPAR(header, 'INSTRUME', count=count) IF count EQ 1 THEN data =miche_read(file, /silent) $ ELSE mask =readfits(file) restore, 'fftstruct.xdr' fftsrruct.mask=mask ;Display array mask = rebin(mask, 640, 480) widget_control, fftstruct.fieldIndex, set_value=' ' WSet, fftstruct.imageIndex Widget_Control, fftstruct.draw1, Draw_Motion_Events=0 fftstruct.imagetype=2 loadct, 13 tvscl, shift(mask, -320, -240), 0 fftstruct.mbox=0 save, fftsruct, filename = 'fftstruct.xdr' endif ENDCASE 'Store Data': BEGIN restore, 'fftstruct.xdr' path=fftstruct.lastpath file =Dialog_Pickfile(FILTER='*.fits',PATH=path, GET_PATH=pathname) if (file NE '') then begin fftstruct.lastpath=pathname fits_write, file, fftstruct.filtered_data, fftstruct.header endif ENDCASE 'Store Mask': BEGIN restore, 'fftstruct.xdr' file =Dialog_Pickfile(FILTER='*.fits') if (file NE '') then begin fits_write, file, fftstruct.datamask ; we need a "miche_wite" routine here endif ENDCASE 'Quit': BEGIN restore, 'fftstruct.xdr' WDelete, fftstruct.pixIndex file_delete, 'fftstruct.xdr', /allow Widget_Control, event.top, /Destroy RETURN ENDCASE ENDCASE Widget_Control, event.top, /No_Copy END ;------------------------------------------------------------------ ;Button events handled here PRO ffits2_event, ev widget_control, ev.id, get_uvalue=uvalue ;define button operations CASE uvalue OF 'Show Data': BEGIN restore, 'fftstruct.xdr' WSet, fftstruct.imageIndex Widget_Control, fftstruct.draw1, Draw_Motion_Events=1 fftstruct.imagetype=1 draw=rebin(fftstruct.data, 640, 480) loadct, 13 tvscl, draw save, fftstruct, filename='fftstruct.xdr' END 'Show Mask': BEGIN ;Restore Mask restore, 'fftstruct.xdr' widget_control, fftstruct.fieldIndex, set_value=' ' WSet, fftstruct.imageIndex Widget_Control, fftstruct.draw1, Draw_Motion_Events=0 fftstruct.imagetype=5 datamask=rebin(fftstruct.datamask, 640, 480) loadct, 3 datamask=shift(datamask, -320, -240) ; the following is to display the constant mask in white not black, ; set one pixel to zero so there is a range of values if max(datamask) EQ min(datamask) then datamask[160,50]=0. tvscl, datamask save, fftstruct, filename='fftstruct.xdr' END 'Show Pwrspec': BEGIN restore,'fftstruct.xdr' psmin=min(fftstruct.fft_image) psmax=max(fftstruct.fft_image) m=239-fftstruct.nshadow if ((psmin eq 0) and (psmax eq 0)) then begin data=fftstruct.data[0:319,0:m] pspec = abs(fft(data)) fftstruct.fft_image=fftstruct.fft_image*0. fftstruct.fft_image[0:319,0:m]=pspec save, fftstruct, filename='fftstruct.xdr' endif pspec=fftstruct.fft_image workarray=pspec[0:319,0:m] workarray=shift(workarray,160,m/2) workarray=alog(workarray) ;log pspec for display widget_control, fftstruct.fieldIndex, set_value=' ' WSet, fftstruct.imageIndex Widget_Control, fftstruct.draw1, Draw_Motion_Events=0 fftstruct.imagetype=2 pspec[0:319,0:m]=workarray pspec=rebin(pspec,640,480) loadct, 13 tvscl, pspec, 0 WSet, fftstruct.pixIndex tvscl, pspec, 0 save, fftstruct, filename='fftstruct.xdr' redraw END 'FFT Filter': BEGIN restore,'fftstruct.xdr' m=239-fftstruct.nshadow data=fftstruct.data[0:319,0:m] mask=fftstruct.datamask[0:319,0:m] newdata = fftclean(data, mask) fftstruct.filtered_data[0:319,0:m]=newdata newdata=fftstruct.filtered_data newdata=rebin(newdata,640,480) WSet, fftstruct.imageIndex Widget_Control, fftstruct.draw1, Draw_Motion_Events=1 fftstruct.imagetype=3 save, fftstruct, filename='fftstruct.xdr' loadct, 13 tvscl, newdata END 'Show Diff': BEGIN restore,'fftstruct.xdr' data=fftstruct.filtered_data-fftstruct.data data=rebin(data,640,480) WSet, fftstruct.imageIndex Widget_Control, fftstruct.draw1, Draw_Motion_Events=1 fftstruct.imagetype=4 save, fftstruct, filename='fftstruct.xdr' loadct, 13 tvscl, data END 'Show Raw Plot': BEGIN restore,'fftstruct.xdr' data=fftstruct.data column = data[100,*] peak = where(column EQ max(column)) spect = miche_xtract(data, center = peak, width = 4) fftstruct.spectrum=spect save, fftstruct, filename='fftstruct.xdr' WSet, fftstruct.drawIndex plot, spect, /xst, title = 'Raw Spectrum', ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') oplot, spect, color=FSC_Color('Green') END 'Show Filtered Plot': BEGIN restore,'fftstruct.xdr' data=fftstruct.filtered_data column = data[100,*] peak = where(column EQ max(column)) spect = miche_xtract(data, center = peak, width = 4) fftstruct.spectrum=spect save, fftstruct, filename='fftstruct.xdr' WSet, fftstruct.drawIndex plot, spect, /xst, title = 'Filtered Spectrum', ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') oplot, spect, color=FSC_Color('Blue') END 'Show Both Plots': BEGIN restore,'fftstruct.xdr' data=fftstruct.filtered_data column = data[100,*] peak = where(column EQ max(column)) spect = miche_xtract(data, center = peak, width = 4) fftstruct.spectrum=spect save, fftstruct, filename='fftstruct.xdr' rawspect = miche_xtract(fftstruct.data, center = peak, width = 4) WSet, fftstruct.drawIndex plot, spect, /xst, title = 'Both Spectra', ytitle = 'Counts', xtitle='Pixel Number', Color=FSC_Color('White') oplot, spect, color=FSC_Color('Blue') oplot, rawspect, color=FSC_Color('Green') END ; the following are holders for when we have boxes on the plots 'Undo Box': BEGIN restore,'fftstruct.xdr' fftstruct.nbox=fftstruct.nbox-2 if fftstruct.nbox LT 0 then fftstruct.nbox=0 save, fftstruct, filename='fftstruct.xdr' redraw makemask END 'Redo Box': BEGIN restore,'fftstruct.xdr' fftstruct.nbox=fftstruct.nbox+2 if fftstruct.nbox GT fftstruct.mbox then fftstruct.nbox=fftstruct.mbox save, fftstruct, filename='fftstruct.xdr' redraw makemask END ENDCASE END function miche_xtract, array, mode = mode, center = center, width = width, file = file, $ verbose=verbose ;+ ; NAME: ; MICHE_ ; ; ; PURPOSE: ; ; ; ; CALLING SEQUENCE: ; ; ; ; INPUTS: ; ; ; ; OPTIONAL INPUT: ; ; ; OUTPUT: ; ; ; ; REVISION HISTORY: ; Written R.S. Fisher, Gemini Observatory, 2005 May ;- ;Print usage if no keywords given by user if n_params() eq 0 then begin print,'Usage: frame = miche_medxtract(array, [mode, center, width])' print, ' ' print,' array -- name of array (or use: file = "filename" ' print,' mode -- mode of extraction (median, sum, mean) ' print,' center -- central pixel of extraction ' print,' width -- width of extraction to collapse (default = 100 pix) ' endif ;Check to see if it is a file or an array if keyword_set(file) then array = miche_read(file, /silent) ;Set default values for center and width if not keyword_set(center) then center = 120 if not keyword_set(width) then width = 100 if not keyword_set(mode) then mode = 'sum' ;define output vector sz = size(array) output = fltarr(sz(1)) ;Set up loop to extract spectrum for i = 0, sz(1) - 1 do begin CASE mode OF 'median': output[i] = median(array[i, (center - (width/2.0)):(center + (width/2.0))]) 'mean': output[i] = mean(array[i, (center - (width/2.0)):(center + (width/2.0))]) 'sum': output[i] = total(array[i, (center - (width/2.0)):(center + (width/2.0))]) ENDCASE endfor print, ' ' if keyword_set(verbose) then print, 'center= ',strtrim(center,2),' width= ',strtrim(width,2), $ ' mode= ',mode ;Return requested array return, output end function miche_read, filename, frametype = frametype, ext = ext, silent=silent ;+ ; NAME: ; MICHE_READ ; ; ; PURPOSE: ; This function can be used to read a single frame from a Michelle Gemini MEF file. ; ; ; CALLING SEQUENCE: ; frame = miche_read('filename', frametype, [extension = extension]) ; ; ; INPUTS: ; frame - name of returned frame ; filename - full name of file to be read ; frametype - type of frame you want returned. ; 'A' - returns chopebeam A frame ; 'B' - returns chopbeam B frame ; 'DIF' - returns the differenced frame (A-B) from this extension ; ; ; OPTIONAL INPUT: ; ext - numerical value of the extension to extract data from. Default = 1 ; ; ; OUTPUT: ; frame - 320x240 Michelle frame defined by frametype and extension ; ; ; REVISION HISTORY: ; Written R.S. Fisher, Gemini Observatory, 2003 Jan ;- ;Print usage if no keywords given by user if n_params() eq 0 then begin print,'Usage: frame = miche_read(filename, frametype=frametype, extension=extension)' print,' filename -- name of input file ' print,' frametype -- type of frame to return, "DIF", "A", "B", default = "DIF" ' print,' ext -- extension of the file to access, default = 1 ' print,' silent -- keyword to supress diagnostic information display ' print,' polar -- keyword to indicate you are reading an IRAF polarization file ' endif ; Check to see if extension keyword is set if not keyword_set(ext) then ext = 1 if not keyword_set(frametype) then frametype = 'DIF' extension = ext ; Read in Michelle file using MRDFITS data = mrdfits(filename, extension, /SILENT) ; Read in header and print information if keyword_set(silent) then silent = 1 else silent = 0 ; Choose which frame to return CASE frametype OF 'DIF': frame = data[*,*,0] 'A': frame = data[*,*,1] 'B': frame = data[*,*,2] ENDCASE if silent eq 0 then begin ;Get header information from original file tempo = mrdfits(filename, 0, header) ;Old line is below here, crashed in 2007 Jan ;header = HEADFITS(filename) print, ' ' print,'******************************' print,' File Information print,'******************************' print, 'Extension: ', strtrim(ext,2) print, 'Frametype: ', frametype print, 'Filter: ', sxpar(header, 'FILTER') print, 'Frametime: ', sxpar(header, 'EXPOSURE') print, ' ' endif ;Return requested array return, frame end ;+ ; NAME: ; ; ; ; PURPOSE: ; ; ; CALLING SEQUENCE: ; ; ; INPUTS: ; ; OPTIONAL INPUT: ; ; ; OUTPUT: ; frame - frame cleaned by "FFT" ; ; ; REVISION HISTORY: ; Written R.S. Fisher, Gemini Observatory, 2006 March ;- function fftclean, data, mask ;Print usage if no keywords given by user if n_params() eq 0 then begin print,'Usage: frame = fftclean(array, mask) ' endif datafft = fft(data) data_real = real_part(datafft) data_imag = imaginary(datafft) ;Apply mask data_real = data_real * mask data_imag = data_imag * mask datanew = complex(data_real, data_imag) datanew = fft(datanew, /inverse) datanew = real_part(datanew) ;Return requested array return, datanew end ;------------------------------------------------------------------ ;Widget interface created here PRO mdefringe DEVICE, DECOMPOSED = 0 ;create main base main = widget_base (title='FFT Spectrum Analysis', /column, /align_center, MBar=menubase) ;create top display window with a header layer1 = widget_base (main, xsize=643, ysize= 485, /column, /align_center, frame=3) label = widget_label (layer1, value='IMAGE DISPLAY AREA', ysize=20) draw1 = widget_draw (layer1, xsize= 640, ysize= 480, uname = 'draw1', /button_events, Event_Pro='imageevents') ; the folowing is a pixmap, which we will use for rapid ; drawing/redrawing of the power spectrum image.... Window, /Free, XSize=640, YSize=480, /Pixmap pixIndex = !D.Window ;Insert spacer above the row of buttons label = widget_label (main, value=' ', ysize=10) ;create menu bar buttons. fileID = Widget_Button(menubase, Value='File', Event_Pro='ffits2_Events') openID = Widget_Button(fileID, Value='Load Data') maskID = Widget_button(fileID, Value='Load Mask') storeID = Widget_Button(fileID, Value='Store Data') storemaskID = Widget_button(fileID, Value='Store Mask') quitID = Widget_Button(fileID, Value='Quit') optionID=Widget_Button(menubase, Value='Option', Event_Pro='Options_Events') setShadowID=Widget_Button(optionID,Value='Set Shadowed Region') resetShadowID=Widget_Button(optionID,Value='Clear Shadowed Region') toggleLineCutID=Widget_Button(optionID,Value='Line Cuts On/OFF') toggleCutID=Widget_Button(optionID,Value='Toggle Line/Column Cuts') cutPspectID=Widget_Button(optionID,Value='"Select Level" Function') ;Create spectrum analysis buttons layer2 = widget_base (main, xsize=643, ysize= 40, /row, /align_center) button = widget_button (layer2, ysize=20, xsize=100, uvalue='Show Data', value='Show Data') button = widget_button (layer2, ysize=20, xsize=100, uvalue='Show Pwrspec', value='Show Pwrspec') button = widget_button (layer2, ysize=20, xsize=100, uvalue='Show Mask', value='Show Mask') button = widget_button (layer2, ysize=20, xsize=100, uvalue='Show Diff', value='Show Diff') button = widget_button (layer2, ysize=20, xsize=100, uvalue='FFT Filter', value='FFT Filter') ;Create 'plot' analysis buttons layer3 = widget_base (main, xsize=643, ysize= 40, /row, /align_center) button = widget_button (layer3, ysize=20, xsize=100, uvalue='Show Raw Plot', value='Old Spectrum') button = widget_button (layer3, ysize=20, xsize=100, uvalue='Show Filtered Plot', value='New Spectrum') button = widget_button (layer3, ysize=20, xsize=100, uvalue='Show Both Plots', value='Both Spectra') button = widget_button (layer3, ysize=20, xsize=100, uvalue='Undo Box', value='Undo Box') button = widget_button (layer3, ysize=20, xsize=100, uvalue='Redo Box', value='Redo Box') labelfield = widget_label(layer3, ysize=20, xsize=100, value='', uname='positionfield') ;Insert spacer below the row of buttons label = widget_label (main, value=' ', ysize=10) ;create bottom display window with a header layer3 = widget_base (main, xsize=643, ysize= 485, /column, /align_center, frame=3) label = widget_label (layer3, value='SPECTRUM DISPLAY AREA', ysize=20) draw2 = widget_draw (layer3, xsize= 640, ysize= 480, uname = 'draw2', Event_Pro='plotevents') widget_control, main, /realize xmanager, 'ffits2', main, /no_block Widget_Control, draw1, Get_Value=imageIndex Widget_Control, draw2, Get_Value=drawIndex ; drawIndex and imageIndex hold the window indexes for the plot and ; spectrum regions respectively. ; make structure for the images and such to be used in this procedure data=fltarr(320,240) spectrum=fltarr(320) boxes=intarr(4,1000) inbox=intarr(4) datamask=data ; this value is set interactively via the cursor (can also be set from ; the header with some work nshadow=0 header=strarr(57600) loadct, 13 fftstruct = { data:data, $ ; the input image header:header, $ ; the original header filtered_data:data, $ ; the output image, after FFT filtering datamask:datamask, $ ; the mask image fft_image:data, $ ; the frequency domain image setshadow:0, $ ; flag for interactively setting nshadow (0=off, 1=on) nshadow:nshadow, $ ; integer for number of rows have no signal boxes:boxes, $ ; array to hold boxes to mask in the frequency domain nbox:0, $ ; integer pointer for box array mbox:0, $ ; integer pointer for box array "undo/redo" function, holds the maximum nbox value that has been used draw1:draw1, $ ; the full ID for the image window (needed for motion event on/off) imageIndex:imageIndex, $ ; the image window ID drawIndex:drawIndex, $ ; the draw window ID fieldIndex:labelfield, $ ; position label field window ID imagetype:0, $ ; flag for image type (0 for undefined, 1 for image, 2 for FFT) lineplot:1, $ ; toggle flag for plotting line plots beased on cursor position in image window (0 = no, 1 = yes) cutplot:1, $ ; toggle flag for plotting line/column plots (1=line, 0=column) pixIndex:pixIndex, $ ; the window ID for the pixmap inbox:inbox, $ ; array to hold box corner coordinates ((x,y) for starting and ending points) spectrum:spectrum, $ ; stores the spectrum for plotting lastpath:'' $ ; file selection/store pathname } ; nshadow values: 0 for lowN/lowQ, 20 for medN1, 25 for medN2, 80 for echelle if (fftstruct.nshadow GT 0) then begin m=239-fftstruct.nshadow values=shift(fftstruct.datamask, -160, -120) values[0:319,0:m]=1. fftstruct.datamask=shift(values, 160, 120) endif ; Store the structure. All work will be done accessing this ; structure and replacing the save file when changes are made save, fftstruct, filename='fftstruct.xdr' END