subroutine tkund(np,x,px,y,py,z,g,dv,pz,al,rho,kx,eps0,coup)

      implicit none
      logical*4 coup
      integer*4 np,ndiv,i,j
      real*8 x(np),px(np),y(np),py(np),z(np),dv(np),g(np),pz(np)
      real*8 x1(np),px1(np),y1(np),py1(np),pxa(np),pya(np)
      real*8 dp(np),aln,rho1(np),thetan(np)
      real*8 al,rho,kx,eps0
      real*8 kap1(np),kap2(np),kap3(np)
      real*8 ai(np),bee(np),cee(np),dee(np)
      real*8 kei1(np),kei2(np),kei3(np)
      real*8 f1x(np), f2y(np),gxy(np),det(np)
      integer*4 method


c       write(*,*) 'rho=',rho,'kx=',kx
        ndiv = 50

        aln = al/ndiv
c        write(*,*) al, aln, ndiv

c     method = 1 non-linear
c              other linear
      method = 1

      if(method .eq. 1)then

c      write(*,*) 'tkund called',dp(1)
       do 220 i=1,np
          dp(i) = (g(i) + 1.d0)*(g(i) + 1.d0) -1.d0
          rho1(i)    = rho*(1+dp(i))
 220   continue

        do 200 i=1,np
         do 210 j=1, ndiv

            x1(i)  = x(i)
            px1(i) = px(i)
            y1(i)  = y(i)
            py1(i) = py(i)

            x1(i) = x(i) + px(i)*aln/2.d0
            y1(i) = y(i) + py(i)*aln/2.d0

            gxy(i) = 1.d0/rho1(i)/rho1(i)/kx
     $                      *cosh(kx*x1(i))*cosh(kx*y1(i))*aln
            f1x(i) = sinh(2.d0*kx*x1(i))/4.d0/rho1(i)/rho1(i)/kx*aln
            f2y(i) = sinh(2.d0*kx*y1(i))/4.d0/rho1(i)/rho1(i)/kx*aln
            det(i) = 1.d0+gxy(i)*gxy(i)/4.d0

c            gxy(i) = 1.d0/rho1(i)/rho1(i)/kx*aln
c            f1x(i) = x1(i)/2.d0/rho1(i)/rho1(i)*aln
c            f2y(i) = y1(i)/2.d0/rho1(i)/rho1(i)*aln
c            det(i) = 1.d0+gxy(i)*gxy(i)/4.d0

            pxa(i)=px(i)-f1x(i)+gxy(i)*py(i)/2.d0
            pya(i)=py(i)-f2y(i)-gxy(i)*px(i)/2.d0
            px1(i)=1.d0/det(i)*(pxa(i)+gxy(i)/2.d0*pya(i))
            py1(i)=1.d0/det(i)*(pya(i)-gxy(i)/2.d0*pxa(i))

            x1(i) = x1(i) + px1(i)*aln/2.d0
            y1(i) = y1(i) + py1(i)*aln/2.d0

            x(i)  = x1(i)
            px(i) = px1(i)
            y(i)  = y1(i)
            py(i) = py1(i)

 210              continue
 200                   continue




      else

        do 120 i=1,np
          dp(i) = (g(i) + 1.d0)*(g(i) + 1.d0) -1.d0
          rho1(i)    = rho*(1+dp(i))

          kei1(i) = sqrt(1.d0/(rho1(i)**4*kx**2)+2.d0/(rho1(i)**2))
          kei2(i) = 1.d0/(rho1(i)**2*kx)
          kei3(i) = -1.d0/(2.d0*rho1(i)**2)

          kap1(i) =  (kei1(i)-kei2(i))/2.d0
          kap2(i) = -(kei1(i)+kei2(i))/2.d0
          kap3(i) = kap1(i)-kap2(i)

          ai(i)  = (sin(kap1(i)*aln)-sin(kap2(i)*aln))/kap3(i)
          bee(i) = (cos(kap1(i)*aln)-cos(kap2(i)*aln))/kap3(i)
          cee(i) = (kap1(i)*cos(kap2(i)*aln)
     $                    -kap2(i)*cos(kap1(i)*aln))/kap3(i)
          dee(i) = (kap1(i)*sin(kap2(i)*aln)
     $                    -kap2(i)*sin(kap1(i)*aln))/kap3(i)
 120    continue

c      write(*,*) 'tkund called',dp(1)
        do 100 i=1,np
         do 110 j=1, ndiv

            x1(i)  = x(i)
            px1(i) = px(i)
            y1(i)  = y(i)
            py1(i) = py(i)

            x1(i) =cee(i)*x(i)+ai(i)*px(i)-dee(i)*y(i)+bee(i)*py(i)
            px1(i)=kei3(i)*ai(i)*x(i)+(-kei2(i)*bee(i)+cee(i))*px(i)
     $           +kei3(i)*bee(i)*y(i)+(kei2(i)*ai(i)-dee(i))*py(i)
            y1(i) =dee(i)*x(i)-bee(i)*px(i)+cee(i)*y(i)+ai(i)*py(i)
            py1(i)=-kei3(i)*bee(i)*x(i)+(-kei2(i)*ai(i)+dee(i))*px(i)
     $           +kei3(i)*ai(i)*y(i)+(-kei2(i)*bee(i)+cee(i))*py(i)

            x(i)  = x1(i)
            px(i) = px1(i)
            y(i)  = y1(i)
            py(i) = py1(i)

 110      continue  
100     continue
      endif
      return
      end
c