!+RBNRPSF subroutine rbnrpf ! --------------------- DESCRIPTION -------------------------------- ! ! Main task routine for rebinning OGIP standard RPSF file with a specified ! amount of minimum counts in each bin ! Input : infile (RPSF format) ! outfile ! cmin (Minimum counts per pixel) ! bkgd (% default, any other value overides this) If % bkgd in ! infile is used. ! ! --- VARIABLES ------------------------------------------------------- IMPLICIT NONE character(180) infile, outfile, subinfo character(8) telrad, instrad character(16) radunit, thetaunit, energunit, rpsfunit integer npts, chatter, ierr integer maxrad, maxtheta, maxen, chanmin, chanmax real pix_size, c_min, bkgd, bkgd_rad, sumtcts, sumrcts real, allocatable :: p_rad_lo(:), p_rad_hi(:) real, allocatable :: p_cts_arc2(:), p_err_arc2(:), p_area_wgt(:) real, allocatable :: p_energ_lo(:), p_energ_hi(:) real, allocatable :: p_theta_lo(:), p_theta_hi(:) integer ntheta, nenerg integer iget, ineed, status logical qerror, qarea, calc, lbin, negpres, rpsf_bkgd, killit ! -------------------- VARIABLE DIRECTORY ------------------------ ! ! infile char : Name of radial profile file (input) ! outfile char : Name of results file (output) ! ierr int : Error flag, ierr = 0 okay ! telrad char : Telescope name for radial data (from infile) ! instrad char : Instrument name for radial data (from infile) ! maxrad int : Maximum size of arrays used for radial profile data ! maxen int : Maximum size of arrays used for pha data ! pix_size real : pixel size (from rpsf infile), converted to arcmins/pixel ! c_min real : Minimum No. counts/bin (user defined) ! bkgd real : Background count rate in counts/pixel (user defined) ! rad_lo real : Array of lower edge of radial bins ( in arcmins) ! rad_hi real : Array of upper edge of radial bins (in arcmins) ! theta_lo real : Array of theta low values of the observation ! theta_hi real : Array of theta high values of the observation ! energ_lo real : Array of energy low values of the observation ! energ_hi real : Array of energy high values of the observation ! ntheta int : No of theta values ! nenerg int : No of energy values ! qerror logical: True if errors present in rad profile file ! qarea logical: True if Area_weighting factors present ! area_wgt real : Array of area weighting factors ! cts_arc2 real : Array of radial profile in counts/arcmin^2 ! err_arc2 real : Array of errors on cts_arc2 ! npts int : No. of observed psf data points ! lbin logical: true if last bin has less than c_min counts ! negpres logical: true if negative counts present in dataset ! ! ---------------- COMPILATION AND LINKING ---------------- ! ! Link with XPI and with FITSIO ! ! ---------------- CALLED ROUTINES ---------------- ! ! subroutine CPSF_GP : Obtains user defined input ! subroutine CPSF_RINFILE : Reads data from radial profile ! subroutine CPSF_RESULTS : Writes rebinned Observed data ! ! ---------------- AUTHORS/MODIFICATION HISTORY ---------------- ! ! Rehana Yusaf (1993 January 15) : Modularisation and minor modifications ! Rehana Yusaf (1993 Febuary 24) : Changing code to read FITS format radial ! profile input ! ! Rehana Yusaf (1993 March 3) : Minor changes, more info in output ! Rehana Yusaf (1993 May 26) : Minor changes, removing stops etc ! Rehana Yusaf (1993 June 8) : Minor changes, lbin and negpres added ! ! Rehana Yusaf (1993 June 9) 1.0.1: Substitute off-axis function for on-axis ! Rehana Yusaf (1993 July 27) 1.0.2; Minor change to Off-Axis algorithm, ! (comment in function) ! Rehana Yusaf (1993 August 2) 1.0.3; Minor change to OFF-AXIS alg' ! (comment in function) ! Rehana Yusaf (1993 Sept 13 1993) 1.0.4; RPSFVER 1993a format is used, ! instead of 1992a this task is changed ! accordingly. ! Rehana Yusaf (1994 Jan 20) 1.1.0; Stripped out code from PSPCRPSF as a ! basis for this 'new' task that only ! rebins the data, PSPCRPSF now takes as ! input the rebinned RPSF file ! Rehana Yusaf (1994 Feb 3) 1.1.1; rpsf_bkgd is used correctly - that is ! if a user-defined bkgd is entered then ! use that value NOT the value read from ! file ! Rehana Yusaf (1995 Jan 13) 1.1.2; update _gp so that defval is no longer ! read from parfile ! Rehana Yusaf (1995 Feb 28) 1.1.3; Add more diagnostics in _rebin ! bkgd can be calculated ! Rehana Yusaf (1995 April 25) 1.1.4; add clobber and update to _rebin ! ! Banashree Mitra Seifert (1996, Jan) 2.0.0: Modifications made ! . Introduced DMA ! . Replaced by call to MVEXT ! . Replaced filenames by character(180) ! . maxtheta=1, maxen=1 ! . Introduced screen display routines ! wtinfo,wtbegm,wtendm,wtferr ! ! Banashree Mitra Seifert (1996, Feb) 2.1.0: ! . modification done in search for PSF file. ! earlier version (2.0.0) didn't look for ! instr(3) = TOTAL. Instead it looked for ! only instr(3) = NET. ! . introduced EXTNAME before call to MVEXT ! ! Banashree Mitra Seifert (1997, Oct 10) 2.2.0: ! . In the subroutine cpsf_rebin, ! The calculation of area_wgt was wrong. binning was ok, but ! while running pcrpsf/hrirpsf after binning, resulted in a ! mismatch between the psf model normalization and the data ! further down the data-analysis pipeline. This error happened ! when binning requires counts from more than one bin from ! unbinned input data. So replaced ! area_wgt(j,1,1) = sumarea/area ! by ! area_wgt(j,1,1) = sumarea/k ! ! Peter D Wilson (1998 Jul 01) 2.2.1: ! . Updated for new FCPARS behavior ! ! Bryan K Irby (2018 Jan 14) 2.3.0: ! . Replaced udmget/udmfre with allocate/deallocate ! MFC (2020 Apr 16) 2.3.1 ! . f90 version ! ! ------------------------------------------------------------------------ character(5) version parameter (version = '2.3.1') character(8) subname character(40) taskname !cc COMMON /task/taskname subname = 'rbnrpsf' taskname = 'rbnrpsf' subinfo = 'using' // subname // version call wtinfo(chatter, 10, 1, subinfo) ! ----------------------------------------------------------------------- ierr = 0 rpsf_bkgd = .false. call cpsf_gp(infile, outfile, c_min, bkgd, rpsf_bkgd, calc, & bkgd_rad, ierr, chatter, killit) call wtbegm(taskname, version, chatter) IF (ierr.NE.0) THEN go to 100 ENDIF maxrad = 1000 maxtheta = 5 maxen = 50 ! ----------------------- Allocation of DMA ---------------------- ! iget = bytes get added after each call for UDMGET ! (this is the actual count of bytes I am asking for) ! just to keep a count on how much memory is asking for ! ---------------------------------------------------------------- iget = 0 status = 0 allocate(p_rad_lo(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * 4 status = 0 allocate(p_rad_hi(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * 4 status = 0 allocate(p_cts_arc2(maxrad * maxtheta * maxen), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * maxen * 4 status = 0 allocate(p_err_arc2(maxrad * maxtheta * maxen), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * maxen * 4 status = 0 allocate(p_area_wgt(maxrad * maxtheta * maxen), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * maxen * 4 status = 0 allocate(p_energ_lo(maxen), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxen * 4 status = 0 allocate(p_energ_hi(maxen), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxen * 4 status = 0 if(maxtheta .lt. 50)maxtheta = 50 allocate(p_theta_lo(maxtheta), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxtheta * 4 status = 0 allocate(p_theta_hi(maxtheta), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxtheta * 4 50 ineed = 2 * maxrad * 4 + 3 * maxrad * 5 * maxen * 4 + & 2 * maxtheta * 4 + 2 * maxen * 4 write(subinfo, '(a,i10)')'DMAsize required for this task=', ineed call wtinfo(chatter, 10, 1, subinfo) write(subinfo, '(a,i10)')'total bytes of memory obtained=', iget call wtinfo(chatter, 10, 1, subinfo) if (status .ne. 0) then ierr = -1 subinfo = 'failed to allocate dynamic memory ' call wtferr(subname, version, status, subinfo) goto 100 endif ! ---------------------------------------------------------------------+ ! end of DMA allocation ! ---------------------------------------------------------------------+ call cpsf_rinfile(infile, pix_size, c_min, npts, p_rad_lo, & p_rad_hi, radunit, ntheta, & p_theta_lo, p_theta_hi, thetaunit, & nenerg, p_energ_lo, p_energ_hi, & energunit, p_cts_arc2, qerror, & p_err_arc2, rpsfunit, qarea, & p_area_wgt, maxrad, & maxtheta, telrad, instrad, negpres, lbin, & rpsf_bkgd, bkgd, chanmin, chanmax, calc, & bkgd_rad, sumrcts, sumtcts, ierr, chatter) IF (ierr.NE.0) THEN goto 100 ENDIF ! ------------ calling cpsf_results -----------------------------------+ call cpsf_results(outfile, infile, version, npts, & p_rad_lo, p_rad_hi, radunit, & ntheta, p_theta_lo, p_theta_hi, & thetaunit, nenerg, p_energ_lo, & p_energ_hi, energunit, p_cts_arc2, & qerror, p_err_arc2, rpsfunit, qarea, & p_area_wgt, pix_size, maxrad, maxtheta, & bkgd, c_min, telrad, instrad, negpres, lbin, & chanmin, chanmax, sumrcts, sumtcts, & ierr, chatter, killit) if (ierr .ne. 0) then goto 100 endif ! --------- free the dynamic memory ----------------------------- status = 0 deallocate(p_rad_lo, stat = status) status = 0 deallocate(p_rad_hi, stat = status) status = 0 deallocate(p_cts_arc2, stat = status) status = 0 deallocate(p_err_arc2, stat = status) status = 0 deallocate(p_energ_lo, stat = status) status = 0 deallocate(p_energ_hi, stat = status) status = 0 deallocate(p_theta_lo, stat = status) status = 0 deallocate(p_theta_hi, stat = status) status = 0 deallocate(p_area_wgt, stat = status) if (status .ne. 0) then subinfo = ' failed to de-allocate memory ' call wtferr(subname, version, status, subinfo) ierr = 99 endif 100 call wtendm(taskname, version, ierr, chatter) return end ! -------------------------------------------------------------------- ! END OF RBNRPSF ! -------------------------------------------------------------------- !+CPSF_GP subroutine cpsf_gp(infile, outfile, c_min, bkgd, rpsf_bkgd, calc, & bkgd_rad, errflg, chatter, killit) ! --- DESCRIPTION ---------------------------------------------------- ! ! This routine obtains user defined filenames and parameters. ! ! ---------------------- VARIABLES ----------------------------------- IMPLICIT NONE character*(*) infile, outfile character(180) strbkgd, ill_files(5) real c_min, bkgd, bkgd_rad integer chatter, errflg, n_ill logical rpsf_bkgd, calc, killit, valfil ! ------------------ INTERNAL VARIABLES ------------------------------ !character(180) filename character(100) subinfo character(80) extinfo integer ierr, try !integer extnum, try !logical ext ! -------------------- VARIABLE DIRECTORY ---------------------------- ! ! infile char : Name of Radial profile file ! outfile char : Name of results file ! c_min real : Minimum No. counts/bin ! bkgd real : Background count rate in counts/pixel ! chatter int : Chattines flag (<5 quiet,>5normal,>20 noisy) ! ierr int : Error flag ! errflg int : Error flag which returns to main, 0 ok ! errstr char : Error string ! ! ----------------- COMPILATION AND LINKING ------------------------- ! ! Link with FTOOLS ! ! ----------------- AUTHORS/MODIFICATION HISTORY ----------------- ! ! Rehana Yusaf (1993 January) ! Rehana Yusaf (1993 Febuary):Change code to use fcecho for screen output ! Rehana Yusaf (1993 June 9) : Temporary measure, off_axis angle is ! prompted for. ! Rehana Yusaf (1993 Sept 23) : Off Axis angle only prompted for if ! detector file not entered ! Rehana Yusaf (1994 Jan 20) 1.0.3; Adapt code for reading parameters ! for RBNRPSF (previously used for PSPCRPSF) ! Rehana Yusaf (1995 Jan 13) 1.0.4; defval no longer read from parfile ! Rehana Yusaf (1995 Feb 28) 1.0.5; bkgd can be calculated ! Rehana Yusaf (1995 April 25) 1.0.6; read in clobber ! Banashree Mitra Seifert (1996, Jan) 1.1.0; ! . Redefined the filenames to be character(180) ! . Dimensions are carried from the calling subroutine ! . Introduced screen display routines ! wtinfo,wterrm ! Peter D Wilson (1998 Jul 01) 1.1.1: ! . Drop INQUIRE test ! ----------------------------------------------------------------- character(5) version parameter (version = '1.1.1') character(12) subname !- ! ----------------- OBTAINING I/O FILENAMES ----------------------- subname = 'cpsf_getpar ' subinfo = 'using' // subname // version call wtinfo(chatter, 10, 1, subinfo) ! ------------------------------------------------------------------- call uclgst('infile', infile, ierr) IF (ierr.NE.0) THEN extinfo = 'getting infile parameter ' call wterrm(subname, version, extinfo) ENDIF call crmvlbk(infile) IF (infile.EQ.' ') THEN extinfo = 'infile must be entered !' call wterrm(subname, version, extinfo) errflg = 1 return ENDIF ! PDW 7/1/98: Don't bother! Let FTOPEN determine if file exists ! call fcpars(infile,filename,extnum,ierr) ! INQUIRE(FILE=filename,EXIST=ext) ! IF (.NOT.ext) THEN ! extinfo = 'infile does not exist : '//filename ! call wterrm(subname,version,extinfo) ! errflg = 1 ! return ! ENDIF ! ---------------- GET OUTFILE NAME ---------------------------- call uclgst('outfile', outfile, ierr) IF (ierr.NE.0) THEN extinfo = 'getting outfile parameter' call wterrm(subname, version, extinfo) errflg = 1 return ENDIF call crmvlbk(outfile) IF (outfile.EQ.' ') THEN extinfo = 'must enter outfile name !!' call wterrm(subname, version, extinfo) errflg = 1 return ENDIF ! CLOBBER PARAMETER ierr = 0 call uclgsb('clobber', killit, ierr) IF (ierr.NE.0) THEN extinfo = 'getting killit parameter !' call wterrm(subname, version, extinfo) errflg = 1 return ENDIF n_ill = 0 call ck_file(outfile, ill_files, n_ill, valfil, killit, chatter) IF (.NOT.valfil) THEN errflg = 2 return ENDIF ! ------------ OBTAINING USER DEFINED PARAMETERS ----------------------+ call uclgsi('chatter', chatter, ierr) IF (ierr.NE.0) THEN extinfo = 'getting chatter parameter' call wterrm(subname, version, extinfo) ENDIF call uclgsr('c_min', c_min, ierr) IF (ierr.NE.0) THEN extinfo = 'getting minimum no. counts/bin,c_min parameter' call wterrm(subname, version, extinfo) ENDIF calc = .false. try = 0 90 call uclgst('bkgd', strbkgd, ierr) try = try + 1 IF (ierr.NE.0) THEN extinfo = 'getting background count rate, bkgd parameter' call wterrm(subname, version, extinfo) ENDIF ierr = 0 call uclpst('bkgd', '%', ierr) IF (ierr.NE.0) THEN extinfo = ' putting default bkgd parameter' call wterrm(subname, version, extinfo) ENDIF call crmvlbk(strbkgd) call ftupch(strbkgd) IF (strbkgd(1:1).EQ.'%') THEN rpsf_bkgd = .true. ELSEIF (strbkgd(1:1).EQ.'C') THEN calc = .true. call uclgsr('bkgd_rad', bkgd_rad, ierr) IF (ierr.NE.0) THEN extinfo = 'getting background radius parameter' call wterrm(subname, version, extinfo) ENDIF ELSE read(strbkgd, *, IOSTAT = ierr) bkgd If (ierr.NE.0) THEN extinfo = ' error in bkgd value ' call wterrm(subname, version, extinfo) errflg = 1 return ENDIF rpsf_bkgd = .false. IF (bkgd.LT.0) THEN IF (try.LT.3) THEN extinfo = 'negative bkgd entered, try again' call wterrm(subname, version, extinfo) goto 90 ELSE errflg = 1 return ENDIF ENDIF ENDIF return end ! -------------------------------------------------------------------- ! END OF SUBROUTINE CPSF_GP ! -------------------------------------------------------------------- !+CPSF_RINFILE ! ---------------------------------------------------------------------+ subroutine cpsf_rinfile(infile, pix_size, c_min, npts, & rad_lo, rad_hi, radunit, ntheta, theta_lo, & theta_hi, thetaunit, nenerg, energ_lo, & energ_hi, energunit, cts_arc2, qerror, & err_arc2, rpsfunit, qarea, area_wgt, maxrad, & maxtheta, telrad, instrad, negpres, lbin, & rpsf_bkgd, bkgd, chanmin, chanmax, calc, & bkgd_rad, sumrcts, sumtcts, ierr, chatter) ! --- DESCRIPTION --------------------------------------------------------- ! ! Reads in Radial profile data from user defined input file. The input file ! is in FITS format. ! ! --- VARIABLES ----------------------------------------------------------- IMPLICIT NONE character*(*) infile character*(*) telrad, instrad character(16) radunit, thetaunit, energunit, rpsfunit integer maxrad, ierr, maxtheta integer npts, chatter, ntheta, nenerg, chanmin, chanmax real pix_size, c_min, sumrcts, sumtcts, bkgd, bkgd_rad real rad_lo(*), rad_hi(*), cts_arc2(maxrad, maxtheta, 1) real err_arc2(maxrad, maxtheta, 1), area_wgt(maxrad, maxtheta, 1) real energ_lo(*), energ_hi(*), theta_lo(*), theta_hi(*) logical negpres, lbin, qerror, qarea, rpsf_bkgd, calc ! ----------------- INTERNAL VARIABLES ---------------------------------- character(180) subinfo, errinfo character(40) comm character(20) instr(4), outhdu(9, 50), extnames(50), outver(9, 50) character(20) extname character(16) hduclas3 integer maxin, iunit, status, ninstr, nsearch, next(50) integer extnum parameter (maxin = 1000) real incts(maxin, 1, 1), inrad_lo(maxin), inrad_hi(maxin) real inerr(maxin, 1, 1), pix ! ----------------- VARIABLE DIRECTORY -------------------------------- ! ! Arguments ... ! ! infile char : Name of radial profile data file (user defined) ! pix real : Pixelsize in degrees (from infile) ! pix_size real : pixel size in arcmins per pixel ! c_min real : Minimum No. counts/bin (user defined) ! maxrad int : Maximum size of arrays for radial profile data ! rad_lo real : Array of lower edge of radial bins ( in arcmins) ! rad_hi real : Array of upper edge of radial bins (in arcmins) ! theta_lo real : Array of lower theta values ! theta_hi real ; Array of upper theta values ! energ_lo real : Lower Energy bounds ! energ_hi real : Upper Energy bounds ! area_wgt real : Array of area weighting factors ! cts_arc2 real : Array of radial profile in counts/arcmin^2 ! err_arc2 real : Array of errors on cts_arc2 ! npts int : Counter for No. of psf data points ! telrad char : Telescope name ! instrad char : Instrument name ! negpres logical: true if negative counts present in data ! lbin logical: true if last bin < c_min counts ! ! ! maxin int : Array dimensions, this should be same as maxrad ! incts real : Array of radial profile in counts/arcmin^2, ! used to read from file ! inerr int : Array of errors on incts, used to read from file ! inrad_lo real : Array of minimum radial profile in arcmin ! used to read from file ! inrad_hi real : Array of max radial profile in arcmin ! subinfo char : User information ! errstr char : Error string for this routine ! iunit int : Fortran i/o unit ! data_pres logical: true if observational data is present ! theo_pres logical: true if predicted data is present ! ! ---------------------- CALLED ROUTINES ---------------------------- ! ! subroutine RD_RPSF1993a : (CALLIB) Reads OGIP FITS format radial profile ! subroutine CPSF_REBIN : Rebins radial profile data if necessary. That ! is, each bin is checked ,and if it is less than ! the Minimum No. counts/bin defined by the user ! then the data is rebinned. ! ! -------------------- AUTHORS/MODIFICATION HISTORY -------------------- ! ! Rehana Yusaf (1993 January) ! Rehana Yusaf (1993 Febuary) : Change code to read FITS format RPSF data ! Rehana Yusaf (1993 Sept 13) 1.0.1; Change code to read FITS RPSFVER 1993a ! (previous read was 1992a) ! Rehana Yusaf (1994 Jan 20) 1.0.2; RD_RPSF1993a has been renamed RDRPF1, ! also bkgd is read from RPSF file if the user ! does not define a value. bkgd is passed to the ! cpsf_rebin routine. ! Rehana Yusaf (1995 Mar 7) 1.0.3; Calculate SUMRCTS and SUMTCTS in _rebin ! Banashree Mitra Seifert (December 1995) 2.0.0: modifications are ! . In the variable declaration so that the dimensions ! are carried over from the parent subroutine. ! . Replaced by call to MVEXT ! . Modified some error messages ! . Introduced screen display routines ! wtinfo,wtferr,wterrm ! -------------------------------------------------------------------------- character(5) version parameter (version = '2.0.0') character(14) subname integer rwmode !- ! -------------------- USER INFO ------------------------------ subname = 'cpsf_rinfile' subinfo = 'using cpsf_rinfile Ver ' // version call wtinfo(chatter, 10, 1, subinfo) ! ---------------------- FIND RPSF EXTENSION --------------------------+ ierr = 0 ninstr = 3 instr(1) = 'RESPONSE' instr(2) = 'RPRF' instr(3) = 'NET' nsearch = 50 rwmode = 0 call mvext (rwmode, infile, iunit, ninstr, instr, nsearch, next, & outhdu, extnames, outver, extname, ierr, chatter) if (ierr .ne. 0) then ierr = 0 ninstr = 3 instr(1) = 'RESPONSE' instr(2) = 'RPRF' instr(3) = 'TOTAL' nsearch = 50 call mver (iunit, extnum, ninstr, instr, nsearch, next, & outhdu, extnames, outver, extname, ierr, chatter) if (ierr .ne. 0) then subinfo = 'error in searching for observed PSF extension' call wterrm(subname, version, subinfo) subinfo = 'RPSF observed data is not present' call wterrm(subname, version, subinfo) call ftclos(iunit, status) subinfo = ' closing PSF obs file' call wtferr(subname, version, status, subinfo) return else subinfo = ' extension with observed PSF data found' call wtinfo(chatter, 10, 3, subinfo) endif else subinfo = ' extension with observed PSF data found' call wtinfo(chatter, 10, 3, subinfo) endif ! --------- STARTS if block for READ OBSERVED DATA --------------------+ call rdrpf1(iunit, hduclas3, npts, inrad_lo, inrad_hi, radunit, & ntheta, theta_lo, theta_hi, thetaunit, nenerg, & energ_lo, energ_hi, energunit, incts, qerror, inerr, & rpsfunit, qarea, area_wgt, telrad, instrad, maxrad, & maxtheta, ierr, chatter) IF (ierr.NE.0) THEN errinfo = ' reading RPSF file' call wterrm(subname, version, errinfo) return ENDIF ! ------------------- READ PIXEL SIZE---------------------------------- status = 0 call ftgkye(iunit, 'PIXSIZE', pix, comm, status) errinfo = ' reading PIXSIZE ' call wtferr(subname, version, status, errinfo) IF (pix.EQ.(0.0)) THEN errinfo = & 'The PIXSIZE keyword is set to 0.0 in radial profile' call wterrm(subname, version, errinfo) ierr = 1 return ENDIF ! --------- Converting pixsize from degrees to arcmin/pixel ------------ pix_size = pix * 60 ! ---------- Read CHANMIN and CHANMAX values if present --------------- status = 0 call ftgkyj(iunit, 'CHANMIN', chanmin, comm, status) status = 0 call ftgkyj(iunit, 'CHANMAX', chanmax, comm, status) status = 0 ! ------------ Read BACKGROUND if RPSF_BKGD is true ------------ IF (rpsf_bkgd) THEN call ftgkye(iunit, 'BACKGRND', bkgd, comm, status) errinfo = 'reading BACKGRND ' call wtferr(subname, version, status, errinfo) ENDIF status = 0 call ftclos(iunit, status) errinfo = 'closing RPSF file' call wtferr(subname, version, status, errinfo) ! ------------ CALLING REBINNING ROUTINE --------------------------- call cpsf_rebin(inrad_lo, inrad_hi, area_wgt, incts, inerr, npts, & pix_size, c_min, rad_lo, rad_hi, cts_arc2, & err_arc2, maxin, maxrad, maxtheta, negpres, lbin, & chatter, bkgd, sumrcts, sumtcts, calc, bkgd_rad, & ierr) return end ! ---------------------------------------------------------------------+ ! END OF CPSF_RINFILE ! ---------------------------------------------------------------------+ !+CPSF_REBIN subroutine cpsf_rebin(inrad_lo, inrad_hi, area_wgt, incts, inerr, & npts, pix_size, c_min, rad_lo, rad_hi, & cts_arc2, err_arc2, maxin, maxrad, maxtheta, & negpres, lbin, chatter, bkgd, sumraw, sumtcts, & calc, bkgd_rad, ierr) ! ------------------ DESCRIPTION ------------------------------------- ! ! This routine rebins the radial profile data by comparing each bin with ! the user defined Minimum counts/bin ,and rebinning if the bin is less ! than it. The arrays inrad_lo,inrad_hi,in_pix,incts,inerr have been used ! to read data from radial file. After rebinning check the arrays rad_lo ! rad_hi,area_wgt,cts_arc2 and err_arc2 are used throughout the program. ! ! --- VARIABLES --------------------------------------------------------- IMPLICIT NONE integer maxin, maxrad, chatter, maxtheta, ierr, npts real inrad_lo(*), inrad_hi(*), incts(maxin, 1, 1), inerr(maxin, 1, 1) real rad_lo(*), rad_hi(*), cts_arc2(maxrad, 1, 1) real err_arc2(maxrad, 1, 1), area_wgt(maxrad, 1, 1), bkgd, bkgd_rad real pix_size, c_min, sumraw, sumtcts logical lbin, negpres, calc ! ! -------------------- INTERNAL VARIABLES -------------------------- ! character(200) subinfo real sumrcts, rawcts, sumpix, absumr, area, sumarea, sumpix_bkgd real per_bkgd, sumcts, pi, totpix, totpix2, sumrcts_bkgd integer i, j, k, onpts, ibin, step, remain, tmaxin parameter (tmaxin = 1000) real in_pix(tmaxin) logical rebin ! ! ----------------------- VARIABLE DIRECTORY -------------------------- ! ! Arguments ... ! ! maxin int : Maximum size of arrays ! pix_size real : pixel size in arcmin per pixel ! c_min real : Minimum counts/bin (user defined) ! after rebin ! inrad_lo real : Before rebin,lower edge of radial bins (in arcmins) ! inrad_hi real : Before rebin,upper edge of radial bins (in arcmins) ! in_pix real : Array of No. of pixels in each bin ! no longer argument 1.0.1; internal array ! incts real : Before rebin, array containing ctc/arcmin^2 ! inerr real : Before rebin, array of errors on incts ! area_wgt real : Area weighting factors ! rad_lo real : After rebin, lower edge of radial bins (in arcmins) ! rad_hi real : After rebin, upper edge of radial bins (in arcmins) ! n_pix real : After rebin, Array of No. of pixels in bins ! cts_arc2 real : After rebin, array of cts/arcmin^2 ! err_arc2 real : After rebin, error on cts/arcmin^2 ! sumraw real : Sum of raw counts ! chatter int : Chattines flag (<5 quiet,>5 normal,>20 noisy) ! negpres logical: true if negative counts present ! lbin logical: true if last bin has < c_min counts ! ! Internals ... ! ! rawcts real : Raw counts ! sumrcts real : Sum of raw counts that are rebinned for current rebin ! sumpix real : Sum of No. of pixels in a bin when rebin occurs ! onpts int : No. radial of data values after rebin ! ibin int : No. of bin which is currently in rebin process ! step int : Counter for No. of bin accumulated for new bin ! remain int : Counter for remaining No. of data values ! i,j,k int : counters for loops ! rebin logical: True if current bin requires rebinning ! absumr real : Absolute value of sumrcts ! tmaxin integer: tmaxin=maxin ! : it was made this way since one cannot assign dimension ! : passed to an unpassed variable ! ! ------------------ AUTHORS/MODIFICATION HISTORY ------------------------ ! ! Rehana Yusaf (1993 January) ! Rehana Yusaf (1993 Febuary) : Change code to rebin PSPCRPSF format data ! Rehana Yusaf (1993 June 8) : Add negpres and lbin ! Rehana Yusaf (1993 Sept 13) 1.0.1; Arguments changed due to reading ! RPSF Ver 1993a in RINFILE instead of 1992a ! ! Rehana Yusaf (1994 Jan 20) 1.0.2; bkgd passed to this routine, and taken away ! from the sum of the raw counts. ! Rehana Yusaf (1994 Feb 3) 1.0.3; sumraw calculated correctly ! ! Rehana Yusaf (1995 Feb 28) 1.0.4; add calc and bkgd_rad parameters ! Rehana Yusaf (1995 May 2) 1.0.5; If calc and calculated bkgd gives ! a -ve sumtcts ( due to excl regions ?) ! then calculate bkgd by setting sumtcts ! to 1, bkgd = (sumcts - 1)/sumpix ! Banashree Mitra Seifert (1996, Jan) 1.1.0: ! . Introduced tmaxin as integer tmaxin=maxin ! It was made this way since one cannot assign dimension ! passed to an unpassed variable ! . Replaced pi = 3.14159.... ! by pi = 4.*atan(1.) ! . Introduced screen display routines ! wtinfo,wterrm,wtwarm ! ! Banashree Mitra Seifert (1997, Oct 10) 1.2.0: ! . The calculation of area_wgt was wrong. binning was ok, but while ! running pcrpsf/hrirpsf after binning, resulted in a mismatch between ! the psf model normalization and the data further down the ! data-analysis pipeline. This error happened when binning requires ! counts from more than one bin from unbinned input data. So replaced ! area_wgt(j,1,1) = sumarea/area ! by ! area_wgt(j,1,1) = sumarea/k ! !------------------------------------------------------------------------- character(5) version parameter (version = '1.2.0') character(11) subname !- ! ------------------- USER INFORMATION ------------------------------ pi = 4. * atan(1.) subname = 'cpsf_rebin' subinfo = 'using ' // subname // ' Ver ' // version call wtinfo(chatter, 10, 1, subinfo) ! ------------------ INITIALISE VARIABLES ---------------------------- IF (calc) THEN IF (bkgd_rad.GE.inrad_hi(npts)) THEN subinfo = 'Inner radius for bkd calculation is > '& // 'Outer radius' call wterrm(subname, version, subinfo) ierr = 1 return ENDIF per_bkgd = bkgd_rad / inrad_hi(npts) IF (per_bkgd.GE.(0.75)) THEN subinfo = 'The inner radius for bkgd calculation'& // ' is > 75% of the total radius' call wtwarm(subname, version, chatter, 1, subinfo) ENDIF ENDIF sumrcts = 0 sumpix_bkgd = 0 sumrcts_bkgd = 0 do i = 1, npts area = pi * (inrad_hi(i)**2 - inrad_lo(i)**2) in_pix(i) = area_wgt(i, 1, 1) * area / pix_size**2 IF (calc) THEN IF (inrad_lo(i).GE.bkgd_rad) THEN sumpix_bkgd = sumpix_bkgd + in_pix(i) rawcts = incts(i, 1, 1) * in_pix(i) * pix_size**2 sumrcts_bkgd = sumrcts_bkgd + rawcts ENDIF ENDIF enddo IF (calc) THEN IF (sumrcts_bkgd.GE.0.0) THEN bkgd = sumrcts_bkgd / sumpix_bkgd ELSE subinfo = 'the sum of the counts in the'& // ' bkgd region is < 0' call wtwarm(subname, version, chatter, 1, subinfo) subinfo = 'the bkgd has been set to 0' call wtinfo(chatter, 10, 1, subinfo) bkgd = 0 ENDIF ENDIF rebin = .false. negpres = .false. lbin = .false. step = 1 onpts = 0 i = 1 j = 1 totpix2 = pi * (inrad_hi(npts)**2 - inrad_lo(1)**2) / pix_size**2 ! --------- A LOOP WHICH CHECKS EACH BIN AND REBINS IF NECESSARY ------ sumraw = 0.0 absumr = 0.0 totpix = 0.0 do WHILE (i.LE.npts) step = 1 rawcts = incts(i, 1, 1) * in_pix(i) * pix_size**2 IF (rawcts.LT.(0.0)) THEN IF (.NOT.negpres) THEN subinfo = 'input data contains -ve counts !' call wtwarm(subname, version, chatter, 1, subinfo) negpres = .true. ENDIF ENDIF IF (rawcts.GE.c_min) THEN rad_lo(j) = inrad_lo(i) rad_hi(j) = inrad_hi(i) cts_arc2(j, 1, 1) = incts(i, 1, 1) err_arc2(j, 1, 1) = inerr(i, 1, 1) area_wgt(j, 1, 1) = area_wgt(i, 1, 1) sumraw = sumraw + rawcts totpix = totpix + in_pix(i) ELSE rebin = .true. ibin = i sumrcts = 0 absumr = 0 sumpix = 0 sumarea = 0 k = 0 remain = npts - ibin do WHILE (rebin.AND.(k.LE.remain)) sumrcts = sumrcts + rawcts absumr = absumr + ABS(rawcts) sumpix = sumpix + in_pix(ibin + k) sumarea = sumarea + area_wgt(ibin + k, 1, 1) IF (sumrcts.GE.c_min) THEN rebin = .false. ENDIF k = k + 1 rawcts = incts(ibin + k, 1, 1) * in_pix(ibin + k) * pix_size**2 IF (rawcts.LT.(0.0)) THEN IF (.NOT.negpres) THEN subinfo = 'input data contains -ve counts!' call wtwarm(subname, version, chatter, 1, subinfo) negpres = .true. ENDIF ENDIF enddo ! ---------- IF LAST BIN < C_MIN PRINT WARNING ---------------- IF (rebin.AND.(k.GT.remain)) THEN subinfo = 'last bin has less than min cnts/bin' call wtwarm(subname, version, chatter, 1, subinfo) lbin = .true. ENDIF rad_lo(j) = inrad_lo(ibin) rad_hi(j) = inrad_hi(ibin + k - 1) cts_arc2(j, 1, 1) = sumrcts / (sumpix * pix_size**2) err_arc2(j, 1, 1) = SQRT(absumr) / (sumpix * pix_size**2) area = pi * (rad_hi(j)**2 - rad_lo(j)**2) !cc this calculation of area_wgt was wrong. binning was ok, but while !cc running pcrpsf/hrirpsf after binning, resulted in a mismatch between !cc the psf model normalization and the data further down the !cc data-analysis pipeline. This error happened when binning requires !cc counts from more than one bin from unbinned input data. So replaced !cc area_wgt(j,1,1) = sumarea/area !cc by !cc area_wgt(j,1,1) = sumarea/k area_wgt(j, 1, 1) = sumarea / k step = k totpix = totpix + sumpix sumraw = sumraw + sumrcts ENDIF j = j + 1 i = i + step enddo sumcts = sumraw sumraw = sumcts - (totpix * bkgd) sumtcts = sumcts - (totpix2 * bkgd) IF (calc) THEN IF (sumtcts.LE.0.0) THEN bkgd = (sumcts - 1) / totpix2 sumtcts = 1 sumraw = sumcts - totpix * bkgd subinfo = 'the background value has been calculated such'& // ' that the sum of the' call wtinfo(chatter, 10, 1, subinfo) subinfo = 'theoretical counts is set to 1. If the background'& // ' is calculated' call wtinfo(chatter, 10, 1, subinfo) subinfo = ' in the usual way then sumtcts is -ve, this may be'& // ' due to excluded regions.' call wtinfo(chatter, 10, 1, subinfo) ENDIF ENDIF IF (calc) THEN write(subinfo, '(a,f20.9)') ' calculated bkgd value :', bkgd call wtinfo(chatter, 10, 2, subinfo) ENDIF write(subinfo, '(a,F20.9)')' sum of pixels calculated using'& // ' area of circle:', totpix2 call wtinfo(chatter, 20, 2, subinfo) write(subinfo, '(a,F20.9)')' actual sum of pixels :', totpix call wtinfo(chatter, 20, 2, subinfo) IF ((totpix / totpix2).LE.(0.9)) THEN subinfo = 'actual sum of pixels is < 90% of sum'& // ' of pixels calculated using area of circle, '& // 'where the outer radius is used' call wtwarm(subname, version, chatter, 9, subinfo) subinfo = 'note: This may be due to regions being excluded' call wtinfo(chatter, 10, 1, subinfo) ENDIF write(subinfo, '(a,F20.9)')'total sum of counts :', sumcts call wtinfo(chatter, 20, 2, subinfo) write(subinfo, '(a,F20.9)')'sum of counts in source :', sumraw call wtinfo(chatter, 20, 2, subinfo) write(subinfo, '(a,a,F20.9)')'theoretical sum of counts'& , '(corrected for any excl regions):', sumtcts call wtinfo(chatter, 20, 2, subinfo) IF (sumrcts.LT.0) THEN subinfo = 'sum of counts is negative this may be due to'& // ' background dominating the source' call wtwarm(subname, version, chatter, 1, subinfo) ierr = 1 return ENDIF IF (sumtcts.LT.0) THEN subinfo = 'theoretical Sum of counts is negative, this'& // ' may be due to background dominating the source' call wtwarm(subname, version, chatter, 1, subinfo) ierr = 1 return ENDIF IF ((sumraw / sumcts).LE.(0.1)) THEN subinfo = 'source is less than 10% of total counts !' call wtwarm(subname, version, chatter, 1, subinfo) ENDIF onpts = j - 1 npts = onpts return end ! ---------------------------------------------------------------------+ ! END OF SUBROUTINE CPSF_REBIN ! ---------------------------------------------------------------------+ !+CPSF_RESULTS subroutine cpsf_results(outfile, infile, mnver, npts, rad_lo, & rad_hi, radunit, ntheta, theta_lo, & theta_hi, thetaunit, nenerg, energ_lo, & energ_hi, energunit, cts_arc2, qerror, & err_arc2, rpsfunit, qarea, area_wgt, & pix_size, maxrad, maxtheta, bkgd, c_min, & telrad, instrad, negpres, lbin, chanmin, & chanmax, sumrcts, sumtcts, & ierr, chatter, killit) ! --------------- DESCRIPTION ---------------------------------------- ! This subroutine writes observed and thoeretical radial psf to output ! file. The output file is in FITS format RPSFVER 1993a ! ! ------------------------ VARIABLES ---------------------------------- IMPLICIT NONE character*(*) outfile, infile character*(*) telrad, instrad character(5) mnver character*(*) rpsfunit, radunit, thetaunit, energunit integer maxrad, maxtheta, ntheta, nenerg, chanmin, chanmax integer npts, chatter, ierr real rad_lo(*), rad_hi(*), cts_arc2(maxrad, maxtheta, 1) real err_arc2(maxrad, maxtheta, 1) real c_min, sumrcts, sumtcts, pix_size, bkgd real theta_lo(*), theta_hi(*), energ_lo(*), energ_hi(*) real area_wgt(maxrad, maxtheta, 1) logical negpres, lbin, qerror, qarea, killit ! -------------------- INTERNAL VARIABLES------------------------------- character(80) subinfo character(16) hduclas3, creator character(8) cc_min, extname integer ounit, ierrstat, maxhist, nk_hist, nk_comm, maxcomm parameter (maxhist = 10, maxcomm = 10) character(80) hist(maxhist), comms(maxcomm) real pix ! ---------------------- VARIABLE DIRECTORY --------------------------- ! ! Arguments ... ! ! outfile char : Name of results file ! infile char : Radial Profile input file ! mnver char : Main program version ! telrad char : Telescope name for radial data ! instrad char : Instrument name for radial data ! maxrad int : Size of arrays for radial profile data ! maxpred int : Size of arrays for predicted radial psf data ! rad_lo real : Array of lower edger of radial bins, in arcmin ! rad_hi real : Array of upper edge of radial bins, in arcmin ! theta_lo real : Lower theta values ! theta_hi real : Upper theta values ! energ_lo real : Lower energy boundary ! energ_hi real : Upper energy boundary ! ntheta int : No. of theta values ! nenerg int : No. of energy values ! qerror logical: True if err_arc2 present ! qarea logical: True if area weighting is present ! cts_arc2 real : Array of radial profile in counts/arcmin^2 ! err_arc2 real : Array of errors on cts_arc2 ! pix_size real : Pixel size in arcmin per pixel ! c_min real : Minimum No. of counts/bin (user defined) ! npts int : No. of observed radial values ! chatter int : Chatter flag (<5 quiet,>5 normal, >20 noisy) ! negpres logical: true if negative counts present in dataset ! lbin logical: true if last has < c_min counts ! ! Internals ... ! ! errtxt char : Error text obtained from ftgerr ! errstr char : Error string ! rad_del real : Radius of predicted bin ! pix real : Pixel size in degrees ! ounit int : Output file number ! ierrstat int : Error flag ! i,j,k int : Counters for loops ! ! -------------------- COMPILATION AND LINKING ------------------------ ! ! Link with FTOOLS, FITSIO and CALLIB ! ! -------------------- CALLED ROUTINES ---------------------------- ! ! subroutine OP_NPA : CALLIB routine which opens a FITS file, writes ! the header, and null primary array ! subroutine WTRPF1 : CALLIB,Writes observed radial psf data to fits file ! subroutine FTCLOS : FITSIO routine which closes a FITS file ! subroutine WT_FERRMSG : CALLIB routine which writes FITSIO and routine ! error message ! ! -------------------- AUTHORS/MODIFICATION HISTORY -------------------- ! ! Rehana Yusaf (1993 January 15) ! Rehana Yusaf (1993 March 4) : Adding more comments ! Rehana Yusaf (1993 June 8) : Adding negpres and lbin ! Rehana Yusaf (1993 Sept 13) 1.0.1; Use WT_RPSF1993a instead of ! WT_RPSF1992a ! Rehana Yusaf (1994 Jan 20) 1.0.2; WT_RPSF1993a has been updated to WTRPF1 ! this write HDUCLASS keywords ! Rehana Yusaf (1995 MAR 7) 1.0.3; sumtcts is also written to outfile ! Rehana Yusaf (1995 April 25) 1.0.4; add killit ! Banashree Mitra Seifert (1996, Jan) 1.1.0; ! . Variables dimensions are carried from the ! calling subroutine ! . implemented screen display routines ! wtinfo,wtferr ! --------------------------------------------------------------------------- character(5) version parameter (version = '1.1.0') character(14) subname !- ! ------------------ USER INFO ----------------------------------------+ subname = 'cpsf_results' subinfo = 'using ' // subname // version call wtinfo(chatter, 10, 1, subinfo) ! --------------- MODIFYING OUTPUT INTO DESIRED FORM ------------------ ! ------------------ WRITING TO FITS FILE ------------------------- ! OBSERVED DATA ... call cgetlun(ounit) call opnpa(outfile, chatter, ounit, killit, ierrstat) nk_hist = 0 write (cc_min, 100) c_min nk_comm = 3 comms(1) = 'DATA OBTAINED FROM RADIAL PROFILE : ' // infile comms(2) = 'THE RADIAL PROFILE DATA IS REBINNED ' comms(3) = 'USING (USER DEFINED) MINIMUM COUNTS/BIN : ' // cc_min IF (lbin) THEN nk_comm = nk_comm + 1 comms(nk_comm) = 'LAST BIN CONTAINS LESS THAN MINUMUM COUNTS/BIN' ENDIF IF (negpres) THEN nk_comm = nk_comm + 1 comms(nk_comm) = 'WARNING:INPUT DATASET CONTAINS NEGATIVE COUNTS' ENDIF extname = 'OBS RPSF' pix = pix_size / float(60) If (bkgd.EQ.(0.0)) THEN hduclas3 = 'NET' ELSE hduclas3 = 'TOTAL' ENDIF ! ------------------ USER INFO ----------------------------------------+ call wtrpf1(ounit, extname, hduclas3, npts, rad_lo, rad_hi, radunit, & ntheta, theta_lo, theta_hi, thetaunit, nenerg, & energ_lo, energ_hi, energunit, cts_arc2, qerror, & err_arc2, rpsfunit, qarea, area_wgt, hist, nk_hist, & comms, nk_comm, telrad, instrad, maxrad, maxtheta, & ierr, chatter) IF (ierr.NE.0) THEN subinfo = 'error in writing RPSF extension ' call wterrm(subname, version, subinfo) return ENDIF ! ---------------- Write Pixelsize ----------------------------- ierrstat = 0 call ftpkye(ounit, 'PIXSIZE', pix, 8, 'pixelsize in deg', ierrstat) subinfo = 'writing PIXSIZE ' call wtferr(subname, version, ierrstat, subinfo) ! ----------------- Write Backgrnd ---------------------------------- ierrstat = 0 call ftpkye(ounit, 'BACKGRND', bkgd, 8, & 'Background count rate in cts/pixel ', ierrstat) subinfo = 'writing BACKGRND ' call wtferr(subname, version, ierrstat, subinfo) ! ------------ Write CHANMIN and CHANMAX if not zero ------------------ ierrstat = 0 IF (chanmin.NE.0) THEN call ftpkyj(ounit, 'CHANMIN', chanmin, & 'Minimum PI channel for image ', ierrstat) subinfo = 'writing CHANMIN' call wtferr(subname, version, ierrstat, subinfo) ENDIF ierrstat = 0 IF (chanmax.NE.0) THEN call ftpkyj(ounit, 'CHANMAX', chanmax, & 'Maximum PI channel for image ', ierrstat) subinfo = 'writing CHANMAX' call wtferr(subname, version, ierrstat, subinfo) ENDIF ! --------------- WRITE SUMRCTS ---------------------------------------+ ierrstat = 0 call ftpkye(ounit, 'SUMRCTS', sumrcts, 8, & 'Sum of source counts under profile', ierrstat) subinfo = 'writing SUMRCTS' call wtferr(subname, version, ierrstat, subinfo) ! --------------- WRITE SUMTCTS ---------------------------------------+ ierrstat = 0 call ftpkye(ounit, 'SUMTCTS', sumtcts, 8, & 'Theoretical Sum source counts(corr for any '& // 'excl regions)', ierrstat) subinfo = 'writing SUMTCTS' call wtferr(subname, version, ierrstat, subinfo) ! --------------- WRITE CREATOR KEYWORD -------------------------------+ ierrstat = 0 creator = 'rbnrpsf ' // mnver call ftpkys(ounit, 'CREATOR', creator, & 's/w task which wrote this dataset', ierrstat) subinfo = 'writing CREATOR keyword' call wtferr(subname, version, ierrstat, subinfo) 100 format(F6.0) ierrstat = 0 call ftclos(ounit, ierrstat) subinfo = 'closing output file' call wtferr(subname, version, ierrstat, subinfo) return end ! --------------------------------------------------------------------- ! END OF SUBROUTINE RESULTS ! ---------------------------------------------------------------------