**==clcrsp.spg  processed by SPAG 4.50J  at 17:44 on 20 Jan 1996
 
      SUBROUTINE CLCRSP(N_energies, Resp_energies, N_channels, 
     &                  Ch_bounds, Rsp_min, Ngroup, Max_tot_groups,
     &                  Max_elements, Disperse, Res_reln, Fwhm, 
     &                  N_r_energies, N_lines, LineEn, LineDt, Ichanb, 
     &                  Ichane, Resp_matrix, Resp_work, Num_elements, 
     &                  Num_tot_groups, Ierr)

      IMPLICIT NONE
 
      INTEGER N_energies , N_channels, N_r_energies, N_lines
      INTEGER Max_elements , Max_tot_groups
 
      REAL Resp_energies(0:N_energies) , Resp_matrix(Max_elements)
      REAL Ch_bounds(N_channels,2), LineEn(N_r_energies)
      REAL Resp_work(N_channels), LineDt(3,N_lines,N_r_energies)
      REAL Rsp_min, Fwhm
 
      INTEGER Ngroup(N_energies)
      INTEGER Ichanb(Max_tot_groups) , Ichane(Max_tot_groups)
      INTEGER Num_elements , Num_tot_groups , Ierr

      CHARACTER res_reln*(*)

      LOGICAL disperse

c This routine actually calculates the response matrix assuming a
c Gaussian response
c Arguments :
c     N_energies       i     i: Number of response energy bins
c     Resp_energies    r     i: Response energy ranges
c     N_channels       i     i: Number of channels
c     Ch_bounds        r     i: Channel energy boundaries
c     Rsp_min          r     i: The minimum response value that will be stored
c     Ngroup           i     r: The number of response groups for each energy
c     Max_tot_groups   i     i: The maximum number of total response groups
c     Ichanb           i     r: Start channel for each response group
c     Ichane           i     r: End channel for each response group
c     Max_elements     i     i: Maximum number of response elements
c     Disperse         c*(*) i: Whether dispersive or non-dispersive
c     Fwhm             r     i: FWHM for fiducial spectral resolution
c     N_r_energies     i     i: Number of resolution points if tabulated
c     N_lines          i     i: Number of lines (photopeak + escape)
c     LineEn           i     i: Energies for tabulated response data
c     LineDt           i     i: Tabulated line (delta energy, 
c                               resolution, norm) data
c     Resp_matrix      r     r: Response matrix elements
c     Resp_work        r     w: Work array for response for current energy
c     Num_elements     i     r: Number of response matrix elements calculated
c     Num_tot_groups   i     r: Total number of response groups calculated
c     Ierr             i     r: Error status    0 == OK

 
      REAL centroid, sigma , norm, energy, effic, gfact
 
      INTEGER igroup, iresp, i, j, k, icen

      CHARACTER wrtstr*255
 
      LOGICAL counting
 
      INTEGER bisearch
      REAL     grspvl
      EXTERNAL grspvl, bisearch

      Ierr = 0

c Initialize the Ichanb and Ichane arrays

      DO i = 1, Max_tot_groups
         Ichanb(i) = 0
         Ichane(i) = 0
      ENDDO

c Loop round the response energies
 
      igroup = 0
      iresp = 0
 
      DO i = 1 , N_energies
 
         energy = (Resp_energies(i-1)+Resp_energies(i))/2

c Initialize the resp_work array

         DO j = 1, N_channels
            resp_work(j) = 0.
         ENDDO

c Loop over the lines (photopeak + escapes)

         DO k = 1, N_lines

c Get the centroid, sigma, and norm for this line

            CALL GETSIG(energy, k, Disperse, Res_reln, Fwhm, 
     &                  N_r_energies, N_lines, lineEn, lineDt,
     &                  centroid, sigma, norm)
            IF ( sigma .LT. 0 ) THEN
               Ierr = 3
               RETURN
            ENDIF

            IF ( i/100. .EQ. FLOAT(i/100) ) THEN
               WRITE(wrtstr, '(2(2x,i8), 4(2x,1pe12.5))') i, k, energy, 
     &            centroid, sigma, norm
               CALL fcecho(wrtstr)
            ENDIF

c Find the channel containing centroid

            icen = bisearch(N_channels, Ch_bounds, centroid)

