C============================================================================ C S U B R O U T I N E F I T S I M R E A D C---------------------------------------------------------------------------- C by Mike Eracleous (12/4/97) C Based on the simpler subroutine FITSREAD C---------------------------------------------------------------------------- C Subroutine to read a FITS file assuming that it consists of a number of C contiguous header records (containing all the header information) followed C by a number of contiguous data records (containing an image). C The entire data and header arrays as passed to the calling routine, but C at the same time a number of specific header keywords are interpreted C and their values are passed to the calling routine, for convenience. C C While the file is being read the subroutine displays a run-time summary of C the file structure for verification. The strategy is to make a first pass C over the entire file to verify that the structure is as assumed. Then a C second pass is made to read the header and the data separately and store C them in apropriate arrays. C C During the second pass, once the header is read some of its keywords are C interpreted and displayed. The number of pixels in the spectrum is probably C the most important header parameter, which is passed to the calling routine. C---------------------------------------------------------------------------- C This subroutine is specific to FITS files from the STScI digitized sky C survey, downloaded from their WWW tool. It expects the image pixels to be C described by 2-byte integers and it also expects to find certain information C in the FITS header. C---------------------------------------------------------------------------- C THIS ROUTINE MUST BE COMPILED WITH: f77 -c -O3 (-autopar) C ---------------------------------- C IF THE "-xl" OPTION IS USED THE COMPILER WILL OVERRIDE THE RECORD LENGTH C SPECIFIED IN FORTRAN "OPEN" STATEMENTS AND ASSIGN ITS OWN DEFAULT RECORD C LENGTH. THIS WILL MAKE A MESS OF EVERYTHING SINCE THE FITS FORMAT REQUIRES C RECORDS OF A FIXED LENGTH OF 2880 BYTES. C---------------------------------------------------------------------------- SUBROUTINE FITSCHARTREAD & (FNAME,VERBOSE, & HLINE,HEAD,KEYWORD,NHEAD, & FLUX,NROWS,NPIX,XSIZE,YSIZE, & XPIXELSZ,YPIXELSZ,PLTSCALE, & RASTRING,DECSTRING,EPOCH) PARAMETER (MAXREC=1000) REAL*4 FLUX(*) INTEGER*2 IARRAY(1440) INTEGER IHEAD(MAXREC),IDATA(MAXREC) CHARACTER*80 BUFF(36),HEAD(*),HLINE(*) CHARACTER*8 FIRSTKEY(MAXREC),RECTYPE(MAXREC),TYPE,KEYWORD(*) CHARACTER FNAME*100,DIMS*11,RASTRING*20,DECSTRING*20 LOGICAL VERBOSE C--> Open the file for direct access and with a fixed record length of C 2880 bytes. OPEN(UNIT=31,FILE=FNAME,STATUS='OLD', & RECL=2880,FORM='UNFORMATTED',ACCESS='DIRECT') IF (VERBOSE) THEN WRITE(*,510) FNAME 510 FORMAT(/,5X,45('-'), & /,5X,'VERIFICATION PASS: ',A30, & /,5X,45('-'), & /,5X,'Rec.#',4x,'Type ',4x,'1st Keywd.', & /,5X,45('-')) ENDIF C--> Make the first pass through the file and find the header and data C records. Report them on the screen for the user's reference and keep C and internal record of them as well for future use. Search the first C header record for the NAXIS keywords which contain the dimensions of C the image and use them to compute the number of image rows. IREC=1 100 READ(31,REC=IREC,END=999) BUFF FIRSTKEY(IREC)=BUFF(1)(1:8) IF (FIRSTKEY(IREC).EQ.'SIMPLE'.OR. & FIRSTKEY(IREC).EQ.'XTENSION') TYPE='HEADER' RECTYPE(IREC)=TYPE DO I=1,36 IF (BUFF(I).EQ.'END') THEN TYPE='DATA' NHEAD=I+(IREC-1)*36 ENDIF IF (IREC.EQ.1) THEN IF (BUFF(I)(1:6).EQ.'NAXIS ') THEN READ(BUFF(I)(10:80),*) NAXIS ENDIF IF (BUFF(I)(1:6).EQ.'NAXIS1') THEN READ(BUFF(I)(10:80),*) NPIX NAXIS1=NPIX ENDIF IF (BUFF(I)(1:6).EQ.'NAXIS2') THEN READ(BUFF(I)(10:80),*) NROWS NAXIS2=NROWS ENDIF ENDIF ENDDO C--> Compute the total number of image rows and the total number of pixels C that the data array must hold to accomodate the entire image. Also C write the exact image dimensions in a label and print it out later C for the user. IF (IREC.EQ.1) THEN IF (NAXIS.EQ.1) THEN WRITE(DIMS,610) NAXIS1 610 FORMAT(I5,' x 1') ELSEIF (NAXIS.EQ.2) THEN NDIG1=1+INT(LOG10(FLOAT(NAXIS1))) NDIG2=1+INT(LOG10(FLOAT(NAXIS2))) WRITE(DIMS,620) NAXIS1,NAXIS2 620 FORMAT(I,' x ',I) ENDIF ENDIF C--> Report that this header record was read IF (VERBOSE) THEN IF (RECTYPE(IREC).EQ.'DATA') THEN FIRSTKEY(IREC)='' IF (RECTYPE(IREC-1).EQ.'HEADER') THEN WRITE(*,520) IREC,RECTYPE(IREC),FIRSTKEY(IREC) ENDIF ELSE WRITE(*,520) IREC,RECTYPE(IREC),FIRSTKEY(IREC) ENDIF 520 FORMAT(5X,2X,I2,5X,A6,4X,A8) ENDIF C--> Go back to the top to read another record. IREC=IREC+1 GOTO 100 C--> This is the end of the FITS file. Now go through all fo the accumulated C record information and organize it: C - make a list of the header and data record numbers. C - display the information on the screen for verification 999 NHREC=0 NDREC=0 DO I=1,IREC-1 IF (RECTYPE(I).EQ.'DATA') THEN NDREC=NDREC+1 IDATA(NDREC)=I ENDIF IF (RECTYPE(I).EQ.'HEADER') THEN NHREC=NHREC+1 IHEAD(NHREC)=I ENDIF ENDDO IF (VERBOSE) THEN WRITE(*,530) IREC-1,NHREC,NDREC,NHEAD,DIMS 530 FORMAT(5X,45('-'), & /,5x,'Summary', & /,5X,45('-'), & /,5x,'Total Records : ',I3, & /,5x,'Header Records : ',I3, & /,5x,'Data Records : ',I3, & /,5x,'Header Lines : ',I3, & /,5x,'Image dimensions: ',A11, & /,5X,45('-')) ENDIF C--> This is the second pass over the file during which the header and C data are loaded into their main storage arrays. Since we know the C structure of the file, we can read it more easily. When the data C are being read in, make sure to read only the first row/column C if there are more than one. This is accomplished by loading one C pixel at a time from the buffer array to the data storage array and C keeping count of the number of pixels: once the total number of pixels C specified in the header is reached (=NAXIS1), stop there. IF (VERBOSE) THEN WRITE(*,540) FNAME 540 FORMAT(5X,'LOADING PASS: ',A30, & /,5X,45('-')) ENDIF NK=0 IF (NHREC.GT.1) THEN DO J=1,NHREC-1 READ(31,REC=IHEAD(J)) BUFF DO I=1,36 NK=NK+1 KEYWORD(NK)=BUFF(I)(1:8) HEAD(NK)=BUFF(I)(10:80) HLINE(NK)=BUFF(I) ENDDO ENDDO ENDIF READ(31,REC=IHEAD(NHREC)) BUFF DO I=1,36 NK=NK+1 KEYWORD(NK)=BUFF(I)(1:8) HEAD(NK)=BUFF(I)(10:80) HLINE(NK)=BUFF(I) ENDDO NP=0 IF (NDREC.GT.1) THEN DO J=1,NDREC-1 READ(31,REC=IDATA(J)) IARRAY DO I=1,1440 NP=NP+1 IF (NP.GT.NPIX*NROWS) GOTO 200 FLUX(NP)=IARRAY(I) ENDDO ENDDO ENDIF READ(31,REC=IDATA(NDREC)) IARRAY DO I=1,720 NP=NP+1 IF (NP.GT.NPIX*NROW) GOTO 200 FLUX(NP)=IARRAY(I) ENDDO C--> The header and data should be loaded now. Go through the header and C interpret some of the keywords for reference. The exposure time is C passed to the calling routine for future use. The existence of a C wavelength calibration is also checked by looking for the apropriate C keywords, and the value of the logical variable WSCALE is set C accordingly. Finally, the presence of the AIRMASS keyword is also C checked, and the value is passed to the calling routine if it exists. 200 DO I=1,NHEAD IF (KEYWORD(I)(1:7).EQ.'OBJCTRA') THEN CALL TRIM(HEAD(I),RASTRING,NCH) ENDIF IF (KEYWORD(I).EQ.'OBJCTDEC') THEN CALL TRIM(HEAD(I),DECSTRING,NCH) ENDIF IF (KEYWORD(I)(1:7).EQ.'EQUINOX') THEN READ(HEAD(I),*) EPOCH ENDIF IF (KEYWORD(I).EQ.'PLTSCALE') THEN READ(HEAD(I),*) PLTSCALE ENDIF IF (KEYWORD(I).EQ.'XPIXELSZ') THEN READ(HEAD(I),*) XPIXELSZ ENDIF IF (KEYWORD(I).EQ.'YPIXELSZ') THEN READ(HEAD(I),*) YPIXELSZ ENDIF ENDDO C--> Compute the angular dimensions of the chart using the plate scale C from the header. Note that the pixel size is in microns while the C plate scale is in "/mm. The final units are arcminutes. XSIZE=PLTSCALE*XPIXELSZ*1.E-3*NPIX/60. YSIZE=PLTSCALE*YPIXELSZ*1.E-3*NROWS/60. C--> Write out a final report for the user. Use the parameters values C extracted from the header as much as possible so that the user can C verify that these parameters have been extracted correctly. IF (VERBOSE) THEN WRITE(*,550) EPOCH,RASTRING,EPOCH,DECSTRING, & NINT(XSIZE),NINT(YSIZE) 550 FORMAT(5X,'Center: R.A. (',F6.1,') = ',A, & /,5X,' Dec. (',F6.1,') = ',A, & /,5X,'Projected Image Size: ',I3,' x ',I3,' arcmin', & /,5X,45('-')) ELSE WRITE(*,560) FNAME,NPIX,NROWS,NINT(XSIZE),NINT(YSIZE) 560 FORMAT(/,5X,A20,5X, & I4,' x ',I4,' pix',10X,I3,' x ',I3,' arcmin') ENDIF RETURN END C============================================================================ C S U B R O U T I N E T R I M C---------------------------------------------------------------------------- C This subroutine looks at a string variable as stored in a FITS header and C finds the range of non blank characters. This is useful for trimming off C leading and trailing blanks from a character variable. The method is based C on the fact that the string variable is enclosed in single quotes, so the C strategy is as follows: C - read the string int a buffer variable using an unformatted read so that C the string value between quotes is passed to the buffer C - scan the buffer varible starting from both ends until the first non C blank characters are found C - the positions of the first non blank characters (from either direction) C in the buffer are used to transfer the characters from the buffer to an C output sting (the trimmed string), which is then passed to the calling C routine C - the number of non blank characters is also passed to the calling routine C for reference C---------------------------------------------------------------------------- C Parameters C ---------- C STRING (C*80, input) : the character string to be scanned and trimmed C VALUE (C*80, output): the trimmed output string C NCHAR (I*4, output) : the number of non-black characters in the output C trimmed string C---------------------------------------------------------------------------- SUBROUTINE TRIM (STRING,VALUE,NCHAR) CHARACTER STRING*80,BUFFER*80,VALUE*80 C--> Read the STRINg into the BUFFER with the an internal unformatted READ READ(STRING,*) BUFFER C--> Go through the buffer (starting at the two ends) and find the C positions of the first non-blank characters J1=1 DO I=1,80 IF (BUFFER(I:I).NE.' '.AND.BUFFER(I+1:I+1).EQ.' ') GOTO 100 IF (BUFFER(I:I).EQ.' '.AND.BUFFER(I+1:I+1).NE.' ') THEN J1=I+1 GOTO 100 ENDIF ENDDO 100 J2=80 DO I=1,80 K=81-I IF (BUFFER(K:K).NE.' '.AND.BUFFER(K-1:K-1).EQ.' ') GOTO 110 IF (BUFFER(K:K).EQ.' '.AND.BUFFER(K-1:K-1).NE.' ') THEN J2=K-1 GOTO 110 ENDIF ENDDO C--> Transfer the right character range to the output string 110 NCHAR=J2-J1+1 WRITE(VALUE,600) BUFFER(J1:J2) 600 FORMAT(A) RETURN END