C***********************************************************************
      SUBROUTINE LOGANA( X, Y ,NY )

C     Calculus of the logarithm by unwrapping the phase

      INTEGER*4 NY, NPI, K
      COMPLEX*16 X( NY ), Y( NY )
      REAL*8 YR( NY ), YI( NY ), PI
      PARAMETER( PI = 3.14159265358979323846264D0 )


      YR( 1 ) = DLOG( CDABS( X( 1 ) ) )
      YI( 1 ) = DATAN2( DIMAG( X( 1 ) ), DREAL( X( 1 ) ) )

      Y( 1 ) = DCMPLX( YR( 1 ), YI( 1 ) ) 

      NPI = 0

      DO K = 2, NY

         YR( K ) = DLOG( CDABS( X( K ) ) )
         YI( K ) = ATAN2( DIMAG( X( K ) ), DREAL( X( K ) ) )

         IF ((YI(K) + NPI * 2.0D0 * PI  - YI(K-1)) .GT. PI) THEN
            NPI = NPI - 1
         ENDIF 

         IF ((YI(K) + NPI * 2.0D0 * PI - YI(K-1)) .LT. - PI) THEN
            NPI = NPI + 1
         ENDIF 

         YI( K ) = YI( K ) + NPI * 2.0D0 * PI

         Y( K ) = DCMPLX( YR( K ), YI( K ) ) 

      END DO
c           write(*,*) 'NPI :',NPI

      RETURN
      END
C*********************************************************************
      SUBROUTINE FNBZEROS( NBZ, FY, NRECT)

C     Calculus of the number of zeros inside the contour

      INTEGER K
      COMPLEX*16 FY( NRECT ), NBZR, Y( NRECT ), CI
      INTEGER*4 NBZ, NRECT
      REAL*8 PI
      PARAMETER( PI = 3.14159265358979323846264D0 )
      PARAMETER( CI = ( 0.0D0 , 1.0D0 ) )

      CALL LOGANA( FY, Y, NRECT )

      NBZR = ( Y( NRECT  ) - Y( 1 ) )/ 2.0D0 / PI / CI
c      write(*,*) 'Nreel:', NBZR
      NBZ = INT( REAL( NBZR ) )
 
      RETURN 
      END
C*********************************************************************
      SUBROUTINE FNBZEROSC( NBZ, FY, NCER )

C     Calculus of the number of zeros inside the contour

      INTEGER K
      COMPLEX*16 FY( NCER ), FYY(50000), NBZP, Y(50000), CI
      INTEGER*4 NBZ, NCER, NCER1
      REAL*8 PI
      PARAMETER( PI = 3.14159265358979323846264D0 )
      PARAMETER( CI = ( 0.0D0 , 1.0D0 ) )

      FYY( 1 ) = FY( NCER )

      DO K = 1, NCER
        FYY( K+1) = FY( K )
      END DO

      NCER1 = NCER + 1

      CALL LOGANA( FYY, Y, NCER1 )

      NBZP = ( Y( NCER1  ) - Y( 1 ) )/ 2.0D0 / PI / CI

      NBZ = INT( REAL( NBZP ) )
 
      RETURN 
      END
C**********************************************************************
      SUBROUTINE RECT( T, REP, DR, CENT, U, V, NCER )

C     Design of the rectangle used as contour
C      U is the horizontal size   V the vertical size

      COMPLEX*16 T( 20000 ), CENT, CI
      REAL*8 RAY, A1, A2, U, V, DR( 2 ), DR1, DR2, DRM
      INTEGER*4  NCER, NCER1, NCER2, REP( 5 ), K
      PARAMETER( CI = ( 0.0D0 , 1.0D0 ) )

      NCER1 = NCER

      A1   = DREAL( CENT ) - U / 2.0D0
      A2   = DIMAG( CENT ) + V / 2.0D0
      
c      write(*,*)'A1 :',A1,A1+U
c      write(*,*)'A2 :',A2,A2-V

c      DR1   = V / DBLE( NCER1 )
c      DR2   = U / DBLE( NCER1 )

