!+MARFRMF SUBROUTINE MARFRF IMPLICIT NONE ! ! Description: ! Program to multiply/divide a detector redistribution matrix through by an ! Ancillary Response file or scalar ! NOTE - The par file is assumed to have been opened. ! The RMF extension will contain a variable length array for the ! "MATRIX" column if (and only if) this leads to a significant ! saving in disk space. ! ! Passed Parameters ! None ! ! User i/ps required (prompted for): ! None here, isolated in GP_MARFRMF (see below) ! ! Include/Common blocks files ! common TASK : (FTOOLS) standard fatal error message thingy ! ! Called routines ! subroutine GP_MARFRMF : (below) Gets parameters from XPI par file ! subroutine DO_MARFRMF : (below) Performs the conversion ! ! Compilation: ! subroutines require CALLIB, FTOOLS, FITSIO ! ! Origin: ! Original ! ! Authors/Modification History: ! Ian M George (1.0.0:1992 Nov 15), quick & v.dirty ! Ian M George (1.0.1:1994 Jan 10), reduced max array sizes ! Ian M George (2.0.0:1994 Mar 07), added Dynamic Memory Allocation ! Ian M George (2.1.0:1994 Jun 27), DMA requirements read from i/p file ! Rehana Yusaf (2.1.1:1994 Sep 13), minor changes, clobber read ! ck_file has additional argument ! Ian M George (2.1.2:1995 Mar 14), bugettes fixed in gp_marfrmf ! Ian M George (2.2.0:1995 Apr 08), added scalar & qdivide parameters ! Ian M George (2.3.0:1995 Jun 21), DO_MARFRMF calls wtebd2 not wtebd1 ! Ian M George (2.3.1:1995 Jun 22), reads chantyp from RMF & writes to o/p ! Ian M George (2.3.2:1995 Aug 21), minor additional error checking ! Ian M George (2.3.3:1995 Nov 22), DO_MARFRMF calls rdrmf2 & wtrmf2 now ! Ian M George (2.3.4:1995 Nov 29) Added wtinfo, wterrm, wtwarm etc ! Ian M George (2.3.5:1995 Dec 18) Added wtbegm & wtendm changes ! Ian M George (3.0.0:96 Oct 04) updated the depths to use rdrmf3,rdebd3, ! wtrmf3 ! & wtebd3 subroutines, which try their best to ! not to get the indexing of the matrix incorrect ! for detectors whose first channel is numbered ! channel zero ! Fixed bug whereby EBOUNDS extension always read ! from rmfil irrespective of value of ebfil param ! Ian M George (3.0.1:96 Nov 05) tweak so rmfversn=1.2.0 if not read from par ! Banashree M Seifert (3.1.0) ! . only in gp_marfrmf subroutine, when asks for input files ! first it asks for ARF then RMF file name. Before (v3.0.1) ! it was otherway round. User didnot like it ! ! Banashree M Seifert (3.2.0, Mar 12, 1997) ! . initialised ierr in do_mathpha ! Peter D Wilson (3.2.1, Feb 20, 1998) ! . Copy HISTORY keywords from input RMF and EBOUNDS extensions ! to the new output file ! Peter D Wilson (3.2.2, Aug 11, 1999) ! . Call wterrm instead of wtferr when rdrmf4 fails ! Ning Gan (3.2.3, Aug 30, 1999) ! . The imaxarf was set to the wrong value ( read from the ! wrong extension if the arf extension is not first one). ! . Correct the bug in rdrmf4. ! Ning Gan (3.2.4, Oct 14, 1999) ! . Switched the order of reading parameter for filename of arf ! and rmf. ! M Tripicco (3.2.5, Oct 31, 2000) ! . upped history/comment string lengths to 72 to match CFITSIO ! default and fcrhky() which starts copying from col 9 now. ! kaa (3.2.6, Jul 18, 2005) ! added (hidden) arfcol option so ARF can be read from a column ! other than SPECRESP. ! Bryan Irby (3.3.0, Dec 04, 2018): ! . Replaced udmget/udmfre with allocate/deallocate ! MFC (3.3.1 Apr 16 2020): ! . f90 version; lengthened some truncated character strings ! ------------------------------------------------------------------ character(7) version parameter (version = '3.3.1') !- ! Internals character(40) taskname integer chatter, ierr real scalar character(5) rmfversn character(20) telescop, instrume, detnam, filter, arfcol character(80) rmfil, arfil, outfil, ebfil logical qoverride, killit, qdivide ! Initialize taskname = 'marfrmf' ierr = 0 ! Get Parameters from the par file call gp_marfrmf(rmfil, ebfil, arfil, outfil, chatter, rmfversn, & qoverride, telescop, instrume, detnam, filter, killit, & scalar, qdivide, arfcol, ierr) if(ierr.NE.0) goto 148 ! Start-up MAIN call wtbegm(taskname, version, chatter) ! Do the nasty deed call do_marfrmf(taskname, version, & rmfil, ebfil, arfil, outfil, chatter, rmfversn, & qoverride, telescop, instrume, detnam, filter, killit, & scalar, qdivide, arfcol, ierr) ! Finish-Off 148 continue call wtendm(taskname, version, ierr, chatter) return end ! ------------------------------------------------------------------------- !+GP_MARF subroutine gp_marfrmf(rmfil, ebfil, arfil, outfil, chatter, & rmfversn, qoverride, telescop, instrume, detnam, filter, killit, & scalar, qdivide, arfcol, ierr) IMPLICIT NONE integer chatter, ierr real scalar character*(*) rmfversn character*(*) telescop, instrume, detnam, filter, arfcol character*(*) rmfil, arfil, outfil, ebfil logical qoverride, killit, qdivide ! ! Description: ! Gets the parameters required by MARFRMF from the parameter file ! NOTE - The par file is assumed to have been opened. ! ! User i/ps required (prompted for): ! RMFIL - file name of i/p RSP_MATRIX extn (redistribution matrix) ! EBFIL - file name of i/p EBOUNDS extn (nom chan-energy boundaries) ! ARFIL - file name of i/p SPECRESP extn (ancillary response dataset) ! OUTFIL - name of o/p RMF (redistribution matrix) to be created ! CHATTER - chattiness flag for o/p (5 quite,10 normal,15 high,>20 silly) ! RMFVERSN - OGIP version of RMF file required ! QOVERRIDE - Whether tele,instr,filt etc keywords from par file are ! to be used, rather than those derived from RSP file ! TELESCOP - if(QOVERRIDE) Telescope/Mission name to be used ! INSTRUME - if(QOVERRIDE) Instrument/Detector name to be used ! DETNAM - if(QOVERRIDE) Supplimentary Detector name to be used (if reqd) ! FILTER - if(QOVERRIDE) Filter name to be used ! KILLIT - Something FTOOLS Grp made me add ! SCALAR - if(ARFIL=NONE) Scalar value by which RSP_MATRIX multiplied ! QDIVIDE - Flag whether * or / performed (qdivide = T for /) ! ARFCOL - Option to use a column other than SPECRESP for the ARF data ! ! Origin: ! Original ! ! Called Routines ! subroutine FCECHO : (FTOOLS) writes to standard o/p device ! subroutine WT_FERRMSG : (CALLIB) Writes standard FITSIO message etc ! ! Compilation: ! requires XPI/Host interface etc and CALLIB ! ! Authors/Modification History: ! Ian M George (1.0.0: 1992 Dec 23), Original ! Ian M George (1.0.1: 1992 Dec 23), Fixed error handling ! Rehana yusaf (1.0.2: 1994 Sep 13), clobber read ! Ian M George (1.1.0: 1995 Mar 14), clobber read as 'b' rather than 's' ! Ian M George (2.0.0: 1995 Apr 08), added scalar & qdivide parameters ! Ian M George (2.0.1: 1995 Nov 29) Added wtinfo, wterrm, wtwarm etc ! ! Banashree M Seifert (2.1.0) ! . when asks for input files first it asks for ARF then RMF ! file name. Before (v2.0.1) it was otherway round. ! User didnot like it ! kaa (2.2.0) added arfcol parameter ! ----------------------------------------------------------------------- character(7) version parameter (version = '2.2.0') !- ! Internals character(10) subname parameter (subname = 'gp_marfrmf') character(80) message character(80) tmpfile ! Initialize ierr = 0 telescop = 'UNKNOWN' instrume = 'UNKNOWN' detnam = ' ' filter = 'NONE' scalar = 0.0 ! Get the name of the RMF call uclgst('rmfil', rmfil, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting RMFIL parameter') goto 999 endif ! Get the name of the ARF call uclgst('arfil', arfil, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting ARFIL parameter') goto 999 endif call crmvlbk(arfil) ! Check for a null filename, implying a scalar will be entered tmpfile = arfil call ftupch(tmpfile) if((tmpfile(1:1).EQ.' ').OR.& (tmpfile(1:4).EQ.'NONE')) then message = 'Assuming you wish to simply operate using ' // & 'a scalar' call wtinfo(chatter, 1, 1, message) call uclgsr('scalar', scalar, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting SCALAR parameter') goto 999 else arfil = 'NONE' endif if(scalar.EQ.1.0) then call wtwarm(subname, version, chatter, 1, 'SCALAR = 1.0') call wtinfo(chatter, 1, 1, & ' The requested output file identical to ' // & 'input file') elseif(scalar.EQ.0.0) then call wterrm(subname, version, & 'SCALAR = 0.0 is meaningless') goto 999 elseif(scalar.LT.0.0) then call wtwarm(subname, version, chatter, 1, 'SCALAR < 0.0') call wtinfo(chatter, 1, 1, 'I will try, ' // & 'but have no idea what it means') endif endif ! Get the name of the o/p file call uclgst('outfil', outfil, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting OUTFIL parameter') goto 999 endif ! Get the name of file containing EBOUNDS extn call uclgst('ebfil', ebfil, ierr) if(ierr.NE.0) then call wtwarm(subname, version, chatter, 1, & 'Problem getting EBFIL parameter') call wtinfo(chatter, 1, 1, 'setting EBFIL = RMFIL') ebfil = rmfil ierr = 0 elseif(ebfil(:1).EQ.'%')then ebfil = rmfil endif ! Get the OGIP version number of the RMF file format to be created call uclgst('rmfversn', rmfversn, ierr) if(ierr.ne.0) then call wtwarm(subname, version, chatter, 1, & 'Problem getting RMFVERSN parameter') rmfversn = '1.2.0' message = 'setting RMFVERSN = ' // rmfversn call wtinfo(chatter, 1, 1, message) ierr = 0 elseif(rmfversn.EQ.'1992a') then rmfversn = '1.2.0' endif ! Get the column to be used to read the ARF call uclgst('arfcol', arfcol, ierr) if(ierr.ne.0) then call wtwarm(subname, version, chatter, 1, & 'Problem getting ARFCOL parameter, assuming SPECRESP') ierr = 0 arfcol = 'SPECRESP' endif ! Get the chattiness flag call uclgsi('chatter', chatter, ierr) if(ierr.NE.0) then call wtwarm(subname, version, chatter, 1, & 'Problem getting CHATTER parameter') ierr = 0 call wtinfo(chatter, 1, 1, 'setting CHATTER = 10') chatter = 10 endif ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 20, 1, message) ! Check whether the RSP_MATRIX extension is to be multiplied or ! divided by the SPECRESP dataset (or scalar) ! ... qdivide = T for DIVIDE call uclgsb('qdivide', qdivide, ierr) if(ierr.NE.0) then call wterrm(subname, version, & 'Problem getting QDIVIDE parameter') call wtinfo(chatter, 1, 1, ' Too dangerous to continue') goto 999 endif ! Get the Override flag as to whether Telecop,Instr etc names derived ! from the RMF file are to be ignored, and those listed below used ! instead call uclgsb('qoverride', qoverride, ierr) if(ierr.NE.0) then call wtwarm(subname, version, chatter, 1, & 'Problem getting QOVERRIDE parameter') ierr = 0 call wtinfo(chatter, 1, 1, ' setting QOVERRIDE = FALSE') qoverride = .false. endif if(qoverride) then ! ... OK, we're gonna override, so get the necessary parameters call uclgst('telescop', telescop, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting TELESCOP parameter') goto 999 endif call uclgst('instrume', instrume, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting INSTRUME parameter') goto 999 endif call uclgst('detnam', detnam, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting DETNAM parameter') goto 999 endif call uclgst('filter', filter, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting FILTER parameter') goto 999 endif endif call uclgsb('clobber', killit, ierr) if(ierr.ne.0) then call wtwarm(subname, version, chatter, 1, & 'Problem getting CLOBBER parameter') ierr = 0 message = ' setting CLOBBER = FALSE' call wtinfo(chatter, 1, 1, message) killit = .false. endif ! Check for errors 999 if(ierr.ne.0) then call wterrm(subname, version, ' unable to continue') endif return end ! ------------------------------------------------------------------------- !+DO_MARFRMF subroutine do_marfrmf(tname, tvers, & rmexp, ebexp, arexp, outfil, chatter, & rmfversn, qoverride, telescop, instrume, detnam, filter, & killit, scalar, qdivide, arfcol, ierr) IMPLICIT NONE ! - Alex: added new variables numgrp, order, qorder, isorder ! integer chatter, ierr !integer numgrp real scalar character*(*) rmfversn, tname, tvers character*(*) telescop, instrume, detnam, filter, arfcol character*(*) rmexp, arexp, outfil, ebexp logical qoverride, killit, qdivide, qorder, isorder ! ! Description: ! Program to multiply an RMF redistribution matrix through by an Ancillary ! response file (or scalar). ! ! User i/ps required (prompted for): ! None ! ! Passed parameters ! RMEXP - filename[ext] of i/p RSP_MATRIX extn (redistribution matrix) ! EBEXP - filename[ext] of i/p EBOUNDS extn (nom chan-energy boundaries) ! AREXP - filename[ext] of i/p SPECRESP extn (ancillary response dataset) ! OUTFIL - name of o/p RMF (redistribution matrix) to be created ! CHATTER - chattiness flag for o/p (5 quite,10 normal,15 high,>20 silly) ! RMFVERSN - OGIP version of RMF file required ! QOVERRIDE - Whether tele,instr,filt etc keywords from par file are ! to be used, rather than those derived from RSP file ! TELESCOP - if(QOVERRIDE) Telescope/Mission name to be used ! INSTRUME - if(QOVERRIDE) Instrument/Detector name to be used ! DETNAM - if(QOVERRIDE) Supplimentary Detector name to be used (if reqd) ! FILTER - if(QOVERRIDE) Filter name to be used ! SCALAR - if(ARFIL=NONE) Scalar value by which RSP_MATRIX multiplied ! QDIVIDE - Flag whether * or / performed (qdivide = T for /) ! ARFCOL - Column to use to get ARF data (standard is SPECRESP) ! IERR o : Error Flag (zero if all OK) ! ! Called Routines: ! subroutine CGETLUN : (CALLIB) Gets a free logical unit ! subroutine CK_FILE : (CALLIB) Checks if file exists, deletes if reqd ! subroutine CLNSTR : (CALLIB) Cleans strings for fitsio ! subroutine FCECHO : (FTOOLS) Writes to standard o/p ! subroutine FNDEXT : (CALLIB) Finds xtens based on EXTNAME value ! subroutine FNDHDU : (CALLIB) Finds xtens based on HDUCLASn values ! subroutine FTCLOS : (FITSIO) Closes a FITS file ! subroutine FTMAHD : (FITSIO) Move to an absolute xtens no. ! subroutine FTMRHD : (FITSIO) Move a relative no. xtens rel to CDU ! subroutine FTOPEN : (FITSIO) Open a FITS file ! subroutine FTPKYx : (FITSIO) writes a keyword of type "x" ! subroutine FTVERS : (FITSIO) Returns version of FITSIO ! subroutine OP_NPA : (CALLIB) Opens & writes a null P.Header ! subroutine RDARF1 : (CALLIB) Reads an SPECRESP xtens ! subroutine RDEBD1 : (CALLIB) Reads an EBOUNDS xtens ! subroutine RDRMF1 : (CALLIB) Reads an RSP_MATRIX xtens ! subroutine REMAP : (CALLIB) Remaps 1-d array onto new grid ! subroutine WT_FERRMSG : (CALLIB) Writes standard FITSIO error message ! subroutine WTRMF1 : (CALLIB) Writes the RMF Xtensn (RMFVERSN = 1992a) ! subroutine WTEBD2 : (CALLIB) Writes the EBOUNDS Xtensn (RMFVERSN = 1992a) ! subroutine RDRMF4 : (CALLIB) Reads an RSP_MATRIX xtens(RMFVERSN = 1999a) ! subroutine WTRMF4 : (CALLIB) Writes the RMF Xtensn (RMFVERSN = 1999a) ! ! Origin: ! Original ! ! Authors/Modification History: ! Ian M George (1.0.0:1993 Oct 20), original ! Ian M George (1.0.1:1994 Jan 10), reduced max array sizes ! Ian M George (2.0.0:1994 Mar 01), added Dynamic Memory Allocation ! Rehana Yusaf (2.0.1:1994 sep 13), added killit argument ! Ian M George (3.0.0:1995 Apr 08), added scalar & qdivide parameters ! Ian M George (3.0.1:1995 Jun 21), call wtebd2 instead of wtebd1 ! Ian M George (3.0.2:1995 Jun 22), reads chantyp from RMF & writes to o/p ! Ian M George (3.0.3:1995 Aug 21), minor additional error checking ! Ian M George (3.0.4:1995 Nov 22), rdrmf1->rdrmf2 & wtrmf1->wtrmf2 ! Ian M George (4.0.0:96 Oct 04) updated to use rdrmf3,rdebd3,wtrmf3 ! & wtebd3 subroutines, which try their best to ! not to get the indexing of the matrix incorrect ! for detectors whose first channel is numbered ! channel zero ! Fixed bug whereby EBOUNDS extension always read ! from rmfil irrespective of value of ebfil param ! M Tripicco (4.0.1:97 Jan 03) added initialization of status ! ! Banashree M Seifert (4.1.0, Mar 12, 1997) ! . error status initialised (line #1203) ! ( Mar 18, 1997) maxgrp made 10 instead of 5 (for HXTE data) ! Banashree M Seifert (4.2.0, June6, 1997) ! . warning messages when looking for new HDUCLAS are made ! vissible at high chatter since otherwise it confuses ! users ! Peter D Wilson (4.2.1, Feb 20, 1998) ! . Copy HISTORY keywords from RMF and EBOUNDS extensions ! Alex M. (4.2.2, Feb 8, 1999) Incorporated new subroutines RDRMF4 and ! WTRMF4 written by Dr. Keith A. Arnaud. ! Modified subroutines multmatrix and ! divmatrix. Fixed bug in 'call fndhdu()' ! associated with passing a wrong parameter ! (iunit1 instead of iunix3). ! Ning Gan (4.2.3, Aug 30, 1999) ! . The imaxarf was set to the wrong value ( read from the ! wrong extension if the arf extension is not first one). ! kaa (4.3.0) Added arfcol option to read ARF from column other than SPECRESP ! ------------------------------------------------------------------------ character(7) version parameter (version = '4.4.0') !- character(10) subname parameter (subname = 'do_marfrmf') ! Maximum Array sizes ( Alex -- added maxelt, numelt ) ! --------------------------------------------------------------------------- integer maxchan, maxne, maxgrp, maxelt, numelt integer maxhist, maxcomm, maxextn parameter (maxhist = 100, maxcomm = 20, maxextn = 99) ! --------------------------------------------------------------------------- ! Commons character(40) taskname COMMON/task/taskname ! Internals ! ... parameters & statis arrays integer block, imove, i, status, jj, htype integer ninstr, nsearch, nfound, clenact integer rmextn, ebextn, arextn, fchan integer next(maxextn) ! character(8) dummy8 character(20) dummy20 character(20) instr(9) character(20) outhdu(9, maxextn), outver(9, maxextn) character(20) extnam(maxextn) character(20) rmfhdu(9), outhdu3 character(20) atelescop, ainstrume, adetnam, afilter character(30) rchantyp, echantyp character(80) rmfil, arfil, ebfil character(80) rmf_rmfversn, ebd_rmfversn integer nk, nk_prhist, nk_rmfhist, nk_ebhist, nk_comm integer ounit, iunit1, iunit2, iunit3 integer imaxgrp, imsg, rflchan integer ichan, iebound, icol integer rien, aien real ftsver real area, lo_thresh character(5) arfversn character(8) key(4) character(20) old(4), new(4) character(20) old_tele, old_inst, old_detn, old_filt ! upped hist/comm strings to 72 since CFITSIO starts them in col 9 now character(110) prhist(maxhist), rmfhist(maxhist), ebhist(maxhist) character(110) comment(maxcomm), newhist character(30) comm character(100) message, dummystr(1) character(80) rmf_extname, rmf_extcomm character(80) ebd_extname, ebd_extcomm logical qokfil, anyflg integer idma, igot integer imaxne, imaxchan, imaxelt, imxarf ! ... pointers to "arrays" to be dynamically allocated ! --- Alex added p_order integer, allocatable :: p_ngrp(:), p_F_chan(:), p_N_chan(:) integer, allocatable :: p_chan(:), p_order(:) real, allocatable :: p_matrix(:), p_sprsp(:), p_factor(:) real, allocatable :: p_re_lo(:), p_re_hi(:) real, allocatable :: p_ae_lo(:), p_ae_hi(:) real, allocatable :: p_e_min(:), p_e_max(:) ! ... "arrays" to be dynamically allocated ! integer ngrp(imaxne) integer F_chan(imaxne,maxgrp) ! integer N_chan(imaxne,maxgrp) real matrix(imaxelt) ! real sprsp(imxarf) real factor(imaxne) ! real re_lo(imaxne) real re_hi(imaxne) ! real ae_lo(imxarf) real ae_hi(imxarf) ! real e_min(imaxchan) real e_max(imaxchan) ! integer chan(imaxchan) integer order(maxgrp) ! ... Other logical qscalar ! Initialize status = 0 echantyp = 'UNKNOWN' rchantyp = 'UNKNOWN' qscalar = .false. ierr = 0 dummystr(1) = ' ' ! -- Alex --------------- isorder = .false. qorder = .false. ! ------------------------ ! Check out whether an ARF or a scalar has been requested if((arexp.EQ.'NONE').AND.(scalar.NE.0.0)) then qscalar = .true. else qscalar = .false. endif ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 20, 1, message) call ftvers(ftsver) write(message, '(a,f6.3)')& ' using FITSIO Version ', ftsver call wtinfo(chatter, 20, 2, message) ! Check that the o/p file doesn't already exist or is illegal call ck_file(outfil, dummystr, 1, qokfil, killit, chatter) if(.NOT.qokfil) then message = 'OUTFIL parameter found offensive to ' // & subname // ' ' // version call wtinfo(chatter, 1, 1, message) ierr = -1 goto 482 endif ! Store the old keyword values (required below if override requested) old_tele = telescop old_inst = instrume old_detn = detnam old_filt = filter ! Parse the supplied filenames, stripping off incld extension numbers call fcpars(rmexp, rmfil, rmextn, status) if(status.NE.0) then message = ' Problem passing rmfil expression' call wtfwrn(subname, version, chatter, 1, status, message) call wtinfo(chatter, 1, 1, ' will search all extensions') rmextn = -99 endif call fcpars(ebexp, ebfil, ebextn, status) if(status.NE.0) then message = ' Problem passing ebfil expression' call wtfwrn(subname, version, chatter, 1, status, message) call wtinfo(chatter, 1, 1, 'will search all extensions') ebextn = -99 endif if(ebexp.EQ.rmexp) ebextn = -99 if(.not.qscalar) then ! --- Alex: initialized arextn arextn = 0 call fcpars(arexp, arfil, arextn, status) if(status.NE.0) then message = ' Problem passing arfil expression' call wtfwrn(subname, version, chatter, 1, status, message) call wtinfo(chatter, 1, 1, 'will search all extensions') arextn = -99 endif endif ! Open i/p RSP_MATRIX file status = 0 call cgetlun(iunit1) call ftopen(iunit1, rmfil, 0, block, status) if(status.ne.0) then message = ' opening RMF file: ' // rmfil(:20) call wtferr(subname, version, status, message) ierr = 1 goto 482 endif ! Read HISTORY keywords from primary header of RSP_MATRIX file nk_prhist = 0 call fcrhky(iunit1, maxhist - 1, nk_prhist, prhist, status) if(status.ne.0) then message = ' reading RMF HISTORY keywords' call wtferr(subname, version, status, message) ierr = 1 goto 482 endif ! Open i/p EBOUNDS file, if necessary if(ebfil.NE.rmfil) then status = 0 call cgetlun(iunit2) call ftopen(iunit2, ebfil, 0, block, status) if(status.ne.0) then message = ' opening EBOUNDS file: ' // ebfil(:20) call wtferr(subname, version, status, message) ierr = 1 goto 482 endif else iunit2 = iunit1 endif ! Open i/p ARF file and get the number of energy bins if(.not.qscalar) then call cgetlun(iunit3) call ftopen(iunit3, arfil, 0, block, status) if(status.ne.0) then message = ' opening ARF file: ' // arfil(:20) call wtferr(subname, version, status, message) ierr = 1 goto 482 endif endif ! Fill in a few history & comment records newhist = 'RSP_MATRIX dataset convolved with SPECRESP array' // & ' by MARFRMF ' // tvers nk_prhist = nk_prhist + 1 prhist(nk_prhist) = newhist comment(1) = 'MARFRMF ' // tvers // 'Summary:' comment(2) = ' I/p RMF file: ' // rmexp if(ebexp.EQ.rmexp) then comment(3) = ' I/p EBOUNDS file: ' // ebfil else comment(3) = ' I/p EBOUNDS file: ' // ebexp endif comment(4) = ' I/p ARF file: ' // arexp nk_comm = 4 if(qscalar) then nk_comm = nk_comm + 1 write(comment(nk_comm), '(a,g12.6)')& ' Scalar factor: ', scalar endif nk_comm = nk_comm + 1 if(qdivide) then if(qscalar) then comment(nk_comm) = & ' ... I/p RMF dataset DIVIDED by scalar' else comment(nk_comm) = & ' ... I/p RMF dataset DIVIDED by ARF dataset' endif else if(qscalar) then comment(nk_comm) = & ' ... I/p RMF dataset MULTIPLIED by scalar' else comment(nk_comm) = & ' ... I/p RMF dataset MULTIPLIED by ARF dataset' endif endif ! -------------------------- RSP_MATRIX extension -------------------- ! Find the RSP_MATRIX extension in the RMF file ! - Extension number NOT given as part of rmexp (search for HDUCLAS/EXTNAM) if(rmextn.LT.0) then ninstr = 2 instr(1) = 'RESPONSE' instr(2) = 'RSP_MATRIX' nsearch = maxextn call fndhdu(chatter, iunit1, ninstr, instr, & nsearch, nfound, next, outhdu, outver, extnam, ierr) ! ... check for old-style EXTNAME values if no OK HDUCLASn values found if(nfound.LE.0) then message = ' Ext w/ allowed HDUCLASn keywrds not found' call wtwarm(subname, version, chatter, 20, message) message = ' offending file: ' // rmfil call wtinfo(chatter, 20, 1, message) message = ' searching for extnsion with EXTNAME = MATRIX' call wtinfo(chatter, 20, 2, message) call fndext(chatter, iunit1, 'MATRIX', & nsearch, nfound, next, outhdu, outver, extnam, ierr) endif ! - Extension number IS given as part of rmfil else call ftmahd(iunit1, rmextn + 1, htype, status) message = ' Problem moving to specified xtens' call wtferr(subname, version, status, message) ! ... grab the HDUCLAS values for reference ninstr = 1 instr(1) = '*' nsearch = 1 call fndhdu(chatter, iunit1, ninstr, instr, & nsearch, nfound, next, outhdu, outver, extnam, ierr) nfound = 1 next(1) = 0 endif ! - sort out what we've got if(nfound.GT.1) then call wterrm(subname, version, & 'RMFIL contains >1 RSP_MATRIX datasets') write(message, '(i12,a)') nfound, ' extensions found:' call wtinfo(chatter, 1, 1, message) do i = 1, nfound write(message, '(a,i12,a)') 'Ext ', next(i), ':' call wtinfo(chatter, 1, 2, message) write(message, '(a,a)') 'EXTNAME = ', extnam(i) call wtinfo(chatter, 1, 3, message) do jj = 1, 4 write(message, '(a,i2,2a)')& 'HDUCLAS', jj, ' = ', outhdu(i, jj) call wtinfo(chatter, 1, 3, message) enddo enddo message = & ' ... Extension number must be specified via rmfil param' call wterrm(subname, version, message) ierr = 2 goto 482 elseif(nfound.LE.0) then message = ' Unable to locate an RSP_MATRIX extension' call wterrm(subname, version, message) ierr = 2 goto 482 else do i = 1, 9 rmfhdu(i) = outhdu(i, 1) enddo outhdu3 = 'FULL' endif ! Move to the Extension if not already there if(next(1).GT.0) then imove = next(1) status = 0 call ftmrhd(iunit1, imove, htype, status) message = ' Problem moving to RSP-MATRIX xtens' call wtferr(subname, version, status, message) endif ! Alex: Get the values of RMF arrays sizes maxne, maxgrp, and maxelt maxne = 0 maxgrp = 0 maxelt = 0 call rmfsz(iunit1, chatter, maxne, maxgrp, maxelt, status) if(status.ne.0) then message = ' reading RMF array sizes maxne,maxgrp, maxelt' call wtferr(subname, version, status, message) ierr = 1 goto 482 endif ! Set max sizes of arrays imaxne = maxne imaxgrp = maxgrp imaxelt = maxelt ! Read the dimensions of the RMF file maxchan = 0 call ftgkyj(iunit1, 'DETCHANS', maxchan, comm, status) if(status.NE.0) then call wtferr(subname, version, status, & ' reading DETCHANS') ierr = 1 goto 482 endif imaxchan = maxchan ! Calc the DMA requirements idma = (4 * imaxne + imaxne * imaxchan + 2 * imaxne * maxgrp + & 3 * imaxchan) * 4 IF (.NOT.qscalar) idma = idma + 12 * imxarf write(message, '(a,i12,a)')& 'DMA requirements: ', idma, ' bytes' call wtinfo(chatter, 20, 1, message) ! ***** ! Allocate dynamic memory ! -- Alex -- added memory allocatio for order array igot = 0 allocate(p_ngrp(imaxne), stat = status) if(status.NE.0) goto 645 igot = igot + imaxne * 4 allocate(p_matrix(imaxelt), stat = status) if(status.NE.0) goto 645 igot = igot + imaxelt * 4 allocate(p_F_chan(imaxgrp), stat = status) if(status.NE.0) goto 645 igot = igot + imaxgrp * 4 allocate(p_N_chan(imaxgrp), stat = status) if(status.NE.0) goto 645 igot = igot + imaxgrp * 4 allocate(p_order(imaxgrp), stat = status) if(status.NE.0) goto 645 igot = igot + imaxgrp * 4 allocate(p_re_lo(imaxne), stat = status) if(status.NE.0) goto 645 igot = igot + imaxne * 4 allocate(p_re_hi(imaxne), stat = status) if(status.NE.0) goto 645 igot = igot + imaxne * 4 allocate(p_e_min(imaxchan), stat = status) if(status.NE.0) goto 645 igot = igot + imaxchan * 4 allocate(p_e_max(imaxchan), stat = status) if(status.NE.0) goto 645 igot = igot + imaxchan * 4 allocate(p_factor(imaxne), stat = status) if(status.NE.0) goto 645 igot = igot + imaxne * 4 allocate(p_chan(imaxchan), stat = status) if(status.NE.0) goto 645 igot = igot + imaxchan * 4 645 if(status.NE.0) then call wterrm(subname, version, & ' Failed to allocate Dynamic Memory') write(message, '(a,i12,a)')& ' an additional ', idma - igot, & ' bytes of swap space is required' call wtinfo(chatter, 1, 1, message) ierr = -1 goto 482 endif ! ***** ! Read in the RMF history keywords nk_rmfhist = 0 call fcrhky(iunit1, maxhist - 1, nk_rmfhist, rmfhist, status) nk_rmfhist = nk_rmfhist + 1 rmfhist(nk_rmfhist) = newhist ! Read in the RMF data ( Alex -- replaced rdrmf3 version by rdrmf4 ) ! Read in the RMF data ( Ziqin -- replaced rdrmf4 version by rdrmf5 ) rmf_extname = ' ' rmf_extcomm = ' ' call rdrmf5(& iunit1, chatter, rmf_extname, rmf_extcomm, & qorder, imaxne, maxgrp, maxelt, rmf_rmfversn, & dummy20, telescop, instrume, detnam, filter, area, rchantyp, & rflchan, ichan, rien, imaxgrp, numelt, p_re_lo, & p_re_hi, p_ngrp, p_F_chan, p_N_chan, & isorder, p_order, p_matrix, lo_thresh, ierr) if(ierr.NE.0) then message = ' Problem reading RSP_MATRIX xtens' call wterrm(subname, version, message) goto 482 endif ! ------------------------------------------------------------------------- ! Checks & Warnings if(qscalar) then message = ' Setting HDUCLAS3 of o/p same as ' // & ' that for i/p RSP_MATRIX' call wtinfo(chatter, 1, 1, message) outhdu3 = rmfhdu(3) message = ' (ie ' // outhdu3(:clenact(outhdu3)) // ')' call wtinfo(chatter, 1, 2, message) else if(rmfhdu(3).EQ.'FULL') then if(qdivide) then message = & ' HDUCLAS3 in o/p RSP_MATRIX set to REDIST' call wtinfo(chatter, 1, 1, message) outhdu3 = 'REDIST' else message = ' i/p RSP_MATRIX is a FULL matrix ' call wtwarm(subname, version, chatter, 1, message) message = & ' ...... Unsure what multiplication by an ARF gives' call wtinfo(chatter, 1, 2, message) message = & ' ... HDUCLAS3 in o/p RSP_MATRIX set to UNKNOWN' call wtinfo(chatter, 1, 1, message) outhdu3 = 'UNKNOWN' endif elseif(rmfhdu(3).EQ.'REDIST') then if(qdivide) then message = ' i/p RSP_MATRIX is a REDIST matrix ' call wtwarm(subname, version, chatter, 1, message) message = & 'Unsure what division by an ARF gives' call wtinfo(chatter, 1, 2, message) message = & ' HDUCLAS3 in o/p RSP_MATRIX set to UNKNOWN' call wtinfo(chatter, 1, 1, message) outhdu3 = 'UNKNOWN' else message = & ' HDUCLAS3 in o/p RSP_MATRIX set to FULL' call wtinfo(chatter, 1, 1, message) outhdu3 = 'FULL' endif elseif(rmfhdu(3).EQ.'DETECTOR') then if(qdivide) then message = ' i/p RSP_MATRIX is a DETECTOR matrix ' call wtwarm(subname, version, chatter, 1, message) message = & ' Assuming division by an ARF ' // & 'gives a REDIST matrix' call wtinfo(chatter, 1, 2, message) message = & ' HDUCLAS3 in o/p RSP_MATRIX set to REDIST' call wtinfo(chatter, 1, 1, message) outhdu3 = 'REDIST' else message = & ' HDUCLAS3 in o/p RSP_MATRIX set to FULL' call wtinfo(chatter, 1, 1, message) outhdu3 = 'FULL' endif else message = ' i/p RSP_MATRIX of unknown type' call wtwarm(subname, version, chatter, 1, message) message = 'HDUCLAS3 = ' // rmfhdu(3) call wtinfo(chatter, 1, 2, message) message = ' HDUCLAS3 in o/p RSP_MATRIX set to UNKNOWN' call wtinfo(chatter, 1, 1, message) outhdu3 = 'UNKNOWN' endif endif ! Rewind the RMF file status = 0 call ftmahd(iunit1, 1, htype, status) if(status.NE.0) then message = ' Problem rewinding RMF file' call wtferr(subname, version, status, message) endif ! -------------------------- EBOUNDS extension -------------------- ! Find the EBOUNDS extension in the EBFIL file ! - Extension number NOT given as part of ebfil (search for HDUCLAS/EXTNAM) if(ebextn.LT.0) then ninstr = 2 instr(1) = 'RESPONSE' instr(2) = 'EBOUNDS' nsearch = maxextn call fndhdu(chatter, iunit2, ninstr, instr, & nsearch, nfound, next, outhdu, outver, extnam, ierr) ! - check for old-style EXTNAME values if no acceptable HDUCLASn values found if(nfound.LE.0) then message = ' Ext w/ allowed HDUCLASn keywrds not found' call wtwarm(subname, version, chatter, 20, message) message = 'offending file: ' // ebfil call wtinfo(chatter, 20, 1, message) message = 'searching for extnsion with EXTNAME = EBOUNDS' call wtinfo(chatter, 20, 1, message) call fndext(chatter, iunit2, 'EBOUNDS', & nsearch, nfound, next, outhdu, outver, extnam, ierr) endif ! - Extension number IS given as part of ebexp else call ftmahd(iunit2, ebextn + 1, htype, status) message = ' Problem moving to specified xtens' call wtferr(subname, version, status, message) ! ... grab the HDUCLAS values for reference ninstr = 1 instr(1) = '*' nsearch = 1 call fndhdu(chatter, iunit1, ninstr, instr, & nsearch, nfound, next, outhdu, outver, extnam, ierr) nfound = 1 next(1) = 0 endif ! - sort out what we've got if(nfound.GT.1) then call wterrm(subname, version, & ' EBFIL contains >1 EBOUNDS dataset') write(message, '(i12,a)') nfound, ' extensions found:' call wtinfo(chatter, 1, 1, message) do i = 1, nfound write(message, '(a,i12,a)') ' Ext ', next(i), ':' call wtinfo(chatter, 1, 2, message) write(message, '(a,a)') 'EXTNAME = ', extnam(i) call wtinfo(chatter, 1, 2, message) do jj = 1, 4 write(message, '(a,i2,2a)')& 'HDUCLAS', jj, ' = ', outhdu(i, jj) call wtinfo(chatter, 1, 3, message) enddo enddo message = & ' ... Extension number must be specified via ebfil param' call wterrm(subname, version, message) ierr = 2 goto 482 elseif(nfound.LE.0) then message = ' Unable to locate an EBOUNDS extension' call wterrm(subname, version, message) ierr = 2 goto 482 endif ! Move to the Extension if necessary if(next(1).GT.0) then imove = next(1) status = 0 call ftmrhd(iunit2, imove, htype, status) message = ' Problem moving to EBOUNDS xtens' call wtferr(subname, version, status, message) endif ! ***** ! Read in the EBOUNDS history keywords nk_ebhist = 0 call fcrhky(iunit1, maxhist - 1, nk_ebhist, ebhist, status) nk_ebhist = nk_ebhist + 1 ebhist(nk_ebhist) = newhist ! Read in the EBOUNDS data ebd_extname = ' ' ebd_extcomm = ' ' call rdebd4(iunit2, chatter, ebd_extname, ebd_extcomm, & imaxchan, & telescop, instrume, detnam, filter, area, & echantyp, fchan, & iebound, p_chan, p_e_min, p_e_max, & ebd_rmfversn, ierr) if(ierr.NE.0) then message = ' Problem reading EBOUNDS xtens' call wterrm(subname, version, message) goto 482 endif ! Determine what the first channel number is in the EBOUNDS extension ! call fixebd(chatter, MEMI(p_chan), fchan, ierr) ! Read additional info required below ! call ftgkys(iunit1,'CHANTYPE',echantyp,comm,status) ! IF (chatter.GE.30) THEN ! message = ' reading CHANTYPE from EBOUNDS file ' ! call wtferr(subname, version, status, message) ! ENDIF ! Perform checks if(rchantyp.NE.echantyp) then message = ' CHANTYPE mis-match' call wtwarm(subname, version, chatter, 1, message) message = ' CHANTYPE from RMF extn: ' // rchantyp call wtinfo(chatter, 1, 1, message) message = ' CHANTYPE from EBOUNDS extn: ' // echantyp call wtinfo(chatter, 1, 1, message) message = ' continuing assuming you know what ' // & 'you are doing' call wtinfo(chatter, 1, 2, message) message = ' o/p file will have CHANTYPE: ' // & rchantyp call wtinfo(chatter, 1, 2, message) endif if(rflchan.ne.fchan) then message = ' Apparent discrepancy in channel numbering between'& // ' input RMF & ARF extns' call wtwarm(subname, version, chatter, 1, message) write(message, '(a,i12)') 'RMF appears to start with channel: ', & rflchan call wtinfo(chatter, 1, 2, message) write(message, '(a,i12)')& 'EBOUNDS appears to start with channel: ', & fchan call wtinfo(chatter, 1, 2, message) message = subname // version // ' unsure what to do' call wtinfo(chatter, 1, 3, message) ierr = 10 goto 482 endif ! Close the RSP_MATRIX & EBOUNDS files status = 0 call ftclos(iunit1, status) if(iunit2.NE.iunit1) call ftclos(iunit2, status) if(.NOT.qscalar) then ! -------------------------- SPECRESP extension -------------------- ! Find the SPECRESP extension in the ARF file ! - Extension number NOT given as part of arexp (search for HDUCLAS/EXTNAM) if(arextn.LT.0) then ninstr = 2 instr(1) = 'RESPONSE' instr(2) = 'SPECRESP' nsearch = maxextn call fndhdu(chatter, iunit3, ninstr, instr, & nsearch, nfound, next, outhdu, outver, extnam, ierr) ! - check for old-style EXTNAME values if no acceptable HDUCLASn values found if(nfound.LE.0) then message = ' Ext w/ allowed HDUCLASn keywrds not found' call wtwarm(subname, version, chatter, 20, message) message = ' offending file: ' // arfil call wtinfo(chatter, 20, 1, message) message = ' searching for extnsion with EXTNAME = SPECRESP' call wtinfo(chatter, 20, 1, message) call fndext(chatter, iunit3, 'SPECRESP', & nsearch, nfound, next, outhdu, outver, extnam, ierr) endif ! - Extension number IS given as part of arexp else call ftmahd(iunit3, arextn + 1, htype, status) message = 'Problem moving to specified xtens' call wtferr(subname, version, status, message) ! ... grab the HDUCLAS values for reference ninstr = 1 instr(1) = '*' nsearch = 1 ierr = 0 call fndhdu(chatter, iunit3, ninstr, instr, & nsearch, nfound, next, outhdu, outver, extnam, ierr) nfound = 1 next(1) = 0 endif ! - sort out what we've got if(nfound.GT.1) then call wterrm(subname, version, & ' ARFIL contains >1 SPECRESP datasets') write(message, '(i12,a)') nfound, ' extensions found:' call wtinfo(chatter, 1, 1, message) do i = 1, nfound write(message, '(a,i12,a)') 'Ext ', next(i), ':' call wtinfo(chatter, 1, 2, message) write(message, '(a,a)') 'EXTNAME = ', extnam(i) call wtinfo(chatter, 1, 2, message) do jj = 1, 4 write(message, '(a,i2,2a)')& 'HDUCLAS', jj, ' = ', outhdu(i, jj) call wtinfo(chatter, 1, 2, message) enddo enddo message = & ' ... Extension number must be specified via arfil param' call wterrm(subname, version, message) ierr = 2 goto 482 elseif(nfound.LE.0) then message = ' Unable to locate an SPECRESP extension' call wterrm(subname, version, message) ierr = 2 goto 482 endif ! Move to the Extension if necessary if(next(1).GT.0) then imove = next(1) status = 0 call ftmrhd(iunit3, imove, htype, status) message = ' Problem moving to SPECRESP xtens' call wtferr(subname, version, status, message) endif ! read the imxarf. CALL ftgkyj(iunit3, 'NAXIS2', imxarf, comm, status) IF (status .NE. 0) THEN call wtferr(subname, version, status, & ' reading NAXIS2 for SPECRESP') ierr = 1 goto 482 ENDIF ! allocate the p_ae_lo, p_ae_hi, p_sprsp if(.NOT.qscalar) then allocate(p_ae_lo(imxarf), stat = status) if(status.NE.0) goto 645 igot = igot + imxarf * 4 allocate(p_ae_hi(imxarf), stat = status) if(status.NE.0) goto 645 igot = igot + imxarf * 4 allocate(p_sprsp(imxarf), stat = status) if(status.NE.0) goto 645 igot = igot + imxarf * 4 endif ! Read in the SPECRESP data call rdarf1(iunit3, chatter, & atelescop, ainstrume, adetnam, afilter, & aien, p_ae_lo, p_ae_hi, p_sprsp, & arfversn, ierr) if(ierr.NE.0) then message = ' Problem reading SPECRESP xtens' call wterrm(subname, version, message) goto 482 endif ! if arfcol is not SPECRESP then try to read the column IF (arfcol .NE. 'SPECRESP') THEN CALL ftgcno(iunit3, .false., arfcol, icol, ierr) IF (ierr .NE. 0) THEN message = ' Cannot find ' // arfcol call wterrm(subname, version, message) goto 482 ENDIF CALL ftgcve(iunit3, icol, 1, 1, aien, 0.0, p_sprsp, & anyflg, ierr) IF (ierr .NE. 0) THEN message = ' Cannot read ' // arfcol call wterrm(subname, version, message) goto 482 ENDIF ENDIF ! Close the ARF file call ftclos(iunit3, status) endif ! ------- End of data collection ------------------- ! ------- Check RSP_MATRIX vs SPECRESP when FILE entered ------- if(.NOT.qscalar) then ! Check/Warn if there could be trouble ! ... Trouble with the energy grids call chkgrid(chatter, rien, aien, p_re_lo, p_re_hi, & p_ae_lo, p_ae_hi, comment, nk_comm) ! ... Trouble with the mission/instrument etc imsg = 0 if(telescop.NE.atelescop) then imsg = imsg + 1 key(imsg) = 'TELESCOP' old(imsg) = telescop new(imsg) = atelescop endif if(instrume.NE.ainstrume) then imsg = imsg + 1 key(imsg) = 'INSTRUME' old(imsg) = instrume new(imsg) = ainstrume endif if(detnam.NE.adetnam) then imsg = imsg + 1 key(imsg) = 'DETNAM' old(imsg) = detnam new(imsg) = adetnam endif if(filter.NE.afilter) then imsg = imsg + 1 key(imsg) = 'FILTER' old(imsg) = filter new(imsg) = afilter endif if(imsg.NE.0) then message = ' Following conflicts found:' call wtwarm(subname, version, chatter, 1, message) nk_comm = nk_comm + 1 comment(nk_comm) = message do i = 1, imsg message = key(i) // ' from RMF = ' // old(i) call wtinfo(chatter, 1, 1, message) nk_comm = nk_comm + 1 comment(nk_comm) = message message = key(i) // ' from ARF = ' // new(i) call wtinfo(chatter, 1, 1, message) nk_comm = nk_comm + 1 comment(nk_comm) = message enddo if(.NOT.qoverride) then message = 'values from RMF will be used in o/p' call wtinfo(chatter, 1, 2, message) nk_comm = nk_comm + 1 comment(nk_comm) = message endif endif ! Remap the SPECRESP dataset onto the energy grid supplied via RSP_MATRIX call remap(chatter, aien, p_ae_lo, p_ae_hi, p_sprsp, & rien, p_re_lo, p_re_hi, p_factor, ierr) if(ierr.GT.0) then message = ' Problem with the Remapping' call wterrm(subname, version, message) goto 482 elseif(ierr.LT.0) then message = ' Carrying on regardless ' call wtwarm(subname, version, chatter, 1, message) endif endif ! ------- Fix Variables & Pop SPECRESP array when SCALAR entered ------- ierr = 0 if(qscalar) then aien = rien call scalfill(rien, scalar, p_factor, ierr) atelescop = telescop ainstrume = instrume adetnam = detnam afilter = filter if(qdivide) then lo_thresh = lo_thresh / scalar else lo_thresh = lo_thresh * scalar endif endif ! --------- Do The Maths ---------------------------------------------- ! Perfom the multiplication (RSP_MATRIX * SPECRESP) ! Alex - modified the calls of divmatrix and multmatrix subroutines if(qdivide) then call divmatrix(imaxelt, imaxgrp, imaxne, p_ngrp, & p_N_chan, p_matrix, p_factor) else call multmatrix(imaxelt, imaxgrp, imaxne, p_ngrp, & p_N_chan, p_matrix, p_factor) endif ! Override if necessary, informing the debugging user of conflicts if(qoverride) then imsg = 0 if(telescop.NE.old_tele) then imsg = imsg + 1 key(imsg) = 'TELESCOP' old(imsg) = telescop new(imsg) = old_tele telescop = old_tele endif if(instrume.NE.old_inst) then imsg = imsg + 1 key(imsg) = 'INSTRUME' old(imsg) = instrume new(imsg) = old_inst instrume = old_inst endif if(detnam.NE.old_detn) then imsg = imsg + 1 key(imsg) = 'DETNAM' old(imsg) = detnam new(imsg) = old_detn detnam = old_detn endif if(filter.NE.old_filt) then imsg = imsg + 1 key(imsg) = 'FILTER' old(imsg) = filter new(imsg) = old_filt filter = old_filt endif if(chatter.GE.20) then if(imsg.NE.0) then message = ' Following conflicts found:' call wtwarm(subname, version, chatter, 1, message) nk_comm = nk_comm + 1 comment(nk_comm) = message do i = 1, imsg message = key(i) // ' from RSP = ' // old(i) call wtinfo(chatter, 1, 1, message) nk_comm = nk_comm + 1 comment(nk_comm) = message message = ' overwritten with ' // new(i) call wtinfo(chatter, 1, 2, message) nk_comm = nk_comm + 1 comment(nk_comm) = message enddo endif endif endif ! Open & Write the FITS file ! ------------------------ PRIMARY HEADER --------------------------- if(rmfversn(1:1).EQ.'1')then ! ... Open the FITS file and write a null primary header call opnpa(outfil, chatter, ounit, killit, status) if(status.ne.0) then ierr = 5 goto 482 endif ! ... Add additional keywords to Primary Header call FTPKYS(ounit, 'CREATOR', & taskname, & 's/w task which wrote this dataset', & status) call FTPKYS(ounit, 'CONTENT', & 'RESPONSE MATRIX', & 'RMF & EBOUNDS xtensions', & status) ! call FTPKYS(ounit,'RMFVERSN', ! & rmfversn , ! & 'OGIP classification of FITS format style', ! & status) do 132 nk = 1, nk_prhist call FTPHIS(ounit, prhist(nk), status) 132 continue else message = 'Unknown format: ' // rmfversn call wterrm(subname, version, message) ierr = 1 goto 482 endif ! ------------------------ finished PRIMARY --------------------------- ! ------------------------ RMF EXTENSION ---------------------------- IF (ebfil.EQ.rmfil) THEN IF(rmextn.gt. ebextn) THEN goto 504 ENDIF ENDIF 502 continue if(rmfversn(1:1).EQ.'1')then ! ... Write the RMF extension within RMF file ! -- Alex: replaced the wtrmf3 version by wtrmf4 ! -- Ziqin: replaced the wtrmf4 version by wtrmf5 call wtrmf5(ounit, chatter, rmf_extname, rmf_extcomm, & nk_rmfhist, rmfhist, & nk_comm, comment, rmfversn, outhdu3, telescop, instrume, & detnam, filter, area, rchantyp, rflchan, numelt, ichan, rien, & imaxgrp, p_re_lo, p_re_hi, p_ngrp, & p_F_chan, p_N_chan, qorder, p_order, & p_matrix, lo_thresh, status) if(status.NE.0) goto 876 call FTPKYS(ounit, 'CREATOR', & taskname, & 's/w task which wrote this dataset', & status) else message = 'Unknown format: ' // rmfversn call wterrm(subname, version, message) ierr = 1 goto 482 endif ! ------------------------ finished RMF EXTENSION ---------------------------- IF (ebfil.EQ.rmfil) THEN IF(rmextn.gt. ebextn) THEN goto 505 ENDIF ENDIF 876 continue 504 continue ! ----------------------------- EBOUNDS EXTENSION ---------------------------- if(rmfversn(1:1).EQ.'1')then ! ... Write the EBOUNDS extension within the RMF file call wtebd4(ounit, chatter, ebd_extname, ebd_extcomm, & nk_ebhist, ebhist, & nk_comm, comment, rmfversn, & telescop, instrume, detnam, filter, area, & echantyp, fchan, & iebound, p_e_min, p_e_max, status) call FTPKYS(ounit, 'CREATOR', & taskname, & 's/w task which wrote this dataset', & status) call ftpkys(ounit, 'CHANTYPE', & rchantyp, & 'Detector Channel Type in use (PHA or PI)', & status) call wtferr(subname, version, status, & ' Problem writing CHANTYPE keyword') else message = 'Unknown format: ' // rmfversn call wterrm(subname, version, message) ierr = 1 goto 482 endif IF (ebfil.EQ.rmfil) THEN IF(rmextn.gt. ebextn) THEN goto 502 ENDIF ENDIF 505 continue ! ------------------------ finished EBOUNDS EXTENSION ---------------------- ! Close the FITS file status = 0 call ftclos(ounit, status) if(status.ne.0) then call wtferr(subname, version, status, & ' Problem closing o/p file') ierr = 6 goto 482 endif ! Check for errors 482 if(ierr.ne.0) then call wterrm(subname, version, ' Fatal - aborting') endif ! Free the allocated memory: deallocate(p_ngrp, stat = status) if(status.NE.0) goto 646 deallocate(p_matrix, stat = status) if(status.NE.0) goto 646 deallocate(p_F_chan, stat = status) if(status.NE.0) goto 646 deallocate(p_N_chan, stat = status) if(status.NE.0) goto 646 deallocate(p_order, stat = status) if(status.NE.0) goto 646 deallocate(p_re_lo, stat = status) if(status.NE.0) goto 646 deallocate(p_re_hi, stat = status) if(status.NE.0) goto 646 deallocate(p_e_min, stat = status) if(status.NE.0) goto 646 deallocate(p_e_max, stat = status) if(status.NE.0) goto 646 deallocate(p_factor, stat = status) if(status.NE.0) goto 646 deallocate(p_chan, stat = status) if(status.NE.0) goto 646 if(.not.qscalar) then deallocate(p_ae_lo, stat = status) if(status.NE.0) goto 646 deallocate(p_ae_hi, stat = status) if(status.NE.0) goto 646 deallocate(p_sprsp, stat = status) if(status.NE.0) goto 646 endif 646 if(status.NE.0) then message = ' Failed to deallocate Dynamic Memory' call wtwarm(subname, version, chatter, 1, message) message = ' (only a potential problem ' // & 'if your NOT running the host version)' call wtinfo(chatter, 1, 1, message) endif ! -------------------------------------------------------------------- return end ! ----------------------------------------------------------- !+ MULTMATRIX ! - old: subroutine multmatrix(maxchan,rien,ichan,matrix,factor) ! - old: integer rien, ichan, maxchan, numelt ! - old: real matrix(maxchan, *), factor(*) SUBROUTINE MULTMATRIX(IMAXELT, IMAXGRP, IMAXEN, NGRP, & NCHAN, MATRIX, FACTOR) IMPLICIT NONE INTEGER IMAXELT, IMAXGRP, IMAXEN, NGRP(IMAXEN), NCHAN(IMAXGRP) REAL MATRIX(IMAXELT), FACTOR(IMAXEN) ! Description ! Multiplies the matrix by the energy-dependent vector factor ! Required to make use of dynamic memory alloaction routines ! ! Author/Modification history ! Ian M George ! Alex modified the matrix multiplication to take into account ! that the matrix is now 1-D array !- ! Internals INTEGER I, J, IGRP, IMTRX, ICH !INTEGER jj IGRP = 0 IMTRX = 0 DO I = 1, IMAXEN DO J = 1, NGRP(I) IGRP = IGRP + 1 DO ICH = 1, NCHAN(IGRP) IMTRX = IMTRX + 1 if(matrix(imtrx).lt.1.0e-20) then MATRIX(IMTRX) = 0.0 else MATRIX(IMTRX) = MATRIX(IMTRX) * FACTOR(I) endif ENDDO ENDDO ENDDO !C ------------------------- old version --------------------- !C do i = 1, rien !C do jj = 1, ichan !C matrix(jj,i) = matrix(jj,i) * factor(i) !C enddo !C enddo !C ----------------------------------------------------------- return end ! ----------------------------------------------------------- !+ DIVMATRIXC ! - old: subroutine divmatrix(maxchan, rien, ichan, matrix, factor) ! - old: integer rien, ichan, maxchan ! - old: real matrix(maxchan, *), factor(*) SUBROUTINE DIVMATRIX(IMAXELT, IMAXGRP, IMAXEN, NGRP, & NCHAN, MATRIX, FACTOR) IMPLICIT NONE INTEGER IMAXELT, IMAXGRP, IMAXEN, NGRP(IMAXEN), NCHAN(IMAXGRP) REAL MATRIX(IMAXELT), FACTOR(IMAXEN) ! Description ! Divides the matrix by the energy-dependent vector factor ! Required to make use of dynamic memory alloaction routines ! ! Author/Modification history ! Ian M George ! Alex modified the matrix division to take into account ! that the matrix is now 1-D array !- ! Internals INTEGER I, J, IGRP, IMTRX, ICH IGRP = 0 IMTRX = 0 DO I = 1, IMAXEN DO J = 1, NGRP(I) IGRP = IGRP + 1 DO ICH = 1, NCHAN(IGRP) IMTRX = IMTRX + 1 MATRIX(IMTRX) = MATRIX(IMTRX) / FACTOR(I) ENDDO ENDDO ENDDO ! ----------------- old version --------------------------- ! do i = 1, rien ! do jj = 1, ichan ! matrix(jj,i) = matrix(jj,i) / factor(i) ! enddo ! enddo ! ---------------------------------------------------------- return end ! ----------------------------------------------------------- !+ CHKGRID subroutine chkgrid(chatter, rien, aien, re_lo, re_hi, & ae_lo, ae_hi, comment, nk_comm) IMPLICIT NONE integer rien, aien, nk_comm integer chatter real re_lo(rien), re_hi(rien) real ae_lo(aien), ae_hi(aien) character*(*) comment(*) ! Description ! Checks the energy grid of the RMF & ARF files ! Required to make use of dynamic memory alloaction routines ! ! Author/Modification history ! Ian M George character(7) version parameter (version = '1.0.0') !- ! Internals character(10) subname parameter (subname = 'chkgrid') character(80) message if(rien.NE.aien .OR. re_lo(1).NE.ae_lo(1) .OR.& re_hi(rien).NE.ae_hi(aien)) then message = ' RMF/ARF energy grid mismatch' nk_comm = nk_comm + 1 comment(nk_comm) = message call wtwarm(subname, version, chatter, 1, message) write(message, '(a,i12,a,f10.5,a,f10.5,a)') 'RMF: ', & rien, ' bins between ', re_lo(1), ' & ', re_hi(rien),' keV' call wtinfo(chatter, 1, 1, message) nk_comm = nk_comm + 1 comment(nk_comm) = message write(message, '(a,i12,a,f10.5,a,f10.5)') 'ARF: ', & aien, ' bins between ', ae_lo(1), ' & ', ae_hi(aien) call wtinfo(chatter, 1, 1, message) nk_comm = nk_comm + 1 comment(nk_comm) = message message = 'ARF dataset will be remapped onto RMF grid' call wtinfo(chatter, 1, 1, message) nk_comm = nk_comm + 1 comment(nk_comm) = message endif return end ! ----------------------------------------------------------- !+ SCALFILL subroutine scalfill(rien, scalar, factor, ierr) IMPLICIT NONE integer rien, ierr real scalar real factor(*) ! Description ! Fills the factor array with the value of scalar ! ! Author/Modification history ! Ian M George (1.0.0: 1995 Apr 09) original !- ! Internals integer i do i = 1, rien factor(i) = scalar enddo ierr = 0 return end ! -----------------------------------------------------------