C*****************************************************************************
C TASK NAME: speconvrmf
C
C FILE NAME: speconvrmf.f
C
C DESCRIPTION: FTOOL to convert an EGRET .resp (or resp.) file into a FITS
C     file in XSPEC format.  The .resp file is generated by Patrick Nolan's
C     SPECTRAL code, and contains the response matrix of the EGRET detector.
C
C AUTHOR/DATE: Peter J.T. Leonard and Jeffrey M.S. Silvis, HSTX, 03/97
C
C NOTES:
C
C USAGE:
C     Host: speconvrmf
C     IRAF: speconvrmf
C
C ROUTINES IMPLEMENTED IN THIS FILE:
C     subroutine speconvrmf  - top level subroutine, called by IRAF or host
C                            C wrapper in heconvrmf.c
C     subroutine spread_resp - reads in EGRET file from SPECTRAL
C     subroutine spset_bins  - set up the incident energy bins
C     subroutine spwrite_rsp - writes out FITS file in XSPEC format
C     subroutine spcmprsp    - compress a response matrix to XSPEC format
C
C MODIFICATION HISTORY:
C*****************************************************************************

C*****************************************************************************
C SUBROUTINE: speconvrmf
C
C DESCRIPTION: Program to convert an EGRET .resp (or resp.) file into a FITS
C     file in XSPEC format.  The .resp file is generated by Patrick Nolan's
C     SPECTRAL code, and contains the response matrix of the EGRET detector.
C
C AUTHOR/DATE: John R. Mattox, 08/96
C
C NOTES:
C
C ARGUMENTS:
C     rspfil - SPECTRAL .resp filename
C     rmffil - XSPEC .rmf filename
C
C PRIMARY LOCAL VARIABLES:
C
C CALLED ROUTINES:
C     subroutine spread_resp - reads in EGRET file from SPECTRAL
C     subroutine spwrite_rsp - writes out FITS file in XSPEC format
C
C MODIFICATION HISTORY:
C     03/97, cosmetic changes, PJTL and JMSS, HSTX
C*****************************************************************************
 
      subroutine speconvrmf (rspfil, rmffil, evclass)

      implicit none

      character*(*)  rspfil, rmffil
      integer*4 ierr, j, NENER, NCHAN, status, nk_hist, nk_comm
      integer*4 MAX_ENERGIES, MAX_CHANNELS, MAX_ELEMENTS
      integer*4 MAX_GROUPS, MAX_TOT_GROUPS, evclass
      parameter (MAX_ENERGIES=300, MAX_CHANNELS=64, MAX_GROUPS=6)
      parameter (MAX_ELEMENTS=MAX_ENERGIES*MAX_CHANNELS)
      parameter (MAX_TOT_GROUPS=MAX_ENERGIES*MAX_GROUPS)
      integer*4 ngroup(MAX_ENERGIES), ichanb(MAX_TOT_GROUPS)
      integer*4 ichane(MAX_TOT_GROUPS)
 
      real*4 e_min(MAX_CHANNELS), e_max(MAX_CHANNELS), RSP_MIN
      real*4 matrix(64,300), rspmatrix(64,300), rspmat(MAX_ELEMENTS)
      real*4 resp_energies(0:MAX_ENERGIES), ch_bounds(2,MAX_CHANNELS)
      parameter (RSP_MIN=1.e-6)
 
      character*80 hist
      character*20 tlscpe, instrm, detnam, filter, hduclas3

      data tlscpe/'EGRET'/
      data instrm/'EGRET'/
      data detnam/'spark chamber'/
      data filter/'EGRET'/
      data hduclas3/'FULL'/
 
C Setup common block for taskname - necessary for fcerr.
      character*80 msg
      character*40 taskname
      common /task/taskname
      taskname = 'speconvrmf'
 
      status = 0

      hist =
     &      'XSPEC RMF file produced from EGRET response file:'//rspfil
      nk_hist = 1
      nk_comm = 0

