      PROGRAM COVAR

C     Takes replica shade files and computes synthesized
C     data vectors with Gaussian random noise added
C     Averages data vector to from covariance matrix

C     Data vectors are produced for the phone positions in
C     the shade file (not those of the mode file).

      INTEGER PRTFIL, SHDFIL, DATFIL, COVFIL

      PARAMETER ( INPFIL = 5, PRTFIL = 6, SHDFIL = 20, DATFIL = 21,
     &   COVFIL = 22, MAXM = 1500,
     &   MAXNSD = 4094, MAXNRD = 4094, MAXNR = 4094,
     &   MXPHONES = 101, PI = 3.141592 )

      INTEGER   ISEED( 4 )
      REAL      R( MAXNR ), SD( MAXNSD ), RD( MAXNRD ),
     &          COVISO( MXPHONES, MXPHONES )
      COMPLEX   P( MAXNR ), D( MAXNRD ), DNOISY( MAXNRD ),
     &          PHIS( MAXM ), PHIR( MAXM, MXPHONES ), CK( MAXM ),
     &          NOISE( MAXNRD ), COV( MXPHONES, MXPHONES )
      CHARACTER TITLE*80, TITLE2*80, PLTTYP*10, NOISETYPE*1

      SAVE ISEED
      DATA ISEED /2301, 4562, 3219, 1823/

C     *** Read data file header ***

      CALL RDHEAD( SHDFIL, ' ', MAXNR, MAXNSD, MAXNRD,
     &   FREQ, R, NR, DELTAR, PLTTYP, XS, YS, THETA,
     &   SD, NSD, RD, NRD, TITLE )

      OPEN ( FILE = 'DATFIL', UNIT = DATFIL, STATUS = 'NEW' )
      OPEN ( FILE = 'COVFIL', UNIT = COVFIL, STATUS = 'NEW' )

C     Read user specification for source position and compute indices

      READ( INPFIL, * ) TITLE
      READ( INPFIL, * ) SDEPTH, SRANGE
      READ( INPFIL, * ) NOISETYPE, SNRDB
      READ( INPFIL, * ) NTIMES, ISEED( 4 )
C     NTIMES = 30000  ! **********************************
C     SNRDB  = 200.0   ! **********************************
      SRANGE = 1000.0 * SRANGE          ! Convert km to meters
      SNR    = 10.0 ** ( SNRDB / 10.0 )

C     *** Calculate source depth index ***

      ISDEP = NEIGH( SDEPTH, SD, NSD )

C     --- Check to see if within tolerance

      IF ( FREQ .EQ. 0.0 ) THEN
         TOL = 0.001
      ELSE
         TOL = 1500.0 / ( 6.0 * FREQ )  ! one wavelength
      END IF

      IF ( ABS( SD( ISDEP ) - SDEPTH ) .GT. TOL ) THEN
         WRITE( PRTFIL, * ) 'Sources exist at the following depths:'
         WRITE( PRTFIL, * ) ( SD( I ), I = 1, NSD )
         CALL ERROUT( PRTFIL, 'F', 'COVAR',
     &                'No source at specified depth' )
      END IF

C     *** Calculate source range index ***

      ISRAN = NEIGH( SRANGE, R, NR )

C     --- Check to see if within tolerance

      IF ( ABS( R( ISRAN ) - SRANGE ) .GT. TOL ) THEN
         WRITE( PRTFIL, * ) 'Sources exist at the following ranges:'
         WRITE( PRTFIL, * ) ( R( IR ), IR = 1, NR )
         CALL ERROUT( PRTFIL, 'F', 'COVAR',
     &                'No source at specified range' )
      END IF

C     --- extract data vector

      WRITE( PRTFIL, * ) 'Data vector being extracted for SDEP, SRAN',
     &   SD( ISDEP), R( ISRAN )

      IREC = ( ISDEP - 1 ) * NRD + 6
      DO IRD = 1, NRD
         IREC = IREC + 1
         READ ( SHDFIL, REC = IREC ) ( P( IR ), IR = 1, NR )
         D( IRD ) = P( ISRAN )
      END DO