c Load the line into the resp_work array

            j = icen
            gfact = 1.
            DO WHILE ( gfact .GT. 0. .AND. j .GE. 1 )
               gfact = 0.
               IF ( sigma .GT. 0. ) THEN 
                  gfact = GRSPVL(centroid,Ch_bounds(j,1),
     &                           Ch_bounds(j,2),sigma)
               ELSE
                  IF ( centroid .GT. Ch_bounds(j,1) .AND. 
     &                 centroid .LE. Ch_bounds(j,2) ) THEN
                     gfact = 1.
                  ENDIF
               ENDIF
               Resp_work(j) = Resp_work(j) + norm*gfact
               j = j - 1
            ENDDO

            j = icen + 1
            gfact = 1.
            DO WHILE ( gfact .GT. 0. .AND. j .LE. N_channels )
               gfact = 0.
               IF ( sigma .GT. 0. ) THEN 
                  gfact = GRSPVL(centroid,Ch_bounds(j,1),
     &                           Ch_bounds(j,2),sigma)
               ELSE
                  IF ( centroid .GT. Ch_bounds(j,1) .AND. 
     &                 centroid .LE. Ch_bounds(j,2) ) THEN
                     gfact = 1.
                  ENDIF
               ENDIF
               Resp_work(j) = Resp_work(j) + norm*gfact
               j = j + 1
            ENDDO

c End loop over lines

         ENDDO

c Transfer the work array into the standard compressed response structure

         Ngroup(i) = 0
         counting = .FALSE.
         effic = 0.
         DO j = 1 , N_channels

c  if response greater than minimum acceptable then include
 
            IF ( Resp_work(j).GT.Rsp_min ) THEN
 
               iresp = iresp + 1
               IF ( iresp.GT.Max_elements ) THEN
                  wrtstr = 'Too many response elements'//
     &                     ' - increase max_elements'
                  CALL fcecho(wrtstr)
                  Ierr = 1
                  RETURN
               ENDIF
               Resp_matrix(iresp) = Resp_work(j)
               effic = effic + Resp_work(j)
 
c  if not currently in a group then start one
 
               IF ( .NOT.counting ) THEN
                  igroup = igroup + 1
                  IF ( igroup.GT.Max_tot_groups ) THEN
                     wrtstr = 'Too many response groups'//
     &                        ' - increase max_tot_groups'
                     CALL fcecho(wrtstr)
                     Ierr = 2
                     RETURN
                  ENDIF
                  Ichanb(igroup) = j
                  counting = .TRUE.
                  Ngroup(i) = Ngroup(i) + 1
               ENDIF
 
               Ichane(igroup) = j
 
c  if response not greater than acceptable minimum and in a group 
c  then end that group
 
 
            ELSEIF ( counting ) THEN
               counting = .FALSE.
 
            ENDIF
 
         ENDDO

c  end loop over energies

      ENDDO
 
      Num_elements = iresp
      Num_tot_groups = igroup
 
      RETURN
      END
 
 
c ******************************************************************

      FUNCTION bisearch(N_channels, Ch_bounds, centroid)


      INTEGER bisearch, N_channels
      REAL Ch_bounds(N_channels,2)
      REAL centroid

c Function to do a bisection search for which element of the Ch_bounds array
c the input centroid lies in.

      INTEGER low, high

      LOGICAL increase

      IF ( Ch_bounds(1,1) .LT. Ch_bounds(N_channels, 2) ) THEN
         increase = .TRUE.
      ELSE
         increase = .FALSE.
      ENDIF

      IF ( increase ) THEN
         IF ( centroid .LT. Ch_bounds(1,1) .OR. 
     &        centroid .GT. Ch_bounds(N_channels,2) ) THEN
            bisearch = -1
            RETURN
         ENDIF
      ELSE
         IF ( centroid .GT. Ch_bounds(1,1) .OR. 
     &        centroid .LT. Ch_bounds(N_channels,2) ) THEN
            bisearch = -1
            RETURN
         ENDIF
      ENDIF

      low = 0
      high = N_channels

      DO WHILE ( (high-low) .GT. 1 )

         bisearch = (low + high) / 2
         IF ( (increase .AND. centroid .GT. Ch_bounds(bisearch,1)) .OR.
     &        (.NOT.increase .AND. centroid .LT. Ch_bounds(bisearch,1))
     &           ) THEN
            low = bisearch
         ELSE
            high = bisearch
         ENDIF

      ENDDO

      bisearch = low

      RETURN
      END

