SUBROUTINE DOGFT (LU,LAST,ARRAY,SYSTEM,INDEX)
************************************************************************
*     set up Gaussian fit simulating the beam
* inputs :
*     LU      integer*2       i/o lu
*     LAST    integer*2       last valid point in the array
*     ARRAY   real*4          data array
*     SYSTEM  character*4     array of system identifiers
*     INDEX   integer*2       index identifying system in use
*
      PARAMETER  (MAXI = 6)
      LOGICAL*2  DUALBEAM
      REAL*4     ARRAY(1), FLUXESTA(11), FLUXESTB(11)
      CHARACTER  SYSTEM(MAXI)*4
*
*     subtract linear fit to baseline of array
      CALL LNOFF (LAST,ARRAY,RMS)
      WRITE (lu,'(9x,''rms noise = '',1pg12.6,'' K'')') RMS
*
*     define the mid beam points in  the data array
      IF (SYSTEM(INDEX) .EQ. '  6D') THEN
          MIDBEAMA = LAST / 2 - 8
          MIDBEAMB = LAST / 2 + 8
          DUALBEAM = .TRUE.
      ELSE IF (SYSTEM(INDEX) .EQ. '3.6D') THEN
          MIDBEAMA = LAST / 2 - 13
          MIDBEAMB = LAST / 2 + 13
          DUALBEAM = .TRUE.
      ELSE
*         for single feeds
          MIDBEAM = LAST / 2
          DUALBEAM = .FALSE.
      END IF
*
      IF (DUALBEAM) THEN
*         make 11 estimates of flux for each beam
*         error in flux is rms/2 per beam owing to gaussian fit
          J = 0
          DO I = MIDBEAMA - 5, MIDBEAMA + 5
              J = J + 1
              CALL GFLUX (ARRAY,I,FLUXESTA(J))
          END DO
          J = 0
          DO I = MIDBEAMB - 5, MIDBEAMB + 5
              J = J + 1
              CALL GFLUX (ARRAY,I,FLUXESTB(J))
          END DO
*
*         find (signed) maximum flux in each beam
          FLUXA = FLUXESTA(1)
          FLUXB = FLUXESTB(1)
          DO J = 2, 11
              FLUXA = MAX (FLUXA, FLUXESTA(J))
              FLUXB = MIN (FLUXB, FLUXESTB(J))
          END DO
          FLUX = (FLUXA - FLUXB) / 2
          FLUXERROR = RMS / ( 2 * 1.414 )
          WRITE (lu,'(9x,''A beam flux = '',1pg12.6,''K'',
     &                3x,''B beam flux = '',1pg12.6,''K''/
     &                9x,''average flux= '',1pg12.6,'' +- _'')')
     &                FLUXA, FLUXB, FLUX
*
      ELSE
*         single beam flux fit
          J = 0
          DO I = MIDBEAM - 5, MIDBEAM + 5
              J = J + 1
              CALL GFLUX (ARRAY,I,FLUXESTA(J))
          END DO
          FLUX = FLUXESTA(1)
          DO J = 2, 11
              FLUX = MAX (FLUX,FLUXESTA(J))
          END DO
          FLUXERROR = RMS / 2
          WRITE (lu,'(9x,''flux = '',1pg12.6,'' +- _'')') FLUX
      END IF
*
      WRITE (lu,'(1pg12.6,'' K'')') FLUXERROR
      IF (FLUX .LT. 5 * FLUXERROR) THEN
          WRITE (lu,'(9x,''flux less than 5 * flux error :''/
     &        3x,''A beam flux estimates (K) :''/
     &        3x,6(1pg10.3,2x) / 3x,5(1pg10.3,2x))')
     &        (FLUXESTA(I), I = 1,11)
          IF (DUALBEAM) THEN
              WRITE (lu,'(3x,''B beam flux estimates (K) :''/
     &            3x,6(1pg10.3,2x) / 3x,5(1pg10.3,2x))')
     &            (FLUXESTB(I), I = 1,11)
          END IF
      END IF
      RETURN
      END