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