SSSSS CCCCC PPPPPP AAA CCCCC KK KK SC HW AR ZC HR IS TO FF EL PA CK AGE SS CC PP PP AA AA CC KK KK SS CC PP PP AA AA CC KK KKK SSS CC PPPPBY AAAAAAA CC KKKKK SS CC PP AA AA CC KK KKK SS CC PP AA AA CC KK KK LL OY DN IC HO LA ST RE FE TH EN SSSSS CCCCC PP AA AA CCCCC KK KK Version 2 USER'S GUIDE Description of a collection of programs for computing the Schwarz-Christoffel conformal map between a disk and a polygon Lloyd N. Trefethen Courant Institute of Mathematical Sciences New York University 251 Mercer St. New York, NY 10012 July 1983 CONTENTS INTRODUCTION... 1 STRUCTURE OF THE PACKAGE... 2 Capabilities and speed Fortran environment Differences between Versions 1 and 2 Routines included Machine dependent constants Common blocks USE OF THE PACKAGE... 6 Outline QINIT: initialization for Gauss-Jacobi quadrature SCSOLV: solution of parameter problem WSC: evaluation of the forward map ZSC: evaluation of the inverse map EXAMPLES... 11 TEST1: check WSC and ZSC against exact solution TEST2: call RESIST to compute a resistance TEST3: generate a rectilinear grid PROGRAM LISTING (excluding library routines)... 18 INTRODUCTION The SCPACK collection contains Fortran routines for computing the Schwarz-Christoffel transformation, a conformal map that maps the interior of the unit disk in the complex plane onto the interior of a polygon P with vertices W(k), k=1,...,N: w = w ( z ) polygon disk The image polygon may have vertices at infinity, so that it is really any simply-connected region in the plane whose boundary consists of a finite number of straight lines. The conformal map can be applied to a variety of problems involving a polyg- onal problem domain: for example, solving Laplace's equation with Dirichlet, Neumann, or mixed boundary conditions; solving Poisson's equation; finding eigenvalues of the Laplace opera- tor; hodograph computations for ideal free-streamline flow; grid generation. The design of SCPACK is described along with some elementary applications in L.N. Trefethen, "Numerical Computation of the Schwarz-Christoffel Transformation," SIAM J. Sci. Stat. Comput. 1 (1080), 82-102. To obtain copies of this User's Guide, or to arrange to obtain a machine-readable copy of the code, contact me at the address given on the title page. For quick response, send a 9 track, odd parity, unalabeled tape of density 1600 bpi along with specifications of blocksize, record length, and ASCII or EBCDIC. Use of the package is unrestricted, except that a reference should be given in any publications that make use of it. I would be pleased to receive reprints of any such publica- tions, or less formal descriptions of applications of SCPACK. Any suggestions for improvements in the program or its document- ation are eagerly solicited. STRUCTURE OF THE PACKAGE Capabilities and speed SCPACK contains routines to solve the "parameter problem" associated with the Schwarz-Christoffel map (subroutine SCSOLV), to evaluate the resulting S-C map (WSC), and to evaluate its inverse (ZSC). At present the package has been tested on VAX, IBM 370, and CYBER 173,175 computers. A typical computation requires on the order of minutes on a Vax 750 and on the order of seconds on the other machines. See the paper referenced in the Introduction for a discussion of execution speed. The bottleneck in all SCPACK computations is a calculation of a complex logarithm within subroutine ZPROD. Some optimization here might speed up the package at particular installations. For high-accuracy conformal mapping of a polygon, the Schwarz- Christoffel approach is extremely satisfactory because it handles the singularities at corners exactly and reduces the map to a finite number of parameters. On the other hand it is not recommended that SCPACK be used for mapping curved domains by means of polygonal approximations. Many superior methods based on integral equations exist for such problems; for references see the recent numerical analysis literature or one of the following books: - D. Gaier, Konstruktive Methoden in der Konformen Abbildung, Springer, 1964. - P. Henrici, Applied and Computational Complex Analysis III, Wiley-Interscience, to appear around 1985. Fortran environment SCPACK is written in a subset of Fortran 77 that should be acceptable to most Fortran 66 compilers. All real numbers are single precision and all complex numbers are complex single precision. Computations involving complex quantities are always performed in complex arithmetic, and variables whose names begin with C, W, or Z are always com- plex. Each routine (other than those in the library routines described below) begins with the statement IMPLICIT COMPLEX(C,W,Z) to make these type conventions automatic. Differences between Versions 1 and 2 SCPACK Version 1 (= 1.1) was developed on an IBM 370 installation at the Stanford Linear Accelerator Center. Version 2 was developed on a VAX 750 / Unix machine at the Institute for Computer Applic- ations in Science and Engineering (ICASE), NASA Langley Research Center. The main difference is that Version 2 is written in single rather than double precision and attempts to be more portable. Here are some specific changes from Version 1: - All real and complex variables are single precision. This change was made mainly because complex double precision is not standard in Fortran 77. If you are using SCPACK on a machine with a short word length (e.g. 32 bits), it is recommended that you compile the code with an automatic precision doubling feature, if available. On the IBM 370 this is the AUTODBL feature. - Generic functions MAX, MIN, ABS, LOG, EXP, ACOS, SIN, COS, REAL, CMPLX, SIGN, and MOD are used instead of specific functions AMIN1, CABS, ALOG, etc. It is assumed that the compiler auto- matically provides results of the same type as the argument. - PLTCON, PLTSCL, FINITE, and ZFNEWT have been removed. A routine NEARW analogous to NEARZ has been added. - The calling sequences of routines QINIT, WSC, ZSC, and NEARZ have been changed. - A gamma function routine GAMMA has been included in GAUSSJ, since some compilers do not include this as an intrinsic function. Routines included SCPACK consists of 19 subroutines, 4 of which the user will control directly. These are: QINIT computes Gauss-Jacobi quadrature nodes and weights SCSOLV solves the mapping problem for a given polygon WSC evaluates the forward map w(z) ZSC evaluates the inverse map z(w) Depending on the application, the user may also wish to make use of ANGLES computes angles, given vertices RPROD computes |w'(z)|**2 at a point z NEARZ returns the number of the nearest prevertex, etc. NEARW returns the number of the nearest vertex, etc. ZQUAD computes the S-C integral between two points Here are the other internal routines of SCPACK, which the user will not normally need direct access to: CHECK SCFUN SCTEST ZQUAD1 ZQSUM YZTRAN SCOUTP ZFODE DIST ZPROD In addition, three sets of library routines are made use of. NS01A, by M.J.D. Powell, solves a nonlinear system of equa- tions whose solution is required for the parameter problem. NS01A includes a few linear algebra routines from LINPACK to invert a matrix. All together: NS01A SGEFA SGEDI SSWAP SAXPY SSCAL ISAMAX QINIT requires GAUSSJ and its subordinate routines, by G.H. Golub and J.H. Welsch, to compute nodes and weights for Gauss-Jacobi quadrature: GAUSSJ CLASS IMTQL2* GAMMA* Finally (and least essentially) ZSC requires ODE, by Shampine and Gordon, to solve an ordinary differential equation: ODE DE* STEP* INTRP The implementer is free to experiment with other comparable routines to replace these three packages; SCPACK should still work satisfactorily after such substitutions. However, all of the above routines are normally included on an SCPACK tape. They are all in the public domain. (Subroutine GAMMA is an adaptation from the proprietary IMSL library, but it may be freely called as part of the GAUSSJ group provided it is not separated for general use.) *Machine dependent constants The routines marked with asterisks in the list above contain machine dependent constants, as follows: IMTQL2: MACHEP = U DE: FOURU = 4*U STEP: FOURU = 4*U , TWOU = 2*U GAMMA: various quantities where U is the smallest floating point number such that 1+U > U. In the standard SCPACK tape these numbers are set to generous values that should give satisfactory although not ideal perfor- mance in single precision on most machines. The user compiling SCPACK in double precision should replace the default values with correct ones if possible. However, SCPACK will perform more reliably in double than in single precision even if the machine- dependent constants are not changed. Subroutine GAMMA contains many machine dependent constants and should be replaced by a locally available gamma function routine if the user wants more than 5 or 6 places of accuracy. Common blocks The package makes use of two labeled common blocks, /PARAM1/ and /PARAM2/. These are needed for getting informa- tion to the functions SCFUN and ZFODE called by NS01A and ODE, respectively, whose calling sequences are fixed. USE OF THE PACKAGE Outline To compute the Schwarz-Christoffel transformation from the unit disk to a given polygon, here is what to do. First, set N (number of vertices), W(k) (the vertices), NPTSQ (number of quadrature points per subinterval), and BETAM(k) (external angles divided by minus pi). (Each angle BETAM(k) can be com- puted by subroutine ANGLES from the vertices W(k-1), W(k), and W(k+1) if these are all finite.) Now call QINIT to compute nodes and weights for Gauss-Jacobi quadrature. Second, set WC (a point interior to the polygon) and the other input parame- ters to SCSOLV, and call SCSOLV to solve the parameter problem for the constant C and the prevertices Z(k). Third, map indi- vidual points as desired from the disk to the polygon with routine WSC and from the polygon to the disk with routine ZSC. More complicated problems may require other sequences of computations. For example, two or more S-C maps can be com- posed in arbitrary ways by calling QINIT and SCSOLV twice, using of course distinct variables in the calling sequences to them. This is necessary whenever one polygon is mapped onto another with the disk as an intermediate domain. See the example TEST3 below for an example of composition of maps. QINIT: initialization for Gauss-Jacobi quadrature QINIT must be called before SCSOLV, WSC, or ZSC. It com- putes Gauss-Jacobi quadrature nodes and weights which make integration of the Schwarz-Christoffel formula possible. SUBROUTINE QINIT(N,BETAM,NPTSQ,QWORK) N, BETAM, and NPTSQ are input parameters. All four argu- ments are described in the SCSOLV documentation below. NPTSQ is the number of quadrature points which will be used on each quadrature subinterval. This is roughly equal to the number of digits of accuracy to which integrals will be calculated. Thus typical values might be 2, 4, or 6. Total SCPACK compu- tation time increases roughly linearly with NPTSQ. QWORK is a work array that must be dimensioned at least NPTSQ*(2N+3), which on output will contain quadrature nodes and weights corresponding to the various vertices of the poly- gon. The same array with the same NPTSQ must be given later as input to SCSOLV, WSC, or ZSC. If more than one S-C map is to be composed, then more than one copy of QWORK must be filled by QINIT and kept for SCSOLV, WSC, or ZSC. SCSOLV: solution of the parameter problem SCSOLV solves the parameter problem for the constants C and the prevertices Z(k) of the S-C map. This map is of the form Z N W = WC + C * INT PROD (1-Z/Z(K))**BETAM(K) DZ . 0 K=1 If the problem involves a map of the disk onto a given poly- gon, SCSOLV must be called before this map is evaluated with WSC or ZSC. (In rarer applications where the prevertices Z(k) and angles BETA(k) are known at the outset rather than the vertices W(k), SCSOLV is not needed. This occurs below in example TEST3.) SUBROUTINE SCSOLV(IPRINT,IGUESS,TOL,ERREST,N,C,Z,WC, & W,BETAM,NPTSQ,QWORK) IPRINT -2,-1,0, OR 1 FOR INCREASING AMOUNTS OF OUTPUT (INPUT) IGUESS 1 IF AN INITIAL GUESS FOR Z IS SUPPLIED, OTHERWISE 0 (INPUT) TOL DESIRED ACCURACY IN SOLUTION OF NONLINEAR SYSTEM (INPUT). RECOMMENDED VALUE: 10.**(-NPTSQ-1) * TYPICAL SIZE OF VERTICES W(K) ERREST ESTIMTATED ERROR IN SOLUTION (OUTPUT). MORE PRECISELY, ERREST IS AN APPROXIMATE BOUND FOR HOW FAR THE TRUE VERTICES OF THE IMAGE POLYGON MAY BE FROM THOSE COMPUTED BY NUMERICAL INTEGRATION USING THE NUMERICALLY DETERMINED PREVERTICES Z(K). N NUMBER OF VERTICES OF THE IMAGE POLYGON (INPUT). MUST BE .LE. 20 C COMPLEX SCALE FACTOR IN FORMULA ABOVE (OUTPUT) Z COMPLEX ARRAY OF PREVERTICES ON THE UNIT CIRCLE. DIMENSION AT LEAST N. IF AN INITIAL GUESS IS BEING SUPPLIED IT SHOULD BE IN Z ON INPUT, WITH Z(N)=1. IN ANY CASE THE CORRECT PREVERTICES WILL BE IN Z ON OUTPUT. WC COMPLEX IMAGE OF 0 IN THE POLYGON, AS IN ABOVE FORMULA (INPUT). IT IS SAFEST TO PICK WC TO BE AS CENTRAL AS POSSIBLE IN THE POLYGON IN THE SENSE THAT AS FEW PARTS OF THE POLYGON AS POSSIBLE ARE SHIELDED FROM WC BY REENTRANT EDGES. W COMPLEX ARRAY OF VERTICES OF THE IMAGE POLYGON (INPUT). DIMENSION AT LEAST N. IT IS A GOOD IDEA TO KEEP THE W(K) ROUGHLY ON THE SCALE OF UNITY. W(K) WILL BE IGNORED WHEN THE VERTEX LIES AT INFINITY; SEE BETAM, BELOW. EACH CONNECTED BOUNDARY COMPONENT MUST INCLUDE AT LEAST ONE VERTEX W(K), EVEN IF IT HAS TO BE A DEGENERATE VERTEX WITH BETAM(K) = 0. W(N) AND W(1) MUST BE FINITE. BETAM REAL ARRAY WITH BETAM(K) THE EXTERNAL ANGLE IN THE POLYGON AT VERTEX K DIVIDED BY MINUS PI (INPUT). DIMENSION AT LEAST N. PERMITTED VALUES LIE IN THE RANGE -3.LE.BETAM(K).LE.1. (EXAMPLES: EACH BETAM(K) IS -1/2 FOR A RECTANGLE, -2/3 FOR AN EQUI- LATERAL TRIANGLE, +1 AT THE END OF A SLIT.) THE SUM OF THE BETAM(K) WILL BE -2 IF THEY HAVE BEEN SET CORRECTLY. BETAM(N-1) SHOULD NOT BE 0 OR 1. W(K) LIES AT INFINITY IF AND ONLY IF BETAM(K).LE.-1. NPTSQ THE NUMBER OF POINTS TO BE USED PER SUBINTERVAL IN GAUSS-JACOBI QUADRATURE (INPUT). RECOMMENDED VALUE: EQUAL TO ONE MORE THAN THE NUMBER OF DIGITS OF ACCURACY DESIRED IN THE ANSWER. MUST BE THE SAME AS IN THE CALL TO QINIT WHICH FILLED THE VECTOR QWORK. QWORK REAL QUADRATURE WORK ARRAY (INPUT). DIMENSION AT LEAST NPTSQ * (2N+3) BUT NO GREATER THAN 460. BEFORE CALLING SCSOLV QWORK MUST HAVE BEEN FILLED BY SUBROUTINE QINIT. WSC: evaluation of the forward map Once the parameter problem has been solved, the function WSC is available to solve the forward map: w = w(z). WSC is a complex function with calling sequence FUNCTION WSC(ZZ,KZZ,Z0,W0,K0,N,C,Z,BETAM,NPTSQ,QWORK) The value returned will be the image of ZZ in the polygon. ZZ POINT IN THE DISK AT WHICH W(ZZ) IS DESIRED (INPUT) KZZ K IF ZZ = Z(K) FOR SOME K, OTHERWISE 0 (INPUT) Z0 NEARBY POINT IN THE DISK AT WHICH W(Z0) IS KNOWN AND FINITE (INPUT) W0 W(Z0) (INPUT) K0 K IF Z0 = Z(K) FOR SOME K, OTHERWISE 0 (INPUT) N,C,Z,BETAM,NPTSQ,QWORK AS IN SCSOLV (INPUT) Z0 should not be too far from ZZ. A simple and adequate choice is the nearest prevertex Z(k) with W(k) finite, or 0 if 0 is closer than any Z(k). This choice of Z0 can be deter- mined automatically by the function NEARZ, which is documented in its comment cards. See TEST3 below for an example. ZSC: evaluation of the inverse map ZSC is parallel to WSC in usage but evaluates the inverse map z = z(w). The calling sequence is FUNCTION ZSC(WW,IGUESS,ZINIT,Z0,W0,K0,EPS,IER,N,C, & Z,WC,W,BETAM,NPTSQ,QWORK) The value returned will be the preimage of WW in the disk. WW POINT IN THE POLYGON AT WHICH Z(WW) IS DESIRED (INPUT) IGUESS (INPUT) .EQ.1 - INITIAL GUESS IS SUPPLIED AS PARAMETER ZINIT .NE.1 - GET INITIAL GUESS FROM PROGRAM ODE (SLOWER). FOR THIS THE SEGMENT WC-WW MUST LIE WITHIN THE POLYGON. ZINIT INITIAL GUESS IF IGUESS.EQ.1, OTHERWISE IGNORED (INPUT). MAY NOT BE A PREVERTEX Z(K) Z0 POINT IN THE DISK NEAR Z(WW) AT WHICH W(Z0) IS KNOWN AND FINITE (INPUT). W0 W(Z0) (INPUT). THE LINE SEGMENT FROM W0 TO WW MUST LIE ENTIRELY WITHIN THE POLYGON. K0 K IF Z0 = Z(K) FOR SOME K, OTHERWISE 0 (INPUT) EPS DESIRED ACCURACY IN ANSWER Z(WW) (INPUT) IER ERROR FLAG (INPUT AND OUTPUT). ON INPUT, GIVE IER.NE.0 TO SUPPRESS ERROR MESSAGES. ON OUTPUT, IER.NE.0 INDICATES UNSUCCESSFUL COMPUTATION -- TRY AGAIN WITH A BETTER INITIAL GUESS. N,C,Z,WC,W,BETAM,NPTSQ,QWORK AS IN SCSOLV (INPUT) As with WSC, here W0 should not be too far from WW. The routine NEARW selects WC or the nearest vertex W(k) for W0, a choice that should be suitable for many applications. NEARW is documented in its comment cards. Note that the segment W0-WW must lie entirely within the image polygon. In treating non-convex polygons the user must bear this in mind when calling ZSC -- it may occasionally be necessary to first find the inverse image of an intermediate point W0'. EXAMPLES TEST1: check WSC and ZSC against exact solution The first program computes the S-C map in a simple case where it is known analytically. Results for both the forward map WSC and the inverse map ZSC are then compared with the exact values. This is an example of a polygon with a vertex at infinity. 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 When I ran TEST1 on the VAX 750 under the UNIX F77 compiler, this was the output. (See p. 94 of the paper referenced in the Introduction for an explanation of the oscillatory error figures.) THE SUM-OF-SQUARES ERROR AT STEP 1 IS .4628e+00 THE SUM-OF-SQUARES ERROR AT STEP 2 IS .4647e+00 THE SUM-OF-SQUARES ERROR AT STEP 3 IS .4616e+00 THE SUM-OF-SQUARES ERROR AT STEP 4 IS .4656e+00 THE SUM-OF-SQUARES ERROR AT STEP 5 IS .3467e+00 THE SUM-OF-SQUARES ERROR AT STEP 6 IS .2426e-01 THE SUM-OF-SQUARES ERROR AT STEP 7 IS .1397e+00 THE SUM-OF-SQUARES ERROR AT STEP 8 IS .1097e-02 THE SUM-OF-SQUARES ERROR AT STEP 9 IS .1231e-04 THE SUM-OF-SQUARES ERROR AT STEP 10 IS .2705e-05 THE SUM-OF-SQUARES ERROR AT STEP 11 IS .1665e-06 THE SUM-OF-SQUARES ERROR AT STEP 12 IS .3947e-04 THE SUM-OF-SQUARES ERROR AT STEP 13 IS .6966e-08 THE SUM-OF-SQUARES ERROR AT STEP 14 IS .1194e-04 THE SUM-OF-SQUARES ERROR AT STEP 15 IS .6152e-12 PARAMETERS DEFINING MAP: (N = 4) (NPTSQ = 5) K W(K) TH(K)/PI BETAM(K) Z(K) --- ---- -------- -------- ---- 1 ( .000, 1.000) .49999988 1.00000 ( .00000031, 1.00000000) 2 ( .000, .000) .99999976 -.50000 (-1.00000000, .00000063) 3 INFINITY 1.50000000 -2.00000 ( .00000001,-1.00000000) 4 ( .000, .000) 2.00000000 -.50000 ( 1.00000000, .00000000) WC = ( .00000000e+00, .14142135e+01) C = ( -.14142140e+01, .00000000e+00) ERREST: .5235e-06 Z,W,WEX,ERR: ( -.300, .210) ( .196, 1.095) ( .196, 1.095) .1229e-06 Z,W,WEX,ERR: ( .000, .270) ( .000, 1.153) ( .000, 1.153) .1199e-06 Z,W,WEX,ERR: ( .300, .330) ( -.133, 1.048) ( -.133, 1.048) .1406e-06 Z,W,WEX,ERR: ( .600, .390) ( -.126, .887) ( -.126, .887) .1490e-06 W,Z,ZEX,ERR: (-1.000, 1.414) ( .383, -.295) ( .383, -.295) .2980e-07 W,Z,ZEX,ERR: ( .000, 1.732) ( .000, -.172) ( .000, -.172) .4479e-07 W,Z,ZEX,ERR: ( 1.000, 2.000) ( -.245, -.383) ( -.245, -.383) .0000e+00 W,Z,ZEX,ERR: ( 2.000, 2.236) ( -.298, -.560) ( -.298, -.560) .5960e-07 W,Z,ZEX,ERR: ( 3.000, 2.449) ( -.296, -.679) ( -.296, -.679) .1333e-06 W,Z,ZEX,ERR: ( 4.000, 2.646) ( -.276, -.757) ( -.276, -.757) .6957e-06 TEST2: call RESIST to compute a resistance This program computes the resistance (conformal module) of an L-shaped hexagon, assuming that a voltage difference is applied between the top and the bottom edges. The exact answer is SQRT(3) = 1.7320508. The subroutine RESIST is used for the computation. RESIST is a driver for SCSOLV, etc. that composes two Schwarz-Christoffel maps in such a way as to map a polygon conformally onto a rectangle and thereby determine the dimen- sions of a rectangle which is electrically equivalent to a given polygonal resistor. RESIST is not documented here, but the program is normally included in the SCPACK tape and the comments there should be sufficient documentation. 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 Here is the output from a run of the above program: NDIG = 1: R = 1.732014775 ERREST = .145e-02 NDIG = 3: R = 1.732097626 ERREST = .457e-03 NDIG = 5: R = 1.732051611 ERREST = .176e-04 TEST3: generate a rectilinear grid The third test program, like RESIST, composes two Schwarz- Christoffel maps to find the conformal map from a polygon with four distinguished vertices onto a rectangle. A rectilinear grid of points in the rectangle is then transferred by ZSC and WSC to the polygon. 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 For a test I ran TEST3 on the VAX 750 and gave the following input data at run time: vertices? (terminate with Re w(k) > 100 ) 0 1 0 0 2 0 3 -1 3.5 -.5 2.5 1 999 0 image of zero in the polygon? 2 .5 four distinguished vertices? 1 2 4 5 no. of streamlines, equipotential lines? 7 20 Here is the terminal output. The slow convergence in this example is caused by the very uneven distribution of the prevertices on the unit circle, which is visible in the values of TH(K) in the table. On the Cyber 175 this run took about 15 secs., most of it spent mapping points from rectangle to polygon after the parameter problem had been solved. THE SUM-OF-SQUARES ERROR AT STEP 1 IS .6139e+01 THE SUM-OF-SQUARES ERROR AT STEP 2 IS .6130e+01 THE SUM-OF-SQUARES ERROR AT STEP 3 IS .6135e+01 THE SUM-OF-SQUARES ERROR AT STEP 4 IS .6135e+01 THE SUM-OF-SQUARES ERROR AT STEP 5 IS .6137e+01 THE SUM-OF-SQUARES ERROR AT STEP 6 IS .6140e+01 THE SUM-OF-SQUARES ERROR AT STEP 7 IS .6703e+01 THE SUM-OF-SQUARES ERROR AT STEP 8 IS .5241e+01 THE SUM-OF-SQUARES ERROR AT STEP 9 IS .1980e+01 THE SUM-OF-SQUARES ERROR AT STEP 10 IS .1048e+01 THE SUM-OF-SQUARES ERROR AT STEP 11 IS .1729e+01 THE SUM-OF-SQUARES ERROR AT STEP 12 IS .6516e+00 THE SUM-OF-SQUARES ERROR AT STEP 13 IS .6521e+00 THE SUM-OF-SQUARES ERROR AT STEP 14 IS .6360e+00 THE SUM-OF-SQUARES ERROR AT STEP 15 IS .2172e+00 THE SUM-OF-SQUARES ERROR AT STEP 16 IS .2987e+00 THE SUM-OF-SQUARES ERROR AT STEP 17 IS .1774e+00 THE SUM-OF-SQUARES ERROR AT STEP 18 IS .1560e+00 THE SUM-OF-SQUARES ERROR AT STEP 19 IS .1109e+00 THE SUM-OF-SQUARES ERROR AT STEP 20 IS .8883e-01 THE SUM-OF-SQUARES ERROR AT STEP 21 IS .8904e-01 THE SUM-OF-SQUARES ERROR AT STEP 22 IS .5549e-01 THE SUM-OF-SQUARES ERROR AT STEP 23 IS .5532e-01 THE SUM-OF-SQUARES ERROR AT STEP 24 IS .5532e-01 THE SUM-OF-SQUARES ERROR AT STEP 25 IS .3392e-01 THE SUM-OF-SQUARES ERROR AT STEP 26 IS .9875e-02 THE SUM-OF-SQUARES ERROR AT STEP 27 IS .4590e-03 THE SUM-OF-SQUARES ERROR AT STEP 28 IS .1225e-04 THE SUM-OF-SQUARES ERROR AT STEP 29 IS .8053e-06 THE SUM-OF-SQUARES ERROR AT STEP 30 IS .3812e-07 THE SUM-OF-SQUARES ERROR AT STEP 31 IS .9499e-04 THE SUM-OF-SQUARES ERROR AT STEP 32 IS .1182e-08 PARAMETERS DEFINING MAP: (N = 6) (NPTSQ = 3) K W(K) TH(K)/PI BETAM(K) Z(K) --- ---- -------- -------- ---- 1 ( .000, 1.000) .84163296 -.50000 ( -.87876660, .47725174) 2 ( .000, .000) .84644878 -.50000 ( -.88588625, .46390247) 3 ( 2.000, .000) 1.38581097 .25000 ( -.35109055, -.93634152) 4 ( 3.000,-1.000) 1.65438819 -.50000 ( .46623042, -.88466334) 5 ( 3.500, -.500) 1.65855527 -.43717 ( .47777134, -.87848425) 6 ( 2.500, 1.000) 2.00000000 -.31283 ( 1.00000000, .00000000) WC = ( .20000000e+01, .50000000e+00) C = ( .56718796e+00, .31701493e+00) ERREST: .9236e-03 VERTICES OF IMAGE RECTANGLE: ( .000000e+00, .100000e+01) ( .000000e+00, .000000e+00) ( .400924e+01, -.154972e-05) ( .400924e+01, .999998e+00) FINISHED ROW 0/ 8 OF GRID POINTS FINISHED ROW 1/ 8 OF GRID POINTS FINISHED ROW 2/ 8 OF GRID POINTS FINISHED ROW 3/ 8 OF GRID POINTS FINISHED ROW 4/ 8 OF GRID POINTS FINISHED ROW 5/ 8 OF GRID POINTS FINISHED ROW 6/ 8 OF GRID POINTS FINISHED ROW 7/ 8 OF GRID POINTS FINISHED ROW 8/ 8 OF GRID POINTS And here is the resulting plot: