      PROGRAM BART

C     Takes covariance matrix and replica shade file and computes a Bartlett
C     surface

      INTEGER REPFIL, COVFIL, AMSFIL
      PARAMETER ( REPFIL = 20, COVFIL = 21, AMSFIL = 25,
     &   MAXNSD = 4094, MAXNRD = 4094, MAXNR = 4094,
     &   MXPHONES = 50, TOL = 0.0001 )

      INTEGER IRDREP( MAXNRD )
      REAL    R( MAXNR ), SD( MAXNSD ), RD( MAXNRD ), SURF( 1 ),
     &        RDREP( MAXNRD )
      COMPLEX P( MXPHONES, MAXNR ), PB( MAXNR ),
     &        COV( MXPHONES, MXPHONES )
      CHARACTER TITLE*80, TITLEREP*80, PLTTYP*10

      POWERMAX = 0.0
      ISDMAX   = 0
      IRMAX    = 0

C     *** Read covariance matrix ***

      OPEN ( FILE = 'COVFIL', UNIT = COVFIL, STATUS = 'OLD' )

      READ( COVFIL, * ) TITLE
      READ( COVFIL, * ) FREQ
      READ( COVFIL, * ) NRD

      IF ( NRD .GT. MXPHONES ) STOP "NRD > MXPHONES"

      DO IRD = 1, NRD
         READ( COVFIL, * ) RD( IRD )
      END DO

      DO I = 1, NRD
         DO J = 1, NRD
            READ( COVFIL, * ) I1, J1, COV( I, J )
         END DO
      END DO

C     *** Form power field (SURF = 1.0) ***

C     *** Read replica file header ***

      CALL RDHEAD( REPFIL, 'REPFIL', MAXNR, MAXNSD, MAXNRD,
     &   FREQREP, R, NR, DELTAR, PLTTYP, XS, YS, THETA,
     &   SD, NSD, RDREP, NRDREP, TITLEREP )

      IF ( FREQ .NE. FREQREP ) THEN
         WRITE( *, * ) 'Mismatch in frequencies'
         STOP
      ENDIF

C     --- Find pointers to receiver depths in replica file

      DO IRD = 1, NRD
         IRDREP( IRD ) = NEIGH( RD( IRD ), RDREP, NRDREP )
         IF ( ABS( RDREP( IRDREP( IRD ) ) - RD( IRD ) ) .GT. TOL )
     &      STOP "No matching depth in replica file"
      END DO

C     *** Form header for power field (SURF = 1.0)

      SURF( 1 )  = 1.0
      NSURF = 1

      CALL HEADER( TITLE, SURF, NSURF, SD, NSD,
     &   R, NR, FREQ, ATTEN, PLTTYP, XS, YS, THETA )

      IRECOUT = 6

C     *** Loop over trial source depths ***

      DO ISD = 1, NSD

         DO IRD = 1, NRD
C           --- read in replica vector
            IRECIN = 6 + NRD * ( ISD - 1 ) + IRDREP( IRD )
            READ ( REPFIL, REC = IRECIN ) ( P( IRD, IR ), IR = 1, NR )
         END DO

C        *** Loop over source ranges ***

         DO IR = 1, NR

C           *** form Pb = e* R e taking advantage of symmetry ***

            POWER  = 0.0
            SUMREP = 0.0
            DO I = 1, NRD

               IF ( I .LT. NRD ) THEN
               DO J = I+1, NRD
                  POWER = POWER + 2.0 * REAL( 
     &            CONJG( P( I, IR ) ) * COV( I, J ) * P( J, IR ) )
               END DO
               END IF

C              --- diagonal terms
               POWER = POWER +
     &            CONJG( P( I, IR ) ) * COV( I, I ) * P( I, IR )

C              --- sum for norm of replica vector
               SUMREP = SUMREP + CONJG( P( I, IR ) ) * P( I, IR )
            END DO

C           --- Normalize
            IF( SUMREP .NE. 0.0 ) THEN
C              --- take sqrt to go from power back to amplitude
C              --- take reciprocal to get proper sign in plot

               IF ( TITLE(1:12) .NE. 'Unnormalized' ) 
     &            POWER = POWER / SUMREP

C              --- make note of peak location
               IF ( POWER .GT. POWERMAX ) THEN
                  POWERMAX = POWER
                  ISDMAX   = ISD
                  IRMAX    = IR
               ENDIF

               PB( IR ) = 1.0 / SQRT( POWER )
            ELSE
               PB( IR ) = 0.0
            ENDIF
         END DO        ! next source range

C        --- Write out Bartlett power vector
         IRECOUT = IRECOUT + 1
         WRITE( AMSFIL, REC = IRECOUT ) ( PB( IR ), IR = 1, NR )

      END DO        ! next source depth

C      WRITE( *, * ) 'Peak power   = ', POWERMAX
C      WRITE( *, * ) 'Source depth = ', SD( ISDMAX )
C      WRITE( *, * ) 'Source range = ', R(   IRMAX )

      WRITE( *, * ) R( IRMAX ), SD( ISDMAX )

      STOP
      END
