
      subroutine cadef(alpha,beta,neta1,neta2,
     .                  sigma,etazp,etazm,pzetap,pzetam)
      implicit none
      real*8 alpha,beta,sigma,etazp,etazm,pzetap,pzetam
      integer neta1,neta2
      integer NPTS,NCOEFFS,NC
      real*8 ACC,ZERO,HALF,ONE,TWO,FOUR,AXISFAC,PI
c
c    NPTS: Number of sample points around each hole.
c NCOEFFS: Number of Cadez coord coefficients
c          C^+_n , C^-_n : n=0..NCOEFFS-1.
c     ACC: Accuracy for finding singular point location.
c AXISFAC: 1/2 if Cadez singularity is to be zone centered in the
c          eta direction and 0 if it is to be on an eta zone face.
c
      parameter(NPTS=1000,
     .          NCOEFFS=12,
     .          NC=2*NCOEFFS,
     .          ACC = 1.D-14,
     .          ZERO=0.D0,
     .          HALF=0.5D0,
     .          ONE=1.D0,
     .          TWO=2.D0,
     .          FOUR=4.D0,
     .          AXISFAC = ZERO,
     .          PI=3.14159265358979323846264338327950D0)
      real*8 x(2,2*NPTS),y(2*NPTS),sig(2*NPTS),coeffs(NC)
      real*8 w(NC),u(2*NPTS,NC),V(NC,NC)
      real*8 zetap,zetam
      real*8 etap,etam,dtheta,amo,tc,ts,change,chisq,scale
      real*8 etasing,xising,zsing,rhosing
      integer i
      external etavector
      common /svdvc/ zetap,zetam
c
c find the bispherical coordinate surfaces which satisfy the
c desired alpha and beta ratios
      etap = beta*beta
      etam = ONE/(alpha*alpha)
      etap = etap + (ONE + etam*etam)/etap -
     .             TWO*(ONE + ONE/etap)*etam
      etap = log(sqrt((etap + sqrt(etap*etap - FOUR))/TWO))
      etam = alpha*sinh(etap)
      etam = -log(etam + sqrt(etam*etam + ONE))
c
c set the scale of the coordinate system relative to the radius
c of the plus hole 
      sigma = sinh(etap)
c
c find the z coordinates of the centers of the holes
      zetap = sigma/tanh(etap)
      zetam = sigma/tanh(etam)
      pzetap = zetap
      pzetam = zetam
c
c make initial guesses for the eta0 coordinate values defining
c the minimal surfaces
      etazp = abs(ONE/sinh(etap))
      etazm = abs(ONE/sinh(etam))
c
c calculate npts coordinate values on the surface of each hole
      dtheta = PI/(NPTS-1)
      amo = ONE/alpha
c*ksr* ptile(i)
      do 10 i=1,NPTS
        tc = cos(i*dtheta)
        ts = sin(i*dtheta)
        x(1,2*i-1) = zetap + tc
        x(2,2*i-1) = ts
        x(1,2*i) = zetam + amo*tc
        x(2,2*i) = amo*ts
        y(2*i-1) = etazp
        y(2*i) = etazm
        sig(2*i-1) = ONE
        sig(2*i) = ONE
10    continue
c
c find the least squares fit by Singular value decomposition
      call svdfitc(x,y,sig,2*NPTS,coeffs,NC,u,v,w,
     .                        2*NPTS,NC,chisq,etavector)
c
c locate the singular point in the coordinate transformation
      call findchis(etasing,xising,zsing,rhosing,ACC,zetap,zetam,
     .             coeffs,NC)
c
      change = (etasing - etazp)/(etasing - etazm)
c
c if the ratio of the difference in the eta coordinate between the
c minimal surfaces and the singular point is not the given rational
c ratio, then recalculate the coefficients with new initial eta0 values
20    if((change.lt.((neta1 + AXISFAC)/(neta2 + AXISFAC) + ACC)).and.
     .   (change.gt.((neta1 + AXISFAC)/(neta2 + AXISFAC) - ACC)))
     .                                                       goto 40
        change = (neta2 + AXISFAC)*(etasing - etazp) -
     .                       (neta1 + AXISFAC)*(etasing - etazm)
        change = change/(neta1 + neta2 + TWO*AXISFAC)
        etazp = etazp + change
        etazm = etazm - change
