!+RPSFQDP subroutine rpsfqp ! ! -------------------- DESCRIPTION -------------------------------- ! ! RPSFQDP reads the input file and looks for whether the file has RPSF data ! or the REEF data (observation and/or theoretical) ! If found to be RPSF data then it calls CQDP_RPSF to read input data ! If found to be REEF data then it calls CQDP_REEF to read input data ! CQDP_RPSF reads OBS RPSF and PRED PSF extensions in CALRPSF format FITS ! file. Then it subsequently writes the RPSF/REEF data into an output file ! in QDP format. ! ! -------------------- VARIABLES -------------------------------- ! implicit none character(180) fits_psf, qdp_psf, pred_psf character(150) subinfo integer maxrad, maxpred, maxtheta, maxenerg, nrad, npred integer chatter, ierr, status, iget, ineed integer nset, i character(20) hdu_sav(9, 50) real, allocatable :: p_rad_lo(:), p_rad_hi(:) real, allocatable :: p_cts(:), p_err(:) real, allocatable :: p_prad_lo(:), p_prad_hi(:) real, allocatable :: p_pred(:), p_pred_reef(:) real, allocatable :: p_del_rad(:), p_pdel_rad(:) real, allocatable :: p_rad_mean(:), p_prad_mean(:) real, allocatable :: p_reef(:), p_reef_err(:) real rescale real pix_size, backgrnd logical theo_pres, data_pres, killit logical rpsf_pres, reef_pres ! ! ------------------- VARIABLE DIRECTORY ------------------------------ ! ! fits_psf char : Name of PSF file (user defined) ! qdp_psf char : Name of QDP file (user defined) ! subinfo char : Subroutine information ! chatter int : Chatter flag (<5 quiet,>5normal,>20 noisy) ! maxrad int : Maximum array size for observed psf data ! maxpred int : Maximum array size for predicted psf data ! nrad int : Counter for observed psf data ! npred int : Counter for predicted psf data ! rad_lo real : Array of lower edge of observed radial bins ! rad_hi real : Array of upper edge of observed radial bins ! cts real : Array of observed radial profile, in counts ! err real : Array of errors on cts ! prad_lo real : Array of lower edge of predicted model bins ! prad_hi real : Array of upper edge of predicted model bins ! pred real : Array of theoretical PSF ! rad_mean real : Array of mean of observed radial bins ! del_rad real : Array of half-width of observed radial bins ! prad_mean real : Array of mean predicted model bins ! pdel_rad real : Array of half-width predicted model bins ! rescale real : Scalar rescaling factor for pred ONLY ! ierr int : error flag, 0 is okay ! thoe_pres logical: true if theoretical extension is present ! ! --------------------- CALLED ROUTINES ------------------------------ ! ! subroutine CQDP_GP : Reads user defined parameters using XPI ! subroutine CQDP_RPSF : Reads FITS format RPSF file ! subroutine CQDP_REEF : Reads FITS format REEF file ! subroutine CQDP_CONV : Converts data into desired format ! subroutine CQDP_WT : Writes data into QDP format output file ! ! ------------------ AUTHORS/MODifICATION ---------------------------- ! ! Rehana Yusaf (1993 Feb) ! Rehana Yusaf (1993 May 27) : If theoretical ext not present,only ! OBS RPSF written ! Rehana Yusaf (1993 SEpt 21) : Change to read RPSFVER 1993a instead of ! 1992a ! Rehana Yusaf (1994 Jan 31) 1.0.2; Additional paramters ... ! . pred_psf, filename for pred data ! . data_pres, logical,true if obs data ! present ! Rehana Yusaf (1994 Feb3) 1.0.3; CQDP_RPSF updated so that the program ! is NOT terminated if pix_size not present ! Rehana Yusaf (1994 Sept 12) 1.0.4; Minor cosmetics to cqdp_rpsf and cqdp_wt ! Rehana Yusaf (1995 Jan 13) 1.0.5; update _gp to no longer read defval from parfile ! Ian M George (1995 Apr 05) 2.0.0; added rescale parameter to allow user to ! rescale the theoretical ONLY ! Rehana Yusaf (1995 April 25) 2.0.1; add clobber ! Banashree Mitra Seifert (1995 December) 3.0.0; ! . replaced calls to FNDHDU, FNDEXT, FTMRHD, FTMAHD by ! CALL MVEXT ! . added reader for REEF data RDEEF1 ! . added Dynamic Memory Allocation ! . screen display subroutines used ! wtbegm,wtendm,wterrm,wtinfo ! Rehana Yusaf (1996 Feb 22) 3.0.1; bugfix ! Peter D Wilson (1998 Jul 01) 3.0.2: Updated for new FCPARS behavior ! ! Bryan K Irby (2018 Jan 14) 3.1.0: ! . Replaced udmget/udmfre with allocate/deallocate ! MFC (2020 Apr 16) 3.1.1 ! . f90 version ! ! ----------------------------------------------------------------------- character(5) version parameter (version = '3.1.1') !- ! ---------------------------------------------------------------------- character(40) taskname character(8) subname !ccc common/task/taskname subname = 'rpsfqdp' taskname = 'rpsfqdp' ! ----------------- Get parameter file --------------------------------+ call cqdp_gp(fits_psf, pred_psf, qdp_psf, rescale, ierr, & killit, chatter) call wtbegm(taskname, version, chatter) if (ierr .ne. 0) then goto 100 endif maxrad = 1000 maxpred = 1000 maxtheta = 5 maxenerg = 5 ! ----------------------- 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_rad_mean(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * 4 status = 0 allocate(p_del_rad(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * 4 status = 0 allocate(p_pdel_rad(maxpred), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxpred * 4 status = 0 allocate(p_prad_mean(maxpred * maxtheta * maxenerg), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxpred * maxtheta * maxenerg * 4 status = 0 allocate(p_cts(maxrad * maxtheta * maxenerg), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * maxenerg * 4 status = 0 allocate(p_err(maxrad * maxtheta * maxenerg), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * maxenerg * 4 status = 0 allocate(p_prad_lo(maxpred), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxpred * 4 status = 0 allocate(p_prad_hi(maxpred), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxpred * 4 status = 0 allocate(p_pred(maxpred * maxtheta * maxenerg), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxpred * maxtheta * maxenerg * 4 status = 0 allocate(p_reef(maxrad * maxtheta * maxenerg), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * maxenerg * 4 status = 0 allocate(p_reef_err(maxrad * maxtheta * maxenerg), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * maxenerg * 4 status = 0 allocate(p_pred_reef(maxpred * maxtheta * maxenerg), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxpred * maxtheta * maxenerg * 4 50 ineed = 4 * maxrad * 4 + 3 * maxpred * 4 + 4 * maxrad * maxtheta * maxenerg * 4 + & 3 * maxpred * maxtheta * maxenerg * 4 write(subinfo, '(a,i10)')'dmasize required for this task=', ineed call wtinfo(chatter, 15, 2, subinfo) write(subinfo, '(a,i10)')'total bytes of memory =', iget call wtinfo(chatter, 15, 2, subinfo) if (status .ne. 0) then ierr = -1 subinfo = ' failed to allocate dynamic memory ' call wterrm(subname, version, subinfo) goto 100 endif ! -------------------------------------------------------------------- ! Calling the subroutine decide to search for the type of input file ! e.g., type RPSF --> predicted or observed ! type REEF --> predicted or observed ! nset gives the no of such sets found ! --------------------- Call Decide ---------------------------------- nset = 0 call decide(fits_psf, nset, hdu_sav, chatter) Do i = 1, nset if(hdu_sav(i, 2) .eq. 'RPRF') then rpsf_pres = .true. reef_pres = .false. call cqdp_rpsf(fits_psf, pred_psf, p_rad_lo, & p_rad_hi, p_cts, p_err, & p_prad_lo, p_prad_hi, p_pred, & maxrad, maxpred, maxtheta, nrad, npred, & theo_pres, data_pres, pix_size, backgrnd, & ierr, chatter) if (ierr .ne. 0) then goto 100 endif call cqdp_conv(p_rad_lo, p_rad_hi, p_prad_lo, & p_prad_hi, p_rad_mean, p_del_rad, & p_prad_mean, p_pdel_rad, maxrad, & maxpred, nrad, npred, theo_pres, data_pres, ierr, & chatter) if (ierr .ne. 0) then goto 100 endif call cqdp_wt(fits_psf, qdp_psf, p_rad_mean, & p_del_rad, p_cts, p_err, & p_prad_mean, p_pdel_rad, p_pred, & maxrad, maxpred, maxtheta, nrad, npred, theo_pres, & data_pres, rescale, pix_size, backgrnd, rpsf_pres, & reef_pres, ierr, chatter, killit) if (ierr .ne. 0) then goto 100 endif elseif(hdu_sav(i, 2) .eq. 'REEF') then rpsf_pres = .false. reef_pres = .true. call cqdp_reef(fits_psf, pred_psf, p_rad_lo, & p_rad_hi, p_reef, p_reef_err, & p_prad_lo, p_prad_hi, & p_pred_reef, maxrad, maxpred, maxtheta, & nrad, npred, theo_pres, data_pres, pix_size, & backgrnd, ierr, chatter) if (ierr .ne. 0) then goto 100 endif call cqdp_conv(p_rad_lo, p_rad_hi, p_prad_lo, & p_prad_hi, p_rad_mean, p_del_rad, & p_prad_mean, p_pdel_rad, maxrad, & maxpred, nrad, npred, theo_pres, data_pres, ierr, & chatter) if (ierr .ne. 0) then goto 100 endif call cqdp_wt(fits_psf, qdp_psf, p_rad_mean, & p_del_rad, p_reef, p_reef_err, & p_prad_mean, p_pdel_rad, & p_pred_reef, & maxrad, maxpred, maxtheta, nrad, npred, theo_pres, & data_pres, 1., pix_size, 0., rpsf_pres, reef_pres, & ierr, chatter, killit) if (ierr .ne. 0) then goto 100 endif else subinfo = 'input file is not an RPSF nor a REEF file' call wterrm(subname, version, subinfo) goto 100 endif ! ---------------------------------------------------------------------+ ENDDO ! --------- free the dynamic memory ----------------------------- status = 0 deallocate(p_rad_lo, stat = status) status = 0 deallocate(p_rad_hi, stat = status) status = 0 deallocate(p_rad_mean, stat = status) status = 0 deallocate(p_del_rad, stat = status) status = 0 deallocate(p_pdel_rad, stat = status) status = 0 deallocate(p_prad_mean, stat = status) status = 0 deallocate(p_cts, stat = status) status = 0 deallocate(p_err, stat = status) status = 0 deallocate(p_prad_lo, stat = status) status = 0 deallocate(p_prad_hi, stat = status) status = 0 deallocate(p_pred, stat = status) status = 0 deallocate(p_reef, stat = status) status = 0 deallocate(p_reef_err, stat = status) status = 0 deallocate(p_pred_reef, stat = status) if (status .ne. 0) then subinfo = ' failed to de-allocate memory ' call wterrm(subname, version, subinfo) ierr = 99 endif 100 call wtendm(taskname, version, ierr, chatter) return end ! ---------------------------------------------------------------------- ! end OF RPSFQDP ! ------------------------------------------------------------------ !+CQDP_GP subroutine cqdp_gp(infile, pred_psf, outfile, rescale, ierr, & killit, chatter) ! ------------------------------------------------------------ ! --- DESCRIPTION ------------------------------------------------------ ! ! This routine gets the user defined parameters, that is the input and ! output file names. As well as chatter flag. ! ! --- VARIABLES -------------------------------------------------------- ! IMPLICIT NONE character(180) infile, outfile, ill_files(3) character(180) filename, pred_psf, tmpfile integer err, chatter, ierr, extnum, n_ill character(80) errmess, subinfo real rescale logical killit, valfil !logical ext ! ! --- VARIABLE DIRECTORY --- ! ! infile char : Name of FITS PSF file (user defined) ! outfile char : Name of QDP output file (user defined) ! chatter int : Chattiness flag (<5quiet,>5normal,>20 noisy) ! rescale real : Scalar rescaling factor for pred ONLY ! err int : Error flag ! ierr int : Error flag, returned to main ! version char : Version of subroutine ! errstr char : Error string for this routine ! errmess char : Error message ! ! --- COMPILATION AND LINKING --- ! ! Link with FTOOLS ! ! --- AUTHORS/MODifICATION HISTORY --- ! ! Rehana Yusaf (1993 Feb) ! Rehana Yusaf (1993 May 27) : INQUIRE added ! Rehana Yusaf (1994 Jan 13) : remove defval ! Ian M George (2.0.0:95 Apr 05) added rescale parameter, removed header ! Rehana Yusaf (2.0.1:95 Apr 25); add clobber ! ! Banashree Mitra Seifert (1996, Jan) 2.1.0: ! . Screen display subroutines used ! wtinfo,wterrm ! Peter D Wilson (1998 Jul 01) 2.1.1: ! . Drop INQUIRE tests. ! -------------------------------------------------------------------- character(5) version parameter (version = '2.1.1') character(8) subname !- ! -------------------------------------------------------------------- subname = 'cqdp_gp' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! -------------- READING INPUT FILE PARAMETER ----------------------- err = 0 call uclgst('datafile', infile, err) if (err.NE.0) then errmess = 'getting datafile parameter' call wterrm(subname, version, errmess) endif call crmvlbk(infile) if (infile .eq. ' ') then errmess = ' Input file must be entered ,or NONE !' call wterrm(subname, version, errmess) ierr = 1 return endif ! PDW 7/1/98: Don't bother! Let FTOPEN determine if file exists ! call fcpars(infile,filename,extnum,err) ! call crmvlbk(filename) ! tmpfile = filename ! call ftupch(tmpfile) ! if (tmpfile.NE.'NONE') then ! INQUIRE(FILE=filename,EXIST=ext) ! if (.NOT.ext) then ! errmess = ' input file does not exist !' ! call wterrm(subname,version,errmess) ! ierr = 1 ! return ! endif ! endif ! ! -------------------- READ PRED_PSF --------------------- ! call uclgst('predfile', pred_psf, err) if (err.NE.0) then errmess = ' getting predfile parameter' call wterrm(subname, version, errmess) endif call uclpst('predfile', '%', err) if (err.NE.0) then errmess = ' Putting default predfile parameter' call wterrm(subname, version, errmess) endif call crmvlbk(pred_psf) if ((pred_psf(1:1).EQ.'%').OR.(pred_psf.EQ.' ')) then if (infile.EQ.'NONE') then errmess = ' Either predfile or infile have to be entered' call wterrm(subname, version, errmess) ierr = 2 return ELSE pred_psf = infile endif ELSE ierr = 0 call fcpars(pred_psf, filename, extnum, ierr) tmpfile = filename call ftupch(tmpfile) if (tmpfile.EQ.'NONE') then if (infile.EQ.'NONE') then errmess = ' Either predfile or infile have to be entered' call wterrm(subname, version, errmess) ierr = 2 return endif ! PDW 7/1/98: Don't bother! Let FTOPEN determine if file exists ! ELSE ! INQUIRE(FILE=filename,EXIST=ext) ! if (.NOT.ext) then ! errmess = ' input predfile does not exist !' ! call wterrm(subname,version,errmess) ! ierr = 1 ! return ! endif endif endif ! ! ------------- READING OUTPUT FILE PARAMETER ----------------------- ! call uclgst('outfile', outfile, err) if (err.NE.0) then errmess = ' getting outfile parameter' call wterrm(subname, version, errmess) endif ! ! <<<--- TERMINATE if USER DOES NOT ENTER OUTFILE--->>> ! call crmvlbk(outfile) if (outfile.EQ.' ') then errmess = ' must enter outfile name !' call wterrm(subname, version, errmess) ierr = 1 return endif ! GET CLOBBER call uclgsb('clobber', killit, err) if (err.NE.0) then errmess = ' getting clobber' killit = .false. call wterrm(subname, version, errmess) endif n_ill = 0 call ck_file(outfile, ill_files, n_ill, valfil, & killit, chatter) if (.NOT.valfil) then errmess = ' Invalid outfile name' call wterrm(subname, version, errmess) ierr = 1 return endif ! ! ----------------- READING CHATTER FLAG ------------------------ ! call uclgsi('chatter', chatter, err) if (err.NE.0) then errmess = ' getting chatter parameter' call wterrm(subname, version, errmess) endif ! ! <<<--- READING RESCALE PARAMETER --->>> ! call uclgsr('rescale', rescale, err) if (err.NE.0) then errmess = ' getting rescale parameter' call wterrm(subname, version, errmess) errmess = 'setting Rescale = 1.0' call wterrm(subname, version, errmess) err = 0 rescale = 1.0 endif return end ! ------------------------------------------------------------------- ! end OF SUBROUTINE CQDP_GP ! ------------------------------------------------------------------- !+CQDP_RPSF subroutine cqdp_rpsf(fits_psf, pred_psf, rad_lo, rad_hi, cts, & err, prad_lo, prad_hi, pred, maxrad, maxpred, maxtheta, nrad, & npred, theo_pres, data_pres, pix_size, backgrnd, ierr, chatter) ! ! --- DESCRIPTION --------------------------------------------------- ! ! Reading RPSF FITS file using FITSIO. ! ! --- VARIABLES ----------------------------------------------------- ! IMPLICIT NONE character(180) fits_psf, pred_psf integer maxrad, maxpred, nrad, npred, chatter, ierr, maxtheta real rad_lo(maxrad), rad_hi(maxrad) real cts(maxrad, maxtheta, 1), err(maxrad, maxtheta, 1) real prad_lo(maxpred), prad_hi(maxpred), pred(maxpred, maxtheta, 1) logical theo_pres, data_pres ! ! --- INTERNAL VARIABLES --- ! character(40) errstr, comm character(80) subinfo character(200) errinfo character(8) telescop, instrume real pix_size, backgrnd, tot_cts integer status, iunit character(16) radunit, thetaunit, energunit, rpsfunit character(16) hduclas3 integer ntheta, nenerg real theta_lo(1), theta_hi(1), energ_lo(1), energ_hi(1) real perr(9092, 1, 1), parea_wgt(9092, 1, 1), area_wgt(300, 1, 1) logical qerror, qarea integer nsearch, ninstr, next(50) integer errflg character(20) instr(50), outhdu(9, 50), outver(9, 50), extname character(20) extnames(9, 50) ! ! --- VARIABLE DIRECTORY --- ! ! Arguments ... ! ! fits_psf char : Name of FITS format PSF file ! maxen int : Maximum array size for observed psf data ! maxpred int : Maximum array size for predicted psf data ! nrad int : Counter for observed psf data ! npred int : Counter for predicted psf data ! rad_lo real : Array of lower edge of observed radial bins ! rad_hi real : Array of upper edge of observed radial bins ! cts real : Array of observed radial profile, in counts ! err real : Array of statistical errors ! prad_lo real : Array of lower edge of predicted model bins ! prad_hi real : Array of upper edge of predicted model bins ! pred real : Array of theoretical PSF ! ! Internals ... ! ! errmess char : Error message text ! errstr char : Error text for this routine ! version char : Subroutine version ! status int : Error flag for FITSIO call ! iunit int : Fortran unit number for file ! ! --- CALLED ROUTINES --- ! ! subroutine FTOPEN : FITSIO routine to open file ! subroutine FTCLOS : FITSIO routine to close file ! subroutine RD_RPSF1992a : (CALLIB) Reads observed PSF extension in FITS file ! subroutine CQDP_RTHEO : Reads theoretical PSF extension in FITS file ! subroutine WT_FERRMSG : Writes FITSIO error text if neccesary ! ! --- COMPILATION AND LINKING --- ! ! Link with FITSIO and FTOOLS ! ! --- AUTHORS/MODifICATION HISTORY --- ! ! Rehana Yusaf (1993 Feb) ! Rehana Yusaf (1994 Jan 31) 1.0.1; data_pres and pred_psf paramters added. ! rd_rpsf1993a renamed tp rdrpf1 ! Rehana Yusaf (1994 Sept 12) 1.0.2; only warn about theo not present at ! chatter > 10 ! Banashree Mitra Seifert (1996, Jan) 1.1.0: ! . Replaced by MVEXT ! . Screen display subroutines used ! wtinfo,wtferr ! Banashree Mitra Seifert (1996, Feb) 1.2.0: ! . Corrected call for MVEXT for both obs and theo ! RPSF file ! ! ----------------------------------------------------------------------- character(5) version parameter (version = '1.2.0') character(10) subname !- ! ----------------------------------------------------------------------- subname = 'cqdp_rpsf' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! ! ----------- OPENING PSF FILE ------------- ! status = 0 ninstr = 3 instr(1) = 'RESPONSE' instr(2) = 'RPRF' instr(3) = 'TOTAL' extname = 'OBS RPSF' nsearch = 50 call mvext (0, fits_psf, iunit, ninstr, instr, nsearch, next, & outhdu, extnames, outver, extname, errflg, chatter) if (errflg .ne. 0) then data_pres = .false. call ftclos(iunit, status) subinfo = ' closing PSF file' call wtferr(subname, version, status, subinfo) else data_pres = .true. endif ! ----------- READING OBSERVED PSF EXTENSION ------------- if(data_pres) then call rdrpf1(iunit, hduclas3, nrad, rad_lo, rad_hi, radunit, & ntheta, theta_lo, theta_hi, thetaunit, nenerg, & energ_lo, energ_hi, energunit, cts, qerror, err, & rpsfunit, qarea, area_wgt, telescop, instrume, & maxrad, maxtheta, ierr, chatter) if (ierr.NE.0) then call ftclos(iunit, status) subinfo = ' closing obs_psf file' call wtferr(subname, version, status, subinfo) return endif ! ! ------------------- READ PIXSIZE ------------------------- ! status = 0 call ftgkye(iunit, 'PIXSIZE', pix_size, comm, status) errinfo = ' reading PIXSIZE !' call wtferr(subname, version, status, errinfo) pix_size = pix_size * 60 ! ------------------- READ BACKGROUND -------------- status = 0 call ftgkye(iunit, 'BACKGRND', backgrnd, comm, status) errinfo = errstr // ' reading BACKGRND' call wtferr(subname, version, status, errinfo) if (status.NE.0) then errinfo = 'If rescaling theo curve, background required' call wtinfo(chatter, 1, 1, errinfo) endif ! ------------- Closing the input file ----------------------- status = 0 call ftclos(iunit, status) errinfo = ' closing' // fits_psf call wtferr(subname, version, status, errinfo) endif ! ----------------- closed the if data present statement ! --------------- READING THEORETICAL PSF EXTENSION -------- errflg = 0 ninstr = 3 instr(1) = 'RESPONSE' instr(2) = 'RPRF' instr(3) = 'PREDICTED' extname = 'THEO RPSF' nsearch = 50 call mvext (0, pred_psf, iunit, ninstr, instr, nsearch, & next, outhdu, extnames, outver, extname, & errflg, chatter) if (errflg .ne. 0) then theo_pres = .false. status = 0 call ftclos(iunit, status) subinfo = ' closing PSF file' call wtferr(subname, version, status, subinfo) return else theo_pres = .true. endif ! --------------------- READ THEO DATA -------------------------------+ if(theo_pres) then ierr = 0 call rdrpf1(iunit, hduclas3, npred, prad_lo, prad_hi, radunit, & ntheta, theta_lo, theta_hi, thetaunit, nenerg, & energ_lo, energ_hi, energunit, pred, qerror, perr, & rpsfunit, qarea, parea_wgt, telescop, instrume, & maxpred, maxtheta, ierr, chatter) if (ierr.EQ.1) then theo_pres = .false. ierr = 0 endif ! ------------------- READ PIXSIZE ------------------------------------+ status = 0 call ftgkye(iunit, 'PIXSIZE', pix_size, comm, status) subinfo = ' reading PIXSIZE !' call wtferr(subname, version, status, subinfo) pix_size = pix_size * 60 ! ------------- READ THEORETICAL TOTAL PSF COUNTS ----------------------+ status = 0 call ftgkye(iunit, 'SUMTCTS', tot_cts, comm, status) subinfo = ' reading SUMTCTS !' call wtferr(subname, version, status, subinfo) if (status .ne. 0) then tot_cts = 0 errinfo = 'SUMTCTS is not found ' call wtferr(subname, version, status, errinfo) endif ! ------------------- READ BACKGROUND FOR PSF DATA ------------------- status = 0 call ftgkye(iunit, 'BACKGRND', backgrnd, comm, status) errinfo = ' reading BACKGRND' call wtferr(subname, version, status, errinfo) if (status .ne. 0) then backgrnd = 0 errinfo = 'background is not supplied' call wtinfo(chatter, 1, 1, errinfo) endif backgrnd = backgrnd * pix_size status = 0 endif ! -------------------- CLOSING FILE ---------------------------------- status = 0 call ftclos(iunit, status) errinfo = ' closing ' // pred_psf call wtferr(subname, version, status, errinfo) 200 if ((.NOT.theo_pres).AND.(.NOT.data_pres)) then ierr = 5 errinfo = ' predicted and observed data not present' call wtinfo(chatter, 1, 1, errinfo) endif return end ! ---------------------------------------------------------------- ! end OF SUBROUTINE CQDP_RPSF ! ---------------------------------------------------------------- !+CQDP_CONV subroutine cqdp_conv(rad_lo, rad_hi, prad_lo, prad_hi, rad_mean, & del_rad, prad_mean, pdel_rad, maxrad, maxpred, nrad, & npred, theo_pres, data_pres, ierr, chatter) ! ! --- DESCRIPTION -------------------------------------------------------- ! ! Converts data read from FITS file into desired format for QDP file ! ! --- VARIABLES ---------------------------------------------------------- ! IMPLICIT NONE integer maxrad, maxpred, nrad, npred, i, j, chatter, ierr character(40) subinfo real rad_lo(maxrad), rad_hi(maxrad), prad_lo(maxpred) real rad_mean(maxrad), del_rad(maxrad), prad_hi(maxpred) real prad_mean(maxpred), pdel_rad(maxpred) logical theo_pres, data_pres ! ! --- VARIABLE DIRECTORY --- ! ! maxrad int : Maximum array size for observed psf data ! maxpred int : Maximum array size for predicted psf data ! nrad int : Counter for observed psf data ! npred int : Counter for predicted psf data ! rad_lo real : Array of lower edge of observed radial bins ! rad_hi real : Array of upper edge of observed radial bins ! prad_lo real : Array of lower edge of predicted model bins ! prad_hi real : Array of upper edge of predicted model bins ! rad_mean real : Array of mean of observed radial bins ! del_rad real : Array of half-width of observed radial bins ! prad_mean real : Array of mean predicted model bins ! pdel_rad real : Array of half-width of observed radial bins ! ! --- AUTHORS/MODifICATION HISTORY --- ! ! Rehana Yusaf (Feb 1993) ! Rehana Yusaf (1994 Jan 31) 1.0.1; add data_pres ! ! Banashree Mitra Seifert (1996, Jan) 1.1.0: ! . Screen display is used ! wtinfo ! ------------------------------------------------------------------------- character(5) version parameter (version = '1.1.0') character(10) subname !- ! -------------- USER INFO -------------------------------------- subname = 'cqdp_conv' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! --- (OBSERVED) DATA MANIPULATION -- if (data_pres) then do i = 1, nrad rad_mean(i) = (rad_lo(i) + rad_hi(i)) / float(2) del_rad(i) = (rad_hi(i) - rad_lo(i)) / float(2) enddo endif ! ! --- (THEORETICAL) DATA MANIPULATION -- ! if (theo_pres) then do j = 1, npred prad_mean(j) = (prad_lo(j) + prad_hi(j)) / float(2) pdel_rad(j) = (prad_hi(j) - prad_lo(j)) / float(2) enddo endif return end ! --------------------------------------------------------------------- ! END OF SUBROUTINE CQDP_CONV ! --------------------------------------------------------------------- !+CQDP_WT subroutine cqdp_wt(fits_psf, qdp_psf, rad_mean, del_rad, cts, err, & prad_mean, pdel_rad, pred, maxrad, maxpred, maxtheta, nrad, & npred, theo_pres, data_pres, rescale, pix_size, backgrnd, & rpsf_pres, reef_pres, ierr, chatter, killit) ! ! --- DESCRIPTION ---------------------------------------------------- ! ! This routine writes to output file in QDP format ! ! --- VARIABLES ------------------------------------------------------ ! IMPLICIT NONE character(180) qdp_psf, fits_psf integer maxrad, maxpred, maxtheta, nrad, npred, chatter, ierr real rad_mean(maxrad), del_rad(maxrad) real cts(maxrad, maxtheta, 1), err(maxrad, maxtheta, 1) real prad_mean(maxpred), pdel_rad(maxpred) real pred(maxpred, maxtheta, 1) real rescale real pix_size, backgrnd logical theo_pres, data_pres, killit logical rpsf_pres, reef_pres ! ! --- INTERNALS --- ! character(40) subinfo character(80) desc1, desc2, desc3 character(250) header integer j ! ! --- VARIABLE DIRECTORY --- ! ! maxrad int : Maximum array size for observed psf data ! maxpred int : Maximum array size for predicted psf data ! nrad int : Counter for observed psf data ! npred int : Counter for predicted psf data ! cts real : Array of observed PSF in counts ! pred real : Array of theoretical PSF ! rad_mean real : Array of mean of observed radial bins ! del_rad real : Array of half-width of observed radial bins ! prad_mean real : Array of mean predicted model bins ! pdel_rad real : Array of half-width predicted model bins ! rescale real : Scalar rescaling factor for pred ONLY ! chatter int : Chattiness flag (<5 quiet,>5 normal,>20 noisy) ! ! --- AUTHORS/MODifICATION HISTORY --->>> ! ! Rehana Yusaf (1993 Feb) ! Rehana Yusaf (1994 Jan 31) 1.0.1; add data_pres ! Rehana Yusaf (1994 Sept 12) 1.0.2; fix y-axis label ! Ian M George (1995 Apr 05) 2.0.0; added rescale to passed parameters ! Rehana Yusaf (1995 April 25) 2.0.1; add clobber ! Banashree Mitra Seifert(1995 December) 2.1.0; ! . modified to accomodate for REEF data, in which case, ! label for y-axis is dimensionless ! . added logical parameters rpsf_pres, reef_pres ! . Screen display used ! wtinfo, ! --------------------------------------------------------------------- character(5) version parameter (version = '2.1.0') character(8) subname !- ! --- USER INFORMATION --- subname = 'cqdp_wt' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! --- OPENING OUTPUT FILE/WRITING HEADER --- call opasci(10, qdp_psf, 2, 80, killit, chatter, ierr) header = '! QDP FILE : ' // qdp_psf write(10, 100) header header = '! CONVERTED FROM FITS FILE : ' // fits_psf write(10, 100) header desc1 = ' Read serr 1,2' desc2 = ' la x Radius (arcmin)' write(10, 100) desc1 write(10, 100) desc2 if (rpsf_pres) then desc3 = ' la y Counts per sq arcmin' write(10, 100) desc3 endif if (reef_pres) then desc3 = ' la y eef ' write(10, 100) desc3 endif ! ! <<<--- WRITING OBSERVED PSF --->>> ! if (data_pres) then header = '! ----------------------------' write (10, 100) header header = '! Observed Radial Profile Data' write (10, 100) header header = '! ----------------------------' write (10, 100) header header = '! Radial Mean, Delta Rad, Counts & Statistical Error :' write(10, 100) header do j = 1, nrad write (10, *) rad_mean(j), del_rad(j), cts(j, 1, 1), err(j, 1, 1) enddo endif ! ! <<<--- WRITING THEORETICAL PSF --->>> ! if (theo_pres) then header = '! -------------------------------' write (10, 100) header header = '! Theoretical Radial Profile Data' write (10, 100) header header = '! -------------------------------' write (10, 100) header write (10, *)' no no no no' if(rescale.NE.1.0) then write(header, '(a,g12.5)')& 'rescaling predicted curve by factor: ', rescale call wtinfo(chatter, 10, 1, header) write(header, '(a,g12.5,a)')& '! Theoretical PSF curve rescaled by factor: ', & rescale, ' (at request of user)' write (10, 100) header endif header = '! Radial Mean, Delta Rad, Theoretical PSF :' write (10, 100) header do j = 1, npred ! write (10,*) prad_mean(j),pdel_rad(j), ! & pred(j,1,1)*rescale,' 0.0' if (backgrnd.NE.0.0) then pred(j, 1, 1) = pred(j, 1, 1) - backgrnd / pix_size**2 endif pred(j, 1, 1) = pred(j, 1, 1) * rescale if (backgrnd.NE.0.0) then pred(j, 1, 1) = pred(j, 1, 1) + backgrnd / pix_size**2 endif write (10, *) prad_mean(j), pdel_rad(j), & pred(j, 1, 1), ' 0.0' enddo endif write(10, *)'log y' close(unit = 10) 100 format(A80) return end ! ! <<<--- end OF SUBROUTINE CQDP_WT --->>> ! !+WTCOLS subroutine wtcols(ncols, columns, maxcol, chatter) ! --- DESCRIPTION ------------------------------------------------- ! ! This routine prints out a character array of column header names, ! using fcecho. ! ! --- VARIABLES --------------------------------------------------- ! IMPLICIT NONE character(100) subinfo character(80) info integer ncols, maxcol, i, chatter character(8) curcol, columns(maxcol) ! ! --- VARIABLE DIRECTORY --- ! ! info char : Comment string ! ncols int : No. of Columns ! maxcol int : Array dimension ! columns char : Array containing column names ! curcol char : current column name ! ! --- CALLED ROUTINES --- ! ! subroutine FCECHO : FTOOLS library routine to write to screen ! ! --- LINKING AND COMPILATION --- ! ! Link with FTOOLS ! ! --- AUTHORS/MODIFICATIONS --- ! ! Rehana Yusaf ( Feb 1993) 1.0.0: ! ! Banashree Mitra Seifert (Jan 1996) 1.1.0: ! . Introduced screen display subroutine ! wtinfo ! -------------------------------------------------------------------- character(5) version parameter (version = '1.1.0') character(7) subname !- subname = 'wtcols' subinfo = 'using' // subname // version call wtinfo(chatter, 10, 1, subinfo) ! ! <<<------>>> ! info = 'The Following Columns are present in the input file :' call wtinfo(chatter, 10, 1, info) do i = 1, ncols curcol = columns (i) call wtinfo(chatter, 10, 1, info) enddo return end ! ---------------------------------------------------------------------+ ! END OF SUBROUTINE WTCOLS ! ---------------------------------------------------------------------+ !+RPSF_DECIDE ! ---------------------------------------------------------------------+ subroutine decide(infile, nset, hdu_sav, chatter) ! --------------- DESCRIPTION OF RPSF_DECIDE --------------------------+ ! ! This subroutines reads the input file to be converted to QDP format ! and then determines which format is it in (that is, the RPSF format ! or the REEF format. ! ! --------------- ROUTINES CALLED -------------------------------------+ ! MVEXT --> opens input file and moves to the desired extension either ! by EXTNUM or by HDUCLAS/EXTNAME ! ! --------------- DECLARE VARIABLES -----------------------------------+ implicit none character(180) infile character(20) hdu_sav(9, 50) integer nset, chatter ! --------------- INTERNAL VARIABLES ----------------------------------+ integer iunit, ninstr, nsearch, status, next(50), errflg character(20) extname character(20) instr(50), outhdu(9, 50), extnames(50), outver(9, 50) character(120) subinfo ! --------------- AUTHORS/MODifICATIONS ------------------------------+ ! ! Banashree Mitra Seifert (Nov. 24, 1995) ! ! --------------------------------------------------------------------+ character(5) version parameter (version = '1.0.0') character(12) subname ! --------------------------------------------------------------------+ !- subname = 'rpsf_decide' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) status = 0 ninstr = 2 instr(1) = 'RESPONSE' instr(2) = 'RPRF' nsearch = 50 call mvext(0, infile, iunit, ninstr, instr, nsearch, next, & outhdu, extnames, outver, extname, errflg, chatter) if ((errflg .eq. 0) .or. (errflg .eq. 2)) then subinfo = 'has psf extension' call wtinfo(chatter, 20, 4, subinfo) nset = nset + 1 hdu_sav(nset, 1) = outhdu(1, 1) hdu_sav(nset, 2) = outhdu(2, 1) else subinfo = 'does not have psf extension' call wtinfo(chatter, 20, 3, subinfo) call ftclos(iunit, status) subinfo = ' closing PSF file' call wtferr(subname, version, status, subinfo) endif status = 0 ninstr = 2 instr(1) = 'RESPONSE' instr(2) = 'REEF' nsearch = 50 call mvext(0, infile, iunit, ninstr, instr, nsearch, next, & outhdu, extnames, outver, extname, errflg, chatter) if ((errflg .eq. 0) .or. (errflg .eq. 2)) then subinfo = 'has reef extension' call wtinfo(chatter, 1, 4, subinfo) nset = nset + 1 hdu_sav(nset, 1) = outhdu(1, 1) hdu_sav(nset, 2) = outhdu(2, 1) else subinfo = 'does not have reef extension' call wtinfo(chatter, 20, 3, subinfo) call ftclos(iunit, status) subinfo = ' closing REEF file' call wtferr(subname, version, status, subinfo) endif return end ! ---------------------------------------------------------------------- ! END OF DECIDE ! ---------------------------------------------------------------------- !+CQDP_REEF subroutine cqdp_reef(fits_psf, pred_psf, rad_lo, rad_hi, reef, & reef_err, prad_lo, prad_hi, pred_reef, maxrad, & maxpred, maxtheta, nrad, npred, theo_pres, & data_pres, pix_size, backgrnd, ierr, chatter) ! ! -------------------------- DESCRIPTION -------------------------- ! ! Reading REEF FITS file using FITSIO. ! ! --- VARIABLES ----------------------------------------------------- ! IMPLICIT NONE character(180) fits_psf, pred_psf integer maxrad, maxpred, nrad, npred, chatter, ierr, maxtheta real rad_lo(maxrad), rad_hi(maxrad) real reef(maxrad, maxtheta, *), reef_err(maxrad, maxtheta, *) real prad_lo(maxpred), prad_hi(maxpred) real pred_reef(maxpred, maxtheta, *) logical theo_pres, data_pres ! ! ---------------------- INTERNAL VARIABLES -------------------------- ! character(40) comm character(80) subinfo character(200) errinfo character(20) instr(50), outhdu(9, 50), outver(9, 50), extname(50) character(20) extnames(9, 50) character(8) telescop, instrume real pix_size, backgrnd integer status, iunit character(16) radunit, thetaunit, energunit, reefunit character(16) hduclas3 integer ntheta, nenerg real theta_lo(1), theta_hi(1), energ_lo(1), energ_hi(1) real perr(9092, 1, 1), parea_wgt(9092, 1, 1), area_wgt(300, 1, 1) logical qerror, qarea integer nsearch, ninstr, next(50) integer errflg ! ! --- VARIABLE DIRECTORY --- ! ! Arguments ... ! ! fits_psf char : Name of FITS format EEF file ! maxen int : Maximum array size for observed eef data ! maxpred int : Maximum array size for predicted eef data ! nrad int : Counter for observed eef data ! npred int : Counter for predicted eef data ! rad_lo real : Array of lower edge of observed radial bins ! rad_hi real : Array of upper edge of observed radial bins ! cts real : Array of observed radial profile, in counts ! err real : Array of statistical errors ! prad_lo real : Array of lower edge of predicted model bins ! prad_hi real : Array of upper edge of predicted model bins ! pred real : Array of theoretical EEF ! ! Internals ... ! ! errmess char : Error message text ! errstr char : Error text for this routine ! version char : Subroutine version ! status int : Error flag for FITSIO call ! iunit int : Fortran unit number for file ! ! --- CALLED ROUTINES --- ! ! subroutine FTOPEN : FITSIO routine to open file ! subroutine FTCLOS : FITSIO routine to close file ! subroutine MVEXT : moves to desired extension ! subroutine RDEEF1 : Reads observed EEF extension in FITS file ! subroutine WT_FERRMSG : Writes FITSIO error text if neccesary ! ! ---------------- AUTHORS/MODifICATION HISTORY ----------------------- ! ! Banashree Mitra Seifert (December 1995) ! ! ----------------------------------------------------------------------- character(5) version parameter (version = '1.0.0') character(10) subname !- ! ----------------------------------------------------------------------- subname = 'cqdp_reef' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! ----------- OPENING PSF FILE ------------- ! status = 0 ninstr = 3 instr(1) = 'RESPONSE' instr(2) = 'REEF' instr(3) = 'OBSERVED' nsearch = 50 call mvext (0, fits_psf, iunit, ninstr, instr, nsearch, next, & outhdu, extnames, outver, extname, errflg, chatter) if (errflg .ne. 0) then data_pres = .false. call ftclos(iunit, status) subinfo = ' closing EEF file' call wtferr(subname, version, status, subinfo) else data_pres = .true. endif ! ----------- READING OBSERVED PSF EXTENSION ------------- if(data_pres) then call rdeef1(iunit, hduclas3, nrad, rad_lo, rad_hi, radunit, & ntheta, theta_lo, theta_hi, thetaunit, nenerg, & energ_lo, energ_hi, energunit, reef, qerror, & reef_err, & reefunit, qarea, area_wgt, telescop, instrume, & maxrad, maxtheta, ierr, chatter) if (ierr.NE.0) then call ftclos(iunit, status) subinfo = ' closing obs_eef file' call wtferr(subname, version, status, subinfo) return endif ! ! ------------------- READ PIXSIZE ------------------------- ! if(pix_size .ne. 0.) then status = 0 call ftgkye(iunit, 'PIXSIZE', pix_size, comm, status) errinfo = ' reading PIXSIZE !' call wtferr(subname, version, status, errinfo) pix_size = pix_size * 60 endif ! ------------------- READ BACKGROUND -------------- if(backgrnd .ne. 0.0) then status = 0 call ftgkye(iunit, 'BACKGRND', backgrnd, comm, status) errinfo = ' reading BACKGRND' call wtferr(subname, version, status, errinfo) if (status.NE.0) then errinfo = 'If rescaling theo curve, background required' call wtferr(subname, version, status, errinfo) endif endif ! ------------- Closing the input file ----------------------- status = 0 call ftclos(iunit, status) errinfo = ' closing ' // fits_psf call wtferr(subname, version, status, errinfo) endif ! ----------------- closed the if data present statement ! --------------- READING THEORETICAL PSF EXTENSION -------- 100 ierr = 0 ninstr = 3 instr(1) = 'RESPONSE' instr(2) = 'REEF' instr(3) = 'PREDICTED' nsearch = 50 call mvext (0, pred_psf, iunit, ninstr, instr, nsearch, & next, outhdu, extnames, outver, extname, & errflg, chatter) if (errflg .ne. 0) then theo_pres = .false. call ftclos(iunit, status) subinfo = ' closing EEF file' call wtferr(subname, version, status, subinfo) return endif ! --------------------- READ THEO DATA -------------------------------+ theo_pres = .true. status = 0 call rdeef1(iunit, hduclas3, npred, prad_lo, prad_hi, radunit, & ntheta, theta_lo, theta_hi, thetaunit, nenerg, & energ_lo, energ_hi, energunit, pred_reef, & qerror, perr, & reefunit, qarea, parea_wgt, telescop, instrume, & maxpred, maxtheta, ierr, chatter) if (ierr.EQ.1) then theo_pres = .false. ierr = 0 endif ! <<<--- CLOSING FILE --->>> status = 0 call ftclos(iunit, status) errinfo = ' closing ' // pred_psf call wtferr(subname, version, status, errinfo) 200 if ((.NOT.theo_pres).AND.(.NOT.data_pres)) then ierr = 5 errinfo = ' predicted and observed data not present' call wtinfo(chatter, 1, 1, errinfo) endif return end ! ---------------------------------------------------------------- ! end OF SUBROUTINE CQDP_REEF ! ----------------------------------------------------------------