C     --- Normalize data vector so that array power = 1

      SUM = 0.0
      DO IRD = 1, NRD
         SUM = SUM + CONJG( D( IRD ) ) * D( IRD )
      END DO

      DO IRD = 1, NRD
         D( IRD ) = D( IRD ) / SQRT( SUM )
      END DO

      PHNPOWER = 1.0 / NRD

C     --- Zero out covariance matrix

      DO I = 1, NRD
         NOISE( IRD ) = 0.0
         DO J = 1, NRD
            COV( I, J ) = 0.0
         END DO
      END DO

C     *** If making colored noise, need modal data ***

      IF ( NOISETYPE .EQ. 'C' ) THEN
         IPROF   = 1
         SD( 1 ) = 1.0       ! noise sources are at 1 m
         NSD     = 1
      
         CALL GETMOD( IPROF, ' ', MAXM, SD, NSD, 'N',
     &      CK, PHIS, M, FREQ, TITLE2 )
         CALL GETMOD( IPROF, ' ', MAXM, RD, NRD, 'N',
     &      CK, PHIR, M, FREQ, TITLE2 )

C        --- Construct PHI and POWER
         POWER = 0.0
         DO MODE = 1, M
            DO IRD = 1, NRD
               IF ( AIMAG( CK( MODE ) ) .GE. 0.0 ) THEN
                  CALL ERROUT( PRTFIL, 'F', 'COVAR',
     &       'Modes present with no loss: Noise field is infinite ' )
               END IF

               PHIR( MODE, IRD ) = PHIS( MODE ) * PHIR( MODE, IRD ) /
     &            SQRT( -REAL( CK( MODE ) ) * AIMAG( CK( MODE ) ) )
               POWER = POWER +
     &            PHIR( MODE, IRD ) * CONJG( PHIR( MODE, IRD ) )
            END DO
         END DO
      END IF

C     --- header for data vector file

      WRITE( DATFIL, * ) '''', TITLE(1:60), ''''
      WRITE( DATFIL, * ) FREQ
      WRITE( DATFIL, * ) NTIMES, NRD

      DO IRD = 1, NRD
         WRITE( DATFIL, * ) RD( IRD )
      END DO

C     *** Main loop: generate realizations of noise field ***

      AVGSIGNAL = 0.0
      AVGNOISE  = 0.0

      SIGMA = SQRT( PHNPOWER / SNR )   ! Noise level on a phone

      DO ITIME = 1, NTIMES

         IF ( NOISETYPE .EQ. 'W' ) THEN
            CALL WHITE( ISEED, SIGMA, NOISE, NRD )
         ELSE IF ( NOISETYPE .EQ. 'C' ) THEN
            CALL COLOR( ISEED, PHIR, M, NRD, SIGMA, POWER, NOISE )
         ELSE
            STOP 'Unknown noise type'
         END IF

         DO IRD = 1, NRD    ! loop over receiver depths

            DNOISY( IRD ) = D( IRD ) + NOISE( IRD )

            AVGSIGNAL = AVGSIGNAL + ABS(     D( IRD ) ) ** 2
            AVGNOISE  = AVGNOISE  + ABS( NOISE( IRD ) ) ** 2
         END DO             ! next receiver depth

C        --- write noise vector to DATFIL
C         WRITE( DATFIL, 1000 )
C     &      ( ITIME, IRD, DNOISY( IRD ) , IRD = 1, NRD )
 1000    FORMAT( I3, I3, ' (', F11.7, ',', F11.7, ' )' )