c*ksr* ptile(i)
        do 30 i=1,NPTS
          y(2*i-1) = etazp
          y(2*i) = etazm
30      continue
        call svdfitc(x,y,sig,2*NPTS,coeffs,NC,u,v,w,
     .                          2*NPTS,NC,chisq,etavector)
        call findchis(etasing,xising,zsing,rhosing,ACC,zetap,zetam,
     .               coeffs,NC)
        change = (etasing - etazp)/(etasing - etazm)
      goto 20
40    continue
c
c renormalize the coefficients and coordinate values
      scale = ONE/(coeffs(1) + coeffs(2))
      do 50 i=1,NC
        coeffs(i) = coeffs(i)*scale
50    continue
      etazp = etazp*scale
      etazm = etazm*scale
      etasing = etasing*scale
      xising = xising*scale
c
c initialize the Cadez coordinate routines with ichicord
      call ichicord(NCOEFFS,coeffs,alpha,zetap,zetam,zsing,rhosing)
c
c print out Cadez coordinate information
      write(*,1001) chisq
1001  format('chisq = ',g15.8,
     .       '           /* Chi squared value for Cadez coord fit */')
      return
      end


      function pythagdumb(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.
c
      implicit none
      real*8 pythagdumb,a,b
      real*8 at,bt,ct,ZERO,ONE
      parameter(ZERO=0.D0,ONE=1.D0)
      at = abs(a)
      bt = abs(b)
      if(at.gt.bt) then
         ct = bt/at
         pythagdumb = at*sqrt(ONE + ct*ct)
      else if(bt.gt.ZERO) then
         ct = at/bt
         pythagdumb = bt*sqrt(ONE + ct*ct)
      else
         pythagdumb = ZERO
      end if
      return
      end


      subroutine etavector(z,rho,p,np)
      implicit none
      real*8 z,rho
      integer np
      real*8 p(np)
c      ***!!!*** WARNING***!!!***
c NCOEFFS here must be >= NCOEFFS in cadef
c      ***!!!*** WARNING***!!!***
      real*8 ZERO,ONE
      integer NCOEFFS
      parameter(NCOEFFS=12,ZERO=0.D0,ONE=1.D0)
      real*8 ciap(NCOEFFS),ciam(NCOEFFS),siap(NCOEFFS),siam(NCOEFFS)
      real*8 imp(NCOEFFS),imm(NCOEFFS),zmzp,zmzm,mzmzp,mzmzm,azmzp,azmzm
      real*8 pythagdumb,zetap,zetam
      integer cn,i
      common /svdvc/ zetap,zetam
c
      cn = np/2 - 1
      zmzp = z - zetap
      zmzm = z - zetam
      mzmzp = pythagdumb(zmzp,rho)
      mzmzm = pythagdumb(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
         ciap(i) = ciap(1)*ciap(i-1) - siap(1)*siap(i-1)
         siap(i) = siap(1)*ciap(i-1) + ciap(1)*siap(i-1)
         ciam(i) = ciam(1)*ciam(i-1) - siam(1)*siam(i-1)
         siam(i) = siam(1)*ciam(i-1) + ciam(1)*siam(i-1)
20    continue
      p(1) = mzmzp
      p(2) = mzmzm
      do 30 i=1,cn
        p(2*i+1) = imp(i)*ciap(i)
        p(2*i+2) = imm(i)*ciam(i)
30    continue
      return
      end


      subroutine findchis(eta,xi,z,rho,acc,zetap,zetam,p,np)
      implicit none
      real*8 eta,xi,z,rho,acc,zetap,zetam
      integer np
c      ***!!!*** WARNING***!!!***
c NCOEFFS here must be >= NCOEFFS in cadef
c      ***!!!*** WARNING***!!!***
      real*8 ZERO,QRTR,ONE
      integer NCOEFFS
      parameter (NCOEFFS=12,ZERO=0.D0,QRTR=0.25D0,ONE=1.D0)
      real*8 p(2,0:NCOEFFS-1)
      real*8 imp(NCOEFFS+2),imm(NCOEFFS+2)
      real*8 tr(NCOEFFS+2),ti(NCOEFFS+2)
      real*8 ciap(NCOEFFS+2),ciam(NCOEFFS+2)
      real*8 siap(NCOEFFS+2),siam(NCOEFFS+2)
      integer cn,iters,i
      real*8 pythagdumb,zmzp,zmzm,mzmzp,mzmzm,azmzp,azmzm
      real*8 temp,dchir,dchii,d2chir,d2chii,deltar,deltai

c initialize quess for z and rho
      z = ZERO
      rho = ZERO
c
      cn = np/2 - 1

c intitialize iteration counter.
      iters = 0
c Newton loop
1     continue
      iters = iters + 1
c Allow only 100 Newton iterations, them HALT
      if(iters .gt. 100) then
         write(*,*) 'sngzeta root not converged in 100 iterations'
         stop
      end if

      zmzp = z-zetap
      zmzm = z-zetam
      mzmzp = pythagdumb(zmzp,rho)
      mzmzm = pythagdumb(zmzm,rho)
      azmzp = atan2(rho,zmzp)
      azmzm = atan2(rho,zmzm)
      imp(1) = ONE/mzmzp
      imm(1) = ONE/mzmzm
      do 10 i=2,cn+2
         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+2
         ciap(i) = ciap(1)*ciap(i-1) - siap(1)*siap(i-1)
         siap(i) = siap(1)*ciap(i-1) + ciap(1)*siap(i-1)
         ciam(i) = ciam(1)*ciam(i-1) - siam(1)*siam(i-1)
         siam(i) = siam(1)*ciam(i-1) + ciam(1)*siam(i-1)
20    continue
c compute individual terms in first derivative of chi series
      do 30 i=1,cn
         tr(i) = p(1,i)*imp(i+1)*ciap(i+1) + p(2,i)*imm(i+1)*ciam(i+1)
         ti(i) = p(1,i)*imp(i+1)*siap(i+1) + p(2,i)*imm(i+1)*siam(i+1)
30    continue
      dchir = ZERO
      dchii = ZERO
c sum first derivative of chi series in reverse order
c to minimize roundoff error
      do 40 i=cn,1,-1
         dchir = dchir - i*tr(i)
         dchii = dchii + i*ti(i)
40    continue
      dchir = dchir + p(1,0)*imp(1)*ciap(1) + p(2,0)*imm(1)*ciam(1)
      dchii = dchii - p(1,0)*imp(1)*siap(1) - p(2,0)*imm(1)*siam(1)
c compute individual terms in second derivative of chi series
      do 50 i=1,cn
         tr(i) = p(1,i)*imp(i+2)*ciap(i+2) + p(2,i)*imm(i+2)*ciam(i+2)
         ti(i) = p(1,i)*imp(i+2)*siap(i+2) + p(2,i)*imm(i+2)*siam(i+2)
50    continue
      d2chir = ZERO
      d2chii = ZERO
c sum second derivative of chi series in reverse order
c to minimize roundoff error
      do 60 i=cn,1,-1
         d2chir = d2chir + i*(i+1)*tr(i)
         d2chii = d2chii - i*(i+1)*ti(i)
60    continue
      d2chir = d2chir - p(1,0)*imp(2)*ciap(2) - p(2,0)*imm(2)*ciam(2)
      d2chii = d2chii + p(1,0)*imp(2)*siap(2) + p(2,0)*imm(2)*siam(2)
c update Newton quess
      temp = d2chir*d2chir + d2chii*d2chii
      if(temp.ne.ZERO) then
         temp = ONE/temp
         deltar = (dchir*d2chir + dchii*d2chii)*temp
         deltai = (dchii*d2chir - dchir*d2chii)*temp
      else
         deltar = acc**(QRTR)
         deltai = -acc**(QRTR)
      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(pythagdumb(deltar,deltai) .gt. acc) goto 1
c compute eta and xi singular coordinates values
c
      zmzp = z - zetap
      zmzm = z - zetam
      mzmzp = pythagdumb(zmzp,rho)
      mzmzm = pythagdumb(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 70 i=2,cn
        imp(i) = imp(i-1)*imp(1)
        imm(i) = imm(i-1)*imm(1)
70    continue
      ciap(1) = cos(azmzp)
      ciam(1) = cos(azmzm)
      siap(1) = sin(azmzp)
      siam(1) = sin(azmzm)
      do 80 i=2,cn
         ciap(i) = ciap(1)*ciap(i-1) - siap(1)*siap(i-1)
         siap(i) = siap(1)*ciap(i-1) + ciap(1)*siap(i-1)
         ciam(i) = ciam(1)*ciam(i-1) - siam(1)*siam(i-1)
         siam(i) = siam(1)*ciam(i-1) + ciam(1)*siam(i-1)
80    continue
c compute individual terms in Cadez coordinate series
      do 90 i=1,cn
         tr(i) = p(1,i)*imp(i)*ciap(i) + p(2,i)*imm(i)*ciam(i)
         ti(i) = p(1,i)*imp(i)*siap(i) + p(2,i)*imm(i)*siam(i)
90    continue
      eta = ZERO
      xi = ZERO
c sum the series in reverse order to minimize roundoff error
      do 100 i=cn,1,-1
         eta = eta + tr(i)
         xi = xi -  ti(i)
100   continue
      eta = eta + p(1,0)*mzmzp + p(2,0)*mzmzm
      xi = xi + p(1,0)*azmzp + p(2,0)*azmzm
      return
      end


c
c Numerical Recipes SVDFIT routine modified to allow for 2-D coordinates.
c
      SUBROUTINE svdfitc(x,y,sig,ndata,a,ma,u,v,w,mp,np,chisq,funcs)
      implicit none
      INTEGER ma,mp,ndata,np,NMAX,MMAX
      REAL*8 chisq,a(ma),sig(ndata),u(mp,np),v(np,np),w(np)
      REAL*8 x(2,ndata),y(ndata),TOL,ZERO,ONE
      EXTERNAL funcs
      PARAMETER (NMAX=4000,MMAX=50,TOL=1.D-5,ZERO=0.D0,ONE=1.D0)
CU    USES svbksb,svdcmp
      INTEGER i,j
      REAL*8 sum,thresh,tmp,wmax,afunc(MMAX),b(NMAX)
      do 12 i=1,ndata
        call funcs(x(1,i),x(2,i),afunc,ma)
        tmp=ONE/sig(i)
        do 11 j=1,ma
          u(i,j)=afunc(j)*tmp
11      continue
        b(i)=y(i)*tmp
12    continue
      call svdcmp(u,ndata,ma,mp,np,w,v)
      wmax=ZERO
      do 13 j=1,ma
        if(w(j).gt.wmax)wmax=w(j)
13    continue
      thresh=TOL*wmax
      do 14 j=1,ma
        if(w(j).lt.thresh)w(j)=ZERO
14    continue
      call svbksb(u,w,v,ndata,ma,mp,np,b,a)
      chisq=ZERO
      do 16 i=1,ndata
        call funcs(x(1,i),x(2,i),afunc,ma)
        sum=0.
        do 15 j=1,ma
          sum=sum+a(j)*afunc(j)
15      continue
        chisq=chisq+((y(i)-sum)/sig(i))**2
16    continue
      return
      END


      SUBROUTINE svdcmp(a,m,n,mp,np,w,v)
      implicit none
      INTEGER m,mp,n,np,NMAX
      REAL*8 a(mp,np),v(np,np),w(np)
      real*8 ZERO,ONE,TWO
      PARAMETER (NMAX=100,ZERO=0.D0,ONE=1.D0,TWO=2.D0)
      INTEGER i,its,j,jj,k,l,nm
      REAL*8 anorm,c,f,g,h,s,scale,x,y,z,rv1(NMAX),pythagdumb
      g=ZERO
      scale=ZERO
      anorm=ZERO
      do 25 i=1,n
        l=i+1
        rv1(i)=scale*g
        g=ZERO
        s=ZERO
        scale=ZERO
        if(i.le.m)then
          do 11 k=i,m
            scale=scale+abs(a(k,i))
11        continue
          if(scale.ne.ZERO)then
            do 12 k=i,m
              a(k,i)=a(k,i)/scale
              s=s+a(k,i)*a(k,i)
12          continue
            f=a(i,i)
            g=-sign(sqrt(s),f)
            h=f*g-s
            a(i,i)=f-g
            do 15 j=l,n
              s=ZERO
              do 13 k=i,m
                s=s+a(k,i)*a(k,j)
13            continue
              f=s/h
              do 14 k=i,m
                a(k,j)=a(k,j)+f*a(k,i)
14            continue
15          continue
            do 16 k=i,m
              a(k,i)=scale*a(k,i)
16          continue
          endif
        endif
        w(i)=scale *g
        g=ZERO
        s=ZERO
        scale=ZERO
        if((i.le.m).and.(i.ne.n))then
          do 17 k=l,n
            scale=scale+abs(a(i,k))
17        continue
          if(scale.ne.ZERO)then
            do 18 k=l,n
              a(i,k)=a(i,k)/scale
              s=s+a(i,k)*a(i,k)
18          continue
            f=a(i,l)
            g=-sign(sqrt(s),f)
            h=f*g-s
            a(i,l)=f-g
            do 19 k=l,n
              rv1(k)=a(i,k)/h
19          continue
            do 23 j=l,m
              s=ZERO
              do 21 k=l,n
                s=s+a(j,k)*a(i,k)
21            continue
              do 22 k=l,n
                a(j,k)=a(j,k)+s*rv1(k)
22            continue
23          continue
            do 24 k=l,n
              a(i,k)=scale*a(i,k)
24          continue
          endif
        endif
        anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
25    continue
      do 32 i=n,1,-1
        if(i.lt.n)then
          if(g.ne.ZERO)then
            do 26 j=l,n
              v(j,i)=(a(i,j)/a(i,l))/g
26          continue
            do 29 j=l,n
              s=ZERO
              do 27 k=l,n
                s=s+a(i,k)*v(k,j)
27            continue
              do 28 k=l,n
                v(k,j)=v(k,j)+s*v(k,i)
28            continue
29          continue
          endif
          do 31 j=l,n
            v(i,j)=ZERO
            v(j,i)=ZERO
31        continue
        endif
        v(i,i)=ONE
        g=rv1(i)
        l=i
32    continue
      do 39 i=min(m,n),1,-1
        l=i+1
        g=w(i)
        do 33 j=l,n
          a(i,j)=ZERO
33      continue
        if(g.ne.ZERO)then
          g=ONE/g
          do 36 j=l,n
            s=ZERO
            do 34 k=l,m
              s=s+a(k,i)*a(k,j)
34          continue
            f=(s/a(i,i))*g
            do 35 k=i,m
              a(k,j)=a(k,j)+f*a(k,i)
35          continue
36        continue
          do 37 j=i,m
            a(j,i)=a(j,i)*g
37        continue
        else
          do 38 j= i,m
            a(j,i)=ZERO
38        continue
        endif
        a(i,i)=a(i,i)+ONE
39    continue
      do 49 k=n,1,-1
        do 48 its=1,30
          do 41 l=k,1,-1
            nm=l-1
            if((abs(rv1(l))+anorm).eq.anorm)  goto 2
            if((abs(w(nm))+anorm).eq.anorm)  goto 1
41        continue
1         c=ZERO
          s=ONE
          do 43 i=l,k
            f=s*rv1(i)
            rv1(i)=c*rv1(i)
            if((abs(f)+anorm).eq.anorm) goto 2
            g=w(i)
            h=pythagdumb(f,g)
            w(i)=h
            h=ONE/h
            c= (g*h)
            s=-(f*h)
            do 42 j=1,m
              y=a(j,nm)
              z=a(j,i)
              a(j,nm)=(y*c)+(z*s)
              a(j,i)=-(y*s)+(z*c)
42          continue
43        continue
2         z=w(k)
          if(l.eq.k)then
            if(z.lt.ZERO)then
              w(k)=-z
              do 44 j=1,n
                v(j,k)=-v(j,k)
44            continue
            endif
            goto 3
          endif
          if(its.eq.30) pause 'no convergence in svdcmp'
          x=w(l)
          nm=k-1
          y=w(nm)
          g=rv1(nm)
          h=rv1(k)
          f=((y-z)*(y+z)+(g-h)*(g+h))/(TWO*h*y)
          g=pythagdumb(f,ONE)
          f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
          c=ONE
          s=ONE
          do 47 j=l,nm
            i=j+1
            g=rv1(i)
            y=w(i)
            h=s*g
            g=c*g
            z=pythagdumb(f,h)
            rv1(j)=z
            c=f/z
            s=h/z
            f= (x*c)+(g*s)
            g=-(x*s)+(g*c)
            h=y*s
            y=y*c
            do 45 jj=1,n
              x=v(jj,j)
              z=v(jj,i)
              v(jj,j)= (x*c)+(z*s)
              v(jj,i)=-(x*s)+(z*c)
45          continue
            z=pythagdumb(f,h)
            w(j)=z
            if(z.ne.ZERO)then
              z=ONE/z
              c=f*z
              s=h*z
            endif
            f= (c*g)+(s*y)
            x=-(s*g)+(c*y)
            do 46 jj=1,m
              y=a(jj,j)
              z=a(jj,i)
              a(jj,j)= (y*c)+(z*s)
              a(jj,i)=-(y*s)+(z*c)
46          continue
47        continue
          rv1(l)=ZERO
          rv1(k)=f
          w(k)=x
48      continue
3       continue
49    continue
      return
      END


      SUBROUTINE svbksb(u,w,v,m,n,mp,np,b,x)
      implicit none
      INTEGER m,mp,n,np,NMAX
      REAL*8 b(mp),u(mp,np),v(np,np),w(np),x(np)
      real*8 ZERO
      PARAMETER (NMAX=100,ZERO=0.D0)
      INTEGER i,j,jj
      REAL*8 s,tmp(NMAX)
      do 12 j=1,n
        s=ZERO
        if(w(j).ne.ZERO)then
          do 11 i=1,m
            s=s+u(i,j)*b(i)
11        continue
          s=s/w(j)
        endif
        tmp(j)=s
12    continue
      do 14 j=1,n
        s=ZERO
        do 13 jj=1,n
          s=s+v(j,jj)*tmp(jj)
13      continue
        x(j)=s
14    continue
      return
      END


      SUBROUTINE svdvar(v,ma,np,w,cvm,ncvm)
      implicit none
      INTEGER ma,ncvm,np,MMAX
      REAL*8 cvm(ncvm,ncvm),v(np,np),w(np)
      REAL*8 ZERO,ONE
      PARAMETER (MMAX=20,ZERO=0.D0,ONE=1.D0)
      INTEGER i,j,k
      REAL*8 sum,wti(MMAX)
      do 11 i=1,ma
        wti(i)=ZERO
        if(w(i).ne.ZERO) wti(i)=ONE/(w(i)*w(i))
11    continue
      do 14 i=1,ma
        do 13 j=1,i
          sum=ZERO
          do 12 k=1,ma
            sum=sum+v(i,k)*v(j,k)*wti(k)
12        continue
          cvm(i,j)=sum
          cvm(j,i)=sum
13      continue
14    continue
      return
      END
