C newrex.f C C Simple application to read New Horizons REX Data from PDS data sets C using PDS labels as maps to FITS data C C Sample C Usage: ls data/20070101_003031/rex_003031500?*.lbl \ C | ./newrex > rex.txt C C Input: List of filepaths, one per line, to PDS labels for REX data C from standard input (STDIN) C C Output: For each filepath in the input, write five, and then three*, C columns of data over 1250 rows: C Column 1: Row number C Columns 2 and 3: I and Q (IQ) values from FITS BINTABLE C Columns 4 and 5: Radiometry and Time Tags (RT) from FITS C BINTABLE C C * I and Q data are continuous, 1250 rows per PDS label, but C RT data have only 10 rows per PDS label C C Assumptions and limitations: C - PDS label filepaths in STDIN are 80 characters or less in length C - FITS filepaths mock PDS label filepaths, ending in 'fit' not 'lbl' C - PDS label statements of interest are single line, ASCII C - PDS label statements of interest are 80 chars or less in length C - Filenames in PDS label pointers are upper case (to look for _ENG_) C - All of the same types of tables have identical schema C - IQ Raw table: 2 INTEGER*2 columns; 1250 rows C - RT Raw table: 2 (INTEGER*8 and INTEGER*4) columns; 10 rows C - IQ Calibrated table: 2 REAL*4 columns; 1250 rows C - RT Calibrated table: 2 REAL*4 columns; 1250 rows C - Therefore, there is no need to parse TABLE OBJECTs in PDS label C - Only need to parse ordinal of 1st 2880-byte record of TABLE data C C Algorithm: C C Each PDS filepath, read from STDIN, is used to open a PDS label file, C the read and parse each label statment until the location of the data C in the corresponding FITS file has been determined. C C As the PDS label statement pointers for the IQ and Radimetry/Time C Tag tables are read and parsed, get the ordinals of the first C 2880-byte records for the corresponding data sections of the C Extension Data Units (EDUs) in the FITS file. C C The EDU pointed to by ^EXTENSION_IQVALS_TABLE is a FITS BINTABLE with C two columns, I and Q, and 1250 rows. C C The EDU pointed to by ^EXTENSION_RAD_TIME_TAGS_TABLE is a FITS C BINTABLE with two columns, Radiometry and Time Tags, and 10 rows. C C For the first 10 lines of output this code prints out the the row C number (1-10) and four values: two values from each of the first 10 C rows of the IQ FITS BINTABLE, and two values from each of the 10 rows C of the RT FITS BINTABLE. C C Then this code prints out the row number (11-1250) and the two values C from the last 1240 rows of the IQ FITS BINTABLE. C C The output looks like this for the files (raw shown; calibrated is C similar; column headers shown here are not in output): C C row# I Q Radiometry Time Tag C ==== ==== ==== ========== ======== C C 1 -414 2382 24560001298 6730 C 2 -2021 1445 2459374881 6731 C 3 -2331 -372 4919752187 6732 C ... C 8 1662 1795 17198708111 6737 C 9 -103 2443 19657874753 6738 C 10 -1836 1682 22111212526 6739 C 11 -2474 -121 C 12 -1670 -1753 C 13 84 -2405 C ... C 1248 -1092 -2341 C 1249 785 -2407 C 1250 2200 -1232 C 1 -414 2382 24560001298 6730* C 2 -2021 1445 2459374881 6731* C 3 -2331 -372 4919752187 6732* C ... C C * Those last three lines show the start of the next file C C The format for the columns: I5 for the row number; I16 or F16.4 C for each of the remaining raw or calibrated column, respectively C C The file is raw if the filename from the PDS label pointer has the C string '_ENG_' in it. C C The file is calibrated if the filename has the string '_SCI_' in it, C but code only tests for '_ENG_' and assumes calibrated if '_ENG_' is C absent. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC program newrex implicitnone C PDS label file variables C ======================== C lblFn, fitFn: paths to PDS label and FITS files. C line: PDS label statment C iEq, kwd: idex of '=' in line; PDS token before '=' in line character*80 lblFn character*80 fitFn character*80 line integer iEq character*80 kwd C FITS file variables C =================== C ioFit: LUN for FITS file; C iqrec0, rtrec0: ordinals of 1st 2880-byte records of IQ and RT C data in FITS file C isRaw: .TRUE. if data are raw; .FALSE. if calibrated integer ioFit/51/ integer iqrec0/0/ integer rtrec0/0/ logical isRaw C i,j: multi-purpose counters integer i,j C TABLE data variables C ==================== C iqRaw: IQ Raw data, 1250 rows C - length: 2 * 2 * 1250 = 5,000 bytes < 2 FITS records integer*2 iqRaw(2,1250) C radRaw, rrgRaw: Radiometry and Time Tags Raw, 10 rows C - length: (8 + 4) * 10 = 120 bytes < 1 FITS record integer*8 radRaw(10) integer*4 ttgRaw(10) C iqCal: IQ Calibrated data, 1250 rows C - length: 4 * 2 * 1250 = 10,000 bytes < 4 FITS records real*4 iqCal(2,1250) C rtcal: Radiometry and Time Tags Calibrated, 10 rows C - length: 4 * 2 * 10 = 80 bytes < 1 FITS record real*4 rtCal(2,10) C Buffer variables for FITS data C ============================== C 2880-byte FITS records will be read into local 2880-byte C character buffer variables, with names starting with c2880. C Binary data (IQ, Radiometry, Time Tags) will overlay the buffers C via EQUIVALENCE statements. C FITS records: 4 * 2880 = 11,520 bytes, which is greater than C largest data quantity (4-byte IQ data, 2 columns, C 1250 rows = 10,000 bytes) character*2880 c2880iq(4) C Special case for Raw RT data: columns have different byte widths C - use two c2880 buffers, one for Radiometry, one for Time Tags character*2880 c2880rt1, c2880rt2 C EQUIVALENCE Raw and Calibrated IQ array to IQ character buffer equivalence (iqRaw,c2880iq) equivalence (iqCal,c2880iq) C EQUIVALENCE Calibrated RT array to first RT character buffer equivalence (rtCal,c2880rt1) C Special equivalence for 8+4-byte Raw Radiometry + Time Tag rows C - use separate buffers for each column C - will duplicate single FITS record and shift bytes to align each C array separately equivalence (radRaw,c2880rt1) equivalence (ttgRaw,c2880rt2) C Utility function declarations (code is at the end of this file) C ================================================================= character*80 lblOpen,lblRead,lblClose external lblOpen,lblRead,lblClose integer getRecno external getRecno C START OF CODE C ============= dowhile ( .true. ) C Loop over STDIN lines, read one PDS label path per line read(*,'(a80)',end=99999) lblFn C Open label file, reset IQ and Radiometry/Time Tag record numbers C to indicate TABLE pointers are not yet found fitFn = lblOpen(lblFn) iqrec0 = 0 rtrec0 = 0 C Loop until both TABLE pointers are found dowhile ( iqrec0.lt.1 .or. rtrec0.lt.1 ) C Read next label line line = lblread() C Skip if no equals sign iEq = index(line,'=') if ( iEq .lt. 2 ) continue C Read the keyword before the equals sign read(line(1:iEq-1),'(a)') kwd C Check kwd for IQ and Radiometry/Time Tag TABLE data pointers C e.g. C ^EXTENSION_IQVALS_TABLE = ("REX_0030315000_0X7B3_SCI_1.FIT", 13) C ^EXTENSION_RAD_TIME_TAGS_TABLE = ("REX_0030315000_0X7B3_SCI_1.FIT", 18) if ( kwd .eq. '^EXTENSION_IQVALS_TABLE') then C This is the pointer to the FITS IQ TABLE data C - get the ordinal to the TABLE's first FITS 2880-byte record C - non-zero iqrec0 clears first dowhile condition above iqrec0 = getRecno(line(iEq+1:80)) C Set isRaw = .TRUE. if _ENG_ is in the filename in the C PDS label statement isRaw = index(line(iEq+1:80),'_ENG_') .gt. 0 elseif ( kwd .eq. '^EXTENSION_RAD_TIME_TAGS_TABLE') then C This is the pointer to the FITS Radiometry/Time Tag TABLE data C - get the ordinal to the TABLE's first FITS 2880-byte record C - non-zero rtrec0 clears second dowhile condition above rtrec0 = getRecno(line(iEq+1:80)) endif C If we get to PDS END statement or to end of file, then we did C not find both TABLE statements if ( line .eq. 'END' .or. line .eq. '***DONE***' ) then print*, 'Failed ...' STOP endif enddo C Close label file lblFn = lblClose() C Done with PDS label; now read the data from the FITS file C ========================================================= C Build FITS filepath: change PDS label .lbl suffix to .fit i = index(fitFn,' ') fitFn(i-3:i) = 'fit ' C Open FITS file for direct access with 2880-byte records open(ioFit,status='OLD',file=fitFn,access='DIRECT',recl=2880) if ( isRaw ) then C Read Raw data: C Raw IQ data: two FITS records read(ioFit,rec=iqrec0) c2880iq(1) read(ioFit,rec=iqrec0+1) c2880iq(2) do i=1,1250 do j=1,2 call swapi2(iqRaw(j,i)) enddo enddo C Raw RT data: 1 FITS record > 120 bytes C - duplicate same FITS record in two places read(ioFit,rec=rtrec0) c2880rt1 c2880rt2 = c2880rt1 C Data are aligned in groups of 12 bytes in both buffers, but C INTEGER*8 and INTEGER*4 values, each such pair being one row C of the FITS BINTABLE, are interleaved: C RRRRRRRRTTTTrrrrrrrrtttt... C - RRRRRRRR = 1st Radiometry INTEGER*8 value C - TTTT = 1st Time Tag INTEGER*4 value C - rrrrrrrr = 2nd Radiometry INTEGER*8 value C - tttt = 2nd Time Tag INTEGER*4 value C Shift Radiometry bytes to be contiguous in buffer c2880rt1, C and shift Time Tag values to be contiguous in buffer c2880rt2: C +--radRaw(1) C | C | +--radRaw(2) C | | C V V C <==+===><==+===> C RRRRRRRRrrrrrrrr... (c2880rt1; INTEGER*8 radRaw(10)) C C TTTTtttt... (c2880rt2; INTEGER*4 ttgRaw(10)) C <+=><+=> C ^ ^ C | | C | +--ttgRaw(2) C | C +--ttgRaw(1) do i=1,10 C Shift bytes in buffers: if (i .gt. 1) c2880rt1(i*8-7:i*8) = c2880rt1(i*12-11:i*12-4) c2880rt2(i*4-3:i*4) = c2880rt2(i*12-3:i*12) C Swap bytes call swapi8(radRaw(i)) call swapi4(ttgRaw(i)) enddo C Raw output: do i=1,1250 if ( i .lt. 11 ) then C Row number, IQ, Radiometry and Time Tags for first 10 rows print'(I5,4I16)', i, iqRaw(1,i), iqRaw(2,i) & , radRaw(i), ttgRaw(i) else C Row number and IQ only for remaining 1250 rows print'(I5,2I16)', i, iqRaw(1,i), iqRaw(2,i) endif enddo else C Read Calibrated data: C Calibrated IQ data: C - Read 4 FITS records into 4 buffers, > 10kbytes C - Swap bytes do i=1,4 read(ioFit,rec=iqrec0+i-1) c2880iq(i) enddo do i=1,1250 do j=1,2 call swapr4(iqCal(j,i)) enddo enddo C Calibrated RT data: C - Read 1 FITS record into 1 buffer > 80 bytes C - Swap bytes C - No special case; Radiometry and Time Tags are REAL*4 read(ioFit,rec=rtrec0) c2880rt1 do i=1,10 do j=1,2 call swapr4(rtCal(j,i)) enddo enddo C Calibrated output: do i=1,1250 if ( i .lt. 11 ) then C Row number, IQ, Radiometry and Time Tags for first 10 rows print'(I5,4F16.4)', i, iqCal(1,i), iqCal(2,i) & , rtCal(1,i), rtCal(2,i) else C Row number and IQ only for remaining 1250 rows print'(I5,2F16.4)', i, iqCal(1,i), iqCal(2,i) endif enddo endif close(ioFit) enddo 99999 continue end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Function that extracts ordinal of first 2880-byte record of a FITS C file EDU section from a PDS label pointer statement C C e.g. get the 12 from this PDS label pointer statement: C C ^EXTENSION_IQVALS_HEADER = ("REX_0030315000_0X7B3_ENG_1.FIT", 12) C integer function getRecno(line) implicitnone character*80 line integer iComma, iParen, iRecno character*80 str C getRecno = 0 str = line iComma = index(str,',') if ( iComma .eq. 0 ) return str = str(iComma+1:80) iParen = index(str,')') if ( iParen .eq. 0 ) return read( str(1:iParen-1), '(i10)') iRecno getRecno = iRecno return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC character*80 function lblOpen(lblFn) implicitnone character*80 lblFn integer i character*80 L character*80 lblRead, lblClose integer ioLbl/50/ character*80 lblFnSave / ' ' / integer lblFnLenSave save lblFnSave, lblFnLenSave, ioLbl C Close if already open: open if saved filename is not blank if (lblFnSave .ne. ' ') close(ioLbl) C Truncate at and after first space i = index(lblFn,' ') - 1 open(ioLbl,file=lblFn(1:i),status='OLD') C Save filename, return filename lblFnLenSave = i lblFnSave = lblFn(1:i) lblopen=lblFnSave return C Read one line ... entry lblRead read(ioLbl,'(A78)',end=99999) L lblread=L return C ... Return ***DONE*** if at end of file 99999 lblread='***DONE***' return C Close file entry lblClose if ( lblFnSave(1:1) .ne. ' ') then close(ioLbl) lblFnSave = ' ' endif lblClose=' ' return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Routines to set up EQUIVALENCE between 8-character array and binary C integers and reals subroutine swappers implicitnone integer*8 argi8 integer*4 argi4 integer*2 argi2 real*4 argr4 real*8 argr8 character*1 cequ(8) integer*8 i8 integer*4 i4 integer*2 i2 real*4 r4 real*8 r8 equivalence (cequ,i8) equivalence (cequ,i4) equivalence (cequ,i2) equivalence (cequ,r4) equivalence (cequ,r8) save return entry swapi8(argi8) i8 = argi8 call doswap(cequ,8) argi8 = i8 return entry swapi4(argi4) i4 = argi4 call doswap(cequ,4) argi4 = i4 return entry swapi2(argi2) i2 = argi2 call doswap(cequ,2) argi2 = i2 return entry swapr4(argr4) r4 = argr4 call doswap(cequ,4) argr4 = r4 return entry swapr8(argr8) r8 = argr8 call doswap(cequ,8) argr8 = r8 return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Subroutine to actually do the swap subroutine doswap(cin,n) implicitnone integer n character*1 cin(8) C character*1 c2(2) integer*2 i2 equivalence (i2,c2) C character*1 c1 integer i save C Do nothing if local machine is already MSByte-first i2 = 256 if ( ichar(c2(1)) .eq. 1 ) return C Swap bytes do i=1,n/2 c1 = cin(i) cin(i) = cin(n+1-i) cin(n+1-i) = c1 enddo return end