c      DRM = MAX( DR1, DR2 )

      DRM = 2.0 * ( U + V ) / DBLE( NCER )
      NCER1 = NINT( V / DRM )
      DR1 = V / DBLE( NCER1 )

      NCER2 = NINT( U / DRM )
      DR2 = U / DBLE( NCER2 )

c      write(*,*) 'CENT  ', CENT
c      write(*,*) 'U  V : ', U, V
c      write(*,*) 'DR1  DR2 : ', DR1, DR2

      REP( 1 ) = 1

      DO K = 1, NCER1 + 1
        T( K ) = A1 + CI * ( A2 - ( K - 1 ) * DR1 )
      END DO
       
      REP( 2 ) = NCER1 + 1

      DO K = 1, NCER2
        T( K + REP( 2 ) ) = A1 + K * DR2 + CI * ( A2 - V )
      END DO

      REP( 3 ) = NCER1 + NCER2 + 1

      DO K = 1, NCER1
        T( K + REP( 3 ) ) = A1 + U + CI * ( A2 - V + K * DR1 )
      END DO

      REP( 4 ) = NCER1 * 2 + NCER2 + 1

      DO K = 1, NCER2
         T( K + REP( 4) ) = A1 + U - K * DR2 + CI * A2
      END DO

      REP( 5 ) = 2 * ( NCER1 + NCER2 ) + 1

      DR( 1 ) = DR1

      DR( 2 ) = DR2

      RETURN
      END
      
C**********************************************************************
      SUBROUTINE FSOLANA( SOL, T, FY, RAY, CENT, NBZ, NCER )

C     Solution in the complex plane

      COMPLEX*16 FY( NCER ), SOL(1:50), A(1:50), AI(1:50), CENT
      COMPLEX*16 FY1(50000), Y(50000), Y1, Y2, SY, S(1:50), SM, CI
      REAL*8 RAY, DT, PI, T( NCER )
      INTEGER*4 NBZ, NCER
      INTEGER I, J, K
      LOGICAL POLISH
      PARAMETER( PI = 3.14159265358979323846264D0 )
      PARAMETER( CI = ( 0.0D0 , 1.0D0 ) )

      DT = 2.0D0 * PI / ( T( 2 ) - T( 1 ) )

      DO I = 1, NCER
         FY1( I ) = EXP( - CI * NBZ * T( I ) ) * FY( I )
      END DO

      CALL LOGANA( FY1, Y, NCER )

      DO J = 1, NBZ
         SY = 0.0D0
         DO K = 1, NCER/2 
            Y1 = Y( 2 * K - 1 ) * EXP( CI * DBLE( J ) * T( 2 * K - 1 ) )
            Y2 = Y( 2 * K     ) * EXP( CI * DBLE( J ) * T( 2 * K     ) )
            SY = SY + 4.0D0 * Y1 + 2.0D0 * Y2
         END DO
         SY = SY / 3.0D0
         S( J ) = - DBLE( J ) / DT * SY  * RAY ** J
      END DO

C     Calculus of the polynomial coefficients

      AI( 1 ) = 1.0D0
      AI( 2 ) = - S( 1 )
      
      DO K = 2, NBZ
         SM = 0.0D0
         DO  J = 1, K - 1
             SM = SM + S( J ) * AI( K - J + 1 )
         END DO
         AI( K + 1 ) = - ( SM + S( K ) ) / DBLE( K )
      END DO

      DO K = 1, NBZ + 1
         A( K ) = AI( NBZ + 2 - K )
      END DO

      POLISH = .FALSE.

      CALL ZROOTS( A, NBZ, SOL, POLISH )

      DO K = 1, NBZ
         SOL( K ) = SOL( K ) + CENT
c         write(*,*) 'Sol C :', SOL(K)
      END DO

      RETURN
      END
C**********************************************************************
      SUBROUTINE FSOLREC( SOL, T, FY, REP, CENT, NBZ, NRECT, DR )

