cc     *** scpack version 2 ***        lloyd nicholas trefethen
cc
cc     history:
cc       scpack 1 - double precision, ibm 370      oct. 1979
cc       scpack 2 - single precision, portable     july 1983
cc
cc     scpack is a fortran package for the numerical implementation
cc     of the schwarz-christoffel conformal map from a disk to
cc     a polygon, which may contain vertices at infinity.  see
cc     trefethen, "numerical computation of the schwarz-christoffel
cc     transformation", siam j. sci. stat. comp. 1 (1980), 82-102.
cc     the package may be freely copied and used, but reference to
cc     the above paper should be given in any publications arising
cc     from its use.
cc
cc     comment cards describe the principal features of the package.
cc     additional documentation is available in the "scpack user's
cc     guide".
cc
cc     the package requires library routines ns01a, gaussj, and ode,
cc     and these contain machine dependent constants.  see the scpack
cc     user's guide.


c*******************************************************************
c* qinit                                      primary subroutine  **
c*******************************************************************
      subroutine qinit(n,betam,nptsq,qwork)
c
c computes nodes and weights for gauss-jacobi quadrature
c
c calling sequence parameters: see comments in scsolv
c
c the work array qwork must be dimensioned at least nptsq * (2n+3).
c it is divided up into 2n+3 vectors of length nptsq: the first
c n+1 contain quadrature nodes on output, the next n+1 contain
c quadrature weights on output, and the final one is a
c scratch vector needed by gaussj.
c
      implicit complex(c,w,z)
      dimension qwork(1),betam(n)
c
c for each finite vertex w(k), compute nodes and weights for
c one-sided gauss-jacobi quadrature along a curve beginning at z(k):
      iscr = nptsq*(2*n+2) + 1
      do 1 k = 1,n
        inodes = nptsq*(k-1) + 1
        iwts = nptsq*(n+k) + 1
    1   if (betam(k).gt.-1.) call gaussj(nptsq,0.,betam(k),
     &    qwork(iscr),qwork(inodes),qwork(iwts))
c
c compute nodes and weights for pure gaussian quadrature:
      inodes = nptsq*n + 1
      iwts = nptsq*(2*n+1) + 1
      call gaussj(nptsq,0.,0.,
     &  qwork(iscr),qwork(inodes),qwork(iwts))
c
      return
      end


c*******************************************************************
c* scsolv                                     primary subroutine  **
c*******************************************************************
c
      subroutine scsolv(iprint,iguess,tol,errest,n,c,z,wc,
     &   w,betam,nptsq,qwork)
