
      SUBROUTINE rgsrmfsmooth

c Program to smooth the RGS rmf for extended sources based on an image (probably from
c Chandra).
c
c  Parameters are
c     rmffil           Input rmf
c     imgfil           File with Image data
c     outfil           Output rmf
c     chatter          Chatter level
c     qclob            If true overwrite output file if it already exists
c
C Local variables

      IMPLICIT NONE

      REAL areascal, lo_thresh

      INTEGER status, chatter, ilun, block, flchan
      INTEGER nenerg, numgrp, numelt, numchn
      INTEGER onumgrp, onumelt, order

      CHARACTER*255 rmffil, imgfil, outfil, contxt
      CHARACTER*72 rmfname, rmfname_comm, rmfversn
      CHARACTER*72 hduclas3, telescop, instrume
      CHARACTER*72 detnam, filter, chantype
      CHARACTER*72 ebdname, ebdname_comm, ebdvers
      LOGICAL qclob, isorder

C Pointers for dynamic arrays

      INTEGER ichannel, ieblow, iebhigh
      INTEGER ienlow, ienhigh, ingroups, igchanb, igchann
      INTEGER iorder, irspelt
      INTEGER iear, irsprow
      INTEGER iongroups, iogchanb, iogchann, iorspelt

      INTEGER fcstln
      EXTERNAL fcstln

c  common block for dynamic memory
 
      INTEGER          MEMI(1)
      REAL             MEMR(1)
      EQUIVALENCE (MEMI, MEMR)
      COMMON /MEM/ MEMR

      status = 0

      CALL fcecho("  rgsrmfsmooth v1.01  (02/22/08)")

C Get the parameters

      CALL gprgs(rmffil, imgfil, order, outfil, chatter, qclob, status)
      contxt = 'Failure in gprgs'
      IF ( status .NE. 0 ) GOTO 999

C Open the input RMF

      CALL getlun(ilun)
      CALL ftopen(ilun, rmffil, 0, block, status)
      contxt = 'Failed to open '//rmffil(:fcstln(rmffil))
      IF ( status .NE. 0 ) GOTO 999

C Move to the RMF extension

      CALL ftmnhd(ilun, 2, 'MATRIX', 0, status)
      IF ( status .NE. 0 ) THEN
         status = 0
         CALL ftmnhd(ilun, 2, 'SPECRESP MATRIX', 0, status)
         contxt = 'Failed to move to MATRIX extension'
         IF ( status .NE. 0 ) GOTO 999
      ENDIF

C Get the size of the arrays required for the RMF

      CALL rmfsz(ilun, chatter, nenerg, numgrp, numelt, status)
      contxt = 'Failed to get size of RMF arrays'
      IF ( status .NE. 0 ) GOTO 999

      WRITE(contxt,'(a,i8,a,i8,a,i8)') 'RMF has nenerg = ', nenerg,
     & ' numgrp = ', numgrp, ' numelt = ', numelt
      CALL fcecho(contxt)

C Grab memory for arrays for RMF

      ienlow = 0
      ienhigh = 0
      ingroups = 0
      igchanb = 0
      igchann = 0
      irspelt = 0
      iorder = 0
      CALL udmget(nenerg, 6, ienlow, status)
      CALL udmget(nenerg, 6, ienhigh, status)
      CALL udmget(nenerg, 4, ingroups, status)
      CALL udmget(numgrp, 4, igchanb, status)
      CALL udmget(numgrp, 4, igchann, status)
      CALL udmget(numelt, 6, irspelt, status)
      CALL udmget(nenerg, 4, iorder, status)
      contxt = 'Failed to get memory for RMF arrays'
      IF ( status .NE. 0 ) GOTO 999

C Read the input RMF

      CALL rdrmf5(ilun, chatter, rmfname, rmfname_comm, .FALSE., nenerg,
     &            numgrp, numelt, rmfversn, hduclas3, telescop, 
     &            instrume, detnam, filter, areascal, chantype, flchan,
     &            numchn, nenerg, numgrp, numelt, MEMR(ienlow), 
     &            MEMR(ienhigh), MEMI(ingroups), MEMI(igchanb), 
     &            MEMI(igchann), isorder, MEMI(iorder), MEMR(irspelt), 
     &            lo_thresh, status)
      contxt = 'Failed to read RMF'
      IF ( status .NE. 0 ) GOTO 999

C Grab memory for arrays for EBOUNDS

      ichannel = 0
      ieblow = 0
      iebhigh = 0
      CALL udmget(numchn, 4, ichannel, status)
      CALL udmget(numchn, 6, ieblow, status)
      CALL udmget(numchn, 6, iebhigh, status)
      contxt = 'Failed to get memory for EBOUNDS arrays'
      IF ( status .NE. 0 ) GOTO 999

C Move to the EBOUNDS extension

      CALL ftmnhd(ilun, 2, 'EBOUNDS', 0, status)
      contxt = 'Failed to move to EBOUNDS extension'
      IF ( status .NE. 0 ) GOTO 999

C Read the input EBOUNDS

      CALL rdebd4(ilun, chatter, ebdname, ebdname_comm, numchn, 
     &           telescop, instrume, detnam, filter, areascal,
     &           chantype, flchan, numchn, MEMI(ichannel),
     &           MEMR(ieblow), MEMR(iebhigh), ebdvers, status)
      contxt = 'Failed to read EBOUNDS extension'
      IF ( status .NE. 0 ) GOTO 999

C Grab the memory for the work arrays

      iear = 0
      irsprow = 0
      CALL udmget(nenerg+1, 6, iear, status)
      CALL udmget(numchn, 6, irsprow, status)
      contxt = 'Failed to get memory for work arrays'
      IF ( status .NE. 0 ) GOTO 999

C and for the output response arrays

      onumgrp = 4 * numgrp
      onumelt = 4 * numelt
      WRITE(contxt,'(a,i8,a,i8,a)') 
     &  'Allocating memory for output response : ', onumgrp, 
     &  ' groups and ', onumelt, ' elements'
      CALL fcecho(contxt)

      iongroups = 0
      iogchanb = 0
      iogchann = 0
      iorspelt = 0
      CALL udmget(nenerg,   4, iongroups, status)
      CALL udmget(onumgrp, 4, iogchanb, status)
      CALL udmget(onumgrp, 4, iogchann, status)
      CALL udmget(onumelt, 6, iorspelt, status)
      contxt = 'Failed to get memory for output RMF arrays'
      IF ( status .NE. 0 ) GOTO 999

C Calculate the new RMF

      CALL dosmooth(imgfil, order, nenerg, numgrp, numelt, numchn, 
     &              MEMI(ingroups), MEMI(igchanb), MEMI(igchann),
     &              MEMR(irspelt), MEMR(ieblow), MEMR(iebhigh),
     &              lo_thresh, MEMR(iear), MEMR(irsprow), onumgrp,
     &              onumelt, MEMI(iongroups), MEMI(iogchanb),
     &              MEMI(iogchann), MEMR(iorspelt), status)

      contxt = 'Failure in dosmooth'
      IF ( status .NE. 0 ) GOTO 999

      WRITE(contxt,'(a,i8,a,i8,a,i8)') 'Output RMF has nenerg = ', 
     &      nenerg, ' numgrp = ', onumgrp, ' numelt = ', onumelt
      CALL fcecho(contxt)
   
C Write the output RMF - copying the header and other extensions from the input file

      CALL wrtrsp(ilun, outfil, chatter, qclob, rmfname, 
     &            rmfname_comm, nenerg, onumgrp, onumelt, rmfversn,
     &            hduclas3, telescop, instrume, detnam, filter, 
     &            areascal, chantype, flchan, numchn, MEMR(ienlow), 
     &            MEMR(ienhigh),MEMI(iongroups), MEMI(iogchanb), 
     &            MEMI(iogchann), isorder, MEMI(iorder), 
     &            MEMR(iorspelt), lo_thresh, status)
      contxt = 'Failure in wrtrsp'
      IF ( status .NE. 0 ) GOTO 999

 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         CALL fcecho(contxt)
         WRITE(contxt,'(a,i8)') 'status = ', status
         CALL fcecho(contxt)
      ENDIF

      RETURN
      END

c ************************************************************************************

      SUBROUTINE gprgs(rmffil, imgfil, order, outfil, chatter, qclob, 
     &                 status)

      IMPLICIT NONE

      CHARACTER*(*) rmffil, imgfil, outfil

      INTEGER order, chatter, status

      LOGICAL qclob

      CHARACTER contxt*72

      status = 0

      CALL UCLGST('rmffil', rmffil, status)
      contxt = 'Failed to get rmffil parameter'
      IF ( status .NE. 0 ) GOTO 999

      CALL UCLGST('imgfil', imgfil, status)
      contxt = 'Failed to get imgfil parameter'
      IF ( status .NE. 0 ) GOTO 999

      CALL UCLGSI('order', order, status)
      contxt = 'Failed to get order parameter'
      IF ( status .NE. 0 ) GOTO 999

      CALL UCLGST('outfil', outfil, status)      
      contxt = 'Failed to get outfil parameter'
      IF ( status .NE. 0 ) GOTO 999

      CALL UCLGSI('chatter', chatter, status)
      contxt = 'Failed to get chatter parameter'
      IF ( status .NE. 0 ) GOTO 999

      CALL UCLGSB('clobber', qclob, status)
      contxt = 'Failed to get clobber parameter'
      IF ( status .NE. 0 ) GOTO 999

 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         CALL fcecho(contxt)
         WRITE(contxt,'(a,i8)') 'status = ', status
         CALL fcecho(contxt)
      ENDIF

      RETURN
      END



c ************************************************************************************

      SUBROUTINE dosmooth(imgfil, order, nenerg, numgrp, numelt, numchn,
     &                    ngroups, gchanb, gchann, rspelt, eblow,
     &                    ebhigh, lo_thresh, ear, rsprow, onumgrp,
     &                    onumelt, ongroups, ogchanb, ogchann, orspelt,
     &                    status)

      IMPLICIT NONE

      INTEGER nenerg, numgrp, numelt, numchn
      INTEGER onumgrp, onumelt

      REAL rspelt(numelt), eblow(numchn), ebhigh(numchn)
      REAL ear(0:numchn), rsprow(numchn) 
      REAL orspelt(onumelt)

      REAL lo_thresh

      INTEGER ngroups(nenerg), gchanb(numgrp), gchann(numgrp)
      INTEGER ongroups(nenerg), ogchanb(onumgrp), ogchann(onumgrp)

      INTEGER order, status

      CHARACTER*(*) imgfil


      INTEGER ie, igrp, irsp, ogrp, orsp, ich, ig
      CHARACTER wrtstr*256
      LOGICAL counting


C Set up xspec energies - note that channels are in order of decreasing energy
C in the XMM RGS response

      DO ich = numchn, 1, -1
         ear(numchn-ich) = eblow(ich)
      ENDDO
      ear(numchn) = ebhigh(1)

C Loop over response energies

      igrp = 0
      irsp = 0
      ogrp = 0
      orsp = 0

      status = 0

      DO ie = 1, nenerg

C Create "spectrum" from the response for this energy

         DO ich = 1, numchn
            rsprow(ich) = 0.0
         ENDDO

         DO ig = 1, ngroups(ie)
            igrp = igrp + 1
            DO ich = gchanb(igrp), gchanb(igrp)+gchann(igrp)-1
               irsp = irsp + 1
               rsprow(numchn-ich+1) = rspelt(irsp)
            ENDDO
         ENDDO  

C Run xspec model to convolve spectrum

         CALL rgsxsrc(ear, numchn, FLOAT(order), 0, rsprow, imgfil)

C Place the results for this energy in the output response

         ongroups(ie) = 0
         counting = .FALSE.

         DO ich = 1, numchn

            IF ( rsprow(numchn-ich+1) .GT. lo_thresh ) THEN

               orsp = orsp + 1
               IF ( orsp .GT. onumelt ) THEN
                  WRITE(wrtstr,'(4(a,i8))') 
     &              'Too many response elements ', orsp, ' > ', onumelt,
     &              '  at ie = ', ie, ' ich = ', ich
                  CALL fcecho(wrtstr)
                  status = -1
                  RETURN
               ENDIF
               orspelt(orsp) = rsprow(numchn-ich+1)
               IF ( .NOT.counting ) THEN
                  ogrp = ogrp + 1
                  IF ( ogrp .GT. onumgrp ) THEN
                     WRITE(wrtstr,'(4(a,i8))') 
     &                'Too many response groups ', ogrp, ' > ', onumgrp,
     &                '  at ie = ', ie, ' ich = ', ich
                     CALL fcecho(wrtstr)
                     status = -1
                     RETURN
                  ENDIF
                  ogchanb(ogrp) = ich
                  ogchann(ogrp) = 0
                  counting = .TRUE.
                  ongroups(ie) = ongroups(ie) + 1
               ENDIF
               ogchann(ogrp) = ogchann(ogrp) + 1

            ELSEIF ( counting ) THEN

               counting = .FALSE.

            ENDIF

         ENDDO