C Read in the EGRET response matrix data.
      call spread_resp (rspfil, rspmatrix, resp_energies, ch_bounds,
     &                NENER, evclass, NCHAN)

C Write channel bounds.
      do j = 1, NCHAN
         e_min(j) = ch_bounds(1,j)
         e_max(j) = ch_bounds(2,j)
      enddo

C Write XSPEC response matrix.
      call spwrite_rsp (NCHAN, NENER, MAX_ELEMENTS, MAX_GROUPS,
     &                MAX_TOT_GROUPS, RSP_MIN, rspmatrix, matrix,
     &                resp_energies, e_min, e_max, tlscpe, instrm,
     &                detnam, filter, hduclas3, rmffil, nk_hist,
     &                hist, nk_comm, hist, rspmat, ngroup, ichanb,
     &                ichane, ierr)

      IF (ierr.NE.0) THEN
         msg = 'Problem writing output file!'
         call fcerr (msg)
         msg = 'Aborting program.'
         call fcerr (msg)
         call exit (1)
      END IF

      return
c      call exit (0)

      end
 
C*****************************************************************************
C SUBROUTINE: spread_resp
C
C DESCRIPTION: Reads in EGRET .resp file produced by SPECTRAL.
C
C AUTHOR/DATE: Keith A. Arnaud, 05/89
C              Mark A. Fardal, 08/90
C              Patrick L. Nolan, 01/91
C              Thomas A. McGlynn, 1993
C              John R. Mattox, 07/94
C
C NOTES:
C
C ARGUMENTS:
c     filnam        - name of file containing response matrix
C     input_resp    - response matrix
C     resp_energies - energy bins for the response matrix
C     ch_bounds     - energy bounds for channels
C     nchan         - number of incoming channels
C     nener         - number of measured energy channels
C
C PRIMARY LOCAL VARIABLES:
C
C CALLED ROUTINES:
C     subroutine spset_bins - constructs incident bins
C  
C MODIFICATION HISTORY:
C     03/97, cosmetic changes, PJTL and JMSS, HSTX
C******************************************************************************

      subroutine spread_resp (filnam, input_resp, resp_energies,
     &                      ch_bounds, nchan, evclass, nener)

      implicit none

C     MAX_ENERGIES    number of energies in response matrix
C     MAX_CHANNELS    number of channels in response matrix
C     MAX_ELEMENTS    maximum number of elements in response matrix
C     MAX_GROUPS      maximum number of groups per energy in response matrix
C     MAX_TOT_GROUPS  maximum number of groups in response matrix
C     RSP_MIN         minimum value of response that will be included
      integer*4 MAX_ENERGIES, MAX_CHANNELS, evclass
      integer*4 MAX_ELEMENTS, MAX_GROUPS, MAX_TOT_GROUPS
      parameter (MAX_ENERGIES=300,MAX_CHANNELS=64)
      parameter (MAX_ELEMENTS=MAX_CHANNELS*50, MAX_GROUPS=6)
      parameter (MAX_TOT_GROUPS=MAX_ENERGIES*MAX_GROUPS)
      real*4 RSP_MIN
      parameter (RSP_MIN=1.e-6)
             
C Detector info to put in file ID record.
C     DETID    detector name
C     NBIN     number of binned channels
C     RSPAREA  area in response
C     NRANGE   number of energies in response
C     E0       starting energy for response
C     AREA     energy-independent effective area
      character*16 DETID
      parameter (DETID='CGRO EGRET      ')
      integer*4 NBIN, NRANGE
      parameter (NBIN=MAX_CHANNELS, NRANGE=MAX_ENERGIES)
      real*4 RSPAREA, E0, AREA
      parameter (RSPAREA=1.0, E0=0.5, AREA=1.0)

C History data for response file.
      integer*4 NHIST
      parameter (NHIST=1)
      character*80 history
      parameter (history='EGRET PHA/RESPONSE v1.0')