c
c this subroutine computes the accessory parameters c and
c z(k) for the schwarz-christoffel transformation
c which sends the unit disk to the interior of the polygon
c w(1),...,w(n).  this mapping is of the form:
c
c                         z    n
c    w   =   wc  +   c * int prod (1-z/z(k))**betam(k) dz .    (1.2)
c                         0   k=1
c
c the image polygon may be unbounded; permitted angles lie in the
c range -3.le.betam(k).le.1.  w(n) and w(1) must be finite.
c we normalize by the conditions:
c
c    z(n) = 1                                                  (2.1)
c    w(0.0) = wc  (a point in the interior of the polygon)     (2.1)
c
c calling sequence:
c
c   iprint  -2,-1,0, or 1 for increasing amounts of output (input)
c
c   iguess  1 if an initial guess for z is supplied, otherwise 0
c           (input)
c
c   tol     desired accuracy in solution of nonlinear system
c           (input).  recommended value: 10.**(-nptsq-1) * typical
c           size of vertices w(k)
c
c   errest  estimtated error in solution (output).  more
c           precisely, errest is an approximate bound for how far
c           the true vertices of the image polygon may be from those
c           computed by numerical integration using the
c           numerically determined prevertices z(k).
c
c   n       number of vertices of the image polygon (input).
c           must be .le. 20
c
c   c       complex scale factor in formula above (output)
c
c   z       complex array of prevertices on the unit circle.
c           dimension at least n.  if an initial guess is
c           being supplied it should be in z on input, with z(n)=1.
c           in any case the correct prevertices will be in z on output.
c
c   wc      complex image of 0 in the polygon, as in above formula
c           (input).  it is safest to pick wc to be as central as
c           possible in the polygon in the sense that as few parts
c           of the polygon as possible are shielded from wc by
c           reentrant edges.
c
c   w       complex array of vertices of the image polygon
c           (input).  dimension at least n.  it is a good idea
c           to keep the w(k) roughly on the scale of unity.
c           w(k) will be ignored when the vertex lies at infinity;
c           see betam, below.  each connected boundary component
c           must include at least one vertex w(k), even if it
c           has to be a degenerate vertex with betam(k) = 0.
c           w(n) and w(1) must be finite.
c
c   betam   real array with betam(k) the external angle in the
c           polygon at vertex k divided by minus pi (input).
c           dimension at least n.  permitted values lie in
c           the range -3.le.betam(k).le.1.  (examples: each
c           betam(k) is -1/2 for a rectangle, -2/3 for an equi-
c           lateral triangle, +1 at the end of a slit.)  the
c           sum of the betam(k) will be -2 if they have been
c           set correctly.  betam(n-1) should not be 0 or 1.
c           w(k) lies at infinity if and only if betam(k).le.-1.
c
c   nptsq   the number of points to be used per subinterval
c           in gauss-jacobi quadrature (input).  recommended
c           value: equal to one more than the number of digits
c           of accuracy desired in the answer.  must be the same
c           as in the call to qinit which filled the vector qwork.
c
c   qwork   real quadrature work array (input).  dimension
c           at least nptsq * (2n+3) but no greater than 460.
c           before calling scsolv qwork must have been filled
c           by subroutine qinit.
c
c the problem is solved by finding the
c solution to a system of n-1 nonlinear equations in the n-1
c unknowns y(1),...,y(n-1), which are related to the points
c z(k) by the formula:
c
c      y(k) = log ((th(k)-th(k-1))/(th(k+1)-th(k)))            (2.7)
c
c where th(k) denotes the argument of z(k).
c subroutine scfun defines this system of equations.
c the original problem is subject to the contraints th(k) < th(k+1),
c but these vanish in the transformation from z to y.
c
c reference:  l. n. trefethen, "numerical computation of the
c   schwarz-christoffel transformation," siam j. sci. stat. comp. 1
c   (1980), 82-102.  equation nos. above are taken from this paper.
c
c lloyd n. trefethen
c courant institute of mathematical sciences, new york university
c 251 mercer st.,   new york, ny 10012
c (212) 460-7224
c october 1979 (version 1); july 1983 (version 2)
c
      implicit complex(c,w,z)
      common /param1/ kfix(20),krat(20),ncomp,nptsq2,c2,
     &  qwork2(460),betam2(20),z2(20),wc2,w2(20)
      dimension z(n),w(n),betam(n),qwork(1)
      dimension ajinv(20,20),scr(900),fval(19),y(19)
      external scfun
      nm = n-1
c
c check input data:
      call check(tol,n,w,betam)
c
c determine number of boundary components, etc.:
c   pass 1: one fixed point for each infinite vertex:
      ncomp = 0
      do 1 k = 2,nm
        if (betam(k).gt.-1.) goto 1
        ncomp = ncomp + 1
        kfix(ncomp) = k - 1
        if (ncomp.eq.1) kfix(ncomp) = 1
    1 continue
      if (ncomp.gt.0) goto 2
      ncomp = 1
      kfix(ncomp) = 1
c   pass 2: one ratio for each line segment:
    2 continue
      neq = 2*ncomp
      do 3 k = 1,nm
        if (neq.eq.n-1) goto 4
        if (betam(k).le.-1..or.betam(k+1).le.-1.) goto 3
        neq = neq + 1
        krat(neq) = k
    3 continue
    4 z(n) = (1.,0.)
c
c initial guess, case iguess.eq.0:
c (vertices equally spaced around circle):
      if (iguess.ne.0) goto 11
      do 5 k = 1,nm
    5   y(k) = 0.
      goto 12
c
c initial guess, case iguess.ne.0:
c (vertices supplied on input):
   11 continue
      do 9 k = 1,nm
        km = k-1
        if (km.eq.0) km = n
        tmp1 = aimag(log(z(k+1)/z(k)))
        if (tmp1.lt.0.) tmp1 = tmp1 + 2. * acos(-1.)
        tmp2 = aimag(log(z(k)/z(km)))
        if (tmp2.lt.0.) tmp2 = tmp2 + 2. * acos(-1.)
    9   y(k) = log(tmp2) - log(tmp1)
   12 continue
