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