      function pythag(a,b)
c
c This function returns the sqrt(a*a + b*b) taking into account roundoff
c effects.  This should be used instead of abs(cmplx(a,b)) whenever roundoff
c effects are important.  This is crutial to the routine chi2cyl.
c
      implicit none
      real*8 pythag,a,b
      real*8 ZERO,ONE
      parameter(ZERO=0.D0,ONE=1.D0)
      real*8 at,bt,ct
      at = abs(a)
      bt = abs(b)
      if(at.gt.bt) then
        ct = bt/at
        pythag = at*sqrt(ONE + ct*ct)
      else if(bt.gt.ZERO) then
        ct = at/bt
        pythag = bt*sqrt(ONE + ct*ct)
      else
        pythag = ZERO
      end if
      return
      end

      subroutine ichicord(n,c,a,zp,zm,z,rho)
c
c initialize the Cadez coordinate common block /chiinf/.
c
c      n : Number of Cadez coefficients
c       -> n-1 is stored in /chiinf/:cn
c      c : Matric of Cadez coefficients c(2,0:n-1)
c          c(1,i) are the C+ coefficients
c       -> c(1,i) is stored in /chiinf/:cp
c          c(2,i) are the C- coefficients
c       -> c(2,i) is stored in /chiinf/:cm
c      a : Cadez coordinate alpha ratio
c          Alpha is the ratio of the radius of
c          the "+" hole to the "-" hole
c          alpha = a+/a-
c       -> a is stored in /chiinf/:alpha
c     zp : scaled z coordinate value of the center of
c          of the "+" hole.  zp = z+/a+
c       -> zp is stored in /chiinf/:zetap
c     zm : scaled z coordinate value of the center of
c          of the "-" hole.  zp = z-/a+
c       -> zm is stored in /chiinf/:zetam
c      z : scaled z coordinate value of the singular point
c          of the coordinate system. z = zsing/a+
c       -> z is stored in /chiinf/:zsing
c    rho : scaled rho coordinate value of the singular
c          point of the coordinate system.  rho = rhosing/a+
c       -> rho is stored in /chiinf/:rhosing
c
      implicit none
      real*8 c(2,0:n-1),a,zp,zm,z,rho
      integer n
      real*8 ZERO,ONE
      parameter(ZERO=0.D0,ONE=1.D0)
      integer i
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      real*8 xitemp,zs,rhos
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
      cn = n-1
      if(cn .gt. 30) then
        write(*,*) 'Too many Cadez coefficints must be <= 30'
        stop
      end if
      do 10 i=0,cn
        cp(i) = c(1,i)
        cm(i) = c(2,i)
10    continue
      alpha = a
      zetap = zp
      zetam = zm
      zsing = z
      rhosing = rho
      call cyl2chi(z,rho,etasing,xising)
      zs = zetap + ONE
      rhos = ZERO
      call cyl2chi(zs,rhos,etazp,xitemp)
      zs = zetam - ONE/alpha
      rhos = ZERO
      call cyl2chi(zs,rhos,etazm,xitemp)
      return
      end

      subroutine cyl2chi(z,rho,eta,xi)
c
c convert cylidrical coordinate values (z,rho) into Cadez coordinate
c values (eta,xi)
c
      implicit none
      real*8 z,rho,eta,xi
      real*8 ZERO,ONE
      parameter(ZERO=0.D0,ONE=1.D0)
      real*8 imp(30),imm(30)
      real*8 teta(30),txi(30)
      real*8 ciap(30),ciam(30),siap(30),siam(30)
      real*8 pythag,zmzp,zmzm,mzmzp,mzmzm,azmzp,azmzm
      integer i
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
      zmzp = z - zetap
      zmzm = z - zetam
      mzmzp = pythag(zmzp,rho)
      mzmzm = pythag(zmzm,rho)
      azmzp = atan2(rho,zmzp)
      azmzm = atan2(rho,zmzm)
      imp(1) = ONE/mzmzp
      imm(1) = ONE/mzmzm
      mzmzp = log(mzmzp)
      mzmzm = log(mzmzm)
      do 10 i=2,cn
        imp(i) = imp(i-1)*imp(1)
        imm(i) = imm(i-1)*imm(1)
10    continue
      ciap(1) = cos(azmzp)
      ciam(1) = cos(azmzm)
      siap(1) = sin(azmzp)
      siam(1) = sin(azmzm)
      do 20 i=2,cn
c compute cos(i*azmzp)
        ciap(i) = ciap(1)*ciap(i-1) - siap(1)*siap(i-1)
c compute sin(i*azmzp)
        siap(i) = siap(1)*ciap(i-1) + ciap(1)*siap(i-1)
c compute cos(i*azmzm)
        ciam(i) = ciam(1)*ciam(i-1) - siam(1)*siam(i-1)
c compute sin(i*azmzm)
        siam(i) = siam(1)*ciam(i-1) + ciam(1)*siam(i-1)
20    continue
c compute individual terms in Cadez coordinate series
      do 30 i=1,cn
        teta(i) = cp(i)*imp(i)*ciap(i) + cm(i)*imm(i)*ciam(i)
        txi(i) = cp(i)*imp(i)*siap(i) + cm(i)*imm(i)*siam(i)
30    continue
      eta = ZERO
      xi = ZERO
c sum the series in reverse order to minimize roundoff error
      do 40 i=cn,1,-1
        eta = eta + teta(i)
        xi = xi -  txi(i)
40    continue
      eta = eta + cp(0)*mzmzp + cm(0)*mzmzm
      xi = xi + cp(0)*azmzp + cm(0)*azmzm
      return
      end

      subroutine diffeta1(z,rho,detadrho,detadz)
c
c Compute the first partial derivatives of eta at the
c cylindrical coordinate locatin (z,rho)
c   detadrho : returns the first derivative of eta(rho,z)
c              with respect to rho 
c     detadz : returns the first derivative of eta(rho,z)
c              with respect to z 
c
      implicit none
      real*8 z,rho,detadrho,detadz
      real*8 ZERO,ONE
      parameter(ZERO=0.D0,ONE=1.D0)
      real*8 imp(0:30),imm(0:30)
      real*8 tdhdz(30),tdhdr(30)
      real*8 ciap(30),ciam(30),siap(30),siam(30)
      real*8 pythag,zmzp,zmzm,mzmzp,mzmzm,azmzp,azmzm
      integer i
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
      zmzp = z-zetap
      zmzm = z-zetam
      mzmzp = ONE/pythag(zmzp,rho)
      mzmzm = ONE/pythag(zmzm,rho)
      azmzp = atan2(rho,zmzp)
      azmzm = atan2(rho,zmzm)
      imp(0) = mzmzp*mzmzp
      imm(0) = mzmzm*mzmzm
      do 10 i=1,cn
        imp(i) = imp(i-1)*mzmzp
        imm(i) = imm(i-1)*mzmzm
10    continue
      ciap(1) = cos(azmzp)
      ciam(1) = cos(azmzm)
      siap(1) = sin(azmzp)
      siam(1) = sin(azmzm)
      do 20 i=2,cn
c compute cos(i*azmzp)
        ciap(i) = ciap(1)*ciap(i-1) - siap(1)*siap(i-1)
c compute sin(i*azmzp)
        siap(i) = siap(1)*ciap(i-1) + ciap(1)*siap(i-1)
c compute cos(i*azmzm)
        ciam(i) = ciam(1)*ciam(i-1) - siam(1)*siam(i-1)
c compute sin(i*azmzm)
        siam(i) = siam(1)*ciam(i-1) + ciam(1)*siam(i-1)
20    continue
c compute the individual terms in the series
      do 30 i=1,cn
        tdhdr(i) = i*(cp(i)*imp(i)*(rho*ciap(i) + zmzp*siap(i)) +
     .      cm(i)*imm(i)*(rho*ciam(i) + zmzm*siam(i)))
        tdhdz(i) = i*(cp(i)*imp(i)*(rho*siap(i) - zmzp*ciap(i)) +
     .      cm(i)*imm(i)*(rho*siam(i) - zmzm*ciam(i)))
30    continue
      detadrho = ZERO
      detadz = ZERO