C Background and response filenames for PHA file.
      character*80 BKGFILE, RSPFILE
      parameter (BKGFILE='none', RSPFILE='dummy.rsp')

C Names of programs that generate rsp and pha files.
      character*80 PHAPROG, RSPPROG
      parameter (PHAPROG='egretpha', RSPPROG='egretrsp')

      integer*4 i, j, nchan, nener
      real*4 input_resp(MAX_CHANNELS,MAX_ENERGIES)
      real*4 resp_energies(0:MAX_ENERGIES), ch_bounds(2,MAX_CHANNELS)
      character*(*) filnam

      integer*4 first, last, spacing, jmax
C First bin in incident spectrum.
      parameter (first=-9)
C Last bin in incident spectrum.
      parameter (last=260)
C Incident bins between 20 MeV and 20 GeV.
      parameter (spacing=200)
C Maximum number of observed bins.
      parameter (jmax=15)
 
      integer*4 lumat, luout, nornmax
      parameter (lumat=4, luout=8, nornmax=40)
 
C Estimate of class C bin effective width.
      real*4 cwidth
      parameter (cwidth=32.0)
      
C     NREGNS:               Number of acceptance cones
C     EOBS*(0:JMAX)         Observed energy bin bounds
C     P(FIRST:LAST,0:JMAX)  Instrument response matrix
C     EREGN(0:JMAX)         Energy intervals for cones
C     SRCLON,SRCLAT         Source coordinates
C     CONETHETA(0:JMAX)     Acceptance cone angles
C     CELEQ:                T = source coordinates are in RA/DEC
      integer*4 nregns, istatus, nfile, norns, nornstr, nornend
      integer*4 norntasc, ihdutype, idum0, idum1, idum2, idum3, jtop, jx
      real*8 ornstart(nornmax), ornend(nornmax), orntime(2,nornmax)
      real*4 eobs1(0:jmax), eobs2(0:jmax), p(first:last,0:jmax)
      real*4 eregn(0:jmax), conetheta(0:jmax), emin, emax, lnemin
      real*4 lnemax, srclon, srclat, cutoff, cutangle(0:jmax)
      character*130 response
      character*30 comment
      character*20 comline
      character*12 version
      character*6 enclass
      character*1 comchar
      equivalence (comchar,comline)
      logical celeq, orntasc(nornmax), anyf

      data nfile/0/
          
C Information used in subroutine ADA.
C     EINC(FIRST-1:LAST)    Incident energy bin bounds
C     EC(FIRST:LAST)        Geometric centers of EINC bins
C     DE(FIRST:LAST)        Widths of EINC bins
C     Q(FIRST:LAST,0:JMAX)  Instrument response matrix
C     IVFIRST:              First bin used in analysis
      common/adablk/einc, ec, de, q, ilconv, enorm, ivfirst, idunit
      integer*4 ilconv(3), ivfirst, idunit
      real*4 einc(first-1:last), ec(first:last), de(first:last)
      real*4 q(first:last,0:jmax), enorm
 
C Quantities related to the number of counts observed.
C     PHS              Phase length of source interval
C     PHB              Phase length of background interval
C     NMEAS*(0:JMAX)   Observed counts
C     NBEXP*(0:JMAX)   Estimate of background counts
C     NSHYP*(0:JMAX)   Source estimate (NMEAS-NBEXP)
C      (* = 1: basic binning; * = 2: rebinned copy)
      common/counts/phs, phb, nmeas1, nmeas2, nbexp1, nbexp2, sbexp1,
     &              sbexp2, nshyp1, nshyp2, sigma1, sigma2, fudge1,
     &              fudge2, sfudge1, sfudge2
      real*4 phs, phb, nmeas1(0:jmax), nbexp1(0:jmax), sbexp1(0:jmax)
      real*4 nshyp1(0:jmax), sigma1(0:jmax), fudge1(0:jmax)
      real*4 sfudge1(0:jmax), nmeas2(0:jmax), nbexp2(0:jmax)
      real*4 sbexp2(0:jmax), nshyp2(0:jmax), sigma2(0:jmax)
      real*4 fudge2(0:jmax), sfudge2(0:jmax)
 
