      subroutine initsym(alpha,zetap,zetam)
      implicit none
      real*8 alpha,zetap,zetam
      real*8 a2
      integer i
      real*8 ZERO,ONE
      integer NMAX
      parameter(ZERO = 0.D0,ONE = 1.D0,NMAX=40)
      real*8 zp(0:NMAX),zm(0:NMAX),wp(0:NMAX),wm(0:NMAX)
      common /symdat/ zp,zm,wp,wm
      zp(0) = zetap
      zm(0) = zetam
      wp(0) = ONE
      wm(0) = ONE/alpha
      a2 = alpha**2
      do i=1,NMAX
        if (mod(i,2) .eq. 1) then 
          zp(i) = ONE/(a2*(zp(i-1) - zm(0))) + zm(0)
          zm(i) = ONE/(zm(i-1) - zp(0)) + zp(0)
          wp(i) = wp(i-1)/(alpha*abs(zp(i-1) - zm(0)))
          wm(i) = wm(i-1)/abs(zm(i-1) - zp(0))
        else
          zp(i) = ONE/(zp(i-1) - zp(0)) + zp(0)
          zm(i) = ONE/(a2*(zm(i-1) - zm(0))) + zm(0)
          wp(i) = wp(i-1)/abs(zp(i-1) - zp(0))
          wm(i) = wm(i-1)/(alpha*abs(zm(i-1) - zm(0)))
        end if
      end do
      return
      end



      subroutine tspsi(z,rho,psi)
      implicit none
      real*8 z,rho,psi
      real*8 r2
      integer i
      real*8 ZERO,ONE
      integer NMAX,NPSI
      parameter(ZERO = 0.D0,ONE = 1.D0,NMAX=40,NPSI=20)
      real*8 zp(0:NMAX),zm(0:NMAX),wp(0:NMAX),wm(0:NMAX)
      common /symdat/ zp,zm,wp,wm
c set for returning psi - 1
      psi = ZERO
c set for returning psi
c     psi = ONE
      r2 = rho*rho
      do i=0,NPSI
        psi = psi + wp(i)/sqrt(r2 + (z - zp(i))**2)
        psi = psi + wm(i)/sqrt(r2 + (z - zm(i))**2)
      end do
      return
      end

      subroutine tsenergy(energy)
      implicit none
      real*8 energy
      integer i
      real*8 ZERO,TWO
      integer NMAX,NENG
      parameter(ZERO = 0.D0,TWO = 2.D0,NMAX=40,NENG=20)
      real*8 zp(0:NMAX),zm(0:NMAX),wp(0:NMAX),wm(0:NMAX)
      common /symdat/ zp,zm,wp,wm
      energy = ZERO
      do i=0,NENG
        energy = energy + wp(i) + wm(i)
      end do
      energy = TWO*energy
      return
      end

      subroutine tsdipole(dipole)
      implicit none
      real*8 dipole
      integer i
      real*8 ZERO,TWO
      integer NMAX,NDIP
      parameter(ZERO = 0.D0,TWO = 2.D0,NMAX=40,NDIP=40)
      real*8 zp(0:NMAX),zm(0:NMAX),wp(0:NMAX),wm(0:NMAX)
      common /symdat/ zp,zm,wp,wm
      dipole = ZERO
      do i=0,NDIP
        dipole = dipole + wp(i)*zp(i) + wm(i)*zm(i)
      end do
      dipole = TWO*dipole
      return
      end

      subroutine exactpsi(RHOzc,Zzc,PsiD,ext,fil)
      implicit none
      integer ext(4),fil(4)
      Real*8 RHOzc(ext(1):ext(3),ext(2):ext(4))
      Real*8 Zzc(ext(1):ext(3),ext(2):ext(4))
      Real*8 PsiD(ext(1):ext(3),ext(2):ext(4))
      integer i,j
      do j = fil(2),fil(4)
        do i = fil(1),fil(3)
          call tspsi(Zzc(i,j),RHOzc(i,j),PsiD(i,j))
        end do
      end do
      return
      end