c sum the series in reverse order to minimize roundoff error
      do 40 i=cn,1,-1
        detadrho = detadrho - tdhdr(i)
        detadz = detadz + tdhdz(i)
40    continue
      detadrho = detadrho + cp(0)*rho*imp(0) + cm(0)*rho*imm(0)
      detadz = detadz + cp(0)*zmzp*imp(0) + cm(0)*zmzm*imm(0)
      return
      end

      subroutine diffeta2(z,rho,d2etadrr,d2etadrz)
c
c Compute the second partial derivatives of eta at the
c cylindrical coordinate locatin (z,rho)
c      d2etadrr : returns the second derivative of eta(rho,z)
c              with respect to rho 
c      d2etadrz : returns the second derivative of eta(rho,z)
c              with respect to rho and z
c
      implicit none
      real*8 z,rho,d2etadrr,d2etadrz
      real*8 ZERO,ONE,TWO
      parameter(ZERO=0.D0,ONE=1.D0,TWO=2.D0)
      real*8 imp(0:30),imm(0:30)
      real*8 td2hdr2(30),td2hdrz(30)
      real*8 ciap(30),ciam(30),siap(30),siam(30)
      real*8 pythag,zmzp,zmzm,mzmzp,mzmzm,azmzp,azmzm
      real*8 sqdiffp,sqdiffm
      integer i
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
      zmzp = z-zetap
      zmzm = z-zetam
      mzmzp = ONE/pythag(zmzp,rho)
      mzmzm = ONE/pythag(zmzm,rho)
      azmzp = atan2(rho,zmzp)
      azmzm = atan2(rho,zmzm)
      imp(0) = mzmzp*mzmzp
      imp(0) = imp(0)*imp(0)
      imm(0) = mzmzm*mzmzm
      imm(0) = imm(0)*imm(0)
      do 10 i=1,cn
        imp(i) = imp(i-1)*mzmzp
        imm(i) = imm(i-1)*mzmzm
10    continue
      ciap(1) = cos(azmzp)
      ciam(1) = cos(azmzm)
      siap(1) = sin(azmzp)
      siam(1) = sin(azmzm)
      do 20 i=2,cn
c compute cos(i*azmzp)
        ciap(i) = ciap(1)*ciap(i-1) - siap(1)*siap(i-1)
c compute sin(i*azmzp)
        siap(i) = siap(1)*ciap(i-1) + ciap(1)*siap(i-1)
c compute cos(i*azmzm)
        ciam(i) = ciam(1)*ciam(i-1) - siam(1)*siam(i-1)
c compute sin(i*azmzm)
        siam(i) = siam(1)*ciam(i-1) + ciam(1)*siam(i-1)
20    continue
	
      sqdiffp = zmzp*zmzp - rho*rho
      sqdiffm = zmzm*zmzm - rho*rho
      zmzp = TWO*zmzp*rho
      zmzm = TWO*zmzm*rho
c compute the individual terms in the series
      do 30 i=1,cn
        td2hdr2(i) = i*(i+1)*(cp(i)*imp(i)*(zmzp*siap(i) -
     .      sqdiffp*ciap(i)) +
     .      cm(i)*imm(i)*(zmzm*siam(i) - sqdiffm*ciam(i)))
        td2hdrz(i) = i*(i+1)*(cp(i)*imp(i)*(zmzp*ciap(i) +
     .      sqdiffp*siap(i)) +
     .      cm(i)*imm(i)*(zmzm*ciam(i) + sqdiffm*siam(i)))
30    continue
      d2etadrr = ZERO
      d2etadrz = ZERO
c sum the series in reverse order to minimize roundoff error
      do 40 i=cn,1,-1
        d2etadrr = d2etadrr + td2hdr2(i)
        d2etadrz = d2etadrz + td2hdrz(i)
40    continue
      d2etadrr = d2etadrr + cp(0)*sqdiffp*imp(0) + cm(0)*sqdiffm*imm(0)
      d2etadrz = d2etadrz - cp(0)*zmzp*imp(0) - cm(0)*zmzm*imm(0)
      return
      end

      subroutine diffeta3(z,rho,d3etadrrr,d3etadrrz)