c
c ns01a control parameters:
      dstep = .005
      dmax = 20.
      maxfun = (n-1) * 15
c
c copy input data to /param1/ for scfun:
c (this is necessary because ns01a requires a fixed calling
c sequence in subroutine scfun.)
      nptsq2 = nptsq
      wc2 = wc
      do 6 k = 1,n
        z2(k) = z(k)
        betam2(k) = betam(k)
    6   w2(k) = w(k)
      nwdim = nptsq * (2*n+3)
      do 7 i = 1,nwdim
    7   qwork2(i) = qwork(i)
c
c solve nonlinear system with ns01a:
      call ns01a(nm,y,fval,ajinv,dstep,dmax,tol,maxfun,
     &  iprint,scr,scfun)
c
c copy output data from /param1/:
      c = c2
      do 8 k = 1,nm
    8   z(k) = z2(k)
c
c print results and test accuracy:
      if (iprint.ge.0) call scoutp(n,c,z,wc,w,betam,nptsq)
      call sctest(errest,n,c,z,wc,w,betam,nptsq,qwork)
      if (iprint.ge.-1) write (6,201) errest
  201 format (' errest:',e12.4/)
      return
c
      end


c*******************************************************************
c* wsc                                        primary subroutine  **
c*******************************************************************
c
      function wsc(zz,kzz,z0,w0,k0,n,c,z,betam,nptsq,qwork)
c
c computes forward map w(zz)
c
c calling sequence:
c
c   zz     point in the disk at which w(zz) is desired (input)
c
c   kzz    k if zz = z(k) for some k, otherwise 0 (input)
c
c   z0     nearby point in the disk at which w(z0) is known and
c          finite (input)
c
c   w0     w(z0)  (input)
c
c   k0     k if z0 = z(k) for some k, otherwise 0 (input)
c
c   n,c,z,betam,nptsq,qwork     as in scsolv (input)
c
c convenient values of z0, w0, and k0 for most applications can be
c supplied by subroutine nearz.
c
      implicit complex(c,w,z)
      dimension z(n),betam(n),qwork(1)
c
      wsc = w0 + c * zquad(z0,k0,zz,kzz,n,z,betam,nptsq,qwork)
c
      return
      end


c*******************************************************************
c* zsc                                        primary subroutine  **
c*******************************************************************
c
      function zsc(ww,iguess,zinit,z0,w0,k0,eps,ier,n,c,
     &  z,wc,w,betam,nptsq,qwork)
c
c computes inverse map z(ww) by newton iteration
c
c calling sequence:
c
c   ww     point in the polygon at which z(ww) is desired (input)
c
c   iguess (input)
c          .eq.1 - initial guess is supplied as parameter zinit
c          .ne.1 - get initial guess from program ode (slower).
c                  for this the segment wc-ww must lie within
c                  the polygon.
c
c   zinit  initial guess if iguess.eq.1, otherwise ignored (input).
c          may not be a prevertex z(k)
c
c   z0     point in the disk near z(ww) at which w(z0) is known and
c          finite (input).
c
c   w0     w(z0)  (input).  the line segment from w0 to ww must
c          lie entirely within the closed polygon.
c
c   k0     k if z0 = z(k) for some k, otherwise 0 (input)
c
c   eps    desired accuracy in answer z(ww)  (input)
c
c   ier    error flag (input and output).
c          on input, give ier.ne.0 to suppress error messages.
c          on output, ier.ne.0 indicates unsuccessful computation --
c          try again with a better initial guess.
c
c   n,c,z,wc,w,betam,nptsq,qwork     as in scsolv (input)
c
c convenient values of z0, w0, and k0 for some applications can be
c supplied by subroutine nearw.
c
      implicit complex(c,w,z)
      dimension scr(142),iscr(5)
      dimension z(n),w(n),betam(n),qwork(1)
      external zfode
      logical odecal
      common /param2/ cdwdt,z2(20),betam2(20),n2
