SUBROUTINE wrtfak(phaobs, varobs, region, areasc, qbckfil) IMPLICIT NONE INCLUDE '../../inc/xspec.inc' REAL phaobs(*), varobs(*), region(*), areasc(*) LOGICAL qbckfil C c wrtfak burton 15 April 94 c XSPEC subroutine to write out the accumulated fake data c routines, using FITS format c 3.0 - rewritten to write out data into FITS files c phaobs r i: The count rate to be written out if counting statistics was set. c else an integer number of actual counts is written out. c varobs r i: The errors on the count rate if Poisson statistics were c not assumed for the background. c region r i: The BACKSCAL values c areasc r i: The AREASCAL values c qbckfil l i: If true then write out the background file c c 7/21/97 - correctly write the keywords in the case of POISSERR = .FALSE. c 4/30/99 - change handling for pre-existing files: now dealt with before this function c is called. B. Dorman c 4/29/02 - added option to write out background files c INTEGER icounts, ierror, ichannl, iqualty, igroup REAL grasp, backscal, areascal INTEGER istat, maxchn, iunit, ifile, i, bitpix, ichan, ibin INTEGER status, nrows, nfield, pcount, block, chan0, icol INTEGER MAXCOL PARAMETER (MAXCOL=10) integer ridum character*20 rhdum character*3 rrchtyp real rrdum CHARACTER*11 cdate INTEGER yr,mo,day CHARACTER filenm*255, name*255 CHARACTER contxt*72, extens*10 CHARACTER cgroup*1 CHARACTER*20 tform(MAXCOL), ttype(MAXCOL), tunit(MAXCOL), extnam LOGICAL simple, extend, qqual, qgroup, qpoiss, qback, qarea LOGICAL qexist REAL dgtime, dgbtme, dgcrnm INTEGER dgndst, lenact, dgungc, dgnbnb, dgnbne, dgchn0, ftgsdt CHARACTER fgflnm*255, dgtlsc*20 CHARACTER fgauxf*255, dginst*20, dggcrd*1 LOGICAL dgqnot, dgqfst, dgdtps EXTERNAL lenact, fgflnm, fgauxf, dgtlsc, dgdtps EXTERNAL dgtime, dgbtme, dgcrnm, dgndst, dgungc, dginst EXTERNAL dgnbnb, dgnbne, dgchn0, dggcrd, dgqnot, ftgsdt status = 0 c Find largest number of channels in a data set. maxchn = 0 DO ifile = 1, dgndst() maxchn = MAX(maxchn, dgungc(ifile)) ENDDO c Allocate memory for workspace arrays. CALL udmget( maxchn, 6, icounts, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) CALL udmget( maxchn, 6, ierror, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) CALL udmget( maxchn, 4, ichannl, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) CALL udmget( maxchn, 4, iqualty, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) CALL udmget( maxchn, 4, igroup, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) c c Define values that are constant c bitpix = 8 block = 2880 status = ftgsdt(day,mo,yr,status) write(cdate,'(i4,"-",i2.2,"-",i2.2)') yr,mo,day status = 0 c c Loop around files to be written c CALL getlun(iunit) DO ifile = 1, dgndst() IF ( qbckfil .AND. dgbtme(ifile) .EQ. 0. ) GOTO 300 qpoiss = .TRUE. IF ( dgbtme(ifile) .GT. 0. .AND. .NOT.dgdtps(ifile,'b') ) & qpoiss = .FALSE. c Create new file. fakeit has already checked whether the user c wishes to overwrite this file, so if it does already exist, delete C it. filenm = fgflnm(ifile) IF ( qbckfil ) THEN extens = filenm(index(filenm,'.')+1:) filenm = filenm(:index(filenm,'.')-1)//'_bkg.'//extens ENDIF INQUIRE(file=filenm, exist=qexist) IF ( qexist ) CALL xdelfil(filenm, status) CALL ftinit(iunit, filenm, block, status) contxt = 'WRTFAK : Failed to create '//filenm(:lenact(filenm)) IF( status.NE.0 )GOTO 100 c c simple = .TRUE. extend = .TRUE. c c Define primary header CALL ftpdef(iunit, bitpix, 0, 0, 0, 1, status) contxt = 'WRTFAK : Failed to define primary header' IF( status.NE.0 )GOTO 100 c c Write primary array keywords c CALL ftphpr(iunit, simple, bitpix, 0, 0, 0, 1, extend, status) contxt = 'WRTFAK : Failed to write primary array keywords' IF( status.NE.0 )GOTO 100 c c Create extension c CALL ftcrhd(iunit, status) contxt = 'WRTFAK : Failed to create extension' IF( status.NE.0 )GOTO 100 c Find out whether the quality and grouping columns need to be c written and if so set the arrays qqual = .FALSE. qgroup = .FALSE. DO 150 i = 1, dgungc(ifile) memi(iqualty+i-1) = 0 memi(igroup+i-1) = 1 cgroup = dggcrd(ifile, i) IF( cgroup.EQ.'-' )THEN memi(igroup+i-1) = -1 qgroup = .TRUE. ELSEIF( cgroup.EQ.'*' )THEN memi(iqualty+i-1) = 2 qqual = .TRUE. ELSEIF( cgroup.EQ.' ' )THEN memi(iqualty+i-1) = 1 qqual = .TRUE. ENDIF 150 CONTINUE c Find out whether the BACKSCAL array needs to be written or whether we can c get away with a scalar keyword qback = .FALSE. backscal = region(dgnbnb(ifile)) DO i = dgnbnb(ifile)+1, dgnbne(ifile) IF ( ABS(backscal-region(i))/backscal .GT. 1.e-6 ) THEN qback = .TRUE. ENDIF ENDDO c Do the same for the AREASCAL qarea = .FALSE. areascal = areasc(dgnbnb(ifile)) DO i = dgnbnb(ifile)+1, dgnbne(ifile) IF ( ABS(areascal-areasc(i))/areascal .GT. 1.e-6 ) THEN qarea = .TRUE. ENDIF ENDDO c Write required header keywords for extension nrows = dgungc(ifile) ttype(1) = 'CHANNEL' tform(1) = 'I' tunit(1) = ' ' IF(qpoiss)THEN ttype(2) = 'COUNTS' tunit(2) = 'COUNTS' tform(2) = 'E' nfield = 2 ELSE ttype(2) = 'RATE' tunit(2) = 'cts/s' tform(2) = 'E' ttype(3) = 'STAT_ERR' tunit(3) = 'cts/s' tform(3) = 'E' nfield = 3 ENDIF IF( qqual )THEN nfield = nfield + 1 ttype(nfield) = 'QUALITY' tform(nfield) = 'I' tunit(nfield) = ' ' ENDIF IF( qgroup )THEN nfield = nfield + 1 ttype(nfield) = 'GROUPING' tform(nfield) = 'I' tunit(nfield) = ' ' ENDIF IF( qback )THEN nfield = nfield + 1 ttype(nfield) = 'BACKSCAL' tform(nfield) = 'E' tunit(nfield) = ' ' ENDIF IF( qarea )THEN nfield = nfield + 1 ttype(nfield) = 'AREASCAL' tform(nfield) = 'E' tunit(nfield) = ' ' ENDIF extnam = 'SPECTRUM' pcount = 0 CALL ftphbn(iunit, nrows, nfield, ttype, tform, tunit, extnam, & pcount, status) contxt = 'WRTFAK : Failed to write extension header keywords' IF( status.NE.0 )GOTO 100 CALL rdrinf(fgauxf(ifile,'r'),ridum,rhdum,ridum,rhdum,rhdum, & rhdum,rrchtyp,ridum,ridum,ridum,rrdum, & ridum,rrdum,ridum,ridum,istat) c Write additional header keywords for extension CALL ftpkys(iunit, 'HDUCLASS', 'OGIP', & 'format conforms to OGIP/GSFC conventions', status) contxt = 'WRTFAK : Failed to write HDUCLASS keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'HDUCLAS1', 'SPECTRUM', & 'extension contains a spectrum', status) contxt = 'WRTFAK : Failed to write HDUCLAS1 keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'HDUCLAS2', 'TOTAL', & 'spectrum is gross PHA spectrum', status) contxt = 'WRTFAK : Failed to write HDUCLAS2 keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'HDUCLAS3', 'COUNT', 'spectrum is counts', & status) contxt = 'WRTFAK : Failed to write HDUCLAS3 keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'HDUVERS1', '1.1.0', & 'version number of the format', status) contxt = 'WRTFAK : Failed to write HDUCLAS3 keyword' IF( status.NE.0 )GOTO 100 name = 'XSPEC '//VERSION//' fakeit' CALL ftpkys(iunit, 'CREATOR', name, & 'program that created this file', status) contxt = 'WRTFAK : Failed to write CREATOR keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'FILENAME', filenm, 'source file for fit', & status) contxt = 'WRTFAK : Failed to write FILENAME keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'DATE ', cdate, & 'FITS creation date (dd/mm/yyyy)', status) contxt = 'WRTFAK : Failed to write DATE keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'TELESCOP', dgtlsc(ifile), & 'telescope (mission) name', status) contxt = 'WRTFAK : Failed to write TELESCOP keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'INSTRUME', dginst(ifile), & 'instrument name', status) contxt = 'WRTFAK : Failed to write INSTRUME keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'FILTER ', 'none ', & 'instrument filter used', status) contxt = 'WRTFAK : Failed to write FILTER keyword' IF( status.NE.0 )GOTO 100 IF ( qbckfil ) THEN CALL ftpkye(iunit, 'EXPOSURE', dgbtme(ifile), 4, & 'Exposure time', status) ELSE CALL ftpkye(iunit, 'EXPOSURE', dgtime(ifile), 4, & 'Exposure time', status) ENDIF contxt = 'WRTFAK : Failed to write EXPOSURE keyword' IF( status.NE.0 )GOTO 100 c If we aren't writing a backscal column then put it in the keyword IF ( .NOT.qback ) THEN CALL ftpkye(iunit, 'BACKSCAL', backscal, 4, & 'Background scale factor', status) contxt = 'WRTFAK : Failed to write BACKSCAL keyword' IF( status.NE.0 )GOTO 100 ENDIF c If we aren't writing a areascal column then put it in the keyword IF ( .NOT.qarea ) THEN CALL ftpkye(iunit, 'AREASCAL', areascal, 4, & 'Area scale factor', status) contxt = 'WRTFAK : Failed to write AREASCAL keyword' IF( status.NE.0 )GOTO 100 ENDIF CALL ftpkye(iunit, 'CORRSCAL', dgcrnm(ifile), 4, & 'Correction scale factor', status) contxt = 'WRTFAK : Failed to write CORRSCAL keyword' IF( status.NE.0 )GOTO 100 name = ' ' CALL ftpkys(iunit, 'OBJECT ', name, 'object observed', status) contxt = 'WRTFAK : Failed to write OBJECT keyword' IF( status.NE.0 )GOTO 100 IF ( qbckfil ) THEN CALL ftpkys(iunit, 'BACKFILE', 'NONE', & 'bgrd FITS file', status) ELSE filenm = fgflnm(ifile) extens = filenm(index(filenm,'.')+1:) IF ( index(filenm,'.') .GT. 0 .AND. & dgbtme(ifile) .NE. 0. ) THEN filenm = filenm(:index(filenm,'.')-1)//'_bkg.'//extens CALL ftpkys(iunit, 'BACKFILE', filenm, & 'bgrd FITS file', status) ELSE CALL ftpkys(iunit, 'BACKFILE', 'NONE', & 'bgrd FITS file', status) ENDIF ENDIF contxt = 'WRTFAK : Failed to write BACKFILE keyword' IF( status.NE.0 )GOTO 100 IF ( qbckfil ) THEN CALL ftpkys(iunit, 'CORRFILE', 'NONE', & 'Correction FITS file', status) ELSE CALL ftpkys(iunit, 'CORRFILE', fgauxf(ifile,'c'), & 'Correction FITS file', status) ENDIF contxt = 'WRTFAK : Failed to write CORFFILE keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'RESPFILE', fgauxf(ifile,'r'), & 'Redistribution matrix file (RMF)', status) contxt = 'WRTFAK : Failed to write RESPFILE keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'ANCRFILE', fgauxf(ifile,'a'), & 'Anc. Resp. file (ARF)', status) contxt = 'WRTFAK : Failed to write ANCRFILE keyword' IF( status.NE.0 )GOTO 100 CALL ftpkys(iunit, 'CHANTYPE', rrchtyp, & 'Chan. assigned by det. electronics', status) contxt = 'WRTFAK : Failed to write CHANTYPE keyword' IF( status.NE.0 )GOTO 100 CALL ftpkyj(iunit, 'DETCHANS', dgungc(ifile), & 'Total no. detector channels available', status) contxt = 'WRTFAK : Failed to write DETCHANS keyword' IF( status.NE.0 )GOTO 100 CALL ftpkyj(iunit, 'TLMIN1', dgchn0(ifile), & 'Minimum value legally allowed in column 1', & status) contxt = 'WRTFAK : Failed to write TLMIN1 keyword' IF( status.NE.0 )GOTO 100 CALL ftpkyj(iunit, 'TLMAX1', dgungc(ifile)-(1-dgchn0(ifile)), & 'Maximum value legally allowed in column 1', & status) contxt = 'WRTFAK : Failed to write TLMIN1 keyword' IF( status.NE.0 )GOTO 100 IF ( qpoiss ) THEN CALL ftpkyj(iunit, 'STAT_ERR', 0, 'No stat error spcf.', & status) contxt = 'WRTFAK : Failed to write STAT_ERR keyword' IF( status.NE.0 )GOTO 100 ENDIF CALL ftpkyl(iunit, 'POISSERR', qpoiss, 'Pois. err. assumed ?', & status) contxt = 'WRTFAK : Failed to write POISSERR keyword' IF( status.NE.0 )GOTO 100 CALL ftpkyj(iunit, 'SYS_ERR', 0, 'No sys. err spcf.', status) contxt = 'WRTFAK : Failed to write SYS_ERR keyword' IF( status.NE.0 )GOTO 100 IF( .NOT.qgroup )THEN CALL ftpkyj(iunit, 'GROUPING', 0, 'No grpng of data spcf.', & status) contxt = 'WRTFAK : Failed to write GROUPING keyword' IF( status.NE.0 )GOTO 100 ENDIF IF( .NOT.qqual )THEN CALL ftpkyj(iunit, 'QUALITY', 0, 'No data qual info spcf.', & status) contxt = 'WRTFAK : Failed to write QUALITY keyword' IF( status.NE.0 )GOTO 100 ENDIF c Define extension data structure CALL ftbdef(iunit, nfield, tform, pcount, nrows, status) contxt = 'WRTFAK : Failed to define extension data structure' IF( status.NE.0 )GOTO 100 c c Write data to extension c chan0 = dgchn0(ifile) DO 200 i = 1, nrows memi(ichannl+i-1) = (i - 1) + chan0 200 CONTINUE CALL ftpclj(iunit, 1, 1, 1, nrows, memi(ichannl), status) contxt = 'WRTFAK : Failed to write CHANNEL data to extension' IF( status.NE.0 )GOTO 100 c Set up the COUNTS/RATE&ERROR arrays with zero values included where a channel c was ignored or grouped up. IF ( qbckfil ) THEN grasp = dgbtme(ifile) ELSE grasp = dgtime(ifile) ENDIF ibin = dgnbnb(ifile) - 1 ichan = 0 DO 250 i = 1, nrows memr(icounts+i-1) = 0. memr(ierror+i-1) = 0. IF( (dggcrd(ifile,i).EQ.'+' .OR. dggcrd(ifile,i).EQ.'*') ) & THEN ichan = ichan + 1 IF( dgqnot(ifile,ichan) )THEN ibin = ibin + 1 IF ( qpoiss ) THEN memr(icounts+i-1) = NINT(grasp*areasc(ibin) & *phaobs(ibin)) ELSE memr(icounts+i-1) = areasc(ibin)*phaobs(ibin) memr(ierror+i-1) = areasc(ibin)*sqrt(varobs(ibin)) ENDIF ENDIF ENDIF 250 CONTINUE CALL ftpcle(iunit, 2, 1, 1, nrows, memr(icounts), status) IF(qpoiss)THEN contxt='WRTFAK : Failed to write COUNTS data to extension' ELSE contxt='WRTFAK : Failed to write RATE data to extension' ENDIF IF( status.NE.0 )GOTO 100 icol = 2 IF ( .NOT. qpoiss ) THEN icol = icol + 1 CALL ftpcle(iunit, icol, 1, 1, nrows, memr(ierror), status) contxt='WRTFAK : Failed to write STAT_ERR data to extension' IF( status.NE.0 )GOTO 100 ENDIF IF( qqual )THEN icol = icol + 1 CALL ftpclj(iunit, icol, 1, 1, nrows, memi(iqualty), status) contxt = & 'WRTFAK : Failed to write QUALITY data to extension' IF( status.NE.0 )GOTO 100 ENDIF IF ( qgroup )THEN icol = icol + 1 CALL ftpclj(iunit, icol, 1, 1, nrows, memi(igroup), status) contxt = & 'WRTFAK : Failed to write GROUPING data to extension' IF( status.NE.0 )GOTO 100 ENDIF c If we need to write BACKSCAL as a vector then do so. As with the actual c data we need to put values in for channels that have been ignored or c grouped up. For these channels we set the BACKSCAL to that for the last c valid bin. This is not perfect but will have to do. IF ( qback ) THEN icol = icol + 1 ibin = dgnbnb(ifile) - 1 ichan = 0 DO i = 1, nrows IF ( ibin .LT. dgnbnb(ifile) ) THEN memr(icounts+i-1) = region(dgnbnb(ifile)) ELSE memr(icounts+i-1) = region(ibin) ENDIF IF( (dggcrd(ifile,i).EQ.'+' .OR. & dggcrd(ifile,i).EQ.'*') ) THEN ichan = ichan + 1 IF( dgqnot(ifile,ichan) )THEN ibin = ibin + 1 memr(icounts+i-1) = region(ibin) ENDIF ENDIF ENDDO CALL ftpcle(iunit, icol, 1, 1, nrows, memr(icounts), & status) contxt = & 'WRTFAK : Failed to write BACKSCAL data to extension' IF( status.NE.0 )GOTO 100 ENDIF c If we need to write AREASCAL as a vector then do so. As with the actual c data we need to put values in for channels that have been ignored or c grouped up. For these channels we set the AREASCAL to that for the last c valid bin. This is not perfect but will have to do. IF ( qarea ) THEN icol = icol + 1 ibin = dgnbnb(ifile) - 1 ichan = 0 DO i = 1, nrows IF ( ibin .LT. dgnbnb(ifile) ) THEN memr(icounts+i-1) = areasc(dgnbnb(ifile)) ELSE memr(icounts+i-1) = areasc(ibin) ENDIF IF( (dggcrd(ifile,i).EQ.'+' .OR. & dggcrd(ifile,i).EQ.'*') ) THEN ichan = ichan + 1 IF( dgqnot(ifile,ichan) )THEN ibin = ibin + 1 memr(icounts+i-1) = areasc(ibin) ENDIF ENDIF ENDDO CALL ftpcle(iunit, icol, 1, 1, nrows, memr(icounts), & status) contxt = & 'WRTFAK : Failed to write AREASCAL data to extension' IF( status.NE.0 )GOTO 100 ENDIF c Write out error message if necessary 100 CONTINUE IF( status.NE.0 )THEN CALL xwrite(contxt, 10) WRITE(contxt, '(a, i4, a, i4)') 'FITSIO error = ', status, & 'for file number ', ifile CALL xwrite(contxt, 10) ENDIF c c Close file c CALL ftclos(iunit, status) 300 CONTINUE ENDDO CALL frelun(iunit) CALL udmfre(icounts, 6, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) CALL udmfre(ierror, 6, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) CALL udmfre(ichannl, 4, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) CALL udmfre(iqualty, 4, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) CALL udmfre(igroup, 4, istat ) IF( istat .NE. 0 ) CALL xwrite( & ' WARNING: Memory allocation failed in WRTFAK.', 2) RETURN END