c
c Compute the third partial derivatives of eta at the
c cylindrical coordinate locatin (z,rho)
c      d3etadrrr : returns the third derivative of eta(rho,z)
c              with respect to rho 
c      d3etadrrz : returns the third derivative of eta(rho,z)
c              with respect to rho, rho and z
c
      implicit none
      real*8 z,rho,d3etadrrr,d3etadrrz
      real*8 ZERO,ONE,TWO,THREE
      parameter(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)
      real*8 imp(0:30),imm(0:30)
      real*8 td3hdr3(30),td3hdrrz(30)
      real*8 ciap(30),ciam(30),siap(30),siam(30)
      real*8 pythag,zmzp,zmzm,mzmzp,mzmzm,azmzp,azmzm
      real*8 sqdiffp,sqdiffm,rr3zp,rr3zm,zz3rp,zz3rm
      integer i
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
      zmzp = z-zetap
      zmzm = z-zetam
      mzmzp = ONE/pythag(zmzp,rho)
      mzmzm = ONE/pythag(zmzm,rho)
      azmzp = atan2(rho,zmzp)
      azmzm = atan2(rho,zmzm)
      imp(0) = mzmzp*mzmzp*mzmzp
      imp(0) = imp(0)*imp(0)
      imm(0) = mzmzm*mzmzm*mzmzm
      imm(0) = imm(0)*imm(0)
      do 10 i=1,cn
        imp(i) = imp(i-1)*mzmzp
        imm(i) = imm(i-1)*mzmzm
10    continue
      ciap(1) = cos(azmzp)
      ciam(1) = cos(azmzm)
      siap(1) = sin(azmzp)
      siam(1) = sin(azmzm)
      do 20 i=2,cn
c compute cos(i*azmzp)
        ciap(i) = ciap(1)*ciap(i-1) - siap(1)*siap(i-1)
c compute sin(i*azmzp)
        siap(i) = siap(1)*ciap(i-1) + ciap(1)*siap(i-1)
c compute cos(i*azmzm)
        ciam(i) = ciam(1)*ciam(i-1) - siam(1)*siam(i-1)
c compute sin(i*azmzm)
        siam(i) = siam(1)*ciam(i-1) + ciam(1)*siam(i-1)
20    continue
	
      rr3zp = rho*(rho*rho - THREE*zmzp*zmzp)
      rr3zm = rho*(rho*rho - THREE*zmzm*zmzm)
      zz3rp = zmzp*(zmzp*zmzp - THREE*rho*rho)
      zz3rm = zmzm*(zmzm*zmzm - THREE*rho*rho)
c compute the individual terms in the series
      do 30 i=1,cn
        td3hdr3(i) = i*(i+1)*(i+2)*(cp(i)*imp(i)*(zz3rp*siap(i) -
     .      rr3zp*ciap(i)) +
     .      cm(i)*imm(i)*(zz3rm*siam(i) - rr3zm*ciam(i)))
        td3hdrrz(i) = i*(i+1)*(i+2)*(cp(i)*imp(i)*(zz3rp*ciap(i) +
     .      rr3zp*siap(i)) +
     .      cm(i)*imm(i)*(zz3rm*ciam(i) + rr3zm*siam(i)))
30    continue
      d3etadrrr = ZERO
      d3etadrrz = ZERO
c sum the series in reverse order to minimize roundoff error
      do 40 i=cn,1,-1
        d3etadrrr = d3etadrrr + td3hdr3(i)
        d3etadrrz = d3etadrrz + td3hdrrz(i)
40    continue
      d3etadrrr = d3etadrrr + cp(0)*TWO*rr3zp*imp(0) +
     .                                     cm(0)*TWO*rr3zm*imm(0)
      d3etadrrz = d3etadrrz - cp(0)*TWO*zz3rp*imp(0) -
     .                                     cm(0)*TWO*zz3rm*imm(0)
      return
      end

      function coefnorm()
c
c Return the Cadez normalization coefficient.
c This value is usually 1 meaning the coordinate have
c already been normalized.  There is no need to normalize
c the coordinates except that it makes the coordinate 
c range for xi 0 <= xi <= PI
c
      implicit none
      real*8 coefnorm
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
      coefnorm = cp(0) + cm(0)
      return
      end

      function chi2cyl(eta,xi,z,rho,acc,quad)