C     Solution in the complex plane

      COMPLEX*16 FY( NRECT ), SOL(1:30), A(1:30), AI(1:30), CENT
      COMPLEX*16 T1( NRECT ), Y1, Y2, S(1:50), SM, CI, Y( NRECT )
      COMPLEX*16 YZ( NRECT ), SY( 4 ), T( NRECT ), SYY
      REAL*8 RAY, PI, DR( 2 )
      INTEGER*4 NBZ, NRECT, REP( 5 ), I, J, K
      LOGICAL POLISH
      PARAMETER( PI = 3.14159265358979323846264D0 )
      PARAMETER( CI = ( 0.0D0 , 1.0D0 ) )

      DO K = 1, NRECT
         T1( K ) = T( K ) - CENT
         YZ( K ) =  T1( K ) **( -NBZ ) * FY( K )
      END DO

      CALL LOGANA( YZ, Y, NRECT )

c      Extended rule

      DO J = 1, NBZ
         DO I = 1, 4
           I1 = REP( I )
           I2 = REP( I + 1 ) 
           SYY = 3.0D0/8.0D0 * Y( I1 ) + 7.0D0/6.0D0 * Y( I1 + 1 )
     &             + 23.0D0/24.0D0 * Y( I1 + 2 )
           DO K = I1 + 3, I2 - 3
              SYY  = SYY + Y( K )
           END DO 
           SYY = SYY + 3.0D0/8.0D0 * Y( I2 ) + 7.0D0/6.0D0 * Y( I2 - 1 )
     &              + 23.0D0/24.0D0 * Y( I2 - 2 )
           SY( I ) = SYY
         END DO
         S( J ) = CI * ( SY( 3 ) - SY( 1 ) ) * DR( 1 )
     &                + ( SY( 2 ) - SY( 4 ) ) * DR( 2 )
         S( J ) =  - DBLE( J ) / 2.0D0 / PI / CI * S( J )
         
         IF ( J .NE. NBZ ) THEN
            DO K = 1, NRECT
              Y( K ) = Y( K ) * T1( K )
            END DO
          ENDIF         
      END DO

C     Calculus of the polynomial coefficients

      AI( 1 ) = 1.0D0
      AI( 2 ) = - S( 1 )

c      write(*,*) 'S : ', S(1:NBZ)

      DO K = 2, NBZ
         SM = 0.0D0
         DO  J = 1, K - 1
             SM = SM + S( J ) * AI( K - J + 1 )
         END DO
         AI( K + 1 ) = - ( SM + S( K ) ) / DBLE( K )
      END DO

      DO K = 1, NBZ + 1
         A( K ) = AI( NBZ + 2 - K )
      END DO

      POLISH = .TRUE.

      CALL ZROOTS( A, NBZ, SOL, POLISH )

      DO K = 1, NBZ
         SOL( K ) = SOL( K ) + CENT
c         write(*,*) 'Sol :', SOL(K)
      END DO

      RETURN
      END
C**********************************************************************
      SUBROUTINE PKSORT(n,arr)
      INTEGER*4 n
      REAL*8 arr(n)
      INTEGER*4 i,j
      REAL*8 a
      do 12 j=2,n
        a=arr(j)
        do 11 i=j-1,1,-1
          if(arr(i).le.a)goto 10
          arr(i+1)=arr(i)
11      continue
        i=0
10      arr(i+1)=a
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software #ksV).
C**********************************************************************
      SUBROUTINE ZPKSORT(n,arr)
      INTEGER*4 n
      COMPLEX*16 arr(n)
      INTEGER*4 i,j
      COMPLEX*16 a
      do 12 j=2,n
        a=arr(j)
        do 11 i=j-1,1,-1
          if(dreal(arr(i)).ge.dreal(a))goto 10
          arr(i+1)=arr(i)
11      continue
        i=0
10      arr(i+1)=a
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software #ksV).
C**********************************************************************C
      SUBROUTINE  GETMINFO( CWAT, CMINI, HEAU, HENS, CIMX ,
     &   CSMIN, CMAXI )

C     Get informations on the structure of the medium

      INCLUDE 'comc1.f'

      LOGICAL ELFLAG
      REAL*8 CWAT, CMINI, HEAU, HENS, AT, OMEGA, CSMIN, CIMX, CMAXI
      COMPLEX*16 CP( MAXN ), CS( MAXN ), CP2, CS2

      COMMON /SPACE/ CP, CS, DUMMY( 4 * MAXN )
      COMMON /CPREV/ ALPHAR, BETAR, RHOR, ALPHAI, BETAI
      
      ELFLAG   = .FALSE.
      CSMIN    = 1.0E6
      CMIN      = 1.0E6
      CMAXI    = 0.0D0
      CIMX      = 0.0D0
      NFACT    = 0
      LOC( 1 ) = 0
      
      OMEGA = 2* PI * FREQ