c
      odecal = .false.
      if (iguess.ne.1) goto 1
      zi = zinit
      goto 3
c
c get initial guess zi from program ode:
    1 n2 = n
      do 2 k = 1,n
        z2(k) = z(k)
    2   betam2(k) = betam(k)
      zi = (0.,0.)
      t = 0.
      iflag = -1
      relerr = 0.
      abserr = 1.e-4
      cdwdt = (ww-wc)/c
      call ode(zfode,2,zi,t,1.,relerr,abserr,iflag,scr,iscr)
      if (iflag.ne.2.and.ier.eq.0) write (6,201) iflag
      odecal = .true.
c
c refine answer by newton iteration:
    3 continue
      do 4 iter = 1,10
        zfnwt = ww - wsc(zi,0,z0,w0,k0,n,c,z,betam,nptsq,qwork)
        zi = zi + zfnwt/(c*zprod(zi,0,n,z,betam))
        if (abs(zi).ge.1.1) zi = .5 * zi/abs(zi)
        if (abs(zfnwt).lt.eps) goto 5
    4   continue
      if (.not.odecal) goto 1
      if (ier.eq.0) write (6,202)
      ier = 1
    5 zsc = zi
c
  201 format (/' *** nonstandard return from ode in zsc: iflag =',i2/)
  202 format (/' *** possible error in zsc: no convergence in 10'/
     &      '     iterations.  may need a better initial guess zinit')
      return
      end


c*******************************************************************
c* zfode                             subordinate(zsc) subroutine  **
c*******************************************************************
c
      subroutine zfode(t,zz,zdzdt)
c
c computes the function zdzdt needed by ode in zsc.
c
      implicit complex(c,w,z)
      common /param2/ cdwdt,z(20),betam(20),n
c
      zdzdt = cdwdt / zprod(zz,0,n,z,betam)
c
      return
      end


c*******************************************************************
c* check                          subordinate(scsolv) subroutine  **
c*******************************************************************
c
      subroutine check(eps,n,w,betam)
c
c checks geometry of the problem to make sure it is a form usable
c by scsolv.
c
      implicit complex(c,w,z)
      dimension w(n),betam(n)
c
      sum = 0.
      do 1 k = 1,n
    1   sum = sum + betam(k)
      if (abs(sum+2.).lt.eps) goto 2
      write (6,301)
    2 if (betam(1).gt.-1.) goto 3
      write (6,302)
      stop
    3 if (betam(n).gt.-1.) goto 4
      write (6,303)
      stop
    4 if (abs(betam(n-1)).gt.eps) goto 5
      write (6,304)
      write (6,306)
    5 if (abs(betam(n-1)-1.).gt.eps) goto 6
      write (6,305)
      write (6,306)
      stop
    6 do 7 k = 2,n
        if (betam(k).le.-1..or.betam(k-1).le.-1.) goto 7
        if (abs(w(k)-w(k-1)).le.eps) goto 8
    7   continue
      if (abs(w(1)-w(n)).gt.eps) goto 9
    8 write (6,307)
      stop
    9 if (n.ge.3) goto 10
      write (6,309)
      stop
   10 if (n.le.20) goto 11
      write (6,310)
      stop
   11 continue
      return
c
  301 format (/' *** error in check: angles do not add up to 2'/)
  302 format (/' *** error in check: w(1) must be finite'/)
  303 format (/' *** error in check: w(n) must be finite'/)
  304 format (/' *** warning in check: w(n-1) not determined'/)
  305 format (/' *** error in check: w(n-1) not determined')
  306 format (/'   renumber vertices so that betam(n-1) is not 0 or 1')
  307 format (/' *** error in check: two adjacent vertices are equal'/)
  309 format (/' *** error in check: n must be no less than 3'/)
  310 format (/' *** error in check: n must be no more than 20'/)
      end


c*******************************************************************
c* yztran                         subordinate(scsolv) subroutine  **
c*******************************************************************
c
      subroutine yztran(n,y,z)
c
c transforms y(k) to z(k).  see comments in subroutine scsolv.
c
      implicit complex(c,w,z)
      dimension y(n),z(n)
      nm = n - 1
      pi = acos(-1.)
