      PROGRAM PLOTGRN

C     Plot the Green's function

      INTEGER PLPFIL, PRTFIL
      COMPLEX CI
      
      PARAMETER ( PLPFIL = 1, PRTFIL = 6, CI = ( 0.0, 1.0 ),
     &            MAXM = 1500, MAXNK = 10000 )

      REAL RK( MAXNK ), GMAG( MAXNK ), KMIN, KMAX
      COMPLEX CK( MAXM ), PHIS( MAXM ), PHIR( MAXM ),
     &        G( MAXNK ), WNO2
      CHARACTER TITLE*80

C     Read in various parameters

      OPEN( FILE = 'PLPFIL', UNIT = PLPFIL, STATUS = 'OLD' )
      READ( PLPFIL, * ) SRCD, RCVRD
      READ( PLPFIL, * ) NKPTS, STABIL
      READ( PLPFIL, * ) KMIN, KMAX

C     Get modal values

      CALL GETMOD( 1, ' ', MAXM, SRCD,  1, 'N',
     &                        CK, PHIS, M, FREQ, TITLE )
      CALL GETMOD( 1, ' ', MAXM, RCVRD, 1, 'N',
     &                        CK, PHIR, M, FREQ, TITLE )

C     Compute vector of k-values

      DELK = ( KMAX - KMIN ) / ( NKPTS - 1 )

      DO J = 1, NKPTS
         RK( J ) = KMIN + ( J - 1 ) * DELK
      END DO

C     *** Evaluate Green's function ***

      DO J = 1, NKPTS
         G( J ) = 0.0
         WNO2 = ( RK( J ) + CI * STABIL ) ** 2

         DO L = 1, M
            G( J ) = G( J ) + PHIS( L ) * PHIR( L ) /
     &                         ( WNO2 - CK( L ) **2 )
         END DO

         G( J )    = 2.0 * G( J )
         GMAG( J ) = ABS( G( J ) )
      END DO ! next k-point

C     *** Write data ***

      LAST = INDEX( TITLE, '    ' ) - 1
      WRITE( PRTFIL, 1000 ) TITLE(1:LAST)
 1000 FORMAT( '@ title "', A, '"' )

      WRITE( PRTFIL, 2000 ) NINT( FREQ ), NINT( SRCD ), NINT( RCVRD )
 2000 FORMAT(  '@ subtitle "',
     &         'F  =', I4, ' Hz, ',
     &         'Sd =', I4, ' m,  ',
     &         'Rd =', I4, ' m"' )
 
      DO IK = 1, NKPTS
         WRITE( PRTFIL, * ) RK( IK ), GMAG( IK )
      END DO

      STOP
      END