C End loop over response energies

      ENDDO

      onumgrp = ogrp
      onumelt = orsp


      RETURN
      END

c ************************************************************************************

      SUBROUTINE wrtrsp(ilun, outfil, chatter, qclob, rmfname, 
     &                  rmfname_comm, nenerg, numgrp, numelt, rmfversn,
     &                  hduclas3, telescop, instrume, detnam, filter, 
     &                  areascal, chantype, flchan, numchn, enlow, 
     &                  enhigh, ngroups, gchanb, gchann, isorder, 
     &                  order, rspelt, lo_thresh, status)

      IMPLICIT NONE

      INTEGER nenerg, numgrp, numelt, numchn

      REAL rspelt(numelt), enlow(nenerg), enhigh(nenerg)
      INTEGER ngroups(nenerg), gchanb(numgrp), gchann(numgrp)
      INTEGER order(*)

      REAL areascal, lo_thresh
      INTEGER ilun, status, flchan, chatter

      CHARACTER*(*) rmfname, rmfname_comm, rmfversn, hduclas3, telescop
      CHARACTER*(*) instrume, detnam, filter, chantype, outfil

      LOGICAL isorder, qclob

      INTEGER olun
      CHARACTER*255 contxt

      INTEGER fcstln
      EXTERNAL fcstln

      EXTERNAL BADDUMMY, GOODDUMMY

      status = 0

C Open the output response matrix

      CALL getlun(olun)
 100  CONTINUE
      CALL ftinit(olun, outfil, 1, status)
      IF ( status .EQ. 105 .AND. qclob ) THEN
         status = 0
         CALL xdelfil(outfil, status)
         GOTO 100
      ENDIF
      contxt = 'Failure creating '//outfil(:fcstln(outfil))
      IF ( status .NE. 0 ) GOTO 999

C Move the input file to the MATRIX extension

      CALL ftmnhd(ilun, 2, rmfname, 0, status)
      contxt = 'Failed to move to MATRIX extension'
      IF ( status .NE. 0 ) GOTO 999

C Copy all HDUs ahead of the MATRIX extension from the input to the output
C file

      CALL ftcpfl(ilun, olun, 1, 0, 0, status)
      contxt = 'Failed to copy extensions before MATRIX'
      IF ( status .NE. 0 ) GOTO 999

C Write the MATRIX extension to the RMF file

      CALL wtrmf5(olun, chatter, rmfname, rmfname_comm, 0, rmfname, 0, 
     &            rmfname, rmfversn, hduclas3, telescop, instrume,
     &            detnam, filter, areascal, chantype, flchan, numelt,
     &            numchn, nenerg, numgrp, enlow, enhigh, ngroups, 
     &            gchanb, gchann, isorder, order, rspelt, lo_thresh,
     &            status)
      contxt = 'Failed to write the MATRIX extension'
      IF ( status .NE. 0 ) GOTO 999

C Copy additional keywords in the MATRIX extension

      CALL copyhead(ilun, olun, .TRUE., .TRUE., BADDUMMY, GOODDUMMY, 
     &              status)
      contxt = 'Failed to copy additional keywords in MATRIX'
      IF ( status .NE. 0 ) GOTO 999

C Copy all HDUs following the MATRIX extension from the input to the output
C file

      CALL ftcpfl(ilun, olun, 0, 0, 1, status)
      contxt = 'Failed to copy extensions after MATRIX'
      IF ( status .NE. 0 ) GOTO 999
      
C Close the output file

      CALL ftclos(olun, status)
      contxt = 'Failed to close output file'
      IF ( status .NE. 0 ) GOTO 999

 999  CONTINUE
      IF ( status .NE. 0 ) THEN
         CALL fcecho(contxt)
         WRITE(contxt,'(a,i8)') 'status = ', status
         CALL fcecho(contxt)
      ENDIF

      RETURN
      END