C Setup common block for taskname - necessary for fcerr.
      character*80 msg
      character*40 taskname
      common /task/taskname
      taskname = 'speconvrmf'

      version = 'Version 2.5'
      istatus = 0


C Read the header of the matrix.
      response = filnam
      call ftopen (lumat, response, 0, idum0, istatus)
      call ftgkyj (lumat, 'first'  , idum1,   comment, istatus)
      call ftgkyj (lumat, 'last'   , idum2,   comment, istatus)
      call ftgkyj (lumat, 'spacing', idum3,   comment, istatus)
      call ftgkyj (lumat, 'nregns',  nregns,  comment, istatus)
      call ftgkyl (lumat, 'celeq',   celeq,   comment, istatus)
      call ftgkye (lumat, 'srclon',  srclon,  comment, istatus)
      call ftgkye (lumat, 'srclat',  srclat,  comment, istatus)
      call ftgkye (lumat, 'cutoff',  cutoff,  comment, istatus)
      call ftgkys (lumat, 'EN-CLAS', enclass, comment, istatus)
 
      jtop = nregns
      if (idum1.ne.first.or.idum2.ne.last.or.idum3.ne.spacing
     &                                       .or.jtop.gt.jmax) then
         msg =
     &       'Parameters in MATRIX do not agree with those in SPECTRAL!'
         call fcerr (msg)
         msg = 'Aborting program.'
         call fcerr (msg)
         call exit (1)
      end if
 
csb-02/11      call getenv ('EVCLASS', evclass)

      if ((evclass.eq.1.and.enclass.ne.'A').or.(evclass.eq.2
     &                                   .and.enclass.ne.'A+B+C')) then
         msg = 'Event energy classes do not match response matrix!'
         call fcerr (msg)
         msg = 'Aborting program.'
         call fcerr (msg)
         call exit (1)
      end if

C Read orientation intervals and TASC coincidence mode.
      call ftgkyj (lumat, 'norns', norns, comment, istatus)
      call ftgknd (lumat, 'ornstr', 1, norns, ornstart,nornstr, istatus)
      call ftgknd (lumat, 'ornend', 1, norns, ornend,  nornend, istatus)
      call ftgknl (lumat, 'ortasc', 1, norns, orntasc, norntasc,istatus)
      do i = 1 , norns
	 orntime(1,i) = ornstart(i)
	 orntime(2,i) = ornend(i)
      end do

C Read matrix values.
      call ftg2de (lumat, 0, 0, last-first+1 ,last-first+1, nregns+1,
     &             p, anyf, istatus)
 
C Read observed energy, selection cones, zenith cuts, and empirical fudge factors.
      call ftmrhd (lumat, 1, ihdutype, istatus)
      call ftgcve (lumat, 1, 1, 1, jtop+1, 0, eregn,     anyf, istatus)
      call ftgcve (lumat, 2, 1, 1, jtop+1, 0, conetheta, anyf, istatus)
      call ftgcve (lumat, 3, 1, 1, jtop+1, 0, cutangle,  anyf, istatus)
      call ftgcve (lumat, 4, 1, 1, jtop+1, 0, fudge1,    anyf, istatus)
      call ftgcve (lumat, 5, 1, 1, jtop+1, 0, sfudge1,   anyf, istatus)
      do i = 0, jtop
	 eobs1(i) = eregn(i)
      end do

C Adjust effective areas according to empirical fudge factors.
      do i = first, last
         do j = 0, nregns
            p(i,j) = p(i,j) / fudge1(j)
         end do
      end do

      emin = eobs1(0)
      emax = eobs1(jtop)
      lnemin = alog(emin)
      lnemax = alog(emax)

