C     This fortran program writes the analytic solution of Verification
C     problem 3.7.1 to a file.

      program tubedrag
      implicit double precision (a-h,o-z)
      complex*16 a,b,c,wavk2,cj, beta,ctemp, v0
      complex*16 el,ex,p
      character*50 filename
C
      one=1.d0
      zero=0.d0
      pi=3.1415926535D0
      cj=cmplx(zero,one)
      cone=cmplx(one,zero)
c
      iin=5
      iout=6
      iout2=8

 10   write(iout,*) ' Please enter freq in cycles/time '
      read (iin,*) freq
      omega= freq * 2 * pi

C This inward velocity corresponds to a prescribed outward acceleration of
C 100+0i

      ainward= -100d0
      v0=  (ainward)/( cj *omega)
      bulk1=.1424d6
      bulk2=.1424d6
      rho1=1.21d0
      rho2=1.21d0
      r2=1d4
      totallen=4.d0
      zlen2=1
      zlen1=totallen-zlen2
      
C     Undamped wavespeeds
      c1=sqrt(bulk1/rho1)
      c2=sqrt(bulk2/rho2)

C     Damped wavenumbers
      wavk1=omega/c1
      wavk2=omega/c2* sqrt( cmplx (one, -r2/(omega*rho2)))

C Calculate parameter Beta
      call tanhab( cj *wavk2*zlen2, ctemp)
      beta= cj*wavk2*rho1*c1/ cmplx(r2, omega*rho2)*ctemp

      write(iout,*) ' Filename on which you want the results written:'
      read (iin, 56 ) filename
 56   format( A50 )
      write(iout,*) ' filename= ',filename
      open( file=filename, unit=iout2, status='unknown')
      el=exp( cj* wavk1 *zlen1)

      write( iout2,*) ' ** x,pmag,ppor'

C     loop for writing out analytic solution:
      do i=1,301
        x=(i-1)*.01d0
        ex=exp( cj* wavk1 *x)

        p=rho1*c1*v0 * ((cone+beta) * el/ex + (cone-beta)/el*ex)/
     $       ( (cone+beta)*el - (cone-beta)/el)
        pmag=  abs(p)
        preal= real(p)
        pimag= -real(cj*p)
        ppor=  180/pi*atan2(pimag, preal)
        write( iout2,99) x,pmag,ppor
 99     format(g10.5,x,g10.5,x,g10.5)
      enddo
      end

      subroutine tanhab(c,zztan)
C  Put complex hyperbolic tangent of c into zztan
      complex*16 zztan,c,x,y
C
      x=exp(c)
      y=exp(-c)
      zztan=(x-y)/(x+y)
      return
      end