c
c Converte the Cadez coordinates (eta,xi) into cylindrical coordinates
c (z,rho).  This routine uses Newton's method for complex functions.
c     acc : the accuracy to which the root finding algorithm will
c           determine the cylindrical coordinate values.
c    quad : the quadrant in which the point (eta,xi) is located
c           Used to initialize values for z and rho.
c              0 - if xi <= xising
c              1 - if xi >= xising
c
      implicit none
      real*8 eta,xi,z,rho,acc
      integer chi2cyl,quad
      real*8 Pi,ZERO,QUARTER,ALMOST,ONE,TWO
      parameter (Pi=3.141592653589793238462643D0,ZERO=0.D0,
     .           QUARTER=0.25D0,ALMOST=0.9D0,ONE=1.D0,TWO=2.D0)
      real*8 imp(31),imm(31),tr(31),ti(31)
      real*8 ciap(31),ciam(31),siap(31),siam(31)
      real*8 coefnorm,pythag,r,theta,chir,chii,dchir,dchii,deltar,deltai
      real*8 temp,smrad,etarad
      real*8 zmzp,zmzm,mzmzp,mzmzm,azmzp,azmzm
      integer i,iters
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn

c
c initialize chi2cyl to valid solution
      chi2cyl = 0

c initialize quess for z and rho
c Return immediately if on the singular point
      if((eta.eq.etasing).and.(xi.eq.xising)) then
        z = zsing
        rho = rhosing
        return
      else if(eta.ge.etasing) then
c exterior to double valued region.
c cp(0) + cp(1) is usually, but not necessarily normalized to 1
        r = exp(eta/coefnorm())
        r = max(r,zetap + ONE,-zetam + ONE/alpha)
        theta = xi/coefnorm()
        z = zsing + r*cos(theta)
        rho = r*sin(theta)
      else if(quad .eq. 0) then
c "+" quadrant
        theta = Pi*xi/xising
        if(eta.ge.etazp) then
          z = zetap + cos(theta)
          rho = abs(sin(theta))
        else if(eta.ge.ZERO) then
          etarad = eta/etazp
          z = zetap + etarad*cos(theta)
          rho = etarad*abs(sin(theta))
        else
          smrad = acc**(QUARTER)
          z = zetap + smrad*cos(theta)
          rho = smrad*abs(sin(theta))
        end if
      else
c "-" quadrant
        theta = Pi*(xi - xising)/(Pi*coefnorm() - xising)
        if(eta.ge.etazm) then
          z = zetam + cos(theta)/alpha
          rho = abs(sin(theta)/alpha)
        else if(eta.ge.ZERO) then
          etarad = eta/(etazm*alpha)
          z = zetam + etarad*cos(theta)
          rho = etarad*abs(sin(theta))
        else
          smrad = acc**(QUARTER)
          z = zetam + smrad*cos(theta)
          rho = smrad*abs(sin(theta))
        end if
      end if
c intitialize iteration counter.
      iters = 0
c Newton loop
1     continue
      zmzp = z-zetap
      zmzm = z-zetam
      mzmzp = pythag(zmzp,rho)
      mzmzm = pythag(zmzm,rho)
      if((mzmzp.lt.acc).or.(mzmzm.lt.acc)) then
        rho = rho + sqrt(TWO*acc)
        goto 1
      end if
      iters = iters + 1
c     if(iters .gt. 90) then
c       write(*,*) iters,'(Rho,Z) : (',rho,',',z,')'
c       write(*,*) '(D_R,D_Z) : (',deltai,',',deltar,')'
c     end if
c Allow only 100 Newton iterations, them HALT
      if(iters .gt. 100) then
        write(*,*) 'chi2cyl root not converged in 100 iterations'
        write(*,*) '(Eta,Xi) : (',eta,',',xi,')'
        chi2cyl = 1
        return
      end if
		
      azmzp = atan2(rho,zmzp)
      azmzm = atan2(rho,zmzm)
      imp(1) = ONE/mzmzp
      imm(1) = ONE/mzmzm
      mzmzp = log(mzmzp)
      mzmzm = log(mzmzm)
      do 10 i=2,cn+1
        imp(i) = imp(i-1)*imp(1)
        imm(i) = imm(i-1)*imm(1)
