!+COL2IMG subroutine col2ig IMPLICIT NONE ! ! Description: ! Task to produce a 2-D image in the FITS Primary Array of an ! input Collimator response dataset. The input dataset must be in ! an OGIP-approved format, and can consist of one of the following: ! - an OBSFACT dataset (fn of X & Y) ! - a COLLRESP dataset (fn of E, X & Y) ! where X & Y represent the spatial coordinates, and E an energy coordinate ! In the case of energy-dependent datasets, the image is constructed ! by adding the collimator responses over a user-defined energy range ! In the case of datasets which are not stored on a uniform Cartesian ! spatial grid, user-defined parameters control whether the image ! is constructed using the stored grid (and hence will be "distorted"), ! or is remapped/rebinned onto a user-defined regular grid. ! ! Passed Parameters ! None ! ! User i/ps required (prompted for): ! None ... isolated GP_COL2IMG (below) ! ! COMMONS/INCLUDES etc ! common TASK : (FTOOLS) standard fatal error message thingy ! ! Called routines ! subroutine GP_COL2IMG : (below) Gets parameters from XPI par file ! subroutine DO_COL2IMG : (below) Performs the conversion ! ! Compilation: ! subroutines require XPI, CALLIB, FTOOLS, FITSIO ! ! Origin: ! Original ! ! Authors/Modification History: ! Ian M George (0.9.0:1994 Mar 10), quick & v.dirty ! Ian M George (1.0.0:1994 Aug 05), Test version for XTE ! Ian M George (2.0.0:1994 Oct 14), Added support for COLLRESP datasets ! Ian M George (2.1.0:1994 Oct 25), Bug fix to COLLRESP energy remapping ! Ian M George (2.1.1:1994 Dec 14), cleaned up for the public ! Rehana Yusaf (2.1.2:1994 Dec 22), fixed variables that are multiply ! defined. Also fixed warnings that ! appeared on legacy, due to ! message = errstr//'...etc ! message and errstr were both defined ! as character(80) ! Ian M George (2.2.0:1995 Jun 27), DMA bug-fix in do_col2img & esquish ! Ian M George (2.2.1;1995 Jun 28), minor stuff in CK_OBSFACT & CK_COLRESP ! Ian M George (2.2.2:1995 Jun 30), added fnd* calls to DO_COL2IMG ! Ian M George (2.2.3:1996 Feb 03), added wtinfo & friends ! Bryan K Irby (2.3.0:2018 Jan 11), replaced udmget/udmfre with allocate/deallocate ! MFC (2.3.1 2020 Apr 16) f90 version character(7) version parameter (version = '2.3.1') !- ! Internals character(40) taskname integer chatter, ierr character(10) emin, emax character(80) infil, outfil integer nxbins, nybins real xmin, xmax, ymin, ymax logical remap, qextrap, clobber ! Initialize ierr = 0 taskname = 'col2img' ! Get Parameters from the par file call gp_col2img(infil, outfil, remap, clobber, chatter, & xmin, xmax, nxbins, ymin, ymax, nybins, qextrap, & emin, emax, ierr) ! Start-up Main call wtbegm(taskname, version, chatter) ! Do the nasty deed call do_col2img(taskname, infil, outfil, remap, clobber, chatter, & xmin, xmax, nxbins, ymin, ymax, nybins, qextrap, & emin, emax, version, ierr) if(ierr.NE.0) goto 148 ! Finish-Off 148 continue call wtendm(taskname, version, ierr, chatter) return end ! ------------------------------------------------------------------------- !+GP_COL2IMG subroutine gp_col2img(infil, outfil, remap, clobber, chatter, & xmin, xmax, nxbins, ymin, ymax, nybins, qextrap, & emin, emax, ierr) IMPLICIT NONE integer chatter, ierr integer nxbins, nybins real xmin, xmax, ymin, ymax character*(*) infil, outfil, emin, emax logical remap, qextrap, clobber ! ! Description: ! Gets the parameters required by COL2IMG from the parameter file ! NOTE - The par file is assumed to have been opened. ! ! User i/ps required (prompted for): ! INFIL - Input File name/expression ! OUTFIL - name of o/p Image file required ! CLOBBER - Clobber flag (required by FTOOLS team) ! CHATTER - chattiness flag for o/p (5 quite,10 normal,15 high,>20 silly) ! IERR - Error Flag (zero = OK) ! ! Origin: ! Original ! ! Called Routines ! subroutine WT* : (CALLIB) CALDB-standard writers ! subroutine UCLG* : (XPI) returns parameters ! ! Compilation: ! requires XPI/Host interface etc and CALLIB ! ! Authors/Modification History: ! Ian M George (0.9.0:1994 Mar 10), Limited capability version ! Ian M George (1.0.0:1994 Sep 22), Clobber read from Par file ! Ian M George (1.1.0:1994 Oct 14), Emin, Emax read from Par file ! Ian M George (1.1.1:1996 Feb 04), added wtinfo & friends ! MFC (1.1.2 2020 Apr 16) f90 version character(7) version parameter (version = '1.1.1') !- ! Internals character(10) subname parameter (subname = 'gp_col2img') character(50) contxt ! Initialize ierr = 0 ! Get the name of the input file call uclgst('infil', infil, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting INFIL parameter') goto 999 endif ! Get the name of the output (Image) file call uclgst('outfil', outfil, ierr) if(ierr.ne.0) then call wterrm(subname, version, & 'Problem getting OUTFIL parameter') goto 999 endif ! Get the min & max energies call uclgst('emin', emin, ierr) if(ierr.ne.0) then call wtwarm(subname, version, 1, 1, & 'Problem getting EMIN parameter') ierr = 0 call wtinfo(1, 5, 1, 'will use lower boundary stored') emin = '%' elseif(emin.EQ.' ') then emin = '%' endif call uclgst('emax', emax, ierr) if(ierr.ne.0) then call wtwarm(subname, version, 1, 1, & 'Problem getting EMAX parameter') ierr = 0 call wtinfo(1, 5, 1, 'will use upper boundary stored') emax = '%' elseif(emax.EQ.' ') then emax = '%' endif ! Read in clobber call uclgsb('clobber', clobber, ierr) if (ierr.NE.0) then clobber = .false. call wtwarm(subname, version, 1, 1, & 'Problem getting CLOBBER parameter') ierr = 0 call wtinfo(1, 5, 1, 'setting CLOBBER = False') endif ! Get the remapping Flag call uclgsb('remap', remap, ierr) if(ierr.NE.0) then call wtwarm(subname, version, 1, 1, & 'Problem getting REMAP parameter') ierr = 0 call wtinfo(1, 5, 1, 'setting REMAP = False') remap = .false. endif ! Get the chattiness flag call uclgsi('chatter', chatter, ierr) if(ierr.NE.0) then call wtwarm(subname, version, 1, 1, & 'Problem getting CHATTER parameter') ierr = 0 call wtinfo(1, 1, 1, 'setting CHATTER = 10') chatter = 10 endif ! Give user info if requested contxt = ' using ' // subname // ' ' // version call wtinfo(chatter, 20, 1, contxt) ! Get the extra stuff required when remap reuested if(remap) then call uclgsb('extrap', qextrap, ierr) if(ierr.NE.0) then call wtwarm(subname, version, 1, 1, & 'Problem getting EXTRAP parameter') ierr = 0 call wtinfo(1, 5, 1, 'setting EXTRAP = False') qextrap = .false. endif call uclgsr('xmin', xmin, ierr) if(ierr.NE.0) then call wterrm(subname, version, & 'Problem getting XMIN parameter') goto 999 endif call uclgsr('xmax', xmax, ierr) if(ierr.NE.0) then call wterrm(subname, version, & 'Problem getting XMAX parameter') goto 999 endif call uclgsi('nxbins', nxbins, ierr) if(ierr.NE.0) then call wterrm(subname, version, & 'Problem getting NXBINS parameter') goto 999 endif call uclgsr('ymin', ymin, ierr) if(ierr.NE.0) then call wterrm(subname, version, & 'Problem getting YMIN parameter') goto 999 endif call uclgsr('ymax', ymax, ierr) if(ierr.NE.0) then call wterrm(subname, version, & 'Problem getting YMAX parameter') goto 999 endif call uclgsi('nybins', nybins, ierr) if(ierr.NE.0) then call wterrm(subname, version, & 'Problem getting NYBINS parameter') goto 999 endif endif 999 if(ierr.ne.0) then call wterrm(subname, version, ' unable to continue') endif return end ! ------------------------------------------------------------------------- !+DO_COL2IMG subroutine do_col2img(taskname, inexp, outfil, remap, clobber, & chatter, xmin, xmax, nxbins, ymin, ymax, nybins, qextrap, & eminstr, emaxstr, mnver, ierr) IMPLICIT NONE integer chatter, ierr integer nxbins, nybins real xmin, xmax, ymin, ymax character*(*) inexp, outfil, eminstr, emaxstr, mnver, taskname logical remap, qextrap, clobber ! ! Description: ! Does the business for COL2IMG ! ! User i/ps required (prompted for): ! INEXP - Input File name/expression ! OUTFIL - name of o/p Image file required ! REMAP - Whether spatial remapping is requested ? ! CLOBBER - Clobber flag ! CHATTER - chattiness flag for o/p (5 quite,10 normal,15 high,>20 silly) ! IERR - Error Flag (zero = OK) ! ! Origin: ! Original ! ! Called Routines ! subroutine FCECHO : (FTOOLS) writes to standard o/p device ! subroutine WT_FERRMSG : (CALLIB) Writes standard FITSIO message etc ! ! Compilation: ! requires XPI/Host interface etc and CALLIB ! ! Authors/Modification History: ! Ian M George (0.9.0:1994 Mar 10), Limited capability version ! Ian M George (1.0.0:1994 Sep 22), Added clobber flag ! Ian M George (2.0.0:1994 Oct 14), Added support for COLLRESP datasets ! Rehana Yusaf (2.0.1:1994 Dec 22), mnver passed to this subroutine, ! so that mnver is used for COL2IMG ! version in HISTORY comment ! Ian M George (2.0.2:1995 Jun 27), fixes with DMA pointers ! Ian M George (2.1.0:1995 Jun 30), inserted fnd* calls ! Ian M George (2.1.1:1996 Feb 04), added wtinfo & friends ! Bryan K Irby (2.2.0:2018 Jan 11), replaced udmget/udmfre with allocate/deallocate character(7) version parameter (version = '2.2.0') !- ! ************************************ ! Max Array sizes integer maxhist, maxcomm integer maxextn, maxen parameter (maxhist = 30, maxcomm = 100) parameter (maxextn = 10, maxen = 100) ! Internals character(10) subname character(160) message parameter (subname = 'do_col2img') integer decimals, clenact integer inextn, status, itemp integer iunit, ounit, block integer i, htype integer nk_hist, nk_comm integer bitpix, naxis, gcount, pcount integer naxes(2), icoord1, icoord2, ienerg integer maxenerg, maxcoord1, maxcoord2 logical simple, extend real emin, emax character(5) obfversn, crspvers character(20) telescop, instrume, detnam, filter character(70) hist(maxhist), comment(maxcomm) character(80) dummystr(1) character(80) infil, string character(20) cunit1, cunit2, eunits character(20) cname1, cname2 integer ninstr, nsearch, nfound integer next(maxextn) character(20) instr(9), ftype character(20) outhdu(9, maxextn), outver(9, maxextn) character(20) hduclas(9), hduvers(9) character(20) extnam(maxextn) logical qokfil, qxlin, qylin real delx, dely, ratio ! ... pointers to "arrays" to be dynamically allocated real, allocatable :: p_x(:), p_y(:) real, allocatable :: p_x_out(:), p_y_out(:) real, allocatable :: p_energ_lo(:), p_energ_hi(:) real, allocatable :: p_ework(:), p_collresp(:) real, allocatable :: p_image(:), p_img_out(:) ! ... "arrays" to be dynamically allocated ! real x(maxpos), real y(maxpos) ! real energ_lo(maxen) real energ_hi(maxen) ! real x_out(maxpos), real y_out(maxpos) ! real image(maxpos,maxpos), real img_out(maxpos,maxpos) ! real collresp(maxen,maxpos,maxpos) ! Initialize ierr = 0 status = 0 decimals = 5 ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 20, 1, message) ! Check that the o/p file doesn't already exist or is illegal call ck_file(outfil, dummystr, 1, qokfil, clobber, chatter) if(.NOT.qokfil) then message = 'Offending file is OUTFIL: ' // outfil call wtinfo(chatter, 1, 2, message) ierr = -1 goto 482 endif ! Fill in a few history & comment records hist(1) = ' IMAGE constructed by COL2IMG ' // mnver nk_hist = 1 comment(1) = ' COL2IMG ' // mnver // 'Summary:' comment(2) = ' I/p file: ' // inexp nk_comm = 2 if(remap) then comment(nk_comm + 1) = ' Spatial remapping requested by user' nk_comm = nk_comm + 1 if(qextrap) then comment(nk_comm + 1) = ' (extrapolation beyond spatial ' // & 'grid allowed by user' nk_comm = nk_comm + 1 endif else comment(nk_comm + 1) = & ' No spatial remapping requested by user' nk_comm = nk_comm + 1 endif ! Translate the name of the input file (taking off any extension) call fcpars(inexp, infil, inextn, status) if(status.NE.0) then call wtwarm(subname, version, chatter, 1, & ' Problem passing infil expression') call wtinfo(chatter, 1, 1, 'will search all extensions') inextn = -99 endif ! Open i/p file status = 0 call cgetlun(iunit) call ftopen(iunit, infil, 0, block, status) IF (status.NE.0) THEN message = ' opening file: ' // infil call wtferr(subname, version, status, message) ierr = 1 goto 482 ENDIF ! GO FIND THE EXTENSION ! - Extension number IS given as part of infil if(inextn.GE.0) then call ftmrhd(iunit, inextn, htype, status) IF (status.NE.0) THEN call wtferr(subname, version, status, & ' Problem moving to specified xtens') ierr = 1 goto 482 ENDIF ! ... grab the HDUCLAS values for reference ninstr = 1 instr(1) = '*' nsearch = 1 call fndhdu(chatter, iunit, ninstr, instr, & nsearch, nfound, next, outhdu, outver, extnam, ierr) nfound = 1 next(1) = 0 ftype = outhdu(2, 1) goto 200 else ! - Extension SEARCHING required ! .... Search for OBFFACTOR nsearch = 99 call fndobf(chatter, iunit, nsearch, & obfversn, hduclas, hduvers, & ierr) if(ierr.EQ.0) then ftype = hduclas(2) goto 200 elseif(ierr.GT.0) then goto 482 else ierr = 0 endif ! .... Search for COLLRESP call fndcol(chatter, iunit, nsearch, & crspvers, hduclas, hduvers, & ierr) if(ierr.EQ.0) then ftype = hduclas(2) else call wterrm(subname, version, & ' Unable to locate acceptable extn') goto 482 endif endif 200 continue ! Check out the input file, getting array sizes if(ftype.EQ.'OBSFACTOR') then call ck_obsfact(iunit, chatter, & obfversn, & icoord1, icoord2, ierr) elseif(ftype.EQ.'COLLRESP') then call ck_colresp(iunit, chatter, & crspvers, & ienerg, icoord1, icoord2, ierr) else call wterrm(subname, version, & ' Unable to locate acceptable extn') ierr = 5 goto 482 endif ! Allocate dynamic memory if(remap) then maxcoord1 = MAX(icoord1, nxbins) maxcoord2 = MAX(icoord2, nybins) maxcoord1 = MAX(maxcoord1, 50) maxcoord2 = MAX(maxcoord2, 50) allocate(p_x_out(maxcoord1), stat = status) if(status.NE.0) goto 765 allocate(p_y_out(maxcoord2), stat = status) if(status.NE.0) goto 765 allocate(p_img_out(maxcoord1 * maxcoord2), stat = status) if(status.NE.0) goto 765 else maxcoord1 = icoord1 maxcoord2 = icoord2 endif allocate(p_x(maxcoord1), stat = status) if(status.NE.0) goto 765 allocate(p_y(maxcoord2), stat = status) if(status.NE.0) goto 765 allocate(p_image(maxcoord1 * maxcoord2), stat = status) if(status.NE.0) goto 765 if(ftype.EQ.'COLLRESP') then maxenerg = MAX(ienerg, 50) allocate(p_energ_lo(maxenerg), stat = status) if(status.NE.0) goto 765 allocate(p_energ_hi(maxenerg), stat = status) if(status.NE.0) goto 765 allocate(p_ework(maxenerg), stat = status) if(status.NE.0) goto 765 allocate(p_collresp(maxenerg * maxcoord1 * maxcoord2), stat = status) if(status.NE.0) goto 765 endif 765 if(status.NE.0) then call wterrm(subname, version, & ' Failed to allocate Dynamic Memory') ierr = -1 goto 482 elseif(chatter.GE.20) then call wtinfo(chatter, 10, 1, & 'Dynamic Memory successfully allocated') endif ! ***** ! --------------------- Start of DATA reading ------------------- ! OBSFACTOR if(ftype.EQ.'OBSFACTOR') then call rdobf1(iunit, chatter, & obfversn, & telescop, instrume, detnam, filter, & cname1, cname2, & maxcoord1, maxcoord2, & icoord1, icoord2, p_X, p_Y, & cunit1, cunit2, p_image, & ierr) if(ierr.NE.0) then call wterrm(subname, version, & ' Reading OBSFACTOR dataset') ierr = 2 goto 482 endif call post_obf(chatter, & nk_comm, comment, & cname1, cname2, & cunit1, cunit2, & icoord1, icoord2, p_x, p_y, & ierr) if(ierr.NE.0) then call wterrm(subname, version, & ' Manipulating Dataset') ierr = 2 goto 482 endif ! COLLRESP elseif(ftype.EQ.'COLLRESP') then call rdcol1(iunit, chatter, & crspvers, & telescop, instrume, detnam, filter, & cname1, cname2, & maxenerg, maxcoord1, maxcoord2, & ienerg, p_energ_lo, p_energ_hi, eunits, & icoord1, icoord2, p_x, p_y, & cunit1, cunit2, p_collresp, ierr) if(ierr.NE.0) then call wterrm(subname, version, & ' Reading COLLRESP dataset') ierr = 2 goto 482 endif call post_col(chatter, & nk_comm, comment, & ienerg, p_energ_lo, p_energ_hi, eunits, & cname1, cname2, & cunit1, cunit2, & icoord1, icoord2, p_x, p_y, & ierr) if(ierr.NE.0) then call wterrm(subname, version, & ' Manipulating Dataset') ierr = 2 goto 482 endif else call wterrm(subname, version, & ' Unrecognized/Unsupported HDUCLAS2 value') ierr = 10 goto 482 endif ! --------------------- End of DATA reading ------------------- ! Close the FITS file call ftclos(iunit, status) if(status.ne.0) then call wtferr(subname, version, status, & 'Problem closing i/p data file') ierr = 6 goto 482 endif ! Spatial Grid Check out ratio = 1.E-03 ! ... Loop around to see whether the spatial grid is regular call ck_sgrid(chatter, & cname1, cname2, & icoord1, icoord2, p_x, p_y, & ratio, qxlin, qylin, & delx, dely, & ierr) if(ierr.NE.0) then goto 482 endif if(qxlin.AND.qylin) then message = ' (linear binning along both axes)' comment(nk_comm + 1) = message nk_comm = nk_comm + 1 call wtinfo(chatter, 9, 2, message) elseif(qxlin) then message = ' (linear binning along axis-1)' comment(nk_comm + 1) = message nk_comm = nk_comm + 1 call wtinfo(chatter, 9, 2, message) message = & ' (non-linear binning along axis-2)' comment(nk_comm + 1) = message nk_comm = nk_comm + 1 call wtinfo(chatter, 9, 2, message) elseif(qylin) then message = & ' (non-linear binning along axis-1)' comment(nk_comm + 1) = message nk_comm = nk_comm + 1 call wtinfo(chatter, 9, 2, message) message = ' (linear binning along axis-2)' comment(nk_comm + 1) = message nk_comm = nk_comm + 1 call wtinfo(chatter, 9, 2, message) else message = & ' (non-linear binning along both axes)' comment(nk_comm + 1) = message nk_comm = nk_comm + 1 call wtinfo(chatter, 9, 2, message) endif ! Energy Grid check out and compression ! Trap out the OBSFACTOR case & warn users if necessary if(ftype.EQ.'OBSFACTOR') then if((eminstr.NE.'%') .AND. (emaxstr.NE.'%')) then call wtwarm(subname, version, chatter, 1, & 'Dataset is an OBSFACTOR') call wtinfo(chatter, 1, 2, '(No energy dependence)') call wtinfo(chatter, 1, 2, & 'User-defining energy range ignored') nk_comm = nk_comm + 1 comment(nk_comm) = message eminstr = '%' emaxstr = '%' endif else call ck_egrid(chatter, & nk_comm, comment, & ienerg, p_energ_lo, p_energ_hi, eunits, & eminstr, emaxstr, emin, emax, & ierr) if(ierr.NE.0) then goto 482 endif call esquish(chatter, & maxenerg, ienerg, p_energ_lo, p_energ_hi, & p_ework, & maxcoord1, maxcoord2, icoord1, icoord2, & p_collresp, & emin, emax, p_image, & ierr) if(ierr.NE.0) then goto 482 endif ! ... Add a few comments nk_comm = nk_comm + 1 comment(nk_comm) = ' I/p dataset was collapsed in ' // & 'energy-space to construct this image' nk_comm = nk_comm + 1 write(comment(nk_comm), '(a,f10.5,a,a)')& ' lower energy boundary: ', emin, & ' ', eunits nk_comm = nk_comm + 1 write(comment(nk_comm), '(a,f10.5,a,a)')& ' upper energy boundary: ', emax, & ' ', eunits ! ... Dump stuff to screen if requested write(message, '(a,f10.5,a,a)')& ' Lower energy boundary in o/p image: ', emin, & ' ', eunits call wtinfo(chatter, 15, 1, message) write(message, '(a,f10.5,a,a)')& ' Upper energy boundary in o/p image: ', emax, & ' ', eunits call wtinfo(chatter, 15, 1, message) endif 655 continue ! -------------------- Start of Remapping -------------------------------- if(remap) then call wtinfo(chatter, 9, 1, & 'Spatial Remapping requested') ! ... set up the desired grids call pre_remap(chatter, & xmin, xmax, nxbins, & ymin, ymax, nybins, & p_x_out, p_y_out, & ierr) if(ierr.NE.0) then goto 482 endif ! ... perform the remapping call sremap(chatter, & maxcoord1, maxcoord2, & icoord1, icoord2, & p_x, p_y, p_image, & maxcoord1, maxcoord2, & nxbins, nybins, & p_x_out, p_y_out, & p_img_out, & qextrap, ierr) ! ... error out if necessary if(ierr.NE.0) then goto 482 endif ! ... Fix up all the values for output call post_remap(chatter, & maxcoord1, maxcoord2, & p_x_out, p_y_out, & p_img_out, & maxcoord1, maxcoord2, & nxbins, nybins, & p_x, p_y, p_image, & ierr) if(ierr.NE.0) then call wterrm(subname, version, & 'Remapping Dataset') ierr = 2 goto 482 endif icoord1 = nxbins icoord2 = nybins qxlin = .true. qylin = .true. delx = (xmax - xmin) / nxbins dely = (ymax - ymin) / nybins ! ... Add a few comments concerning the remapping comment(nk_comm + 1) = & ' Input dataset was spatially remapped ' // & 'to construct this image' nk_comm = nk_comm + 1 write(message, '(a,g15.5)')& ' axis-1 min: ', xmin comment(nk_comm + 1) = message write(message, '(a,g15.5)')& ' axis-1 max: ', xmax comment(nk_comm + 2) = message write(message, '(a,i4)')& ' # (new) axis-1 bins: ', nxbins comment(nk_comm + 3) = message nk_comm = nk_comm + 3 write(message, '(a,g15.5)')& ' axis-2 min: ', ymin comment(nk_comm + 1) = message write(message, '(a,g15.5)')& ' axis-2 max: ', ymax comment(nk_comm + 2) = message write(message, '(a,i4)')& ' # (new) axis-2 bins: ', nybins comment(nk_comm + 3) = message nk_comm = nk_comm + 3 if(qextrap) then comment(nk_comm + 1) = & ' NOTE: extrapolation was enabled' else comment(nk_comm + 1) = & ' NOTE: extrapolation was disabled' endif nk_comm = nk_comm + 1 ! Write user info message = ' output image: ' // & cname1(:MIN(10, clenact(cname1))) // & ' vs ' // & cname2(:MIN(10, clenact(cname2))) call wtinfo(chatter, 9, 1, message) write(string, '(g15.5,a,g15.5)')& Xmin, ':', xmax call crmvblk(string) write(message, '(a,a, i4,a,a,a,a)')& cname1(:MIN(10, clenact(cname1))), ' has ', & nxbins, ' elements (', string(:clenact(string)), & ') ', cunit1 call wtinfo(chatter, 9, 2, message) write(string, '(g15.5,a,g15.5)')& ymin, ':', Ymax call crmvblk(string) write(message, '(a,a, i4,a,a,a,a)')& cname2(:MIN(10, clenact(cname2))), ' has ', & nybins, ' elements (', string(:clenact(string)), & ') ', cunit2 call wtinfo(chatter, 9, 2, message) call wtinfo(chatter, 9, 2, & '(linear binning along both axes)') endif ! -------------------- End of Remapping -------------------------------- ! --------------------- Start OUTPUT File ------------------------------ ! OPEN THE IMAGE FILE block = 2880 call cgetlun(ounit) call opfits(ounit, outfil, clobber, chatter, status) if(status.NE.0) then call wtferr(subname, version, status, & ' Opening outfile') ierr = 9 goto 482 endif call wtinfo(chatter, 20, 1, 'Opened the Output File') simple = .true. bitpix = -32 naxis = 2 naxes(1) = icoord1 naxes(2) = icoord2 pcount = 0 gcount = 1 extend = .true. call ftphpr(ounit, simple, bitpix, naxis, naxes, pcount, & gcount, extend, status) if(status.NE.0) then call wtferr(subname, version, status, & ' With FITSIO/FTPHPR call') ierr = 9 goto 482 endif call wtinfo(chatter, 20, 1, & ' Defined the Primary Header Data structure') ! Write a bunch of keywords status = 0 call ftpkys(ounit, 'CREATOR', taskname, & 's/w task which wrote this dataset', & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing CREATOR keyword') status = 0 endif ! ... The Coordinate Keywords call wt_ckys(chatter, ounit, & cname1, cunit1, qxlin, delx, p_x, & cname2, cunit2, qylin, dely, p_y, & ierr) ! ... The Mission/Instrument Keywords status = 0 call ftpkys(ounit, 'TELESCOP', telescop, & 'Mission/Satellite Name', & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing TELESCOP keyword') status = 0 endif call ftpkys(ounit, 'INSTRUME', instrume, & 'Instrument Name', & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing INSTRUME keyword') status = 0 endif if(detnam.NE.' '.OR.detnam.NE.'NONE') then call ftpkys(ounit, 'DETNAM', detnam, & 'Sub-Instrument Name', & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing DETNAM keyword') status = 0 endif endif ! ... The HDUCLAS/VERS Keywords status = 0 call ftpkys(ounit, 'HDUCLASS', 'OGIP', & 'Format conforms to OGIP Standards', & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing HDUCLASS keyword') status = 0 endif status = 0 call ftpkys(ounit, 'HDUCLAS1', 'IMAGE', & 'Dataset is an Image', & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing HDUCLAS1 keyword') status = 0 endif status = 0 call ftpkys(ounit, 'HDUVERS1', '1.0.0', & 'Version of family of formats', & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing HDUVERS1 keyword') status = 0 endif ! Add the energy bounds if appropriate if(ftype.EQ.'COLLRESP') then status = 0 string = 'Lower energy bound of image (' // & eunits(:clenact(eunits)) // ')' call ftpkye(ounit, 'ENERG_LO', emin, decimals, & string, & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing ENERG_LO keyword') status = 0 endif string = 'Upper energy bound of image (' // & eunits(:clenact(eunits)) // ')' call ftpkye(ounit, 'ENERG_HI', emax, decimals, & string, & status) if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem writing ENERG_HI keyword') status = 0 endif endif ! Add the (passed) history cards itemp = 0 do i = 1, nk_hist call FTPHIS(ounit, hist(i), status) if(status.NE.0) then itemp = status status = 0 call FTPHIS(ounit, & ' - (missing record) fitsio illegal character ?', & status) endif enddo if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem putting at least one HISTORY record') status = 0 endif call wtinfo(chatter, 20, 1, 'written the history keywords') status = 0 ! Add the (passed) comment cards itemp = 0 do i = 1, nk_comm call FTPCOM(ounit, comment(i), status) if(status.NE.0) then itemp = status status = 0 call FTPCOM(ounit, & ' - (missing record) fitsio illegal character ?', & status) endif enddo if(status.ne.0) then call wtwarm(subname, version, chatter, 1, & 'problem putting at least one COMMENT record') status = 0 endif call wtinfo(chatter, 20, 1, 'written the comment keywords') status = 0 ! Define the data structure status = 0 call ftpdef(ounit, bitpix, naxis, naxes, pcount, gcount, status) if(status.NE.0) then call wtferr(subname, version, status, & ' With FITSIO/FTPDEF call') ierr = 9 goto 482 endif call wtinfo(chatter, 20, 1, & ' Data structure has been defined') ! Write the Data status = 0 call ftp2de(ounit, 0, icoord1, naxes(1), naxes(2), & p_image, status) if(status.NE.0) then call wtferr(subname, version, status, & ' Writing the data') ierr = 9 goto 482 endif call wtinfo(chatter, 20, 1, & ' data written successfully') ! Close the FITS file status = 0 call ftclos(ounit, status) if(status.ne.0) then call wtferr(subname, version, status, & 'Problem closing o/p data file') ierr = 6 goto 482 endif ! --------------------- End of OUTPUT File ------------------------------ ! Check for errors 482 if(ierr.ne.0) then call wterrm(subname, version, 'Fatal') endif ! DeAllocate the dynamic memory if(remap) then deallocate(p_x_out, stat = status) if(status.NE.0) goto 766 deallocate(p_y_out, stat = status) if(status.NE.0) goto 766 deallocate(p_img_out, stat = status) if(status.NE.0) goto 766 endif deallocate(p_x, stat = status) if(status.NE.0) goto 766 deallocate(p_y, stat = status) if(status.NE.0) goto 766 deallocate(p_image, stat = status) if(status.NE.0) goto 766 if(ftype.EQ.'COLLRESP') then deallocate(p_energ_lo, stat = status) if(status.NE.0) goto 766 deallocate(p_energ_hi, stat = status) if(status.NE.0) goto 766 deallocate(p_ework, stat = status) if(status.NE.0) goto 766 deallocate(p_collresp, stat = status) if(status.NE.0) goto 766 endif 766 if(status.NE.0) then call wterrm(subname, version, & ' Failed to deallocate Dynamic Memory') ierr = 99 else call wtinfo(chatter, 20, 1, & ' ... Dynamic Memory successfully deallocated') endif ! ***** return end ! ------------------------------------------------------------------------- !+CK_OBSFACT subroutine ck_obsfact(iunit, chatter, & obfversn, & nix, niy, ierr) IMPLICIT NONE integer chatter, ierr integer iunit integer nix, niy character*(*) obfversn ! ! Description: ! Checks out an OBSFACT extension in one of the formats ! conforming to the HDUVERS2='1.*.*' family (for DMA purposes) ! !!! Note !!!! File is left open at the end ! ! Passed parameters ! iunit i : FORTRAN unit number of open RMF file ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! OBFVERSN i : String denoting OGIP HDUVERS2 family ! NIX i : Number of values for coordinate-1 ! NIY i : Number of values for coordinate-2 ! IERR o : Error flag (0 = OK) ! ! Called Routines: ! subroutine FCECHO (FTOOLS) Writes to standard o/p ! subroutine FTGBCL (FITSIO) ! subroutine FTGCNO (FITSIO) ! subroutine FTGHBN (FITSIO) ! subroutine FTGKYS (FITSIO) ! subroutine GT_CSYSNMS (CALLIB) Gets coordinate sytem column names ! subroutine WT_FERRMSG (CALLIB) Standard error string writer ! ! Compilation & Linking ! link with FITSIO & CALLIB & FTOOLS ! ! Origin: ! Original ! ! Authors/Modification History: ! Ian M George (1.0.0;94 Oct 17), original ! Ian M George (1.0.1;95 Jun 28), replaced GT_CSYSNMS w/ GCSYNM ! Ian M George (1.0.2:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '1.0.2') !- ! Internals character(10) subname parameter (subname = 'ck_obsfact') character(20) dummy20 double precision doubdum integer status, decimals, nfields integer irows, ncols, ivar integer colnum, enull, clenact parameter (nfields = 3, decimals = 6) character(8) dummy8 character(16) ttype(nfields), tform(nfields), tunits(nfields) character(20) csystem, col1, col2 character(40) comm character(80) message ! Initialize ierr = 0 status = 0 ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 20, 1, message) ! --- READING KEYWORD VALUES --- status = 0 call ftghbn(iunit, nfields, irows, ncols, ttype, tform, tunits, & comm, ivar, status) IF (status.NE.0) THEN call wtferr(subname, version, status, & ' reading binary header info') ierr = 2 goto 998 ENDIF ! OBFVERSN ... obfversn = ' ' status = 0 call ftgkys(iunit, 'HDUVERS2', obfversn, comm, status) IF(status.ne.0) then call wtferr(subname, version, status, & ' reading HDUVERS2 keyword') status = 0 ENDIF IF (obfversn.EQ.' ') THEN call ftgkys(iunit, 'OBFVERSN', obfversn, comm, status) if(status.ne.0) then call wtferr(subname, version, status, & ' reading OBFVERSN keyword') status = 0 ENDIF ENDIF ! Coordinate system in use status = 0 call ftgkys(iunit, 'CSYSNAME', csystem, comm, status) call wtferr(subname, version, status, & ' reading CSYSNAME keyword') IF (status.EQ.202) THEN csystem = 'XMA_POL' call wtinfo(chatter, 20, 2, 'assuming CSYSNAME = XMA_POL') elseif(status.ne.0) then csystem = 'XMA_POL' call wtwarm(subname, version, chatter, 1, & 'Setting CSYSNAME = XMA_POL') ENDIF ! Sort out the coordinate system & column names call gcsynm(chatter, 0, csystem, & col1, col2, dummy8, dummy8, ierr) IF (ierr.NE.0) goto 998 ! --- READ IN THE DATA ! ... The spatial grid ! ...... Coord-1 status = 0 call ftgcno(iunit, .false., col1, colnum, status) If (status.NE.0) THEN message = col1(:clenact(col1)) // & ' column not present' call wtferr(subname, version, status, message) ierr = 4 goto 998 ENDIF call ftgbcl(iunit, colnum, dummy20, dummy20, dummy20, & nix, doubdum, doubdum, enull, dummy20, status) IF (status.NE.0) THEN call wtferr(subname, version, status, & ' reading Coord-1 dataset (FITSIO/FTGBCL)') ierr = 1 goto 998 ENDIF ! ...... Coord-2 status = 0 call ftgcno(iunit, .false., col2, colnum, status) If (status.NE.0) THEN message = col2(:clenact(col2)) // & ' column not present' call wtferr(subname, version, status, message) ierr = 4 goto 998 ENDIF call ftgbcl(iunit, colnum, dummy20, dummy20, dummy20, & niy, doubdum, doubdum, enull, dummy20, status) IF (status.NE.0) THEN call wtferr(subname, version, status, & ' reading Coord-2 dataset (FITSIO/FTGBCL)') ierr = 1 goto 998 ENDIF ! ----- 998 if(ierr.NE.0) then call wterrm(subname, version, ' Punting') endif return end ! ------------------------------------------------------------------------- !+CK_COLRESP subroutine ck_colresp(iunit, chatter, & crspvers, & nen, nix, niy, ierr) IMPLICIT NONE integer chatter, ierr integer iunit integer nix, niy, nen character*(*) crspvers ! ! Description: ! Checks a COLLRESP extension in one of the formats ! conforming to the HDUVERS2='1.*.*' family (for DMA) ! !!! Note !!!! File is left open at the end ! ! Passed parameters ! iunit i : FORTRAN unit number of open RMF file ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! CRSPVERS o : String denoting OGIP HDUVERS2 family ! NEN o : Number of values in energy arrays ! NIX o : Number of values for coordinate-1 ! NIY o : Number of values for coordinate-2 ! IERR o : Error flag (0 = OK) ! ! Called Routines: ! subroutine FCECHO (FTOOLS) Writes to standard o/p ! subroutine FTGBCL (FITSIO) ! subroutine FTGCNO (FITSIO) ! subroutine FTGHBN (FITSIO) ! subroutine FTGKYS (FITSIO) ! subroutine GT_CSYSNMS (CALLIB) Gets coordinate sytem column names ! subroutine WT_FERRMSG (CALLIB) Standard error string writer ! ! Compilation & Linking ! link with FITSIO & CALLIB & FTOOLS ! ! Origin: ! Original ! ! Authors/Modification History: ! Ian M George (1.0.0;1994 Oct 17), original ! Ian M George (1.0.1;1995 Jun 28), replaced GT_CSYSNMS w/ GCSYNM ! Ian M George (1.0.2:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '1.0.2') !- ! Internals character(10) subname parameter (subname = 'ck_colresp') character(20) dummy20 double precision doubdum integer status, decimals, nfields integer irows, ncols, ivar integer colnum, enull, clenact parameter (nfields = 3, decimals = 6) character(8) dummy8 character(16) ttype(nfields), tform(nfields), tunits(nfields) character(20) csystem, col1, col2 character(40) comm character(80) message ! Initialize ierr = 0 status = 0 ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 20, 1, message) ! --- READING KEYWORD VALUES --- status = 0 call ftghbn(iunit, nfields, irows, ncols, ttype, tform, tunits, & comm, ivar, status) IF (status.NE.0) THEN call wtferr(subname, version, status, & ' reading binary header info') ierr = 2 goto 998 ENDIF IF (irows.GT.1) then call wtwarm(subname, version, 1, 1, & ' Dataset contains more than one row') ENDIF ! CRSPVERS ... crspvers = ' ' status = 0 call ftgkys(iunit, 'HDUVERS2', crspvers, comm, status) IF(status.ne.0) then call wtferr(subname, version, status, & ' reading HDUVERS2 keyword') status = 0 ENDIF IF (crspvers.EQ.' ') THEN call ftgkys(iunit, 'CRSPVERS', crspvers, comm, status) if(status.ne.0) then call wtferr(subname, version, status, & ' reading CRSPVERS keyword') status = 0 ENDIF ENDIF ! Coordinate system in use status = 0 call ftgkys(iunit, 'CSYSNAME', csystem, comm, status) call wtferr(subname, version, status, & ' reading CSYSNAME keyword') IF (status.EQ.202) THEN csystem = 'XMA_POL' call wtinfo(chatter, 20, 2, 'assuming CSYSNAME = XMA_POL') elseif(status.ne.0) then csystem = 'XMA_POL' call wtwarm(subname, version, chatter, 1, & 'Setting CSYSNAME = XMA_POL') ENDIF ! Sort out the coordinate system & column names call gcsynm(chatter, 0, csystem, & col1, col2, dummy8, dummy8, ierr) IF (ierr.NE.0) goto 998 ! --- READ IN THE DATA ! ... The energy grid ! ...... Lower bound status = 0 call ftgcno(iunit, .false., 'ENERG_LO', colnum, status) If (status.NE.0) THEN call wtinfo(chatter, 15, 2, & 'ENERG_LO column not present (will assume its a keyword)') nen = 1 goto 123 ENDIF call ftgbcl(iunit, colnum, dummy20, dummy20, dummy20, & nen, doubdum, doubdum, enull, dummy20, status) enull = 0 IF (status.NE.0) THEN call wterrm(subname, version, & 'problem reading ENERG_LO dataset') ierr = 1 goto 998 ENDIF ! ...... Upper bound 123 status = 0 call ftgcno(iunit, .false., 'ENERG_HI', colnum, status) If (status.NE.0) THEN call wtinfo(chatter, 15, 2, & 'ENERG_HI column not present (will assume its a keyword)') nen = 1 goto 124 ENDIF call ftgbcl(iunit, colnum, dummy20, dummy20, dummy20, & nen, doubdum, doubdum, enull, dummy20, status) enull = 0 IF (status.NE.0) THEN call wterrm(subname, version, & 'problem reading ENERG_HI dataset') ierr = 1 goto 998 ENDIF ! ... The spatial grid ! ...... Coord-1 124 status = 0 call ftgcno(iunit, .false., col1, colnum, status) If (status.NE.0) THEN message = col1(:clenact(col1)) // & ' column not present' call wtferr(subname, version, status, message) ierr = 4 goto 998 ENDIF call ftgbcl(iunit, colnum, dummy20, dummy20, dummy20, & nix, doubdum, doubdum, enull, dummy20, status) enull = 0 IF (status.NE.0) THEN call wtferr(subname, version, status, & ' reading Coord-1 dataset (FITSIO/FTGBCL)') ierr = 1 goto 998 ENDIF ! ...... Coord-2 status = 0 call ftgcno(iunit, .false., col2, colnum, status) If (status.NE.0) THEN message = col2(:clenact(col2)) // & ' column not present' call wtferr(subname, version, status, message) ierr = 4 goto 998 ENDIF call ftgbcl(iunit, colnum, dummy20, dummy20, dummy20, & niy, doubdum, doubdum, enull, dummy20, status) enull = 0 IF (status.NE.0) THEN call wtferr(subname, version, status, & ' reading Coord-2 dataset (FITSIO/FTGBCL)') ierr = 1 goto 998 ENDIF ! ----- 998 if(ierr.NE.0) then call wterrm(subname, version, ' Punting') endif return end ! ------------------------------------------------------------------------- !+SREMAP subroutine sremap(chatter, & msizx_in, msizy_in, & nx_in, ny_in, x_in, y_in, img_in, & msizx_out, msizy_out, & nx_out, ny_out, x_out, y_out, img_out, & qextrap, ierr) IMPLICIT NONE integer chatter, ierr integer msizx_in, msizy_in integer msizx_out, msizy_out integer nx_in, ny_in integer nx_out, ny_out real x_in(*), y_in(*) real img_in(msizx_in, msizy_in) real x_out(*), y_out(*) real img_out(msizx_out, msizy_out) logical qextrap ! Description ! Remaps spatial coordinate grid onto a new (linear) grid - extracted ! from DO_COL2IMG for dynamic memory allocation reasons ! ! Passed parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! MSIZX_IN i : max size of x-dimensions of input arrays ! MSIZY_IN i : max size of y-dimensions of input arrays ! NX_IN i : No elements in x-dimensions of input arrays ! NY_IN i : No elements in y-dimensions of input arrays ! X_IN i : Input X-coord grid points array (for IMG_IN) ! Y_IN i : Input Y-coord grid points array (for IMG_IN) ! IMG_IN i : Input Image ! MSIZX_OUT i : max size of x-dimensions of output arrays ! MSIZY_OUT i : max size of y-dimensions of output arrays ! NX_OUT i : No elements in x-dimensions of output arrays ! NY_OUT i : No elements in y-dimensions of output arrays ! X_OUT i : Requested remapped X-coord grid points array (4 IMG_OUT) ! Y_OUT i : Requested remapped Y-coord grid points array (4 IMG_OUT) ! IMG_OUT o : Remapped Output Image ! QEXTRAP i : Flag as to whether to extrapolate oustide X_IN,Y_IN ! IERR o : Error flag (0 = OK) ! ! Called Routines: ! subroutine BILINT (CALLIB) performs bilinear interpolation ! subroutine CLENACT (CALLIB) returns actual length of a string ! subroutine CRMVBLK (CALLIB) removes blanks from a string ! subroutine FCECHO (FTOOLS) writes to Standard Output ! ! Compilation & Linking ! link with FITSIO & CALLIB & FTOOLS ! ! Origin: ! Original ! ! Authors/Modification History: ! Ian M George (1.0.0:1994 Oct 17), original ! Ian M George (1.0.1:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '1.0.1') !- ! Internals character(6) subname parameter (subname = 'sremap') integer i, j, clenact real value character(20) string character(80) message logical qprob ! Initialize ierr = 0 qprob = .false. ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 20, 1, message) ! Check for sillies ! ... Funny inputs if((nx_in.GT.msizx_in).or.(ny_in.GT.msizy_in))then call wtwarm(subname, version, chatter, 1, & ' Incompatible array sizes') write(string, '(i12,a,i12)') msizx_in, 'x', msizy_in call crmvblk(string) message = ' I/p image dimensions : ' // & string(:clenact(string)) call wtinfo(chatter, 1, 2, message) write(string, '(i12,a,i12)') nx_in, 'x', ny_in call crmvblk(string) message = ' I/p array grid dimensions: ' // & string(:clenact(string)) call wtinfo(chatter, 1, 2, message) call wtinfo(chatter, 1, 1, & ' Attempting to continue using reduced grid') if(nx_in.GT.msizx_in) nx_in = msizx_in if(ny_in.GT.msizy_in) ny_in = msizy_in endif ! ... Funny outputs if((nx_out.GT.msizx_out).or.(ny_out.GT.msizy_out))then call wterrm(subname, version, & ' Requested dimensions too big') write(string, '(i12,a,i12)') msizx_out, 'x', msizy_out call crmvblk(string) message = ' Max allowed O/p image dimensions: ' // & string(:clenact(string)) call wtinfo(chatter, 1, 1, message) write(string, '(i12,a,i12)') nx_out, 'x', ny_out call crmvblk(string) message = ' Requested O/p image dimensions : ' // & string(:clenact(string)) call wtinfo(chatter, 1, 2, message) ierr = 1 goto 999 endif ! ... Perform the re-mapping do j = 1, ny_out do i = 1, nx_out call bilint(chatter, x_out(i), y_out(j), & msizx_in, msizy_in, & nx_in, ny_in, X_in, Y_in, img_in, & qextrap, value, ierr) if(ierr.EQ.0) then img_out(i, j) = value else img_out(i, j) = -99 if(ierr.NE.-1) then qprob = .true. endif endif enddo enddo ! Dump info if desired if(qprob) call wtwarm(subname, version, chatter, 2, & 'At least one problem with remapping') if(.NOT.qprob) then ierr = 0 call wtinfo(chatter, 15, 1, & 'successfully remapped the dataset') endif ! Error checking 999 if(ierr.NE.0) then call wterrm(subname, version, 'Aborting') endif return end ! ------------------------------------------------------------------------- !+PRE_REMAP subroutine pre_remap(chatter, & xmin, xmax, nxbins, & ymin, ymax, nybins, & x_out, y_out, & ierr) IMPLICIT NONE integer chatter, ierr integer nxbins, nybins real xmin, xmax, ymin, ymax real x_out(*), y_out(*) ! Description ! Pathetic little guy to populate the X & Y grids ! (necessary for the Dynamic Memory Allocation in calling routine) ! ! Passed Parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! XMIN i : Min X-coord of remapped grid ! XMAX i : Max X-coord of remapped grid ! NXBINS i : No X-coord elements remapped grid ! YMIN i : Min Y-coord of remapped grid ! YMAX i : Max Y-coord of remapped grid ! NYBINS i : No Y-coord elements remapped grid ! X_OUT o : Requested remapped X-coord grid points array (4 IMG_OUT) ! Y_OUT o : Requested remapped Y-coord grid points array (4 IMG_OUT) ! IERR o : Error flag (0 = OK) ! ! Called Routines ! None ! ! Compilation & Linking ! No special requirements ! ! Origin: ! Original ! ! Author/Modification History ! Ian M George (1.0.0:1994 Oct 18) original ! Ian M George (1.0.1:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '1.0.1') !- ! Internals character(9) subname parameter (subname = 'pre_remap') integer i real delx_tmp, dely_tmp character(80) message ! Initialize ierr = 0 ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 30, 1, message) ! Just do it delx_tmp = (xmax - xmin) / nxbins dely_tmp = (ymax - ymin) / nybins do i = 1, nxbins x_out(i) = xmin + (i - 1) * delx_tmp enddo do i = 1, nybins y_out(i) = ymin + (i - 1) * dely_tmp enddo return end ! ------------------------------------------------------------------------- !+POST_REMAP subroutine post_remap(chatter, & msizx_in, msizy_in, & x_in, y_in, img_in, & msizx_out, msizy_out, & nx_out, ny_out, & x_out, y_out, img_out, & ierr) IMPLICIT NONE integer chatter, ierr integer msizx_in, msizy_in integer msizx_out, msizy_out integer nx_out, ny_out real x_in(*), y_in(*) real img_in(msizx_in, msizy_in) real x_out(*), y_out(*) real img_out(msizx_out, msizy_out) ! Description ! Pathetic little guy to populate the X & Y grids ! (necessary for the Dynamic Memory Allocation in calling routine) ! ! Passed Parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! MSIZX_IN i : max size of x-dimensions of input arrays ! MSIZY_IN i : max size of y-dimensions of input arrays ! X_IN i : Input X-coord grid points array (for IMG_IN) ! Y_IN i : Input Y-coord grid points array (for IMG_IN) ! IMG_IN i : Input Image ! MSIZX_OUT i : max size of x-dimensions of output arrays ! MSIZY_OUT i : max size of y-dimensions of output arrays ! NX_OUT i : No elements in x-dimensions of output arrays ! NY_OUT i : No elements in y-dimensions of output arrays ! X_OUT i : Requested remapped X-coord grid points array (4 IMG_OUT) ! Y_OUT i : Requested remapped Y-coord grid points array (4 IMG_OUT) ! IMG_OUT o : Remapped Output Image ! IERR o : Error flag (0 = OK) ! ! Called Routines ! None ! ! Compilation & Linking ! No special requirements ! ! Origin: ! Original ! ! Author/Modification History ! Ian M George (1.0.0:1994 Oct 18) original ! Ian M George (1.0.1:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '1.0.1') !- ! Internals character(10) subname parameter (subname = 'post_remap') integer i, j character(80) message ! Initialize ierr = 0 ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 30, 1, message) ! Just do it do i = 1, nx_out X_out(i) = x_in(i) enddo do j = 1, nx_out Y_out(j) = y_in(j) enddo do j = 1, ny_out do i = 1, nx_out img_out(i, j) = img_in(i, j) enddo enddo return end ! ------------------------------------------------------------------------- !+POST_OBF subroutine post_obf(chatter, & nk_comm, comment, & cname1, cname2, & cunit1, cunit2, & icoord1, icoord2, x, y, & ierr) IMPLICIT NONE integer chatter, ierr integer nk_comm integer icoord1, icoord2 real x(*), y(*) character*(*) comment(*) character*(*) cname1, cname2 character*(*) cunit1, cunit2 ! Description ! Pathetic little guy to dump some junk to screen & comment cards ! (necessary for the Dynamic Memory Allocation in calling routine) ! ! Passed Parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! {incomplete} ! IERR o : Error flag (0 = OK) ! ! Called Routines ! subroutine CLENACT (CALLIB) returns actual length of a string ! subroutine FCECHO (FTOOLS) writes to Standard Output ! ! Compilation & Linking ! No special requirements ! ! Origin: ! Original ! ! Author/Modification History ! Ian M George (1.0.0: 94 Oct 18) original ! Ian M George (1.0.1:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '1.0.1') !- ! Internals character(8) subname parameter (subname = 'post_obf') integer clenact character(80) message, string ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 30, 1, message) ! Write user info message = 'input image: ' // & cname1(:MIN(10, clenact(cname1))) // & ' vs ' // & cname2(:MIN(10, clenact(cname2))) // & ' (OBSFACTOR dataset)' call wtinfo(chatter, 9, 1, message) write(string, '(g15.5,a,g15.5)')& X(1), ':', X(icoord1) call crmvblk(string) write(message, '(a,a, i4,a,a,a,a)')& cname1(:MIN(10, clenact(cname1))), ' has ', & icoord1, ' elements (', string(:clenact(string)), & ') ', cunit1 call wtinfo(chatter, 9, 2, message) write(string, '(g15.5,a,g15.5)')& y(1), ':', Y(icoord2) call crmvblk(string) write(message, '(a,a, i4,a,a,a,a)')& cname2(:MIN(10, clenact(cname2))), ' has ', & icoord1, ' elements (', string(:clenact(string)), & ') ', cunit2 call wtinfo(chatter, 9, 2, message) ! Add a few comments comment(nk_comm + 1) = & ' COL2IMG identified i/p dataset as an OBSFACTOR' nk_comm = nk_comm + 1 comment(nk_comm + 1) = & ' axis-1 label: ' // cname1 comment(nk_comm + 2) = & ' axis-1 units: ' // cunit1 write(message, '(g15.5,a,g15.5)')& X(1), ':', X(icoord1) call crmvblk(message) comment(nk_comm + 3) = & ' axis-1 range: ' // message nk_comm = nk_comm + 3 comment(nk_comm + 1) = & ' axis-2 label: ' // cname2 comment(nk_comm + 2) = & ' axis-2 units: ' // cunit2 write(message, '(g15.5,a,g15.5)')& Y(1), ':', X(icoord2) call crmvblk(message) comment(nk_comm + 3) = & ' axis-2 range: ' // message nk_comm = nk_comm + 3 return end ! ------------------------------------------------------------------------- !+POST_COL subroutine post_col(chatter, & nk_comm, comment, & ienerg, energ_lo, energ_hi, eunits, & cname1, cname2, & cunit1, cunit2, & icoord1, icoord2, x, y, & ierr) IMPLICIT NONE integer chatter, ierr integer nk_comm integer icoord1, icoord2, ienerg real x(*), y(*) real energ_lo(*), energ_hi(*) character*(*) comment(*) character*(*) cname1, cname2 character*(*) cunit1, cunit2, eunits ! Description ! Pathetic little guy to dump some junk to screen & comment cards ! (necessary for the Dynamic Memory Allocation in calling routine) ! ! Passed Parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! {incomplete} ! IERR o : Error flag (0 = OK) ! ! Called Routines ! subroutine CLENACT (CALLIB) returns actual length of a string ! subroutine FCECHO (FTOOLS) writes to Standard Output ! ! Compilation & Linking ! No special requirements ! ! Origin: ! Original ! ! Author/Modification History ! Ian M George (1.0.0: 94 Oct 18) original ! Ian M George (1.0.1:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '1.0.1') !- ! Internals character(8) subname parameter (subname = 'post_col') integer clenact character(80) message, string ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 30, 1, message) ! Write user info message = ' input image: ' // & cname1(:MIN(10, clenact(cname1))) // & ' vs ' // & cname2(:MIN(10, clenact(cname2))) // & ' (COLLRESP dataset)' call wtinfo(chatter, 9, 1, message) write(string, '(g15.5,a,g15.5)')& X(1), ':', X(icoord1) call crmvblk(string) write(message, '(a,a, i4,a,a,a,a)')& cname1(:MIN(10, clenact(cname1))), ' has ', & icoord1, ' elements (', string(:clenact(string)), & ') ', cunit1 call wtinfo(chatter, 9, 2, message) write(string, '(g15.5,a,g15.5)')& y(1), ':', Y(icoord2) call crmvblk(string) write(message, '(a,a, i4,a,a,a,a)')& cname2(:MIN(10, clenact(cname2))), ' has ', & icoord1, ' elements (', string(:clenact(string)), & ') ', cunit2 call wtinfo(chatter, 9, 2, message) ! Add a few comments comment(nk_comm + 1) = & ' COL2IMG identified i/p dataset as a COLLRESP' nk_comm = nk_comm + 1 write(message, '(g15.5,a,g15.5)')& energ_lo(1), ':', energ_hi(ienerg) call crmvblk(message) comment(nk_comm + 1) = & ' energy range: ' // message // ' ' // & eunits(:clenact(eunits)) nk_comm = nk_comm + 1 comment(nk_comm + 1) = & ' axis-1 label: ' // cname1 comment(nk_comm + 2) = & ' axis-1 units: ' // cunit1 write(message, '(g15.5,a,g15.5)')& X(1), ':', X(icoord1) call crmvblk(message) comment(nk_comm + 3) = & ' axis-1 range: ' // message nk_comm = nk_comm + 3 comment(nk_comm + 1) = & ' axis-2 label: ' // cname2 comment(nk_comm + 2) = & ' axis-2 units: ' // cunit2 write(message, '(g15.5,a,g15.5)')& Y(1), ':', X(icoord2) call crmvblk(message) comment(nk_comm + 3) = & ' axis-2 range: ' // message nk_comm = nk_comm + 3 return end ! ------------------------------------------------------------------------- !+CK_SGRID subroutine ck_sgrid(chatter, & cname1, cname2, & icoord1, icoord2, x, y, & ratio, qxlin, qylin, & delx, dely, & ierr) IMPLICIT NONE integer chatter, ierr integer icoord1, icoord2 real x(*), y(*) real ratio character cname1, cname2 logical qxlin, qylin ! Description ! Pathetic little guy to checkout the spatial grid ! (necessary for the Dynamic Memory Allocation in calling routine) ! ! Passed Parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! {incomplete} ! IERR o : Error flag (0 = OK) ! ! Called Routines ! None ! ! Compilation & Linking ! No special requirements ! ! Origin: ! Original ! ! Author/Modification History ! Ian M George (1.0.0: 94 Oct 18) original ! Ian M George (2.0.0: 94 Dec 14) addede delx & dely as passed parameters ! Ian M George (2.0.1:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '2.0.1') !- ! Internals character(8) subname parameter (subname = 'ck_sgrid') integer i, clenact real delta, delx, dely character(80) message ! Initialization ierr = 0 qxlin = .true. qylin = .true. ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 30, 1, message) ! Do it if(icoord1.LE.1) then message = 'Only 1 element along ' // & cname1(:clenact(cname1)) call wtwarm(subname, version, chatter, 1, message) goto 654 endif delx = X(2) - X(1) do i = 2, icoord1 delta = X(i) - X(i - 1) if(MOD(delta, delx).GT.(delx / ratio)) then qxlin = .false. goto 654 endif enddo 654 if(icoord2.LE.1) then message = 'Only 1 element along ' // & cname2(:clenact(cname2)) call wtwarm(subname, version, chatter, 1, message) goto 655 endif dely = Y(2) - Y(1) do i = 2, icoord2 delta = Y(i) - Y(i - 1) if(MOD(delta, dely).GT.(dely / ratio)) then qylin = .false. goto 655 endif enddo 655 continue return end ! ------------------------------------------------------------------------- !+ WT_CKYS subroutine wt_ckys(chatter, ounit, & cname1, cunit1, qxlin, delx, x, & cname2, cunit2, qylin, dely, y, & ierr) IMPLICIT NONE integer chatter, ierr, ounit real delx, dely real x(*), y(*) character*(*) cname1, cname2, cunit1, cunit2 logical qxlin, qylin ! Description ! Pathetic little guy to write a few keywords to the output file ! (necessary for the Dynamic Memory Allocation in calling routine) ! ! Passed Parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! {incomplete} ! IERR o : Error flag (0 = OK) ! ! Called Routines ! None ! ! Compilation & Linking ! No special requirements ! ! Origin: ! Original ! ! Author/Modification History ! Ian M George (1.0.0: 94 Oct 18) original ! Ian M George (1.0.1:1996 Feb 04), added wtinfo & friends character(7) version parameter (version = '1.0.1') !- ! Internals character(7) subname parameter (subname = 'wt_ckys') integer status, clenact real inreal integer decimals character(80) message, string ! Initialization ierr = 0 decimals = 6 ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 30, 1, message) ! Do it status = 0 call ftpkys(ounit, 'CTYPE1', cname1, & 'Name of 1st coordinate axis', & status) call wtferr(subname, version, status, & ' problem writing CTYPE1 keyword') status = 0 call ftpkys(ounit, 'CUNIT1', cunit1, & 'Units of 1st coordinate axis', & status) call wtferr(subname, version, status, & ' problem writing CUNIT1 keyword') if(qxlin) then status = 0 string = 'Pixel Size of 1st coordinate axis (' // & cunit1(:clenact(cunit1)) // ')' call ftpkye(ounit, 'CDELT1', delx, decimals, & string, & status) call wtferr(subname, version, status, & ' problem writing CDELT1 keyword') status = 0 inreal = 1.0 call ftpkye(ounit, 'CRPIX1', inreal, decimals, & 'Axis-1 coord of Reference Pixel', & status) call wtferr(subname, version, status, & ' problem writing CRPIX1 keyword') status = 0 string = cname1(:clenact(cname1)) // & ' Value at Reference Pixel (' // & cunit1(:clenact(cunit1)) // ')' inreal = X(1) call ftpkye(ounit, 'CRVAL1', inreal, decimals, & string, & status) call wtferr(subname, version, status, & ' problem writing CRVAL1 keyword') endif status = 0 call ftpkys(ounit, 'CTYPE2', cname2, & 'Name of 2nd coordinate axis', & status) call wtferr(subname, version, status, & ' problem writing CTYPE2 keyword') status = 0 call ftpkys(ounit, 'CUNIT2', cunit2, & 'Units of 2nd coordinate axis', & status) call wtferr(subname, version, status, & ' problem writing CUNIT2 keyword') if(qylin) then status = 0 string = 'Pixel Size of 2nd coordinate axis (' // & cunit2(:clenact(cunit2)) // ')' call ftpkye(ounit, 'CDELT2', dely, decimals, & string, & status) call wtferr(subname, version, status, & ' problem writing CDELT2 keyword') status = 0 inreal = 1.0 call ftpkye(ounit, 'CRPIX2', inreal, decimals, & 'Axis-2 coord of Reference Pixel', & status) call wtferr(subname, version, status, & ' problem writing CRPIX2 keyword') status = 0 string = cname2(:clenact(cname2)) // & ' Value at Reference Pixel (' // & cunit2(:clenact(cunit2)) // ')' inreal = Y(1) call ftpkye(ounit, 'CRVAL2', inreal, decimals, & string, & status) call wtferr(subname, version, status, & ' problem writing CRVAL2 keyword') endif return end ! ------------------------------------------------------------------------- !+ CK_EGRID subroutine ck_egrid(chatter, & nk_comm, comment, & ienerg, energ_lo, energ_hi, eunits, & eminstr, emaxstr, emin, emax, & ierr) IMPLICIT NONE integer chatter, ierr integer nk_comm, ienerg integer maxcomm parameter (maxcomm = 100) real emin, emax real energ_lo(*), energ_hi(*) character*(*) eminstr, emaxstr character*(*) comment(maxcomm), eunits ! Description ! Check-out routine for energy grid ! (necessary for the Dynamic Memory Allocation in calling routine) ! ! Passed Parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! {incomplete} ! IERR o : Error flag (0 = OK) ! ! Called Routines ! None ! ! Compilation & Linking ! No special requirements ! ! Origin: ! Original ! ! Author/Modification History ! Ian M George (1.0.0: 94 Oct 19) original ! Ian M George (1.0.1: 1996 Feb 04) add wtinfo & friends character(7) version parameter (version = ' 1.0.1 ') !- ! Internals character(8) subname parameter (subname = 'ck_egrid') integer clenact character(80) message ! Initialization ierr = 0 ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 30, 1, message) ! Translate the user-defined energy ranges ! ... lower bound if(eminstr.EQ.'%') then emin = energ_lo(1) else read(eminstr, *, iostat = ierr) emin if(ierr.NE.0) then call wtferr(subname, version, ierr, & ' Unable to translate energy bound ') message = ' EMIN = ' // eminstr(:clenact(eminstr)) call wtinfo(chatter, 1, 1, message) goto 986 endif endif ! ... upper bound if(emaxstr.EQ.'%') then emax = energ_hi(ienerg) else read(emaxstr, *, iostat = ierr) emax if(ierr.NE.0) then call wtferr(subname, version, ierr, & ' Unable to translate energy bound ') message = ' EMAX = ' // emaxstr(:clenact(emaxstr)) call wtinfo(chatter, 1, 1, message) goto 986 endif endif ! Final error check 986 if(ierr.NE.0) then call wterrm(subname, version, ' Punting') endif 987 continue return end ! ------------------------------------------------------------------------- !+ ESQUISH subroutine esquish(chatter, & maxenerg, ienerg, energ_lo, energ_hi, ework, & maxcoord1, maxcoord2, icoord1, icoord2, collresp, & emin, emax, image, & ierr) IMPLICIT NONE integer chatter, ierr integer maxenerg, maxcoord1, maxcoord2 integer ienerg, icoord1, icoord2 real emin, emax real energ_lo(*), energ_hi(*), ework(*) real collresp(maxenerg, maxcoord1, maxcoord2) real image(maxcoord1, maxcoord2) ! Description ! Compresses a COLLRESP dataset in the energy dimension ! ! Passed Parameters ! CHATTER i : chattiness flag for o/p (5 quite,10 normal,>20 silly) ! {incomplete} ! IERR o : Error flag (0 = OK) ! ! Called Routines ! {incomplete} ! ! Compilation & Linking ! CALLIB, FTOOLS ! ! Origin: ! Original ! ! Author/Modification History ! Ian M George (1.0.0: 94 Oct 19) original ! Ian M George (1.1.0: 94 Oct 25), Bug fix to energy remapping ! Ian M George (2.0.0: 95 Jun 27), Major overhaul ! Ian M George (2.0.1: 1996 Feb 04) add wtinfo & friends character(7) version parameter (version = ' 2.0.1 ') !- ! Internals character(8) subname parameter (subname = 'ck_egrid') integer i, j, k real acc, value(1) character(80) message ! Initialization acc = 1.E-3 ierr = 0 ! Give user info if requested message = ' using ' // subname // ' ' // version call wtinfo(chatter, 30, 1, message) ! Do the stuff do k = 1, icoord2 do j = 1, icoord1 do i = 1, ienerg ework(i) = collresp(i, j, k) enddo call rmap1d(chatter, & ienerg, energ_lo, energ_hi, ework, & 1, emin, emax, value, & 3, acc, ierr) if(ierr.NE.0) then goto 953 endif image(j, k) = value(1) enddo enddo ! Final error check 953 if(ierr.NE.0) then call wterrm(subname, version, 'Aborting') endif return end ! -------------------------------------------------------------------------