!+st2rpsf subroutine st2rpf ! ! --- DESCRIPTION --------------------------------------------------------- ! ! ST2RPSF reads a STW format FITS file and writes required data to an output ! file in CALRPSF format FITS file, that is an OGIP standard format for ! Radial profiles. ! ! ---------------------------- VARIABLES ------------------------------+ IMPLICIT NONE character(180) infil, outfil character(100) source, subinfo character(8) telstw, inststw, telescope, instrume integer maxrad, maxtheta, chanmin, chanmax real, allocatable :: p_cts_arc2(:), p_err_arc2(:) real, allocatable :: p_rad_lo(:), p_rad_hi(:) real, allocatable :: p_rad1(:), p_rad2(:) real, allocatable :: p_net(:), p_neterr(:) real, allocatable :: p_area_wgt(:) integer, allocatable :: p_n_pix(:) real arc_pix, pixsize, bkgd, sumrcts, bkgd_rad, sumtcts integer npts, chatter, extno integer ineed, iget, status, ierr logical cont, calc, killit ! ! --- VARIABLE DIRECTORY --- ! ! infil char : STW file name (user defined) ! outfil char : CALRPSF file name (user defined) ! maxrad int : Array dimensions ! net real : Array of observed radial profile in counts ! neterr real : Array of errors on net ! cts_arc2 real : Array of observed radial profile in counts/arcmin^2 ! err_arc2 real : Array of error on cts_arc2 ! rad_lo real : Array of lower edge of radial bins (in arcmin) ! rad_hi real : Array of upper edge of radial bins (in arcmin) ! rad1 real : Array of lower edge of radial bins (in pixels) ! rad2 real : Array of upper edge of radial bins (in pixels) ! n_pix int : Array of No of pixels in bins ! area_wgt real : Area weighting factor ! arc_pix real : arcmins per pixel ! telstw char : telescope name ! inststw char : Instrument/Detector name ! pixsize real : Pixsize in degrees ! npts int : Counter for observed psf data ! chatter int : Chattiness flag (<5 quiet,>5normal,>20 noisy) ! subinfo char : Routine info for user ! cont logical : IF false program stops ! ! ---------------------- CALLED ROUTINES ------------------------------ ! ! subroutine RPSF_GP : Reads user defined parameters using XPI ! subroutine RPSF_RSTW : Reads FITS format STW file ! subroutine RPSF_CONV : Converts data into desired format ! subroutine RPSF_WT : Writes data into calrpsf format output file ! ! ------------------ AUTHORS/MODIFICATION ----------------------------- ! ! Rehana Yusaf (1993 Feb) ! Rehana Yusaf (1993 March 3) ! Rehana Yusaf (1993 August 17) Ver 1.0.2; update rpsf_gp, rpsf_src ! more error checking ! Rehana Yusaf (1993 September 3) 1.0.3; New more general radial profile ! writing routine is used WT_RPSF1993a ! Rehana Yusaf (1994 Jan 13) 1.0.4; WT_RPSF1993a has been renamed to ! WTRPF1, and HDUCLAS has been added. ! New parameters - CHANMIN, CHANMAX ! TELESCOPE and INSTRUMENT and BKGD ! added ! Rehana Yusaf (1994 Feb 3) 1.0.5; Write SUMRCTS to outfile ! Rehana Yusaf (1995 Jan 13) 1.0.6; update _gp to no longer use defval ! Rehana Yusaf (1995 Feb 21) 1.0.7; update to add option of calculating ! the bkgd ! Rehana Yusaf (1995 April 25) 1.0.8; update to add clobber ! Banashree Mitra Seifert (1995 Dec) 2.0.0; modified to ! . filename = character(180) ! . HDUCLAS introduced ! . replaces calls to FNDHDU/EXT,FTMRHD,FTMAHD by MVEXT ! . DMA introduced ! . error checking modified ! . screen display subroutines used ! wtbegm,wtendm,wtinfo,wterrm,wtwarm ! Peter D Wilson (1998 Jul 01) 2.0.1: ! . Updated for new FCPARS behavior ! Bryan K Irby (Jan 14, 2018) 2.1.0: ! . Replaced udmget/udmfre with allocate/deallocate ! MFC (Apr 16, 2020) 2.1.1 ! . f90 version ! ! ---------------------------------------------------------------------- character(5) version parameter (version = '2.1.1') character(40) taskname character(7) subname !ccc COMMON/task/taskname taskname = 'st2rpsf' subname = 'st2rpsf' cont = .true. ! ! ---------------- GET PARAMETERS ------------------------------------ ! call rpsf_gp(infil, extno, telescope, instrume, chanmin, chanmax, & bkgd, outfil, calc, bkgd_rad, cont, chatter, killit) call wtbegm(taskname, version, chatter) IF (.NOT.cont) THEN go to 100 ENDIF maxrad = 1000 maxtheta = 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_rad1(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * 4 status = 0 allocate(p_rad2(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * 4 status = 0 allocate(p_net(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * 4 status = 0 allocate(p_neterr(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * 4 status = 0 allocate(p_cts_arc2(maxrad * maxtheta * 1), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * 1 * 4 status = 0 allocate(p_err_arc2(maxrad * maxtheta * 1), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * 1 * 4 status = 0 allocate(p_area_wgt(maxrad * maxtheta * 1), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + maxrad * maxtheta * 1 * 4 status = 0 allocate(p_n_pix(maxrad), stat = status) if (status .ne. 0) then goto 50 endif iget = iget + 2 * maxrad * 4 50 ineed = 8 * maxrad * 4 + 3 * maxrad * maxtheta * 1 * 4 write(subinfo, '(a,i10)')'Memory requirements for this task =', & ineed call wtinfo(chatter, 15, 2, subinfo) write(subinfo, '(a,i10)')'Total memory dynamically allocated =', & 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 ! ---------------------------------------------------------------------+ ! End of DMA allocation ! ---------------------------------------------------------------------+ ! ------------------------ READ STW file ------------------------------+ ! call rpsf_rstw(infil, extno, p_n_pix, p_net, & p_neterr, p_rad1, p_rad2, & source, npts, maxrad, telstw, inststw, & cont, chatter) IF (.NOT.cont) THEN go to 100 ENDIF ! ! ------ Check that telescope/instrument values that are read ---------- ! ------------ from the file and defined by user match ----------------- ! IF (telstw.NE.'UNKNOWN') THEN IF (telstw.NE.telescope) THEN subinfo = ' Telescope defined by user : ' // telescope call wtinfo(chatter, 1, 1, subinfo) subinfo = ' Telescope value read from STW file : ' // telstw call wtinfo(chatter, 5, 2, subinfo) subinfo = 'mismatch in defining telescope!' call wtwarm(subname, version, chatter, 5, subinfo) subinfo = 'User defined value used' call wtinfo(chatter, 5, 1, subinfo) ENDIF ENDIF IF (inststw.NE.'UNKNOWN') THEN IF (inststw.NE.instrume) THEN subinfo = ' Instrument defined by user : ' // instrume call wtinfo(chatter, 1, 1, subinfo) subinfo = ' Instrument value read from STW file : ' // inststw call wtinfo(chatter, 1, 1, subinfo) subinfo = 'mismatch in defining instrument!' call wtwarm(subname, version, chatter, 5, subinfo) subinfo = 'User defined value used' call wtinfo(chatter, 5, 1, subinfo) ENDIF ENDIF telstw = telescope inststw = instrume ! ---------------------------------------------------------------------+ ! CONVERT TO OGIP RPSF FORMAT ! ---------------------------------------------------------------------+ call rpsf_conv(p_rad1, p_rad2, p_rad_lo, & p_rad_hi, p_net, p_neterr, & p_cts_arc2, p_err_arc2, source, & arc_pix, pixsize, sumrcts, sumtcts, npts, & maxrad, maxtheta, p_n_pix, p_area_wgt, & bkgd, calc, bkgd_rad, cont, chatter) IF (.NOT.cont) THEN goto 100 ENDIF ! ---------------------------------------------------------------------+ ! WRITE RPSF FILE ! ---------------------------------------------------------------------+ call rpsf_wt(outfil, infil, version, p_rad_lo, & p_rad_hi, p_cts_arc2, p_err_arc2, & p_area_wgt, pixsize, npts, maxrad, maxtheta, & telstw, inststw, bkgd, chanmin, chanmax, & sumrcts, sumtcts, cont, chatter, killit) IF (.NOT.cont) 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_rad1, stat = status) status = 0 deallocate(p_rad2, stat = status) status = 0 deallocate(p_net, stat = status) status = 0 deallocate(p_neterr, stat = status) status = 0 deallocate(p_cts_arc2, stat = status) status = 0 deallocate(p_err_arc2, stat = status) status = 0 deallocate(p_area_wgt, stat = status) status = 0 deallocate(p_n_pix, 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 ST2RPSF ! ---------------------------------------------------------------------- !+RPSF_GP ! ---------------------------------------------------------------------+ subroutine rpsf_gp(infil, extno, telescope, instrume, chanmin, & chanmax, bkgd, outfil, calc, bkgd_rad, cont, & chatter, killit) ! -------------------------- DESCRIPTION ------------------------------ ! ! This routine gets the user defined parameters, that is the input and ! output file names . In addition the Chatter flag is obtained. ! ! -------------------------- VARIABLES ------------------------------- IMPLICIT NONE character*(*) infil, outfil character(180) file, strbkgd, ill_files(5) character*(*) telescope, instrume character(80) errmess, subinfo, rdstr, rdstr2 integer chanmin, chanmax, extno, try, n_ill integer err, chatter real bkgd, bkgd_rad logical cont, calc, killit, valfil !logical ext ! ! ------------------------ VARIABLE DIRECTORY ------------------------ ! ! infil char : Name of FITS STW PSF file (user defined) ! outfil char : Name of PSF output file (user defined) ! err int : Error flag ! chatter int : Chattiness flag (<5 quiet,>5normal,>20 noisy) ! version char : Version of subroutine ! errstr char : Error string ! errmess char : Error message ! extno int : Extension number of RPSF data (default is 1) ! cont logical : False, If infil does not exist ! ext logical : True if infil exists ! ! ---------------------- COMPILATION AND LINKING --------------------- ! ! Link with HOST and FTOOLS ! ! --- AUTHORS/MODIFICATION HISTORY --- ! ! Rehana Yusaf (1993 Feb) ! Rehana Yusaf (1993 August 17) Ver 1.0.1; Remove leading blanks from ! infile and outfile ! Rehana Yusaf (1994 Jan 13) 1.0.2; Add new parameters TELESCOP and ! INSTRUM - also CHANMIN and CHANMAX ! (dependent on telescop and instrume ! values). BKGD is also prompted for. ! Also extname is no longer read ! Rehana Yusaf (1995 Jan 13) 1.0.3; defval is no longer read from parfile ! Rehana Yusaf (1995 Feb 21) 1.0.5; read in calc, if calc then code ! calculates the bkgd ! Rehana Yusaf (1995 April 25) 1.0.6; read in clobber ! ! Banashree Mitra Seifert (1996, Jan) 1.1.0: ! . Screen display routines used ! wterrm,wtinfo,wtwarm ! Peter D Wilson (1998 Jul 01) 1.1.1: ! . Drop INQUIRE test on infil ! ----------------------------------------------------------------------- character(5) version parameter (version = '1.1.1') character(7) subname !- ! ----------------------------------------------------------------------- subname = 'rpsf_gp' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! ! ---------------- READING INPUT FILE PARAMETER --------------------------- err = 0 call uclgst('infil', file, err) IF (err.NE.0) THEN errmess = ' Getting infil parameter' call wterrm(subname, version, errmess) ENDIF call crmvlbk(file) IF (infil.EQ.' ') THEN errmess = ' Must enter infile name !' call wterrm(subname, version, errmess) cont = .false. returN ENDIF err = 0 call fcpars(file, infil, extno, err) IF (extno.LE.0) THEN extno = 1 ENDIF ! PDW 7/1/98: Don't bother! Let FTOPEN determine if file exists ! INQUIRE(FILE=infil,EXIST=ext) ! IF (.NOT.ext) THEN ! errmess = 'Infile does not exist !' ! call wterrm(subname,version,errmess) ! cont = .false. ! return ! ENDIF ! ------------------- READ TELESCOPE NAME ------------------------------- call uclgst('telescope', rdstr, err) IF (err.NE.0) THEN errmess = ' getting telescope parameter' call wterrm(subname, version, errmess) ENDIF call crmvlbk(rdstr) call rmvexsp(rdstr, rdstr2) telescope = rdstr2 ! --------------------- READ INSTRUMENT NAME ------------------------------ call uclgst('instrume', rdstr, err) IF (err.NE.0) THEN errmess = ' getting instrume parameter' call wterrm(subname, version, errmess) ENDIF call crmvlbk(rdstr) call rmvexsp(rdstr, rdstr2) instrume = rdstr2 call ftupch(instrume) call ftupch(telescope) ! ------------ READ CHANMIN and CHANMAX if ROSAT PSPC -------------------- IF ((telescope.EQ.'ROSAT').AND.(instrume(1:4).EQ.'PSPC')) THEN call uclgsi('chanmin', chanmin, err) If (err.NE.0) THEN errmess = ' getting chanmin parameter' call wterrm(subname, version, errmess) ENDIF call uclgsi('chanmax', chanmax, err) IF (err.NE.0) THEN errmess = ' getting chanmax parameter' call wterrm(subname, version, errmess) ENDIF ENDIF ! ------------------ READ BKGD ------------------------------------------ calc = .false. try = 0 90 call uclgst('bkgd', strbkgd, err) try = try + 1 IF (err.NE.0) THEN errmess = ' getting bkgd parameter' call wterrm(subname, version, errmess) ENDIF call crmvlbk(strbkgd) call ftupch(strbkgd) IF (strbkgd(1:1).EQ.'C') THEN calc = .true. call uclgsr('bkgd_rad', bkgd_rad, err) IF (err.NE.0) THEN errmess = ' getting bkgd_rad parameter' call wterrm(subname, version, errmess) ENDIF IF (telescope.EQ.'ROSAT') THEN IF (instrume(1:4).EQ.'PSPC') THEN IF (bkgd_rad.LE.2) THEN errmess = 'The background estimate could be'& // ' contaminated by the PSF' call wtwarm(subname, version, chatter, 5, errmess) errmess = ' using a radius < 2 arcmins for the calc' call wtwarm(subname, version, chatter, 5, errmess) ENDIF ELSEIF (instrume(1:3).EQ.'HRI') THEN IF (bkgd_rad.LE.(0.5)) THEN errmess = 'The background estimate could be'& // ' contaminated by the PSF' call wtwarm(subname, version, chatter, 5, errmess) errmess = ' using a radius < 1/2 arcmin for the calc' call wtwarm(subname, version, chatter, 5, errmess) ENDIF ENDIF ENDIF ELSE read(strbkgd, *, IOSTAT = err) bkgd IF (err.NE.0) THEN errmess = 'parsing bkgd parameter' call wterrm(subname, version, errmess) cont = .false. ENDIF ENDIF ! ---------------- READING OUTPUT FILE PARAMETER ------------------------ call uclgst('outfil', outfil, err) IF (err.NE.0) THEN errmess = ' getting outfile parameter' call wterrm(subname, version, errmess) ENDIF call crmvlbk(outfil) IF (outfil.EQ.' ') THEN errmess = ' must enter outfile name ! ' call wterrm(subname, version, errmess) cont = .false. return ENDIF ! GET CLOBBER call uclgsb('clobber', killit, err) IF (err.NE.0) THEN killit = .false. errmess = 'getting clobber value' call wterrm(subname, version, errmess) ENDIF ! OUTFILE VALIDATION ! ----------------------------------------------------------------------+ n_ill = 0 call ck_file(outfil, ill_files, n_ill, valfil, killit, chatter) IF (.NOT.valfil) THEN cont = .false. return ENDIF ! -------------------- GETTING CHATTER FLAG ---------------------------- call uclgsi('chatter', chatter, err) IF (err.NE.0) THEN errmess = ' getting chatter parameter' call wterrm(subname, version, errmess) ENDIF return end ! ---------------------------------------------------------------------- ! END OF SUBROUTINE RPSF_GP ! ---------------------------------------------------------------------- !+RPSF_RSTW ! ---------------------------------------------------------------------+ subroutine rpsf_rstw(infil, extno, n_pix, net, neterr, rad1, & rad2, source, npts, maxrad, telstw, inststw, & cont, chatter) ! ! --- DESCRIPTION ---------------------------------------------------- ! ! This routine reads STW FITS file. It uses FITSIO to read columns ... ! ! PIXELS : Column of No. of pixels in each bin ! NET : Column of observed radial profile in counts ! NETERR : Column of errors on net ! RAD1 : Column of lower edge of radial bins, in pixels ! RAD2 : Column of upper edge of radial bins, in pixels ! ! In addition to the standard keywords, it also reads ... ! ! TELESCOPE : Mission/Telescope name ! INSTRUME : Instrument/Detector name ! SOU_A : Which defines Source information in string format. ! The inner and outer radi , and No. of steps is extracted ! from this string, in order to later calculate pixelsize. ! ! ----------------------------- VARIABLES ----------------------------------- ! IMPLICIT NONE character*(*) infil !character(180) filename character*(*) source character*(*) telstw, inststw integer maxrad, npts, chatter integer extno real net(*), neterr(*), rad2(*), rad1(*) integer n_pix(*) logical cont ! real net(maxrad),neterr(maxrad) ! real rad2(maxrad),rad1(maxrad) ! integer*4 n_pix(maxrad) ! ! -------------------------- INTERNAL VARIABLE ------------------------------ ! character(70) sou1, subinfo character(40) comm, extinfo, errstr character(30) sou_comm integer maxcol, clenact parameter (maxcol = 14) character(8) columns(maxcol) real enull integer status, block, htype, ncols, colnum integer felem, frow, iunit, find_nstep !integer extnum logical anyflg ! --------------------------- VARIABLE DIRECTORY ----------------------- ! ! maxrad int : Array dimensions ! region int : Array of region Nos ! net real : Array of observed radial profile in counts ! neterr real : Array of errors on net ! rad2 real : Array of upper edge of radial bins (in pixels) ! n_pix int : Array of No of pixels in bins ! npts int : Counter for observed psf data ! telstw char : Used to read Telescope/Mission name ! inststw char : Used to read Instrument/Detector name ! chatter int : Chattiness flag ( >20 noisy) ! extname char : Name of Extname keyword in STW file ! cont : False if column is not found, program returns to main ! ! Internals ... ! ! errmess char : Error message text ! errstr char : Error text for this routine ! extinfo char : More error information ! version char : Subroutine version ! enull real : Used to represent undefined values in FITS file ! status int : Error flag for FITSIO call ! iunit int : Fortran unit number for file ! block int : FITSIO record blocking factor ! htype int : Type of header unit in FITS file ! felem int : First pixel of element array in FITS file ! frow int : Beginning row No. in FITS file ! inull int : Used to represent undefined values in FITS file ! colnum int : No. of column ! i,j int : Counters for loops ! anyflg logical: True if any data undefined in FITS file ! upstw char : Extension name from file in upper case (for comparison) ! upext char : Extname in upper case ! extfind logical : True if required extension is found ! ! --- CALLED ROUTINES --- ! ! subroutine FTOPEN : FITSIO routine to open FITS file ! subroutine FTCLOS : FITSIO routine to close FITS file ! subroutine FTMAHD : FITSIO routine to move to extension header ! subroutine FTGKYs : FITSIO routine to read extension header keyword, ! where the s is for string ! subroutine FTGKYj : FITSIO routine to read extension header keyword, ! where the j, is for an integer ! subroutine FTGKNS : FITSIO routine to read extension header keyword, ! where a rootstring is given, thus an array of ! keywords can be read ! subroutine FTGCVe : FITSIO routine to read columns of data in an ! extension. The e indicates that it is real data. ! ! subroutine FTGERR : FITSIO routine, which gives error text for a ! given error flag. ! ! subroutine FCECHO : FTOOLS library routine to write to screen ! subroutine WT_FERRMSG : Writes Error text if required ! ! --- COMPILATION AND LINKING --- ! ! Link with FITSIO and FTOOLS ! ! --- AUTHORS/MODIFICATION HISTORY --- ! ! Rehana Yusaf (Feb 1993) ! Rehana Yusaf (March 1993) : Add telescope and instrument read ! Rehana Yusaf (22 April 1993) : NET column read from STW file instead ! of cnts/pix, as conversion to ! cnts/arcmin^2 is more accurate using NET. ! This is due to rounding errors in cnts/pix. ! Rehana Yusaf (6 May) : NETERR column read instead of err/pix for ! same reason as above. ! Rehana Yusaf (14 jan 1994) 1.0.2; move to extno as specified by user, extname ! search no longer needed ! Banashree Mitra Seifert (1996, Jan) 1.1.0: ! . variable dimension carried from calling routine ! . Introduced screen display routine ! wtinfo,wterrm,wtferr ! Peter D Wilson (1998 Jul 01) 1.1.1: ! . Drop bogus call to fcpars ! -------------------------------------------------------------------- character(5) version parameter (version = '1.1.1') character(10) subname !- ! ------------------------ USER INFO ---------------------------------- subname = 'rpsf_rstw' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! ------------------------ OPENING STW FILE ------------------------------ ! PDW 7/1/98: No need to call fcpars... done in rpsf_gp ! status=0 ! call fcpars(infil,filename,extnum,status) status = 0 call cgetlun(iunit) call ftopen(iunit, infil, 0, block, status) extinfo = 'opening file : ' // infil IF (status.NE.0) THEN call wtferr(subname, version, status, extinfo) cont = .false. go to 150 ENDIF ! <<<--- MOVING TO DATA EXTENSION --->>> status = 0 call ftmahd(iunit, extno + 1, htype, status) extinfo = ' problem moving to RPSF data extension' call wtferr(subname, version, status, extinfo) IF (status.NE.0) THEN cont = .false. go to 150 ENDIF ! ------------ READING INFORMATION FROM HEADER --------------- call ftgkyj(iunit, 'NAXIS2', npts, comm, status) errstr = 'reading NAXIS2 parameter' call wtferr(subname, version, status, errstr) IF (status.NE.0) THEN cont = .false. go to 150 ENDIF status = 0 call ftgkns(iunit, 'TTYPE', 1, maxcol, columns, ncols, status) errstr = 'reading TTYPE parameter' call wtferr(subname, version, status, errstr) status = 0 call ftgkys(iunit, 'TELESCOP', telstw, comm, status) errstr = 'reading TELESCOP parameter' call wtferr(subname, version, status, errstr) IF (status.EQ.202) THEN telstw = 'UNKNOWN' ENDIF status = 0 call ftgkys(iunit, 'INSTRUME', inststw, comm, status) errstr = 'reading INSTRUME parameter' call wtferr(subname, version, status, errstr) IF (status.EQ.202) THEN inststw = 'UNKNOWN' ENDIF status = 0 call ftgkys(iunit, 'SOU_A', sou1, sou_comm, status) find_nstep = 0 find_nstep = index(sou1(1:), 'n=') IF (find_nstep.EQ.0) THEN source = sou1(:clenact(sou1)) // ''' ' // sou_comm ELSE source = sou1 // sou_comm ENDIF extinfo = 'SOURCE string just read' call wtinfo(chatter, 20, 2, extinfo) call wtinfo(chatter, 20, 2, source) ! ------------------ CHECK TO FIND PIXELS COLUMN -------------------- status = 0 call ftgcno(iunit, .false., 'PIXELS', colnum, status) IF (status.NE.0) THEN extinfo = 'finding PIXELS column ' call wtferr(subname, version, status, extinfo) cont = .false. goto 150 ENDIF ! ---------------------- READING PIXELS COLUMN ------------------------ frow = 1 felem = 1 enull = 0 status = 0 call ftgcvj(iunit, colnum, frow, felem, npts, enull, n_pix, anyflg, & status) extinfo = 'reading pixels column' call wtferr(subname, version, status, extinfo) IF (status.NE.0) THEN cont = .false. go to 150 ENDIF ! ------------------ CHECK TO FIND NET COLUMN ------------------------ status = 0 call ftgcno(iunit, .false., 'NET ', colnum, status) IF (status.NE.0) THEN extinfo = 'finding NET Column ' call wtferr(subname, version, status, extinfo) cont = .false. go to 150 ENDIF ! ----------------- READING NET COLUMN --------------------------------- frow = 1 felem = 1 enull = 0 status = 0 call ftgcve(iunit, colnum, frow, felem, npts, enull, net, anyflg, & status) extinfo = 'reading net column' call wtferr(subname, version, status, extinfo) IF (status.NE.0) THEN cont = .false. go to 150 ENDIF ! --------------- CHECK TO FIND NETERR COLUMN --------------------------- status = 0 call ftgcno(iunit, .false., 'neterr', colnum, status) IF (status.NE.0) THEN extinfo = 'finding NETERR column' call wtferr(subname, version, status, extinfo) cont = .false. goto 150 ENDIF ! ---------------------- READING NETERR COLUMN---------------------------- frow = 1 felem = 1 enull = 0 status = 0 call ftgcve(iunit, colnum, frow, felem, npts, enull, neterr, anyflg, & status) call wtferr(subname, version, status, errstr) IF (status.NE.0) THEN cont = .false. go to 150 ENDIF ! -------------------------- CHECK TO FIND RAD1 COLUMN ------------------- status = 0 call ftgcno(iunit, .false., 'RAD1', colnum, status) IF (status.NE.0) THEN extinfo = 'finding RAD1 column' call wtferr(subname, version, status, extinfo) cont = .false. goto 150 ENDIF ! ---------------- READING RAD1 COLUMN ------------------------------------ frow = 1 felem = 1 enull = 0 status = 0 call ftgcve(iunit, colnum, frow, felem, npts, enull, rad1, anyflg, & status) extinfo = ' reading rad1 column' IF (status.NE.0) THEN extinfo = 'finding RAD1 column' call wtferr(subname, version, status, extinfo) cont = .false. go to 150 ENDIF ! ------------------- CHECK TO FIND RAD2 COLUMN ------------------------ status = 0 call ftgcno(iunit, .false., 'RAD2', colnum, status) IF (status.NE.0) THEN extinfo = ' Error finding RAD2 column' call wtferr(subname, version, status, extinfo) cont = .false. goto 150 ENDIF ! ------------ READING RAD2 COLUMN ------------------------------------- frow = 1 felem = 1 enull = 0 status = 0 call ftgcve(iunit, colnum, frow, felem, npts, enull, rad2, anyflg, & status) extinfo = ' reading rad2 column' IF (status.NE.0) THEN call wtferr(subname, version, status, extinfo) cont = .false. goto 150 ENDIF ! ----------------- CLOSE STW FITS FILE ---------------------------------- 150 status = 0 call ftclos(iunit, status) extinfo = ' closing stw file' IF (status.NE.0) THEN call wtferr(subname, version, status, extinfo) cont = .false. ENDIF call cfrelun(iunit) return end ! ------------------ END OF SUBROUTINE RPSF_RSTW ------------------------ !+WT_COLS ! ---------------------------------------------------------------------+ subroutine wt_cols(ncols, columns, maxcol, chatter) ! --- DESCRIPTION ------------------------------------------------- ! ! This routine prints out a character array of column header names, ! using fcecho. ! ! --- VARIABLES --------------------------------------------------- IMPLICIT NONE 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 ------------------------------------- ! ! Rehana Yusaf ( Feb 1993) ! ! Banashree Mitra Seifert (1995 Dec) 1.1.0: ! . made version no. character(5) instead of *4 (1.0) ! . passed chatter ! . screen display used ! wtinfo ! ------------------------------------------------------------------ character(5) version parameter (version = '1.1.0') !- info = 'The following columns are present in the input file :' call wtinfo(chatter, 20, 1, info) do i = 1, ncols curcol = columns (i) call wtinfo(chatter, 20, 21, curcol) enddo return end ! ---------------------------------------------------------------------- ! END OF SUBROUTINE WT_COLS ! ---------------------------------------------------------------------- !+RPSF_CONV ! ---------------------------------------------------------------------+ subroutine rpsf_conv(rad1, rad2, rad_lo, rad_hi, net, neterr, & cts_arc2, err_arc2, source, arc_pix, pixsize, & sumrcts, sumtcts, npts, maxrad, maxtheta, n_pix, & area_wgt, bkgd, calc, bkgd_rad, cont, chatter) ! ---------------- DESCRIPTION ----------------------------------------+ ! ! This routine converts data into OGIP CALRPSF format ! ! ----------------- VARIABLES ----------------------------------------- IMPLICIT NONE character*(*) source character(80) errinfo, errmess character(200) subinfo integer maxrad, maxtheta real arc_pix, del_rad, pixsize, sumrcts, sumrcts_bkgd, area, pi real bkgd_rad, bkgd, sumpix, sumpix_bkgd, per_bkgd real sumpix2, sumcts, sumtcts real rad1(*), rad2(*), rad_lo(*), rad_hi(*) real net(*), neterr(*), area_wgt(maxrad, maxtheta, 1), theocts real cts_arc2(maxrad, maxtheta, 1), err_arc2(maxrad, maxtheta, 1) integer npts, i, chatter integer n_pix(*) logical findpix, cont, calc ! real rad1(maxrad),rad2(maxrad),arc_pix,del_rad,pixsize ! real rad_lo(maxrad),rad_hi(maxrad),sumrcts,sumrcts_bkgd ! real net(maxrad),area_wgt(maxrad,maxtheta,1),area,pi,bkgd_rad ! real neterr(maxrad),bkgd,sumpix,sumpix_bkgd ! real per_bkgd ! real cts_arc2(maxrad,maxtheta,1),err_arc2(maxrad,maxtheta,1) ! integer*4 n_pix(maxrad) ! ! ----------------------- VARIABLE DIRECTORY ---------------------------- ! ! maxrad int : Array dimensions ! maxtheta int : Array of dimunsions ! rad1 real : Array of lower edge of radial bins (in pixels) ! rad2 real : Array of upper edge of radial bins (in pixels) ! rad_lo real : Array of lower edge of radial bins (in arcmin) ! rad_hi real : Array of upper edge of radial bins (in arcmin) ! net real : Array of counts ! neterr real : Array of errors on net ! cts_arc2 real : Array of counts per arcmin^2 ! err_arc2 real : Array of errors on cts_arc2 ! area_wgt real : Area weighting factor for each bin ! arc_pix real : arcmins per pixel ! del_rad real : del_rad = (outer radius - inner radius)/nbins ! cont logical : False if error occurs ! (in arcmin) ! pixsize real : Pixsize in degrees ! npts int : Counter for No of radial values ! findpix logical: True if still finding arcmins per pixel value ! chatter int : Chattiness flag (>20 noisy) ! subinfo char : Subroutine information ! ! --- CALLED ROUTINES --- ! ! subroutine RPSF_SRC : String manipulation routine, which extracts ! values from source string ! ! --- AUTHORS/MODIFICATION HISTORY --- ! ! Rehana Yusaf (Feb 1993) ! Rehana Yusaf (23 April 1993) : Conversion to counts per arcmin^2 is now ! made using NET and N_PIX arrays. The ! CNT/PIX array was used before (no longer ! read in due to trucation -loss in precision) ! Rehana Yusaf (3 Sept)1.0.3; Area_wgt array added - new rpsf format ! Rehana Yusaf (1994 Feb 3) 1.0.4; Calculate sumrcts ! Rehana Yusaf (1995 Feb 21) 1.0.5; Calculate bkgd if user requests it ! Also add more diagonostics at high chatter ! Banashree Mitra Seifert (1995 Dec)1.1.0: ! . Dimension carried from calling routine ! . filename *180 character ! . Replaced pi=3.141592654 by pi=4*atan(1) ! . screen display used ! wtinfo,wterrm,wtwarm ! Rehana Yusaf (1996 Jan 23) 1.1.1; add theocts, previously ! => sumtcts = sumcts - sumpix * bkgd ! => now sumtcts = theocts - sumpix * bkgd ! ---------------------------------------------------------------- character(5) version parameter (version = '1.1.1') character(10) subname !- ! ------------------------- USER INFO ---------------------------------- pi = 4. * atan(1.) subname = 'rpsf_conv' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! ----------------------- OBTAIN DEL_RAD ----------------------------- call rpsf_src(source, del_rad, cont, chatter) IF (.NOT.cont) THEN errinfo = ' SOU_A string unsuccesfully parsed' call wterrm(subname, version, errinfo) return ENDIF ! --------------- CALCULATE ARCMINS PER PIXEL -------------------------- findpix = .false. i = 1 do WHILE(.NOT.findpix) IF ((rad2(i) - rad1(i)).NE.0) THEN arc_pix = del_rad / (rad2(i) - rad1(i)) findpix = .true. ELSE i = i + 1 ENDIF enddo ! ---------- CALCULATE RAD_LO AND RAD_HI FROM ARC_PIX ------------------ do i = 1, npts rad_lo(i) = rad1(i) * arc_pix rad_hi(i) = rad2(i) * arc_pix enddo IF (bkgd_rad.GE.rad_hi(npts)) THEN errmess = ' inner radius for bkd calculation is > Outer radius' call wterrm(subname, version, errmess) cont = .false. return ENDIF per_bkgd = bkgd_rad / real(rad_hi(npts)) IF (per_bkgd.GE.(0.75)) THEN errmess = 'The inner radius for bkgd calculation is'& // ' > 75% of the total radius' call wtwarm(subname, version, chatter, 9, errmess) ENDIF ! ---------- CONVERT COUNTS AND ERROR INTO ARCMINS**2 UNITS ------------ sumrcts = 0 sumrcts_bkgd = 0 sumpix_bkgd = 0 do i = 1, npts sumcts = sumcts + net(i) IF (calc) THEN IF (rad_lo(i).GE.bkgd_rad) THEN sumpix_bkgd = sumpix_bkgd + n_pix(i) sumrcts_bkgd = sumrcts_bkgd + net(i) ENDIF ENDIF cts_arc2(i, 1, 1) = (net(i) / real(n_pix(i))) * (1 / (arc_pix**2)) err_arc2(i, 1, 1) = (neterr(i) / real(n_pix(i))) * (1 / (arc_pix**2)) enddo sumpix = pi * (rad2(npts)**2 - rad1(1)**2) ! -------------------------------------------------------------------+ ! ------------------- CALCULATE AREA_WGT RATIOS --------------------- sumpix2 = 0 theocts = 0 do i = 1, npts area = pi * ((rad2(i)**2) - (rad1(i)**2)) IF ((area.NE.0).AND.(n_pix(i).NE.0)) THEN area_wgt(i, 1, 1) = real(n_pix(i)) / area ELSE errinfo = ' calculating area_wgt factor !' call wterrm(subname, version, errinfo) ENDIF sumpix2 = sumpix2 + n_pix(i) theocts = theocts + net(i) / area_wgt(i, 1, 1) enddo ! ---------------- CALCULATE PIXSIZE IN DEGREES ------------------------ pixsize = arc_pix / real(60) IF (calc) THEN IF (sumrcts_bkgd.GE.0.0) THEN bkgd = sumrcts_bkgd / sumpix_bkgd ELSE errmess = 'The sum of the counts in the bkgd region is < 0' call wtwarm(subname, version, chatter, 5, errmess) errmess = ' The bkgd has been set to 0' call wtinfo(chatter, 5, 1, errmess) bkgd = 0 ENDIF ENDIF sumtcts = theocts - sumpix * bkgd sumrcts = sumcts - sumpix2 * bkgd IF (calc) THEN IF (sumtcts.LE.0.0) THEN bkgd = (sumcts - 1) / sumpix sumtcts = 1 sumrcts = sumcts - sumpix2 * bkgd subinfo = ' the background value has been calculated such'& // ' that the sum of the' call wtinfo(chatter, 9, 1, subinfo) subinfo = ' theoretical counts is set to 1. If the background'& // ' is calculated in the usual way then sumtcts '& // ' is -ve, this may be due to excluded regions.' call wtinfo(chatter, 9, 1, subinfo) ENDIF ENDIF write(subinfo, '(a,f20.9)') ' calculated bkgd value :', bkgd IF (calc) THEN call wtinfo(chatter, 20, 1, subinfo) ENDIF write(subinfo, '(a,F20.9)')' sum of pixels calculated using'& // ' area of circle:', sumpix call wtinfo(chatter, 20, 2, subinfo) write(subinfo, '(a,F20.9)')' actual sum of pixels :', sumpix2 call wtinfo(chatter, 20, 2, subinfo) IF ((sumpix2 / sumpix).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 wtinfo(chatter, 5, 2, subinfo) subinfo = 'note: this may be due to regions being excluded' call wtinfo(chatter, 5, 2, subinfo) ENDIF write(subinfo, '(a,F20.9)')'total sum of counts :', sumcts call wtinfo(chatter, 20, 2, subinfo) write(subinfo, '(a,F20.9)')'counts in source :', sumrcts 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' call wterrm(subname, version, subinfo) subinfo = ' this may be due to background dominating'& // ' the source' call wtinfo(chatter, 0, 1, subinfo) cont = .false. return ENDIF IF (sumtcts.LT.0) THEN subinfo = 'theoretical Sum of counts is negative' call wterrm(subname, version, subinfo) subinfo = ' this may be due to background dominating'& // ' the source' call wtinfo(chatter, 0, 1, subinfo) cont = .false. return ENDIF IF ((sumrcts / sumcts).LE.(0.1)) THEN subinfo = 'source is less than 10% of total counts !' call wtwarm(subname, version, chatter, 5, subinfo) ENDIF return end ! --------------------------------------------------------------------- ! END OF SUBROUTINE RPSF_CONV ! --------------------------------------------------------------------- !+RPSF_WT ! ---------------------------------------------------------------------+ subroutine rpsf_wt(outfil, infil, stver, rad_lo, & rad_hi, cts_arc2, err_arc2, area_wgt, pixsize, & npts, maxrad, maxtheta, telstw, inststw, & bkgd, chanmin, chanmax, & sumrcts, sumtcts, cont, chatter, killit) ! ----------------- DESCRIPTION ------------------------------------- ! ! This routine writes to FITS output file in CALRPSF format. ! It opens the FITS file and passes the data to WTRPF1 ! which writes the data in the form of a BINARY TABLE extension. ! Finally it closes the FITS file. ! ! ---------------------- VARIABLES -------------------------------- IMPLICIT NONE character*(*) infil, outfil character*(*) telstw, inststw character*(*) stver integer maxrad, npts, chatter, maxtheta, chanmin, chanmax real cts_arc2(maxrad, maxtheta, *), err_arc2(maxrad, maxtheta, *) real bkgd, sumrcts, pixsize, sumtcts real rad_lo(*), rad_hi(*), area_wgt(maxrad, maxtheta, *) logical cont, killit ! -------------------------- INTERNAL VARIABLES ---------------------- character(80) subinfo integer maxhist, maxcomm parameter (maxhist = 10) parameter (maxcomm = 10) character(80) hist(maxhist), comms(maxcomm), extinfo integer ounit, nk_hist, status, nk_comm, ierr real theta_lo(1), theta_hi(1), energ_lo(1), energ_hi(1) integer ntheta, nenerg character(8) extname character(16) hduclas3, taskname character(16) radunit, thetaunit, energunit, rpsfunit logical qarea, qerror ! ------------------------- VARIABLE DIRECTORY ------------------------ ! ! Arguments ... ! ! outfil char : Output filename ! infil char : Input filename ! stver char : Main program version ! maxrad int : Array dimensions for radial profile data ! cts_arc2 real : Array of radial profile in counts per arcmin^2 ! err_arc2 real : Array of errors on cts_arc2 (in arcmin^2) ! rad_lo real : Array of lower edge of radial bins (in arcmin) ! rad_hi real : Array of upper edge of radial bins (in arcmin) ! n_pix int : Array of No of pixels in bin ! Pixsize real : Pixsize in degrees ! npts int : Counter for No. radial profile values ! telstw char : Mission/Telescope name ! inststw char : Instrument/Detector name ! chatter int : Chattines flag (>20 noisy) ! ! Internals ... ! ! ounit int : Fortran file unit number ! maxhist int : Dimension for history comment array ! maxcomm int : Dimension for comments array ! nk_comm int : No. of lines of comments ! nk_hist int : No of lines of history comments ! status int : Error flag ! hist char : History comment array ! comms char : Comment array ! version char : Subroutine version ! errstr char : Error string for this routine ! errmess char : Error message obtained from FITSIO call ! ! --------------------- CALLED ROUTINES --------------------------------- ! ! subroutine OP_NA : (FITSIO) Opens FITS file,writes header and null ! primary array ! subroutine FTGERR : (FITSIO) Obtains appropriate error text ! subroutine WTRPF1 : Writes data extension ! subroutine FTCLOS : (FITSIO) Closes FITS file ! subroutine WT_FERRMSG : Writes Error text if required ! ! ---------------- COMPILATION AND LINKING ---------------- ! ! Link with FITSIO and FTOOLS ! ! ---------------- AUTHORS/MODIFICATION HISTORY ---------------- ! ! Rehana Yusaf (Feb 1993) ! Rehana Yusaf (2 Sept 1993) 1.0.1; Change call to WT_RPSF1993a - new ! more general routine ! Rehana Yusaf (14 Jan 1994) 1.0.2; Renamed WT_RPSF1993a to WTRPF1, for ! CALLIB naming convention ! Also additional keywords are written ! BACKGRND and if ROSAT PSPC then ! CHANMIN and CHANMAX ! Rehana Yusaf (Feb 3 1994) 1.0.3; Write SUMRCTS ! Rehana Yusaf (Mar 7 1995) 1.0.4; Write SUMTCTS, theoretical sum of ! counts, should be same as SUMRCTS ! when no regions are excluded ! Rehana Yusaf (1995 April 25) 1.0.5; add killit (clobber) parameter ! ! Banashree Mitra Seifert (1995 Dec) 1.1.0: ! . Filename character *180 ! . Dimensions to be carried over by calling routine ! . Screen display is used ! wtinfo,wtferr,wtfwrn ! ----------------------------------------------------------------- character(5) version parameter (version = '1.1.0') character(8) subname !- ! -------------------------------------------------------------- subname = 'rpsf_wt' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) ! -------------------- OPENING FITS FILE ------------------------ call cgetlun(ounit) status = 0 call opnpa(outfil, chatter, ounit, killit, status) extinfo = ' setting up primary header' call wtfwrn(subname, version, chatter, 15, status, extinfo) IF (status.NE.0) THEN cont = .false. goto 200 ENDIF ! ----------------- SETTING UP HISTORY COMMENT FIELDS ---------------- nk_hist = 2 hist(1) = 'st2rpsf converts from STW -> RPSF format' hist(2) = 'STW FILE : ' // infil nk_comm = 0 ! ------- SETTING UP UNDEFINED ARRAYS ! THETA AND ENERG VALUES ------- theta_lo(1) = -99 theta_hi(1) = -99 energ_lo(1) = -99 energ_hi(1) = -99 ntheta = 1 nenerg = 1 ! ---------------------- WRITING DATA -------------------------------- radunit = 'arcmin' thetaunit = 'arcmin' energunit = 'KeV' rpsfunit = 'count/arcmin**2' qarea = .true. qerror = .true. ierr = 0 extname = 'OBS RPSF' IF (bkgd.EQ.0.0) THEN hduclas3 = 'NET' ELSE hduclas3 = 'TOTAL' ENDIF 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, telstw, inststw, maxrad, & maxtheta, ierr, chatter) IF (ierr.NE.0) THEN cont = .false. goto 200 ENDIF ! ----------------------- CREATOR KEYWORD ---------------------------- taskname = 'st2rpsf' // ' ' // stver status = 0 call ftpkys(ounit, 'CREATOR', taskname, & 's/w task which wrote this dataset', & status) extinfo = ' writing CREATOR keyword' call wtfwrn(subname, version, chatter, 15, status, extinfo) ! --------------- WRITE PIXELSIZE ---------------------------------- 200 status = 0 call ftpkye(ounit, 'PIXSIZE', pixsize, 8, 'In decimal degrees', & status) extinfo = ' writing PIXSIZE keyword' call wtferr(subname, version, status, extinfo) ! ------------------ WRITE BACKGROUND -------------------- status = 0 call ftpkye(ounit, 'BACKGRND', bkgd, 8, 'In counts per pixel', & status) extinfo = ' writing BACKGRND keyword' call wtferr(subname, version, status, extinfo) ! ------------ WRITE CHANMIN and CHANMAX for ROSAT PSPC -------------- status = 0 call ftpkyj(ounit, 'CHANMIN', chanmin, & 'Minimum PI channel for image', status) extinfo = ' writing CHANMIN keyword' call wtferr(subname, version, status, extinfo) status = 0 call ftpkyj(ounit, 'CHANMAX', chanmax, & 'Maximum PI channel for image', status) extinfo = ' writing CHANMAX keyword' call wtferr(subname, version, status, extinfo) ! --------------------- WRITE SUMRCTS ------------------------------- status = 0 call ftpkye(ounit, 'SUMRCTS', sumrcts, 8, & 'Sum of source counts under profile', status) extinfo = ' writing SUMRCTS' call wtferr(subname, version, status, extinfo) ! ------------------- WRITE SUMTCTS -----------------------------------+ status = 0 call ftpkye(ounit, 'SUMTCTS', sumtcts, 8, & 'Theoretical Sum source counts(corr for '& // 'excl regions)', status) extinfo = ' writing SUMTCTS' call wtferr(subname, version, status, extinfo) ! ----------------- CLOSE FITS FILE --------------------------------- status = 0 call ftclos(ounit, status) extinfo = 'closing output file' call wtferr(subname, version, status, extinfo) IF (status.NE.0) THEN cont = .false. ENDIF call cfrelun(ounit) return end ! ---------------------------------------------------------------------+ ! END OF SUBROUTINE RPSF_WT ! ---------------------------------------------------------------------+ !+RPSF_SRC ! ---------------------------------------------------------------------+ subroutine rpsf_src(sou_a, del_rad, cont, chatter) ! ----------------------- DESCRIPTION ---------------------------------- ! This subroutine is a parser for a particular character string. It ! extracts inner and outer radius values, as well as nsteps from the ! string. The units of inner and outer have to be arcsecs or arcmins ! otherwise this routine returns to main and stops. ! ! ------------------------- VARIABLES ---------------------------------- character*(*) sou_a character(100) source character(70) subinfo, extinfo character(12) charnum, cnum character(1) dpost, spost real inner, outer, nbins, del_rad, cnvfact integer beg, end, pos, skipstr, i integer curpos, chatter, eflg, ck_char logical dunit, cont ! ------------------------ VARIABLE DIRECTORY ------------------------ ! ! Passed Parameters ... ! ! source char : String containing source information ! del_rad real : Del_Rad = (outer rad - inner rad)/Nsteps (in arcmin) ! chatter int : Flag for writing to screen (>20 write) ! cont logical : False if program is to be terminated ! ! Internals ... ! ! charnum char : String used for extracting numbers from source ! dpost char : Used for conversion factor ' or " ! inner real : Inner radial value in arcmins or arcsecs ! outer real : Outer radial value in arcmins or arcsecs ! nbins real : No of divisions of outer - inner ! cnvfact real : Conversion facter, to convert to arcmins ! beg int : Used as position marker when finding beginning of ! number in the source string ! end int : Used as position marker when finding end of ! number in the source string ! pos int : Used as position marker ! skipstr int : Counter for skipped expressions in source ! curpos int : Another position marker, used with function index ! subinfo char : User information about subroutine ! extinfo char : Additional comments ! ! ! --------------------------- CALLED ROUTINES ----------------------- ! ! subroutine FCECHO : (FTOOLS) Library routine to write to screen ! subroutine CRMVBXK : (CALLIB) Compresses two or more blanks to one ! function INDEX : (FORTRAN) Library funnction which searches ! for a substring in specified string ! ! ------------------ COMPILATION AND LINKING -------------------------- ! ! Link with FTOOLS and XANLIB ! ! ---------------- AUTHORS/MODIFICATION HISTORY ------------------------ ! ! Rehana Yusaf ! Rehana Yusaf (18 august 1993) 1.0.1; add more error checking ! ! Banashree Mitra Seifert (1995 Dec) 1.1.0: ! . Filename character *180 ! . Dimensions to be carried over by calling routine ! . Screen display is used ! wtinfo,wterrm ! ----------------------------------------------------------------- character(5) version parameter (version = '1.1.0') character(9) subname !- ! ---------- PROVIDING USER INFORMATION IF DESIRED ------------------- subname = 'rpsf_src' subinfo = 'using ' // subname // ' ' // version call wtinfo(chatter, 10, 1, subinfo) call crmvlbk(sou_a) call rmvexsp(sou_a, source) extinfo = ' Source string : ' call wtinfo(chatter, 20, 2, extinfo) call wtinfo(chatter, 20, 2, source) ! ---------------- SKIPPING FIRST THREE STRINGS --------------------- i = 1 pos = 0 skipstr = 0 do WHILE (skipstr.LT.3) curpos = index(source(i:), ' ') skipstr = skipstr + 1 i = i + curpos pos = pos + curpos enddo beg = pos + 1 ! ---------------- READING INNER RADIAL VALUE --------------------- dunit = .false. dpost = '"' pos = beg curpos = index(source(pos:), ' ') end = beg + curpos - 2 IF (source(end:end).EQ.dpost) THEN dunit = .true. end = end - 1 ENDIF eflg = 0 charnum = source(beg:end) read(charnum, *, IOSTAT = eflg) inner IF (eflg.NE.0) THEN extinfo = 'rpsf_src : reading inner radius !' call wterrm(subname, version, extinfo) extinfo = ' extracting from :' // charnum call wterrm(subname, version, extinfo) cont = .false. return ENDIF extinfo = 'extracting inner radius from :' // charnum call wtinfo(chatter, 30, 3, extinfo) ! ---------- READING OUTER RADIAL VALUE ----------- IF (dunit) THEN pos = end + 3 ELSe pos = end + 2 ENDIF beg = pos curpos = index(source(pos:), ' ') end = beg + curpos - 3 charnum = source(beg:end) read(charnum, *, IOSTAT = eflg) outer IF (eflg.NE.0) THEN extinfo = 'reading outer radius' call wterrm(subname, version, extinfo) extinfo = ' extracting from :' // charnum call wtinfo(chatter, 1, 1, extinfo) cont = .false. return ENDIF ! ------------------ READING CNVFACT ------------------------------- i = end + 1 IF (source(i:i).EQ.dpost) THEN cnvfact = 1 / real(60) ELSE cnvfact = 1 ENDIF extinfo = ' extracting outer radius from :' // charnum call wtinfo(chatter, 30, 3, extinfo) ! ------- UNITS CHECKS (arcmins or arcsecs, otherwise invalid) ------- spost = '''' IF ((source(i:i).NE.dpost).and.(source(i:i).NE.spost)) THEN extinfo = 'invalid units for inner and outer radius !' call wterrm(subname, version, extinfo) extinfo = 'should be arcsecs or arcmins ' call wtinfo(chatter, 0, 1, extinfo) cont = .false. return ENDIF ! ----------------------- READING NBINS ----------------------------- pos = end + 5 beg = pos curpos = 0 curpos = index(source(pos:), '.') IF (curpos.EQ.0) THEN curpos = index(source(pos:), ' ') ENDIF end = beg + curpos charnum = source(beg:end) ck_char = ichar(source(beg:beg)) IF ((ck_char.LT.49).OR.(ck_char.GE.58)) THEN beg = beg + 1 ENDIF ck_char = ichar(source(end:end)) IF ((ck_char.LT.49).OR.(ck_char.GE.58)) THEN end = end + 1 ENDIF eflg = 0 read(charnum, *, IOSTAT = eflg) nbins IF (eflg.NE.0) THEN extinfo = 'reading nbins !' call wterrm(subname, version, extinfo) extinfo = ' extracting nbins from : ' // charnum call wtinfo(chatter, 0, 1, extinfo) cont = .false. return ENDIF extinfo = ' extracting nbins from :' // charnum call wtinfo(chatter, 30, 3, extinfo) ! -------------------- USER INFO --------------------------------- ! --------------- cnvfact ... IF (cnvfact.EQ.1) THEN extinfo = 'data in source already in arcmins,cnvfact=1' ELSE extinfo = 'data in source in arcsecs,cnvfact=1/60' ENDIF call wtinfo(chatter, 30, 3, extinfo) ! --------------- inner radius ... inner = inner * cnvfact eflg = 0 write(cnum, 100, IOSTAT = eflg) inner extinfo = 'inner radius value (in arcmins) : ' // cnum call wtinfo(chatter, 20, 2, extinfo) ! ---------------- outer radius ... outer = outer * cnvfact eflg = 0 write(cnum, 100, IOSTAT = eflg) outer extinfo = ' Outer radius value (in arcmins) : ' // cnum call wtinfo(chatter, 20, 2, extinfo) ! --------------- Nbins ... eflg = 0 write(cnum, 100, IOSTAT = eflg)nbins extinfo = ' N steps : ' // cnum call wtinfo(chatter, 20, 2, extinfo) ! ----------------- CALCULATE DEL_RAD --------------------------- del_rad = ((outer - inner) / nbins) 100 FORMAT(F12.8) return end ! ---------------------------------------------------------------------+ ! END OF SUBROUTINE RPSF_SRC ! ---------------------------------------------------------------------+