C Set up the incident energy bins.
      call ftmrhd (lumat, 1, ihdutype, istatus)
      call ftgcve (lumat, 1, 1, 1, last-first+2, 0, einc, anyf, istatus)
      call spset_bins (einc, ec, de, first, last)
      eobs2(0) = emin

C Read the select output file.
      call ftmrhd (lumat, 1, ihdutype, istatus)
      do j = 0, jtop
         nmeas1(j) = 0.0
         nbexp1(j) = 0.0
      end do

C Number of incoming energy channels.
      nchan = last - first + 1
C Number of measured energy channels.
      nener = nregns

C All energies are in keV from here on.
      resp_energies(0) = einc(first-1) * 1000.
      do i = first, last
            j = i - first + 1
            resp_energies(j) = einc(i) * 1000.
      end do

      do i = 1, nener
         ch_bounds(1,i) = eobs1(i-1) * 1000.
         ch_bounds(2,i) = eobs1(i) * 1000.
      end do

      do i = 1, nener
         do j = first, last
            jx = j - first + 1
            input_resp(i,jx) = p(j,i)
         enddo
      enddo

      return

      end
 
C*****************************************************************************
C SUBROUTINE: spset_bins
C
C DESCRIPTION: Set up the incident energy bins.
C
C AUTHOR/DATE: Mark A. Fardal, 09/90
C              Patrick L. Nolan, 01/91
C
C NOTES: SPECTRAL requires a vector of incident photon energies at which
C     to evaluate the model.  The range must extend beyond the observed
C     energy range on either side.  The bins are constructed to be even
C     on a log scale.  EINC(0) and EINC(SPACING) are set to the observed
C     energy range; then FIRST and LAST determine the extent the overlap.
C     The bins in MATRIX must be identical to those used here.
C
C ARGUMENTS:
C     EINC(FIRST-1:LAST) - incident energy bins
C     EC(FIRST:LAST)     - bin centers
C     DE(FIRST:LAST)     - bin widths
C     FIRST, LAST        - indices of first, last incident bins
C
C PRIMARY LOCAL VARIABLES:
C
C CALLED ROUTINES: None.
C
C MODIFICATION HISTORY:
C     03/97, cosmetic changes, PJTL and JMSS, HSTX
C*****************************************************************************

      subroutine spset_bins (einc, ec, de, first, last)

      implicit none

      integer first, last, i
      real einc(first-1:last), ec(first:last), de(first:last)

C Incident bins.
      do i = first, last
C Center of bin i.
	 ec(i) = sqrt(einc(i-1) * einc(i))
C Width of bin i.
         de(i) = einc(i) - einc(i-1)             
      end do

      return

      end
 