c
      dth = 1.
      thsum = dth
      do 1 k = 1,nm
        dth = dth / exp(y(k))
    1   thsum = thsum + dth
c
      dth = 2. * pi / thsum
      thsum = dth
      z(1) = cmplx(cos(dth),sin(dth))
      do 2 k = 2,nm
        dth = dth / exp(y(k-1))
        thsum = thsum + dth
    2   z(k) = cmplx(cos(thsum),sin(thsum))
c
      return
      end


c*******************************************************************
c* scfun                          subordinate(scsolv) subroutine  **
c*******************************************************************
      subroutine scfun(ndim,y,fval)
c
c this is the function whose zero must be found in scsolv.
c
      implicit complex(c,w,z)
      dimension fval(ndim),y(ndim)
      common /param1/ kfix(20),krat(20),ncomp,nptsq,c,
     &  qwork(460),betam(20),z(20),wc,w(20)
      n = ndim+1
c
c transform y(k) to z(k):
      call yztran(n,y,z)
c
c set up: compute integral from 0 to z(n):
      wdenom = zquad((0.,0.),0,z(n),n,n,z,betam,nptsq,qwork)
      c = (w(n)-wc) / wdenom
c
c case 1: w(k) and w(k+1) finite:
c (compute integral along chord z(k)-z(k+1)):
      nfirst = 2*ncomp + 1
      if (nfirst.gt.n-1) goto 11
      do 10 neq = nfirst,ndim
        kl = krat(neq)
        kr = kl+1
        zint = zquad(z(kl),kl,z(kr),kr,n,z,betam,nptsq,qwork)
        fval(neq) = abs(w(kr)-w(kl)) - abs(c*zint)
   10 continue
c
c case 2: w(k+1) infinite:
c (compute contour integral along radius 0-z(k)):
   11 do 20 nvert = 1,ncomp
        kr = kfix(nvert)
        zint = zquad((0.,0.),0,z(kr),kr,n,z,betam,nptsq,qwork)
        zfval = w(kr) - wc - c*zint
        fval(2*nvert-1) = real(zfval)
        fval(2*nvert) = aimag(zfval)
   20 continue
      return
c
      end


c*******************************************************************
c* scoutp                         subordinate(scsolv) subroutine  **
c*******************************************************************
c
      subroutine scoutp(n,c,z,wc,w,betam,nptsq)
c
c prints a table of k, w(k), th(k), betam(k), and z(k).
c also prints the constants n, nptsq, wc, c.
c
      implicit complex(c,w,z)
      dimension z(n),w(n),betam(n)
c
      write (6,102) n, nptsq
      pi = acos(-1.)
      do 1 k = 1,n
        thdpi = aimag(log(z(k))) / pi
        if (thdpi.le.0.) thdpi = thdpi + 2.
        if (betam(k).gt.-1.) write (6,103) k,w(k),thdpi,betam(k),z(k)
    1   if (betam(k).le.-1.) write (6,104) k,thdpi,betam(k),z(k)
      write (6,105) wc,c
      return
c
  102 format (/' parameters defining map:',15x,'(n =',
     &  i3,')',6x,'(nptsq =',i3,')'//
     &  '  k',10x,'w(k)',10x,'th(k)/pi',5x,'betam(k)',
     &  13x,'z(k)'/
     &  ' ---',9x,'----',10x,'--------',5x,'--------',
     &  13x,'----'/)
  103 format (i3,'    (',f6.3,',',f6.3,')',f14.8,f12.5,
     &  3x,'(',f11.8,',',f11.8,')')
  104 format (i3,'        infinity   ',f14.8,f12.5,
     &  3x,'(',f11.8,',',f11.8,')')
  105 format (/' wc = (',e15.8,',',e15.8,')'/
     &          '  c = (',e15.8,',',e15.8,')'/)
      end


c*******************************************************************
c* sctest                         subordinate(scsolv) subroutine  **
c*******************************************************************
c
      subroutine sctest(errest,n,c,z,wc,w,betam,nptsq,qwork)
c
c tests the computed map for accuracy.
c
      implicit complex(c,w,z)
      dimension z(n),w(n),betam(n),qwork(1)