C     *** Loop over media ***

      DO MED = 1, NMEDIA
         IF ( MED .NE. 1 ) LOC( MED ) = LOC( MED-1 ) + N( MED-1 ) + 1
         N1 = N(   MED ) + 1
         I  = LOC( MED ) + 1

C        *** Load diagonals of the finite-difference equations ***

         IF ( CS( I ) .EQ. ( 0.0, 0.0 ) ) THEN

C           --- Case of an acoustic medium ---

            MATER( MED ) = 'ACOUSTIC'
            IF ( NFACT .EQ. 0 ) NFACT = MED
            NLACT = MED

            DO J = I, I + N( MED )
               CMIN   = MIN( CICR( CP( J ) ), CMIN )
               CMAXI  = MAX( CICR( CP( J ) ), CMAXI )
               CIMX   = MAX( ATCI( CP( J ), OMEGA ) , CIMX )
            END DO
C
            IF ( NFACT .EQ. MED ) THEN
              CWAT  = CMIN
              HEAU  = DEPTH( MED + 1 )
            ENDIF

         ELSE

C           --- Case of an elastic medium ---

            MATER( MED ) = 'ELASTIC'
            ELFLAG = .TRUE.
            TWOH   = 2.0 * H( MED )

            DO J = I, I + N( MED )
               CMIN    = MIN( CICR( CS( J ) ), CMIN   )
               CMAXI   = MAX( CICR( CS( J ) ), CMAXI )
               CMAXI   = MAX( CICR( CP( J ) ), CMAXI )
               CSMIN   = MIN( CICR( CS( J ) ), CSMIN )
               CIMX    = MAX( ATCI( CP( J ), OMEGA ), CIMX )
               CIMX    = MAX( ATCI( CS( J ), OMEGA ), CIMX )
            END DO

         ENDIF
C        ------ Next MED:
      END DO

C     *** Bottom properties ***

      IF ( BOTOPT(1:1) .EQ. 'A' ) THEN
         IF ( CSB .NE. ( 0.0, 0.0 ) ) THEN
C           ------ Elastic bottom:
            ELFLAG = .TRUE.
            CMIN   = MIN( CMIN, CICR( CSB ) )
            CMAXI = MAX( CICR( CSB ), CMAXI )
            CMAXI = MAX( CICR( CPB ), CMAXI )
            CSMIN = MIN(CICR( CSB ), CSMIN )
            CIMX   = MAX( CIMX, ATCI( CSB, OMEGA ) )
            CIMX   = MAX( CIMX, ATCI( CPB, OMEGA ) )
         ELSE
C           ------ Acoustic bottom:
            CMIN   = MIN( CMIN, CICR( CPB ) )
            CMAXI = MAX( CICR( CPB ), CMAXI )
            CIMX   = MAX( CIMX, ATCI( CPB, OMEGA ) )
         ENDIF
      ENDIF

C     *** Top properties ***

      IF ( TOPOPT(2:2) .EQ. 'A' ) THEN
         IF ( CST .NE. ( 0.0, 0.0 ) ) THEN
C           ------ Elastic top:
            ELFLAG = .TRUE.
            CMIN = MIN( CMIN, CICR( CST ) )
            CMAXI = MAX( CICR( CST ), CMAXI )
            CSMIN = MIN(CICR( CST ), CSMIN )
            CIMX = MAX( CIMX, ATCI( CST, OMEGA ) )
            CIMX = MAX( CIMX, ATCI( CPT, OMEGA ) )
         ELSE
C           ------ Acoustic top:
            CMIN = MIN( CMIN, CICR( CPT ) )
            CMAXI = MAX( CICR( CPT ), CMAXI )
            CIMX = MAX( CIMX, ATCI( CPT, OMEGA ) )
         ENDIF
      ENDIF

      CMINI = CMIN
      HENS  = DEPTH (NMEDIA + 1)
      IF ( CSMIN .EQ. 1.0E6 ) CSMIN = 0.0

      IF ( ELFLAG ) CMIN = 0.9 * CMIN

      RETURN
      END