C*****************************************************************************
C SUBROUTINE: spwrite_rsp
C
C DESCRIPTION: Writes out a FITS format response file.  The response matrix
C     is passed as XSPEC compressed format.
C
C AUTHOR/DATE: Keith A. Arnaud, 10/94
C              John R. Mattox, 08/96
C
C NOTES:
C
C ARGUMENTS:
C     NCHAN    - number of channels in response matrix
C     NENERG   - number of energy bins in response matrix
C     MXELEM   - maximum number of non-zero response elements
C     MXGRPS   - maximum number of groups at a given energy
C     MXTGRP   - maximum number of total groups
C     RSP_MIN  - minimum value of response that is stored
C     matrix   - response matrix
C     energies - energy bins for the response
C     e_min    - nominal lower energy bounds for channels
C     e_max    - nominal upper energy bounds for channels
C     tlscpe   - telescope name 
C     instrm   - instrument name 
C     detnam   - detector name 
C     filter   - filter name 
C     hduclas3 - if = REDIST   => photon redistribution matrix (only)
C                   = DETECTOR => convolved with detector effects (only)
C                   = FULL     => convolved with all effects (det+optic)
C     rspfil   - output filename
C     nk_hist  - number of history records
C     hist     - history records (optional)
C     nk_comm  - number of comment records
C     comment  - comment records (optional)
C     rspmat   - non-zero response elements
C     ngroup   - number of contiguous channel sets
C     ichanb   - start channel of a group
C     nchang   - number of channels in a group
C     ierr     - error code (0 = OK)
C
C PRIMARY LOCAL VARIABLES:
C
C CALLED ROUTINES:
C     subroutine spcmprsp - compress a response matrix to XSPEC format
C
C MODIFICATION HISTORY:
C     03/97, cosmetic changes, PJTL and JMSS, HSTX
C*****************************************************************************

      subroutine spwrite_rsp (NCHAN, NENERG, MXELEM, MXGRPS, MXTGRP,
     &                      RSP_MIN, rspmatrix, matrix, energies, e_min,
     &                      e_max, tlscpe, instrm, detnam, filter,
     &                      hduclas3, rspfil, nk_hist, hist, nk_comm,
     &                      comment, rspmat, ngroup, ichanb, nchang,
     &                      ierr)
 
      IMPLICIT NONE

      INTEGER*4 NCHAN, NENERG, MXELEM, MXGRPS, MXTGRP
      INTEGER*4 nk_hist, nk_comm, ierr, i, j
      INTEGER*4 ngroup(NENERG), ichanb(MXTGRP), nchang(MXTGRP)

      REAL*4 rspmatrix(64,300), matrix(NCHAN,NENERG), rspmat(MXELEM)
      REAL*4 RSP_MIN, energies(0:NENERG), e_min(NCHAN), e_max(NCHAN)

      CHARACTER*(*) rspfil, hist(nk_hist), comment(nk_comm)
      CHARACTER*(*) tlscpe, instrm, detnam, filter, hduclas3

      INTEGER*4 ounit, chatter, nelem, ntgrps

      chatter = 10
      ierr = 0

      do i = 1, NCHAN
         do j = 1, NENERG
             matrix(i,j) = rspmatrix(i,j)
         enddo
      enddo

C Calculate the compressed matrix.
      call spcmprsp (NCHAN, NENERG, MXELEM, MXGRPS, MXTGRP, RSP_MIN,
     &             matrix, rspmat, ngroup, ichanb, nchang, nelem, 
     &             ntgrps)

      
C Open the output FITS file.
      ounit = 15
      call ftinit (ounit, rspfil, 1, ierr) 

C Write the primary header.
      call ftpdef (ounit, 8, 0, 0, 0, 1, ierr)

C Write the basic primary array keywords.
      call ftphpr (ounit, .TRUE., 8, 0, 0, 0, 1, .TRUE., ierr)

C Write out the additional keywords about the creation of the FITS file.
      call ftpkys (ounit, 'CONTENT', 'RESPONSE',
     & 'spectral response matrix', ierr)
      call ftpkys (ounit, 'ORIGIN', 'NASA/GSFC',
     & 'origin of FITS file', ierr)

C Write the energy bounds extensions.
      call wtebd2 (ounit, chatter, nk_hist, hist, nk_comm, comment, 
     &             '1.1.0', tlscpe, instrm, detnam, filter, 1.0, 
     &             1, NCHAN, e_min, e_max, ierr) 

C Write the response matrix extensions.  Note that at present WTRMF1 wants
C the uncompressed matrix passed and not the compressed matrix.
      call wtrmf1 (ounit, chatter, nk_hist, hist, nk_comm, comment,
     &             '1.1.0', hduclas3, tlscpe, instrm, detnam, filter, 
     &             1.0, NCHAN, NCHAN, NENERG, NENERG, energies(0), 
     &             energies(1), MXGRPS, ngroup, ichanb, nchang, 
     &             matrix, RSP_MIN, ierr)

C Close the output file.
      call ftclos (ounit, ierr)

      RETURN
 
      END
 