c
c test length of radii:
      errest = 0.
      do 10 k = 2,n
        if (betam(k).gt.-1.) rade = abs(wc -
     &    wsc((0.,0.),0,z(k),w(k),k,n,c,z,betam,nptsq,qwork))
        if (betam(k).le.-1.) rade = abs(wsc((.1,.1),0,
     &    z(k-1),w(k-1),k-1,n,c,z,betam,nptsq,qwork)
     &    - wsc((.1,.1),0,z(k+1),w(k+1),k+1,
     &    n,c,z,betam,nptsq,qwork))
        errest = max(errest,rade)
   10   continue
      return
      end


c*******************************************************************
c* zquad                                      primary subroutine  **
c*******************************************************************
c
      function zquad(za,ka,zb,kb,n,z,betam,nptsq,qwork)
c
c computes the complex line integral of zprod from za to zb along a
c straight line segment within the unit disk.  function zquad1 is
c called twice, once for each half of this integral.
c
      implicit complex(c,w,z)
      dimension z(n),betam(n),qwork(1)
c
      if (abs(za).gt.1.1.or.abs(zb).gt.1.1) write (6,301)
  301 format (/' *** warning in zquad: z outside the disk')
c
      zmid = (za + zb) / 2.
      zquad = zquad1(za,zmid,ka,n,z,betam,nptsq,qwork)
     &  - zquad1(zb,zmid,kb,n,z,betam,nptsq,qwork)
      return
      end


c*******************************************************************
c* zquad1                          subordinate(zquad) subroutine  **
c*******************************************************************
c
      function zquad1(za,zb,ka,n,z,betam,nptsq,qwork)
c
c computes the complex line integral of zprod from za to zb along a
c straight line segment within the unit disk.  compound one-sided
c gauss-jacobi quadrature is used, using function dist to determine
c the distance to the nearest singularity z(k).
c
      implicit complex(c,w,z)
      dimension z(n),betam(n),qwork(1)
      data resprm /2./
c
c check for zero-length integrand:
      if (abs(za-zb).gt.0.) goto 1
      zquad1 = (0.,0.)
      return
c
c step 1: one-sided gauss-jacobi quadrature for left endpoint:
    1 r = min(1.,dist(za,ka,n,z)*resprm/abs(za-zb))
      zaa = za + r*(zb-za)
      zquad1 = zqsum(za,zaa,ka,n,z,betam,nptsq,qwork)
c
c step 2: adjoin intervals of pure gaussian quadrature if necessary:
   10 if (r.eq.1.) return
      r = min(1.,dist(zaa,0,n,z)*resprm/abs(zaa-zb))
      zbb = zaa + r*(zb-zaa)
      zquad1 = zquad1 + zqsum(zaa,zbb,0,n,z,betam,nptsq,qwork)
      zaa = zbb
      goto 10
      end


c*******************************************************************
c* dist                            subordinate(zquad) subroutine  **
c*******************************************************************
c
      function dist(zz,ks,n,z)
c
c determines the distance from zz to the nearest singularity z(k)
c other than z(ks).
c
      implicit complex(c,w,z)
      dimension z(n)
c
      dist = 99.
      do 1 k = 1,n
        if (k.eq.ks) goto 1
        dist = min(dist,abs(zz-z(k)))
    1   continue
      return
      end


c*******************************************************************
c* zqsum                           subordinate(zquad) subroutine  **
c*******************************************************************
c
      function zqsum(za,zb,ka,n,z,betam,nptsq,qwork)
c
c computes the integral of zprod from za to zb by applying a
c one-sided gauss-jacobi formula with possible singularity at za.
c
      implicit complex(c,w,z)
      dimension z(n),betam(n),qwork(1)
c
      zs = (0.,0.)
      zh = (zb-za) / 2.
      zc = (za+zb) / 2.
      k = ka
      if (k.eq.0) k = n+1
      iwt1 = nptsq*(k-1) + 1
      iwt2 = iwt1 + nptsq - 1
      ioffst = nptsq*(n+1)
      do 1 i = iwt1,iwt2
    1   zs = zs + qwork(ioffst+i)*zprod(zc+zh*qwork(i),ka,n,z,betam)
      zqsum = zs*zh
      if (abs(zh).ne.0..and.k.ne.n+1)
     &  zqsum = zqsum*abs(zh)**betam(k)
      return
      end