C**************************************************************
      FUNCTION CICR( C )
      
      REAL*8 CICR
      COMPLEX*16 C
      
      IF ( C .EQ. ( 0.0D0,0.0D0 ) ) THEN
         CICR = 0.0D0
      ELSE
         CICR = C * DCONJG( C ) / DREAL( C )
      ENDIF
      
      RETURN
      END
C**************************************************************
      FUNCTION ATCI( C, OMEGA )
      
      REAL*8 ATCI, OMEGA
      COMPLEX*16 C
      
      IF ( C .EQ. ( 0.0,0.0 ) ) THEN
         ATCI = 0.0
      ELSE
         ATCI = OMEGA * DIMAG( C ) / C / DCONJG( C ) 
      ENDIF
      
      RETURN
      END
C**************************************************************      
      SUBROUTINE zroots(a,m,roots,polish)
      INTEGER m,MAXM
      DOUBLE PRECISION EPS
      COMPLEX*16 a(m+1),roots(m)
      LOGICAL polish
      PARAMETER (EPS=1.d-20,MAXM=101)
CU    USES laguer
      INTEGER i,j,jj,its
      COMPLEX*16 ad(MAXM),x,b,c
      do 11 j=1,m+1
        ad(j)=a(j)
11    continue
      do 13 j=m,1,-1
        x=dcmplx(0.d0,0.d0)
        call laguer(ad,j,x,its)
        if(abs(dimag(x)).le.2.d0*EPS**2*abs(dreal(x))) x=dcmplx(dreal(x)
     *,0.d0)
        roots(j)=x
        b=ad(j+1)
        do 12 jj=j,1,-1
          c=ad(jj)
          ad(jj)=b
          b=x*b+c
12      continue
13    continue
      if (polish) then
        do 14 j=1,m
          call laguer(a,m,roots(j),its)
14      continue
      endif
      do 16 j=2,m
        x=roots(j)
        do 15 i=j-1,1,-1
          if(dreal(roots(i)).le.dreal(x))goto 10
          roots(i+1)=roots(i)
15      continue
        i=0
10      roots(i+1)=x
16    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software #ksV).
C******************************************************
      SUBROUTINE laguer(a,m,x,its)
      INTEGER m,its,MAXIT,MR,MT
      DOUBLE PRECISION EPSS
      COMPLEX*16 a(m+1),x
      PARAMETER (EPSS=3.d-16,MR=8,MT=10,MAXIT=MT*MR)
      INTEGER iter,j
      DOUBLE PRECISION abx,abp,abm,err,frac(MR)
      COMPLEX*16 dx,x1,b,d,f,g,h,sq,gp,gm,g2
      SAVE frac
      DATA frac /.5d0,.25d0,.75d0,.13d0,.38d0,.62d0,.88d0,1.d0/
      do 12 iter=1,MAXIT
        its=iter
        b=a(m+1)
        err=abs(b)
        d=dcmplx(0.d0,0.d0)
        f=dcmplx(0.d0,0.d0)
        abx=abs(x)
        do 11 j=m,1,-1
          f=x*f+d
          d=x*d+b
          b=x*b+a(j)
          err=abs(b)+abx*err
11      continue
        err=EPSS*err
        if(abs(b).le.err) then
          return
        else
          g=d/b
          g2=g*g
          h=g2-2.d0*f/b
          sq=sqrt((m-1)*(m*h-g2))
          gp=g+sq
          gm=g-sq
          abp=abs(gp)
          abm=abs(gm)
          if(abp.lt.abm) gp=gm
          if (max(abp,abm).gt.0.d0) then
            dx=m/gp
          else
            dx=exp(dcmplx(log(1.d0+abx),dble(iter)))
          endif
        endif
        x1=x-dx
        if(x.eq.x1)return
        if (mod(iter,MT).ne.0) then
          x=x1
        else
          x=x-dx*frac(iter/MT)
        endif
12    continue
      pause 'too many iterations in laguer'
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software #ksV).