C*****************************************************************************
C SUBROUTINE: spcmprsp
C
C DESCRIPTION: Routine to compress a response matrix to XSPEC format.
C
C AUTHOR/DATE: Keith A. Arnaud, 03/89
C
C NOTES:
C
C ARGUMENTS:
C     NCHAN   - number of channels in response matrix
C     NENERG  - number of energy bins in response matrix
C     MXELEM  - maximum number of non-zero response elements
C     MXGRPS  - maximum number of groups at a given energy
C     MXTGRP  - maximum number of total groups
C     RSP_MIN - minimum value of response that is stored
C     matrix  - response matrix
C     rspmat  - non-zero response elements
C     ngroup  - number of contiguous channel sets
C     ichanb  - start channel of a group
C     nchang  - number of channels in a group
C     nelem   - number of non-zero response elements
C     ntgrps  - number of groups
C
C PRIMARY LOCAL VARIABLES:
C
C CALLED ROUTINES: None.
C
C MODIFICATION HISTORY:
C     03/97, cosmetic changes, PJTL and JMSS, HSTX
C*****************************************************************************
 
      subroutine spcmprsp (NCHAN, NENERG, MXELEM, MXGRPS, MXTGRP,
     &           RSP_MIN, matrix, rspmat, ngroup, ichanb, nchang,
     &            nelem, ntgrps)
 
      IMPLICIT NONE

      INTEGER NCHAN, NENERG, MXELEM, MXGRPS, MXTGRP
      INTEGER ngroup(NENERG), ichanb(MXTGRP), nchang(MXTGRP)
      INTEGER nelem, ntgrps, i, j, igroup, iresp

      REAL RSP_MIN, effic, rsp
      REAL matrix(NCHAN,NENERG), rspmat(MXELEM)

      LOGICAL counting
 
C Setup common block for taskname - necessary for fcerr.
      character*80 msg
      character*40 taskname
      common /task/taskname
      taskname = 'speconvrmf'

C Create the response matrix - start by looping over energies.
      igroup = 0
      iresp = 0
      DO i = 1, NENERG

C Loop over channels.
         effic = 0
         ngroup(i) = 0
         counting = .FALSE.
         DO j = 1, NCHAN

            rsp = matrix(j, i)

C If response greater than minimum acceptable then include.
            IF (rsp.GT.RSP_MIN) THEN

               iresp = iresp + 1
 
               IF (iresp.GT.MXELEM) THEN
                  msg = 'Too many response elements - increase MXELEM!'
                  call fcerr (msg)
                  msg = 'Aborting program.'
                  call fcerr (msg)
                  call exit (1)
               ENDIF
 
               rspmat(iresp) = rsp
               effic = effic + rsp

C If not currently in a group then start one.
               IF (.NOT.counting) THEN
                  igroup = igroup + 1
 
                  IF (igroup.GT.MXTGRP) THEN
                     msg = 'Too many response groups - increase MXTGRP!'
                     call fcerr (msg)
                     msg = 'Aborting program.'
                     call fcerr (msg)
                     call exit (1)
                  ENDIF
 
                  ichanb(igroup) = j
                  counting = .TRUE.
                  ngroup(i) = ngroup(i) + 1
 
                  IF (ngroup(i).GT.MXGRPS) THEN
                     msg = 'Too many response groups - increase MXGRPS!'
                     call fcerr (msg)
                     msg = 'Aborting program.'
                     call fcerr (msg)
                     call exit (1)
                  ENDIF
 
               ENDIF

C If response not greater than acceptable minimum and in group then end group.
            ELSE

               IF (counting) THEN
                  counting = .FALSE.
                  nchang(igroup) = j - ichanb(igroup) + 1
               ENDIF

            ENDIF

         ENDDO

         IF (counting) THEN
            nchang(igroup) = NCHAN - ichanb(igroup) + 1
         ENDIF

      ENDDO

      nelem = iresp
      ntgrps = igroup
      
      RETURN
 
      END