c*******************************************************************
c* zprod                           subordinate(zquad) subroutine  **
c*******************************************************************
c
      function zprod(zz,ks,n,z,betam)
c
c computes the integrand
c
c           n
c         prod  (1-zz/z(k))**betam(k)  ,
c          k=1
c
c taking argument only (not modulus) for term k = ks.
c
c *** note -- in practice this is the innermost subroutine
c *** in scpack calculations.  the complex log calculation below
c *** may account for as much as half of the total execution time.
c
      implicit complex(c,w,z)
      dimension z(n),betam(n)
c
      zsum = (0.,0.)
      do 1 k = 1,n
        ztmp = (1.,0.) - zz/z(k)
        if (k.eq.ks) ztmp = ztmp / abs(ztmp)
    1   zsum = zsum + betam(k)*log(ztmp)
      zprod = exp(zsum)
      return
      end


c*******************************************************************
c* rprod                                      primary subroutine  **
c*******************************************************************
c
      function rprod(zz,n,z,betam)
c
c computes the absolute value of the integrand
c
c           n
c         prod  (1-zz/z(k))**betam(k)
c          k=1
c
      implicit complex(c,w,z)
      dimension z(n),betam(n)
c
      sum = 0.
      do 1 k = 1,n
        ztmp = (1.,0.) - zz/z(k)
    1   sum = sum + betam(k)*log(abs(ztmp))
      rprod = exp(sum)
      return
      end


c*******************************************************************
c* nearz                                      primary subroutine  **
c*******************************************************************
c
      subroutine nearz(zz,zn,wn,kn,n,z,wc,w,betam)
c
c returns information associated with the nearest prevertex z(k)
c to the point zz, or with 0 if 0 is closer than any z(k).
c zn = prevertex position, wn = w(zn), kn = prevertex no. (0 to n)
c
      implicit complex(c,w,z)
      dimension z(n),w(n),betam(n)
c
      dist = abs(zz)
      kn = 0
      zn = (0.,0.)
      wn = wc
      if (dist.le..5) return
      do 1 k = 1,n
        if (betam(k).le.-1.) goto 1
        distzk = abs(zz-z(k))
        if (distzk.ge.dist) goto 1
        dist = distzk
        kn = k
    1   continue
      if (kn.eq.0) return
      zn = z(kn)
      wn = w(kn)
      return
      end


c*******************************************************************
c* nearw                                      primary subroutine  **
c*******************************************************************
c
      subroutine nearw(ww,zn,wn,kn,n,z,wc,w,betam)
c
c returns information associated with the nearest vertex w(k)
c to the point ww, or with wc if wc is closer than any w(k).
c zn = prevertex position, wn = w(zn), kn = vertex no. (0 to n)
c
      implicit complex(c,w,z)
      dimension z(n),w(n),betam(n)
c
      dist = abs(ww-wc)
      kn = 0
      zn = (0.,0.)
      wn = wc
      do 1 k = 1,n
        if (betam(k).le.-1.) goto 1
        distwk = abs(ww-w(k))
        if (distwk.ge.dist) goto 1
        dist = distwk
        kn = k
    1   continue
      if (kn.eq.0) return
      zn = z(kn)
      wn = w(kn)
      return
      end


c*******************************************************************
c* angles                                     primary subroutine  **
c*******************************************************************
c
      subroutine angles(n,w,betam)
c
c computes external angles -pi*betam(k) from knowledge of
c the vertices w(k).  an angle betam(k) is computed for each
c k for which w(k-1), w(k), and w(k+1) are finite.
c to get this information across any vertices at infinity
c should be signaled by the value w(k) = (99.,99.) on input.
c
      implicit complex(c,w,z)
      dimension w(n),betam(n)
      c9 = (99.,99.)
c
      pi = acos(-1.)
      do 1 k = 1,n
        km = mod(k+n-2,n)+1
        kp = mod(k,n)+1
        if (w(km).eq.c9.or.w(k).eq.c9.or.w(kp).eq.c9) goto 1
        betam(k) = aimag(log((w(km)-w(k))/(w(kp)-w(k))))/pi - 1.
        if (betam(k).le.-1.) betam(k) = betam(k) + 2.
    1   continue
      return
      end