10    continue
      ciap(1) = cos(azmzp)
      ciam(1) = cos(azmzm)
      siap(1) = sin(azmzp)
      siam(1) = sin(azmzm)
      do 20 i=2,cn+1
c compute cos(i*azmzp)
        ciap(i) = ciap(1)*ciap(i-1) - siap(1)*siap(i-1)
c compute sin(i*azmzp)
        siap(i) = siap(1)*ciap(i-1) + ciap(1)*siap(i-1)
c compute cos(i*azmzm)
        ciam(i) = ciam(1)*ciam(i-1) - siam(1)*siam(i-1)
c compute sin(i*azmzm)
        siam(i) = siam(1)*ciam(i-1) + ciam(1)*siam(i-1)
20    continue
c compute individual terms in chi series
      do 30 i=1,cn
        tr(i) = cp(i)*imp(i)*ciap(i) + cm(i)*imm(i)*ciam(i)
        ti(i) = cp(i)*imp(i)*siap(i) + cm(i)*imm(i)*siam(i)
30    continue
      chir = ZERO
      chii = ZERO
c sum chi series in reverse order to minimize roundoff error
      do 40 i=cn,1,-1
        chir = chir + tr(i)
        chii = chii - ti(i)
40    continue
      chir = chir + cp(0)*mzmzp + cm(0)*mzmzm - eta
      chii = chii + cp(0)*azmzp + cm(0)*azmzm - xi
c compute individual terms in first derivative of chi series
      do 50 i=1,cn
        tr(i) = cp(i)*imp(i+1)*ciap(i+1) + cm(i)*imm(i+1)*ciam(i+1)
        ti(i) = cp(i)*imp(i+1)*siap(i+1) + cm(i)*imm(i+1)*siam(i+1)
50    continue
      dchir = ZERO
      dchii = ZERO
c sum first derivative of chi series in reverse order
c to minimize roundoff error
      do 60 i=cn,1,-1
        dchir = dchir - i*tr(i)
        dchii = dchii + i*ti(i)
60    continue
      dchir = dchir + cp(0)*imp(1)*ciap(1) + cm(0)*imm(1)*ciam(1)
      dchii = dchii - cp(0)*imp(1)*siap(1) - cm(0)*imm(1)*siam(1)
c update Newton quess
      temp = dchir*dchir + dchii*dchii
      if(temp.ne.ZERO) then
        temp = ONE/temp
        deltar = (chir*dchir + chii*dchii)*temp
        deltai = (chii*dchir - chir*dchii)*temp
      else if(quad .eq. 0) then
c if we wind up on the singular point, move a little
        deltar = -acc**(QUARTER)
        deltai = -acc**(QUARTER)
      else
        deltar = acc**(QUARTER)
        deltai = -acc**(QUARTER)
      end if
      z = z - deltar
      rho = rho - deltai
c if newton moves into wrong area, correct
      if(rho .lt. ZERO) rho = ZERO
c iterate until magnitude of change is < acc
      if(ABS(z).ge.ONE) deltar = deltar/z
      if(ABS(rho).ge.ONE) deltai = deltai/rho
      if(pythag(deltar,deltai) .gt. acc) goto 1
      return
      end

      subroutine cylsing(z,rho)
c
c return in (z,rho) the cylindrical coordinates of the singular point
c
      implicit none
      real*8 z,rho
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
      z = zsing
      rho = rhosing
      return
      end

      subroutine chising(eta,xi)
c
c return in (eta,xi) the Cadez coordinates of the singular point
c
      implicit none
      real*8 eta,xi
      real*8 cp(0:30),cm(0:30),alpha,zetap,zetam,zsing,rhosing
      real*8 etasing,xising,etazp,etazm
      integer cn
      common /chiinf/cp,cm,alpha,zetap,zetam,
     .   zsing,rhosing,etasing,xising,etazp,etazm,cn
      eta = etasing
      xi = xising
      return
      end

