program test1 implicit complex(c,w,z) dimension z(20),w(20),betam(20),qwork(344) zero = (0.,0.) zi = (0.,1.) c c set up problem: n = 4 wc = cmplx(0.,sqrt(2.)) w(1) = zi w(2) = zero w(3) = (1.e20,1.e20) w(4) = zero betam(1) = 1. betam(2) = -.5 betam(3) = -2. betam(4) = -.5 c c compute nodes and weights for parameter problem: nptsq = 5 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(cmplx(0.,k-4.)) tol = 1.e-6 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 = (.3,0.) * cmplx(i-2.,.2*i+.5) 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.,0.)) 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 = cmplx(i-2.,sqrt(i+1.)) ier = 0 zz = zsc(ww,0,zero,zero,wc,0,tol,ier, & n,c,z,wc,w,betam,nptsq,qwork) wtmp = zi * sqrt((-1.,0.)-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('(',f6.3,',',f6.3,') '),d11.4) 202 format (' w,z,zex,err: ',3('(',f6.3,',',f6.3,') '),d11.4) end program test2 implicit complex (c,w,z) dimension w(10),ibrk(4),qwork(300) n = 6 w(1) = (0.,0.) w(2) = (2.,0.) w(3) = (2.,1.) w(4) = (1.,1.) w(5) = (1.,2.) w(6) = (0.,2.) wc = (.5,.5) ibrk(1) = 1 ibrk(2) = 2 ibrk(3) = 5 ibrk(4) = 6 c c main loop: different accuracy specifications: do 10 ndig = 1,5,2 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 courant institute of mathematical sciences c new york university c november 1979, revised july 1983 c implicit complex(c,w,z) 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.,0.) pi = acos(-1.) 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-5, 10.**(-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.,0.) do 9 k = 1,4 betam2(k) = -.5 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.*eest) / xlen1 xwid1 = abs(w2(2)-w2(1)) xwid2 = abs(w2(4)-w2(3)) errw = max(abs(xwid1-xwid2),2.*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