C        --- Compute cross-sensor correlation matrix R

         DO I = 1, NRD
            DO J = 1, NRD
               COV( I, J ) = COV( I, J ) +
     &              DNOISY( I ) * CONJG( DNOISY( J ) )
            END DO
          END DO

      END DO

      AVGSIGNAL = AVGSIGNAL / ( NTIMES * NRD )
      AVGNOISE  = AVGNOISE  / ( NTIMES * NRD )

      WRITE( PRTFIL, * ) 'AVGSIGNAL = ', AVGSIGNAL
      WRITE( PRTFIL, * ) 'AVGNOISE  = ', AVGNOISE

      WRITE( PRTFIL, * ) 'Average signal (dB) = ',
     &   10.0 * LOG10( AVGSIGNAL )
      WRITE( PRTFIL, * ) 'Average noise  (dB) = ',
     &   10.0 * LOG10( AVGNOISE  )
      WRITE( PRTFIL, * ) 'Average SNR    (dB) = ',
     &   10.0 * LOG10( AVGSIGNAL / AVGNOISE )

C     *** Analytic formula for covariance matrix in isovelocity ocean
      DO I = 1, NRD
         DO J = 1, NRD
            COVISO( I, J ) = 0.0
            DEPTH = 100.0
            DO MODE = 1, M
               GAMMA = ( MODE - .5 ) * PI / DEPTH
               COVISO( I, J ) = COVISO( I, J ) +
     &             SIN( GAMMA * SD( 1 ) ) ** 2 *
     &             SIN( GAMMA * RD( I ) ) * SIN( GAMMA * RD( J ) )
            END DO
         END DO
      END DO

C     *** Write normalized cross-sensor correlation matrix column-wise

      WRITE( COVFIL, * ) '''', TITLE(1:60), ''''
      WRITE( COVFIL, * ) FREQ
      WRITE( COVFIL, * ) NRD
      DO IRD = 1, NRD
         WRITE( COVFIL, * ) RD( IRD )
      END DO

      DO I = 1, NRD
         DO J = 1, NRD
            COV( I, J ) = NRD * COV( I,J ) /
     &              ( ( 1.0 + SIGMA**2 ) * NTIMES )
            WRITE( COVFIL, 1000 ) I, J, COV( I, J )
C     &            ABS( COVISO( I, J ) * COV( 1, 1 ) / COVISO( 1, 1 ) )
         END DO
      END DO

      STOP
      END

C----------------------------------------------------------------------C

      SUBROUTINE WHITE( ISEED, SIGMA, NOISE, NRD )

C     Generates Gaussian random white noise

      COMPLEX CI
      PARAMETER ( CI = (0.0, 1.0), PI = 3.1415926 )

      INTEGER   ISEED( 4 )
      COMPLEX   NOISE( * )

      DO IRD = 1, NRD
         X = RAN( ISEED ) + 0.000001
         Y = RAN( ISEED )

         NOISE( IRD ) = SIGMA *
     &      SQRT( -LOG( X ) ) * EXP( CI * 2.0 * PI * Y )

      END DO

C     --- Note power = NRD * SIGMA ** 2

      RETURN
      END

C----------------------------------------------------------------------C

      SUBROUTINE COLOR( ISEED, PHI, M, NRD, SIGMA, POWER, NOISE )

C     Generates colored noise vector

      COMPLEX CI
      PARAMETER ( MAXM = 1500, CI = (0.0, 1.0), PI = 3.1415926 )

      INTEGER   ISEED( 4 )
      COMPLEX PHI( MAXM, * ), NOISE( * ), COEF

C     *** Zero out noise vector ***

      DO IRD = 1, NRD
         NOISE( IRD ) = 0.0
      END DO

C     --- loop over modes
      DO MODE = 1, M

C        --- compute the random coefficient
         X = RAN( ISEED ) + 0.000001
         Y = RAN( ISEED )
         COEF = SQRT( -LOG( X ) ) * EXP( CI * 2.0 * PI * Y )

C        --- loop over receivers
         DO IRD = 1, NRD
            NOISE( IRD ) = NOISE( IRD ) + COEF * PHI( MODE, IRD)
         END DO
      END DO

C     *** Normalize so that power = NRD * SIGMA**2

      DO IRD = 1, NRD
         NOISE( IRD ) = SIGMA * NOISE( IRD ) * SQRT( NRD / POWER )
      END DO

      RETURN
      END
