program test1
implicit complex*16(c,w,z), real*8(a-b,d-h,o-v,x-y)
dimension z(20),w(20),betam(20),qwork(344)
zero = (0.d0,0.d0)
zi = (0.d0,1.d0)
c
c set up problem:
n = 4
wc = dcmplx(0.d0,sqrt(2.d0))
w(1) = zi
w(2) = zero
w(3) = (1.d20,1.d20)
w(4) = zero
betam(1) = 1.d0
betam(2) = -.5d0
betam(3) = -2.d0
betam(4) = -.5d0
c
c compute nodes and weights for parameter problem:
nptsq =12
call qinit(n,betam,nptsq,qwork)
c
c solve parameter problem:
c (initial guess must be given to avoid accidental exact solution)
iprint = 0
iguess = 1
do 1 k = 1,4
1 z(k) = exp(dcmplx(0.d0,k-4.d0))
tol = 1.d-14
call scsolv(iprint,iguess,tol,errest,n,c,z,
& wc,w,betam,nptsq,qwork)
c
c compare wsc(z) to exact values for various z:
do 10 i = 1,4
zz = (.3d0,0.d0) * dcmplx(i-2.d0,.2d0*i+.5d0)
ww = wsc(zz,0,zero,wc,0,n,c,z,betam,nptsq,qwork)
ztmp = -zi * (zz-zi) / (zz+zi)
wwex = zi * sqrt(-ztmp**2 + (1.d0,0.d0))
err = abs(ww-wwex)
write (6,201) zz,ww,wwex,err
10 continue
write (6,200)
c
c compare zsc(w) to exact values for various w:
do 20 i = 1,6
ww = dcmplx(i-2.d0,sqrt(i+1.d0))
ier = 0
zz = zsc(ww,0,zero,zero,wc,0,tol,ier,
& n,c,z,wc,w,betam,nptsq,qwork)
wtmp = zi * sqrt((-1.d0,0.d0)-ww**2)
zzex = -zi * (wtmp-zi) / (wtmp+zi)
err = abs(zz-zzex)
write (6,202) ww,zz,zzex,err
20 continue
c
stop
200 format (1x)
201 format (' z,w,wex,err: ',3('(',f12.9,',',f12.9,') '),d11.4)
202 format (' w,z,zex,err: ',3('(',f12.9,',',f12.9,') '),d11.4)
end
program test2
implicit complex*16(c,w,z), real*8(a-b,d-h,o-v,x-y)
dimension w(10),ibrk(4),qwork(300)
n = 6
w(1) = (0.d0,0.d0)
w(2) = (2.d0,0.d0)
w(3) = (2.d0,1.d0)
w(4) = (1.d0,1.d0)
w(5) = (1.d0,2.d0)
w(6) = (0.d0,2.d0)
wc = (.5d0,.5d0)
ibrk(1) = 1
ibrk(2) = 2
ibrk(3) = 5
ibrk(4) = 6
c
c main loop: different accuracy specifications:
do 10 ndig = 1,11,5
write (6,202) ndig
r = resist(n,w,wc,ibrk,ndig,errest,qwork)
write (6,201) r,errest
10 continue
stop
c
201 format (' r,errest:',2d23.15)
202 format (/' ndig =',i3,':')
end
function resist(n,w,wc,ibrk,ndig,errest,qwork)
c
c this function returns the resistance of a polygonally
c shaped resistor. computations are based on the schwarz-
c christoffel transformation. we normalize by assuming that
c a square has resistance 1.
c
c input parameters:
c
c n number of vertices of the polygon. n must satisfy
c 4 .le. n .le. 20.
c
c w complex array of dimension at least n containing
c the positions of the vertices, viewed as complex numbers.
c the vertices must be listed in counterclockwise order
c around the polygon. it is a good idea to keep the
c w(k) roughly on the scale of unity.
c
c wc a point in the interior of the polygon. try to
c pick wc as central as possible in the sense that
c as few parts of the polygon as possible are shielded
c from it by reentrant edges.
c
c ibrk array of dimension at least 4 containing indices of
c the vertices which define the breaks
c between constant-voltage and insulated portions
c of the boundary. the program will assume that
c the voltage is applied between side w(ibrk(1))-w(ibrk(2))
c and side w(ibrk(3))-w(ibrk(4)), with the other two
c sides insulated. the break vertices must be numbered
c in counterclockwise order; thus the integers ibrk(i)
c must increase with i, except that they may wrap
c once around n.
c
c ndig input integer giving the desired number of digits
c of accuracy in the result. ndig must be at least 1.
c it should be no more than a couple of digits less
c than full-word precision.
c
c output parameter:
c
c errest rough but conservative estimate of the size of
c the error in the value returned.
c
c work space parameter:
c
c qwork real work array. dimension at least (ndig+1)*(2n+3),
c but no more than 460.
c
c
c some advice:
c
c the program is not infallible. if it yields strange
c messages about not converging, it's probably best to
c try a simpler geometry.
c
c the amount of time required is roughly proportional to ndig
c and also roughly proportional to n**3.
c thus problems with many corners can
c take quite a while -- several minutes of cpu time on
c the ibm 370-168 for a problem with n = 20.
c
c lloyd n. trefethen
c department of mathematics
c massachusetts institute of technology
c november 1979, revised july 1983
c
implicit complex*16(c,w,z), real*8(a-b,d-h,o-v,x-y)
c map from disk to resistor:
dimension z(20),w(1),ibrk(1),betam(20),qwork(1)
c map from disk to rectangle with equal resistance:
dimension z2(4),w2(4),betam2(4)
c
zero = (0.d0,0.d0)
pi = acos(-1.d0)
c
c compute angles and check input parameters:
call angles(n,w,betam)
if (ndig.lt.1) write (6,101)
if (ndig.lt.1) stop
if (n.lt.4) write (6,102)
if (n.lt.4) stop
do 1 k = 1,4
if (ibrk(k).gt.n .or. ibrk(k).lt.1) goto 2
1 continue
goto 3
2 write (6,103)
stop
3 continue
c
c set up first map:
nptsq = ndig + 1
tol = max(1.e-12, 10.d0**(-nptsq))
call qinit(n,betam,nptsq,qwork)
call scsolv(-2,0,tol,eest,n,c,z,wc,w,betam,nptsq,qwork)
c
c set up map to rectangle:
c (qwork is overwritten to save space)
n2 = 4
c2 = (1.d0,0.d0)
do 9 k = 1,4
betam2(k) = -.5d0
9 z2(k) = z(ibrk(k))
nptsq2 = ndig + 1
call qinit(n2,betam2,nptsq2,qwork)
do 12 k = 1,4
12 w2(k)=wsc(z2(k),k,zero,zero,0,n2,c2,z2,betam2,nptsq2,qwork)
c
c compute length, width, resistance, and error estimate:
xlen1 = abs(w2(2)-w2(3))
xlen2 = abs(w2(4)-w2(1))
errl = max(abs(xlen1-xlen2),2.d0*eest) / xlen1
xwid1 = abs(w2(2)-w2(1))
xwid2 = abs(w2(4)-w2(3))
errw = max(abs(xwid1-xwid2),2.d0*eest) / xwid1
resist = (xlen1+xlen2) / (xwid1+xwid2)
errest = resist * (errl + errw)
c
return
c
101 format (' *** error in resist *** ndig should be at',
& ' least 1.'/)
102 format (' *** error in resist *** n must be no',
& ' greater than 20'/' and no smaller than 4')
103 format (' *** error in resist *** each ibrk(i) must',
& ' be in the',/,' range from 1 to n')
end
program test3
c
c computes the conformal map from a polygon (1) onto a rectangle (2)
c and plots a corresponding grid of size nx by ny.
c the corners of the rectangle are (0,i), (0,0), (h,0), (h,i).
c
c nick trefethen, icase, july 1983
c
implicit complex(c,w,z)
dimension wgrid(0:41,0:41)
dimension z1(12),w1(12),betam1(12),qwork1(270),ibrk(4)
dimension z2(4),w2(4),betam2(4),qwork2(110)
data zero /(0.,0.)/
c
c input data
print 90
90 format (' vertices? (terminate with re w(k) > 100 )')
k = 0
91 k = k+1
read *,x,y
w1(k) = cmplx(x,y)
if (x.lt.100.) goto 91
n1 = k-1
print 92
92 format (' image of zero in the polygon?')
read *,x,y
wc1 = cmplx(x,y)
print 93
93 format (' four distinguished vertices?')
read *,(ibrk(k),k=1,4)
nq1 = 3
nq2 = nq1
tol = 10.**(-(nq1+1))
print 94
94 format (' no. of streamlines, equipotential lines?')
read *,ny,nx
nxp = nx + 1
nyp = ny + 1
c
c compute parameters for map from disk to polygon
call angles(n1,w1,betam1)
call qinit(n1,betam1,nq1,qwork1)
call scsolv(0,0,tol,eest,n1,c1,z1,wc1,w1,betam1,nq1,qwork1)
c
c compute parameters for map from disk to rectangle
n2 = 4
c2 = (1.,0.)
do 2 k = 1,4
betam2(k) = -.5
2 z2(k) = z1(ibrk(k))
call qinit(n2,betam2,nq2,qwork2)
do 3 k = 1,4
3 w2(k) = wsc(z2(k),k,zero,zero,0,n2,c2,z2,betam2,nq2,qwork2)
c2 = (0.,1.) / (w2(1)-w2(2))
wc2 = -c2*w2(2)
do 4 k = 1,4
4 w2(k) = c2*w2(k) + wc2
write (6,102) (w2(k),k=1,4)
102 format (' vertices of image rectangle: ',
& 4(/' (',e13.6,',',e13.6,')')/)
h = abs(w2(3)-w2(2))
c
c compute grid points
do 12 j = 0,nyp
i1 = 0
i2 = nxp
if (j.eq.0.or.j.eq.nyp) i1 = 1
if (j.eq.0.or.j.eq.nyp) i2 = nx
do 11 i = i1,i2
ww2 = cmplx((i*h)/nxp,float(j)/nyp)
call nearw(ww2,zn2,wn2,kn2,n2,z2,wc2,w2,betam2)
ier = 0
iguess = 1
if (i.eq.i1) iguess = 0
zz = zsc(ww2,iguess,zz,zn2,wn2,kn2,1.e-3,ier,n2,c2,
& z2,wc2,w2,betam2,nq2,qwork2)
call nearz(zz,zn1,wn1,kn1,n1,z1,wc1,w1,betam1)
wgrid(i,j) = wsc(zz,0,zn1,wn1,kn1,n1,c1,z1,betam1,nq1,qwork1)
11 continue
12 write (6,105) j,nyp
105 format (' finished row',i3,'/',i2,' of grid points')
c
c draw plot
do 5 k = 1,n1
5 write (10,103) w1(k)
write (10,104) w1(1)
103 format (2f10.5)
104 format (2f10.5,'" "')
do 6 j = 1,ny
write (10,103) (wgrid(i,j),i=0,nx)
6 write (10,104) wgrid(nxp,j)
do 7 i = 1,nx
write (10,103) (wgrid(i,j),j=0,ny)
7 write (10,104) wgrid(i,nyp)
c
stop
end
.