# to unbundle, sh this file (in an empty directory) echo readme 1>&2 sed >readme <<'//GO.SYSIN DD readme' 's/^-//' -From ttacs1.ttu.edu!kkent Tue Dec 12 21:14:47 1989 -Date: 1 Dec 89 16:42:00 CST -From: "Pearce, Kent" -To: "ehg" - -Introduction: - -This package GEARLIKE is a collection of programs designed to -solve the accessory parameter problem for constructing the -conformal mapping from the open unit disk D to a gearlike image -domain. In this context a gearlike domain is a simply connected -domain which contains 0 as an interior point and for which the -boundary is composed of arcs of circles centered at 0 and -segments of radial lines emanating from 0. The mapping function -generated by the package is normalized so that 0 is mapped to 0 -and the point 1 on the boundary of D is mapped to a user -specified vertex on the gearlike domain. - -It is known that if G is a gearlike domain with vertices Wk, 1 <= -k <= N, and if F is a conformal mapping which takes D onto G such -that F(0) = 0, then F must satisfy the following differential -equation - - zF'(z) n -Bk - -------- = prod (1 - z/Zk) (1) - F(z) k=1 - -where Bk * pi is the exterior angle at the vertex Wk if Wk is -finite and Bk is 1 if the vertex Wk is infinite. Also, in (1) -the points Zk are the prevertices of G under the map F, i.e., for -each k we have Zk belongs to the boundary of D and F(Zk) = Wk. -Further, the following three conditions must be satisfied - - each Bk is either -1, -1/2, 0, 1/2 or 1 (a) - - N - sum Bk = 0 (b) - k=1 - - N _ - sum Bk * Arg(Zk) = 0 (mod pi). (c) - k=1 - - -Equation (1) can be solved to represent F as - - z - F(z) = C * z * exp [ integral (P(w) - 1)/w dw (2) - 0 - -where the function P in (2) denotes the right hand side of (1). -Thus, from (2) in order to explicitly compute the mapping -function f one has to determine N+1 accessory parameters -- C and -the N prevertices Zk. One of the prevertices can be fixed by -distinguishing one of the vertices as WN and requiring that ZN = -1. The remaining accessory parameters can in general only be -determined numerically. - -GEARLIKE: - -The package GEARLIKE allows a user to specify the number of -vertices (( N )) for a gearlike domain, 2 <= N <= 20. The user -specifies the vertices (( W(k), 1 <= k <= N )) in rectangular -coordinates. The implementation in the program requires that the -last vertex W(N), which we are also calling the distinguished -vertex, be finite. Also, the user specifies the values ((BETA(k) -= - Bk, 1 <= k <= N )), where Bk is as described above, i.e., if -the vertex W(k) is finite, then BETA(k) * pi is the negative of -the exterior angle at W(k) and if the vertex W(k) is infinite, -then BETA(k) = -1. For technical reasons the interior angle at -the vertex W(N-1) can not be a straight angle, i.e., BETA(N-1) -can not be equal to 0. - -If the prevertices (( Z(k), 1 <= k <= N )) are known or if a -reasonable approximation is known, those values can be specified -via their arguments; otherwise a uniform distribution for the -prevertices is taken as an initial guess. Whether or not an -initial estimate of the prevertices is available for use in the -program is passed to the program through a control parameter -IGUESS. The program takes the the constrained arguments of the -prevertices and transforms them to a set of unconstrained -variables for the nonlinear solver NS01A. - -The user specifies the number of nodes (( NPTSQ )) for the Gauss- -Jacobi quadrature approximations to the integral in (2). The -user also specifies an error tolerance (( TOL )) for terminating -the nonlinear solver NS01A used to approximate the unknown -accessory parameters. Finally, the user specifies an input value -IPRINT to control the level of output from the program. - -Generally speaking the parameters and subroutines (with argument -lists) were modeled to be interchangeable with or compatible with -the parameters and subroutines in Trefethen's SCPACK for -computing Schwarz-Christoffel mappings. - - -IMPLEMENTATION: - -The main calling program GLSOLV is contained in the file -GLDRV.F. All of the subsidary files for supporting GLSOLV are -contained in GLDRV.F except for: - - - NS01A: The nonlinear solver called for solving - the accessory parameter problem. - - ODE & Routines: The differential equations package by - Shampine & Gordon. ****** THE ODE - ROUTINES ARE ONLY USED FOR COMPUTING THE - INVERSE MAPPING AND THAT ONLY AFTER THE - FORWARD MAPPING PROBLEM HAS BEEN SOLVED. - ****** - - LINPACK Routines: NS01A & ODE call several LINPACK library - routines. - - GAUSSJ & Routines: Library routines for computing the nodes - and weights for Gauss-Jacobi quadrature. - - - ----------------------------------------------------------------- - - -PRIMARY SUBROUTINE: QINIT Calls the library routine GAUSSJ to - compute the nodes and weights for - Gauss-Jacobi quadrature - - INPUT VARIABLES: - - N Number of vertices of gearlike - image domain, 2 <= N <= 20 - - BETA Real array containing the negatives - of the exterior angles at the - vertices divided by pi. For - vertices W(k) which lie at infinity - BETA(k) = -1. Array should be - dimensioned at least size [0:N] - ****** RETAINED FOR COMPATABILITY - WITH CALLING SEQUENCE IN SCPACK - ****** - - NPTSQ The number of points to be used in - the Gauss-Jacobi quadrature - approximations to the integral in - (2) above - - - - OUTPUT VARIABLES: - - QWORK Real array containing the - quadrature nodes and weights. The - dimension of QWORK should be at - least NPTSQ * 7, but no greater - than 100 - - -PRIMARY SUBROUTINE: GLSOLV Calls the nonlinear equation solver - NS01A to solve the accessory - parameter problem for the gearlike - image domain - - INPUT VARIABLES: - - IPRINT Controls the level of output from - the program -2,-1, 0, 1 - (increasing output) - - IGUESS Specifies whether an initial guess - for the prevertices Z(k) is - available: 1 if an estimate is - supplied; 0 otherwise - - TOL Error tolerance for exiting from - the nonlinear equations solver. - Typically taken to be 10^-(NPTSQ+1) - times the mean modulus of the - vertices. - - N Number of vertices of gearlike - image domain, 2 <= N <= 20 - - WC 0 : the image of 0 under the - gearlike map. ****** RETAINED FOR - COMPATIBILITY WITH THE CALLING - SEQUENCE IN SCPACK ****** - - W Complex array containing the - vertices of the gearlike image - domain, specified in rectangular - coordinate form. For vertices - which lie at infinity, any - coordinates can be specified. - Array should be dimensioned at - least size [0:N] - - BETA Real array containing the negatives - of the exterior angles at the - vertexs divided by pi. For - vertices W(k) which lie at infinity - BETA(k) = -1. Array should be - dimensioned at least size [0:N] - - NPTSQ The number of points to be used in - the Gauss-Jacobi quadrature - approximations to the integral in - (2) above - - QWORK Real array containing the - quadrature nodes and weights. The - dimension of QWORK should be at - least NPTSQ * 7, but no greater - than 100 - - INPUT/OUTPUT VARIABLES: - - Z Complex array containing the - prevertices for the gearlike image - domain. If IGUESS = 1, then the - values in Z will be used as initial - inputs. For all cases the output - values will contain the prevertices - as constructed by the program. - Array should be dimensioned at - least size [0:N] - - - OUPUT VARIABLES: - - ERREST Error estimate of deviation of - boundary of constructed mapping to - boundary of exact mapping. - Estimate constructed by comparing - computed values at finite vertices - with exact values at finite - vertices - - C Complex scale factor in formula (2) - above. |C| is the conformal - mapping radius for the constructed - gearlike mapping - -PRIMARY SUBROUTINE: GLFUN NS01A calls an external routine - GLFUN which calcuates the nonlinear - function values for which the - solution will be the accessory - parameter set. The algorithm in - GLFUN is lengthy because (1) the - generality allowed the user for the - choice of W(N); (2) the multipli- - city of choices of argument for - polar coordinate representation of - complex numbers. - - - - INPUT VARIABLES: - - NDIM The number of unconstrained - parameters in the accessory - parameter problem. NDIM = N-1 - - Y The transformed, unconstrained - parameters. - - OUTPUT VARIABLES: - - FVAL The values of the nonlinear - function constructed in the - algorithm in GLFUN. - - -PRIMARY SUBROUTINE: WMAP Computes the mapping values for WW - = F(ZZ) for the gearlike mapping F - being constructed - - INPUT VARIABLES: - - ZZ Input point at which the mapping - value WW = F(ZZ) is to be computed - - Z0 Nearby point at which the mapping - value W0 = F(Z0) is known (either - as a vertex value, or the value 0 - at 0 or a previously computed - value) - - W0 The known value W0 = F(Z0) - - K0 0 if the input point is not a - prevertex; otherwise the index of - the prevertex Z0, i.e., Z0 = Z(K0) - - OTHER VARIABLES: - - N,C,Z,BETA, As in GLSOLV - NPTSQ,QWORK - - -PRIMARY SUBROUTINE: WMAPP Computes the derivative mapping - values for WW = F'(ZZ) for the - gearlike mapping f being - constructed - - INPUT VARIABLES: - - ZZ Input point at which the mapping - value WW = F(ZZ) is to be computed - - Z0 Nearby point at which the mapping - value W0 = F(Z0) is known (either - as a vertex value, or the value 0 - at 0 or a previously computed - value) - - W0 The known value W0 = F(Z0) - - K0 0 if the input point is not a - prevertex; otherwise the index of - the prevertex Z0, i.e., Z0 = Z(K0) - - OTHER VARIABLES: - - N,C,Z,BETA, As in GLSOLV - NPTSQ,QWORK - - -PRIMARY SUBROUTINE: ZMAP Computes the mapping values for ZZ - = INVERSE[F(WW)] for the inverse of - the gearlike mapping F being - constructed - - INPUT VARIABLES: - - WW Input point at which the mapping - value ZZ = INVERSE[F(WW)] is to be - computed - - IGUESS Specifies whether an initial guess - is available for the value of ZZ = - INVERSE[F(WW)]: IGUESS >< 0 is an - initial guess is supplied; IGUESS - = 0 otherwise - - ZGUESS Initial guess for the value of - INVERSE[F(WW)]. Ignored if IGUESS - = 0 - - Z0 A know value: Z0 = INVERSE[F(W0)] - for some nearby point W0, i.e., a - point at which W0 = F(Z0) is known, - finite and near WW - - W0 The value (finite) at Z0 - - K0 0 if the input point Z0 is not a - prevertex; otherwise the index of - the prevertex Z0, i.e., Z0 = Z(K0) - - TOL Desired accuracy for the solution - ZZ = INVERSE[F(WW)] - - - INPUT/OUTPUT VARIABLES: - - IERR Input code for expressing- - suppressing printed error messages - IERR = 0 print messages - IERR >< 0 suppress messages - Output code denoting whether the - subroutine terminated successfully - IERR = 0 successful run - IERR >< 0 failed to converge - - OTHER VARIABLES: - - N,C,Z,BETA, As in GLSOLV - NPTSQ,QWORK - - -SECONDARY ROUTINES: - - CALLED FROM GLSOLV: - - CHECK Checks the input data for - consistency in conforming to the - problem requirements. ***** CHECK - WILL TERMINATE THE PROGRAM IF A - DATA ERROR IS DETECTED. ***** - - OUTPUT Prints a formated summary of the - the prevertices and vertices as - found by NS01A. - - TEST Calculates an error estimate for - the prevertices found by NS01A. - - - CALLED FROM GLFUN: - - YZTRAN Transforms the unconstrained - parameters back to arguments for - the prevertices Z on the unit - circle - - THETA Calculates the arguments of the - prevertices Z on the unit circle, 0 - <= THETA(k) <= THETA(K+1) <= 2 * pi - - RELARG Calcuates the signed change in - argument for the vertices W - - - - - - - CALLED FROM WMAP: - - ZQUAD Calculates the quadrature - approximation to the integral in - (2) - - - OTHER SUBROUTINES: - - - ZODE External subroutine called by ODE - which represents the differential - equation DZZ/DWW = 1/F'(ZZ) - - ---------------------------------------------------------------- - -Secondary programs: - - GL.F & GL.TXT: The program GL.F with its accompanying - text input file GL.TXT is sample program - for specifying the input parameters for - a given gearlike domain. - - PLTTER & PLTTR2: After solving the accessory parameter - problem, the program GL.F calls a - plotting subroutine PLTTER, contained in - the file DISSPLA.F, for graphing the - gearlike domain and accompanying levels - sets using DISSPLA graphics. The plot - calls in PLTTER to DISSPLA should easily - be changeable to calls for an alternate - graphics language. The routine PLTTR2 - will plot the inverse map. - - - //GO.SYSIN DD readme echo gldrv.f 1>&2 sed >gldrv.f <<'//GO.SYSIN DD gldrv.f' 's/^-//' C********************************************************** C* QINIT GEARLIKE PRIMARY SUBROUTINE QUINT *** C********************************************************** SUBROUTINE QINIT (N,BETA,NPTSQ,QW) C C COMPUTES NODES AND WEIGHTS FOR GAUSS-JACOBI QUADRATURE C C CALLING SEQUENCE PARAMETERS: SEE COMMENTS IN GLSOLV C C THE WORK ARRAY QW MUST BE DIMENSIONED AT LEAST NPTSQ * 7 C IT IS DIVIDED UP INTO 7 VECTORES OF LENGTH NPTSQ: THE FIRST 6 C CONTAIN QUADRATURE NODES AND WIEGHTS ON OUTPUT. THE FINAL C ONE IS SCRATCH VECTOR NEEDED BY GAUSSJ. C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION BETA(0:N),QW(1),BB(3) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA BB/0.D0,-0.5D0,0.5D0/ C C FOR THE VALUES 0, -1/2 & 1/2 COMPUTE NODES AND WEIGHTS FOR C ONE-SIDED GAUSS-JACOBI QUADRATURE ALONG A CURVE BEGINNING AT Z(K): C ISCR = NPTSQ*6 + 1 DO 100 K = 1,3 INODES = NPTSQ*2*(K-1) + 1 IWTS = INODES + NPTSQ CALL GAUSSJ (NPTSQ,D0,BB(K), & QW(ISCR),QW(INODES),QW(IWTS)) 100 CONTINUE C RETURN END C********************************************************** C* GLSOLV GEARLIKE PRIMARY SUBROUTINE GLSOLV *** C********************************************************** SUBROUTINE GLSOLV (IPRINT,IGUESS,TOL,ERREST, & N,C,Z,WC,W,BETA,NPTSQ,QW) C C C THIS SUBROUTINE COMPUTES THE ACCESSORY PARAMENTERS C AND C Z(K) FOR THE GEAR-LIKE TRANSFORMATION WHICH SENDS THE UNIT C DISK TO THE INTERIOR OF THE GEAR-LIKE DOMAIN WITH VERTICES C W(1) . . . W(N). THIS MAPPING IS OF THE FORM: C C Z C W = C * Z * EXP( INT ( (P(Z) - 1) / Z ) DZ ). C 0 C C WHERE P(Z) IS OF THE FORM C C N C P(Z) = PROD (1 - Z / Z(K) ) ** BETA(K) C K=1 C C THE IMAGE GEAR-LIKE DOMAIN MAY BE UNBOUNDED. THE VALUES BETA(K) C MUST ALL BE -1, -1/2, 0, 1/2 OR 1. W(N) MUST BE FINITE. C WE NORMALIZE BY THE CONDITIONS: C C Z(N) = 1 C C IN ORDER FOR THE IMAGE DOMAIN TO BE GEAR-LIKE WE REQUIRE C C N C SUM BETA(K) = 0 C K=1 C C N C SUM BETA(K) * ARG(Z(K)) = 0 (MOD PI) C K=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 C 0 (INPUT) C C TOL DESIRED ACCURACY IN SOLUTION OF NONLINEAR SYSTEM C INPUT. RECOMMENED VALUE: 10.0**(NPTSQ-1) * TYPICAL C SIZE OF VERTICES OF W(K) C C ERREST ESTIMATED ERROR IN SOLUTION (OUTPUT). MORE C PRECISELY, ERREST IS AN APPROZIMATE BOUND FOR HOW FAR C THE TRUE VERTICES OF THE IMAGE DOMAIN MAY BE FROM C THOSE COMPUTED BY NUMERICAL INTEGRATION USING THE C NUMERICALLY DETERMINED PREVERTICES Z(K). C C N NUMBER OF VERTICES OF THE IMAGE GEAR-LIKE DOMAIN C (INPUT). MUST BE <= 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 GUES IS C BEING SUPPLIED IT SHOULD BE IN Z ON INPUT, WITH C Z(N) = 1. IN ANY CASE THE CORRECT PREVERTICES WILL C BE IN Z ON OUTPUT. C C WC RETAINED FOR COMPATIBLITY WITH CALLING SEQUENCE IN C SCHWARZ-CHRISTOFEL DRIVER C C W COMPLEX ARRAY OF VERTICES ON THE IMAGE GEAR-LIKE C DOMAIN (INPUT). DIMENSION AT LEAST N. IT IS A GOOD C IDEA 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 BETA, BELOW. EACH CONNECTED BOUNDARY COMPONENT C MUST INCLUDE AT LEAST ONE VERTEX W(K), EVEN IF IT HAS C TO BE A DEGENERATE VERTEX WITH BETA(K) = 0. W(N) C MUST BE FINITE. C C BETA REAL ARRAY WITH BETA(K) THE EXTERNAL ANGLE IN THE C GEARLIKE DOMAIN AT THE VERTEX K, DIVIDED BY MINUS C PI FOR FINITE VERTICES AND SET TO BE -1 AT INFINITE C VERTICES (INPUT). EACH BETA(K) MUST BE -1, -1/2, 0, C 1/2 OR 1. EXAMPLES: AT THE TIP OF A SLIT (CIRCULAR C OR RADIAL) BETA(K) IS ALWAYS 1. AT THE JUNCITON C OF A RADIAL SEGMENT AND A CIRCULAR ARC BETA(K) IS C EITHER -1/2 OR 1/2. W(K) LIES AT INFINITY IF AND C ONLY IF BETA(K) = -1. BETA(N-1) CAN NOT BE 0, I.E., C THE INTERNAL ANGLE AT W(N-1) CAN NOT BE A STRAIGHT C ANGLE. IF IT IS, THEN W(N-1) CAN BE REMOVED FROM THE C THE SET OF VERTICES OF THE GEARLIKE DOMAIN. 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 FILLS THE VECTOR QW. C C QW REAL QUADRATURE WORK ARRAY (INPUT). DIMENSION AT C LEAST NPTSQ * 7 BUT NO MORE THAN 100. BEFORE CALLING C GLSOLV QW MUST HAVE BEEN FILLED BY THE SUBROUTINE C QINIT. C C THE PROBLEM IS SOLVED BY FINDING THE SOLUTIONS TO A SYSTEM OF N-1 C NONLINEAR EQUATIONS IN THE N-1 UNKNOWNS Y(1) . . . Y(N-1), WHICH ARE C RELATED TO THE POINTS Z(K) BY THE FORMULA: C C Y(K) = LOG ((TH(K)-TH(K-1))/(TH(K+1)-TH(K))) C C WHERE TH(K) DENOTES THE ARGUMENT OF Z(K). SUBROUINTE GLFUN DEFINES C THIS SYSTEM OF EQUATIONS. THE ORIGINAL PROBLEM IS SUBJECT TO C CONSTRAINTS TH(K) < TH(K+1), BUT THESE VANISH IN THE TRANSFORMATION C FROM Z TO Y. C C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) COMMON /PARAM1/ W2(0:20),Z2(0:20),C2,WC2,BETA2(0:20), & QW2(100),TOL2,NPTSQ2,JPOLE2 DIMENSION Z(0:N),W(0:N),BETA(0:N),QW(1) DIMENSION AJINV(20,20),SCR(900),FVAL(20),Y(20) EXTERNAL GLFUN DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ PI = D4*DATAN(D1) C C SET VALUE OF Z(N) C Z(N) = C1 C C IDENTIFY (Z(0),W(0),BETA(0)) WITH (Z(N),W(N),BETA(N)) C Z(0) = Z(N) W(0) = W(N) BETA(0) = BETA(N) C C CHECK INPUT DATA: C CALL CHECK (TOL,N,W,BETA) C C CHECK INITIAL GUESS ON LOCATION OF PARAMETERS Z(K) C IF (IGUESS.EQ.0) THEN C C NO VERTICES SUPPLIED C INITIAL ESTIMATES WILL BE GENERATED C ARG = D2*PI/DFLOAT(N) DO 100 K = 1,N-1 TH = ARG*DFLOAT(K) Z(K) = CDEXP(CI*TH) 100 CONTINUE C ENDIF C C TRANSFORM ARGUMENT Z(K) VALUES TO UNCONSTRAINED Y(K) VALUES C DO 200 K = 1,N-1 TMP1 = DIMAG(CDLOG(Z(K+1)/Z(K))) IF (TMP1.LT.D0) TMP1 = TMP1 + D2*PI TMP2 = DIMAG(CDLOG(Z(K)/Z(K-1))) IF (TMP2.LT.D0) TMP2 = TMP2 + D2*PI Y(K) = DLOG(TMP2) - DLOG(TMP1) 200 CONTINUE C C CHECK TO SEE IF W(N-1) OR W(N) LIES AT THE END OF A RADIAL SEGMENT C AND THE INTERIOR ANGLE IS PI/2 AT THAT END C JPOLE2 = 0 IF ((BETA(N-1).EQ.-DH).OR.(BETA(N).EQ.-DH)) THEN IF (BETA(N-1).EQ.-DH) THEN JJ = 0 JPOLE2 = 1 ELSE JJ = 1 JPOLE2 = 2 ENDIF DO 250 J = JJ,N-2 IF (BETA(J).EQ.-D1) GOTO 275 IF (BETA(J).NE.D0) THEN JPOLE2 = 0 GOTO 275 ENDIF 250 CONTINUE JPOLE2 = 0 275 CONTINUE ENDIF C C NS01A CONTROL PARAMETERS: C DSTEP = 1.D-8 DMAX = 20.D0 MAXFUN = (N-1) * 15 C C COPY INPUT DATA TO /PARAM1/ FOR GLFUN: C (THIS IS NECESSARY BECUASE NS01A REQUIRES A FIXED CALLING C SEQUENCE IN SUBROUTINE GLFUN.) C NPTSQ2 = NPTSQ TOL2 = TOL WC2 = WC Z2(0) = Z(0) Z2(N) = Z(N) C DO 300 K = 0,N W2(K) = W(K) BETA2(K) = BETA(K) 300 CONTINUE C DO 400 K = 1,NPTSQ*7 QW2(K) = QW(K) 400 CONTINUE C C SOLVE NONLINEAR SYSTEM WITH NS01A: C CALL NS01A (N-1,Y,FVAL,AJINV,DSTEP,DMAX, & TOL,MAXFUN,IPRINT,SCR,GLFUN) C C COPY OUTPUT DATA FROM /PARAM1/: C C = C2 DO 500 K = 1,N-1 Z(K) = Z2(K) 500 CONTINUE C C PRINT RESULTS AND TEST ACCURACY C IF (IPRINT.GE.0) CALL OUTPUT (N,C,Z,WC,W,BETA,NPTSQ,QW) CALL TEST (ERREST,N,C,Z,WC,W,BETA,NPTSQ,QW) IF (IPRINT.GE.-1) WRITE(6,620) ERREST IF (IPRINT.GE.-1) WRITE(20,620) ERREST C 620 FORMAT (/' ERREST = ',E12.4/) C RETURN END C********************************************************** C* CHECK GEARLIKE SUBROUTINE CHECK *** C********************************************************** SUBROUTINE CHECK (EPS,N,W,BETA) C C CHECKS GEOMETRY OF PROBLEM TO MAKE SURE IS A FORM USABLE C BY GLSOLV. C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION W(0:N),BETA(0:N),BB(5) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA BB/-1.D0,-0.5D0,0.D0,0.5D0,1.D0/ C DO 100 K = 1,N DO 150 J = 1,5 IF (DABS(BETA(K)-BB(J)).LT.EPS) THEN BETA(K) = BB(J) GOTO 100 ENDIF 150 CONTINUE WRITE(6,120) WRITE(20,120) STOP 100 CONTINUE C SUM = D0 DO 200 K = 1,N SUM = SUM + BETA(K) 200 CONTINUE IF (DABS(SUM).GT.EPS) THEN WRITE(6,220) WRITE(20,220) STOP ENDIF C IF (BETA(N).EQ.-D1) THEN WRITE(6,320) WRITE(20,320) STOP ENDIF C IF (BETA(N-1).EQ.D0) THEN WRITE(6,330) WRITE(20,330) STOP ENDIF C DO 400 K = 1,N IF ((BETA(K).EQ.-D1).OR.(BETA(K-1).EQ.-D1)) GOTO 400 WW = CDLOG(W(K)/W(K-1)) R = DIMAG(WW*WW) IF (DABS(R).GT.EPS) THEN WRITE(6,420) WRITE(20,420) STOP ENDIF 400 CONTINUE C IF (N.LT.3) THEN WRITE(6,510) WRITE(20,510) STOP ENDIF C IF (N.GT.20) THEN WRITE(6,520) WRITE(20,520) STOP ENDIF C 120 FORMAT (/' ERROR IN CHECK: ANGLES MUST BE -1,-1/2,0,1/2 OR 1'/) 220 FORMAT (/' ERROR IN CHECK: ANGLES DO NOT ADD UP TO 0'/) 320 FORMAT (/' ERROR IN CHECK: W(N) MUST BE FINITE'/) 330 FORMAT (/' ERROR IN CHECK: BETA(N-1) MUST NOT BE 0'/) 420 FORMAT (/' ERROR IN CHECK: IMAGE DOMAIN NOT GEAR-LIKE'/) 510 FORMAT (/' ERROR IN CHECK: N MUST BE NO LESS THAN 3'/) 520 FORMAT (/' ERROR IN CHECK: N MUST BE NO MORE THAN 20'/) C RETURN END C********************************************************** C* OUTPUT GEARLIKE SUBROUTINE OUTPUT *** C********************************************************** SUBROUTINE OUTPUT (N,C,Z,WC,W,BETA,NPTSQ,QW) C C PRINTS A TABLE OF K, TH(K)/PI, Z(K), BETA(K), WMAP(K) & LOG(WMAP(K)) C ALSO PRINTS THE CONSTANTS N, NPTSQ & C. C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),W(0:N),BETA(0:N),QW(1),TH(0:20) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ PI = D4*DATAN(D1) C CALL THETA (N,Z,TH) C WRITE(6,120) N, NPTSQ WRITE(20,120) N, NPTSQ DO 100 K = 1,N IF (BETA(K).GT.-D1) THEN WK = WMAP (Z(K),K,CO,WC,0,N,C,Z,BETA,NPTSQ,QW) WRITE(6,121) K,TH(K)/PI,Z(K),BETA(K),WK, & CDABS(WK),DIMAG(CDLOG(WK))/PI WRITE(20,121) K,TH(K)/PI,Z(K),BETA(K),WK, & CDABS(WK),DIMAG(CDLOG(WK))/PI ELSE WRITE(6,122) K,TH(K)/PI,Z(K),BETA(K) WRITE(20,122) K,TH(K)/PI,Z(K),BETA(K) ENDIF 100 CONTINUE WRITE(6,123) C WRITE(20,123) C C 120 FORMAT (/' PARAMETERS DEFINING MAP:',15X,'(N =', & I3,')',6X,'(NPTSQ =',I3,')'// & ' K',5X,'TH(K)/PI',15X,'Z(K)',16X,'BETA(K)', & 13X,'W(K)',21X,'POLAR(W(K))'/ & ' ---',4X,'--------',15X,'----',16X,'-------', & 13X,'----',21X,'-----------'/) 121 FORMAT (I3,F14.8,' (',F11.8,',',F11.8,')',F12.5, & 3X,'(',F11.8,',',F11.8,')',3X,'(',F11.8,',',F11.8,')') 122 FORMAT (I3,F14.8,' (',F11.8,',',F11.8,')',F12.5, & ' INFINITY ') 123 FORMAT (/' C = (',E15.8,',',E15.8,')') C RETURN END C********************************************************** C* TEST GEARLIKE SUBROUTINE TEST *** C********************************************************** SUBROUTINE TEST (ERREST,N,C,Z,WC,W,BETA,NPTSQ,QW) C C TEST THE COMPUTED MAP FOR ACCURACY. C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),W(0:N),BETA(0:N),QW(1) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ PI = D4*DATAN(D1) C C TEST LENGTH OF RADDI: C ERREST = D0 DO 100 K = 1,N IF (BETA(K).GT.-D1) THEN WK = WMAP (Z(K),K,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) RADE = CDABS(W(K)-WK) ELSE IF (K.GT.1) THEN WKM1=WMAP(Z(K-1),K-1,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) ELSE WKM1 = WMAP (Z(N),N,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) ENDIF WKP1 = WMAP (Z(K+1),K+1,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) RAD1 = CDABS(W(K-1)-W(K+1)) RAD2 = CDABS(WKM1 - WKP1) RADE = DABS(RAD1-RAD2) ENDIF ERREST = DMAX1(ERREST,RADE) 100 CONTINUE C RETURN END C********************************************************** C* GLFUN GEARLIKE PRIMARY SUBROUTINE GLFUN *** C********************************************************** SUBROUTINE GLFUN (NDIM,Y,FVAL) C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) COMMON /PARAM1/ W(0:20),Z(0:20),C,WC,BETA(0:20), & QW(100),TOL,NPTSQ,JPOLE DIMENSION TH(0:20) DIMENSION FVAL(NDIM),Y(NDIM) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ PI = D4*DATAN(D1) C C SHIFT DIMENSION LEVEL ON NUMBER OF EQUATIONS C N = NDIM + 1 C C TRANSFORM Y(K) TO Z(K) C CALL YZTRAN (N,Y,Z) C C CALCULATE ARGUMENTS OF PREVERTICES Z(K) C CALL THETA (N,Z,TH) C C SET: COMPUTE INTEGRAL FROM 0 TO Z(N): C ZTEMP = ZQUAD (C0,0,Z(N),N,N,Z,BETA,NPTSQ,QW) C = W(N) / CDEXP(ZTEMP) C IF (BETA(1).NE.-D1) THEN IPOLE = 0 DW = CDABS(W(1)) - CDABS(W(0)) IF (DABS(DW).GT.TOL) THEN IX = 1 ELSE IX = 0 ENDIF ENDIF C I = 0 K = 0 KPOLE = 0 C DO WHILE (K.LT.N-2) I = I+1 IF (BETA(I).EQ.-D1) THEN IPOLE = 1 ELSE IF ((IPOLE.EQ.1).AND.(BETA(I).EQ.D0)) THEN WI = WMAP (Z(I),I,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) K = K+1 FVAL(K) = CDABS(WI) - CDABS(W(I)) ELSE IF (IPOLE.EQ.1) THEN IF (DABS(BETA(I)).EQ.DH) THEN KPOLE = K+1 ELSE KPOLE = 0 ENDIF WI = WMAP (Z(I),I,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) K = K+1 IVMOD = K FVAL(K) = CDABS(WI) - CDABS(W(I)) K = K+1 IVARG = K FVAL(K) = DIMAG(CDLOG(WI))-DIMAG(CDLOG(W(I))) ELSE IF (IX.EQ.1) THEN WI = WMAP (Z(I),I,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) K = K+1 IVMOD = K FVAL(K) = CDABS(WI) - CDABS(W(I)) ELSE CALL RELARG (I,ARGF,ARGW,N,TH,Z,W,BETA,NPTSQ,QW) K = K+1 IVARG = K FVAL(K) = ARGF - ARGW ENDIF C IF (IPOLE.EQ.0) THEN IF (DABS(BETA(I)).EQ.DH) THEN IX = 1-IX ENDIF ELSE IF ((IPOLE.EQ.1).AND.(BETA(I).NE.-D1). & AND.(BETA(I).NE.D0)) THEN IPOLE = 0 IF (DABS(BETA(I)).EQ.DH) THEN IX = 0 ELSE IX = 1 ENDIF ENDIF ENDDO C C CHECK TO SEE IF SOME CONDITION NEEDS TO BE REWRITTEN C IF (JPOLE.GT.0) THEN IF (JPOLE.EQ.1) THEN I = N-1 K = IVARG ELSE I = N IF (IX.EQ.D0) THEN K = IVMOD ELSE K = IVARG ENDIF ENDIF CALL RELARG (I,ARGF,ARGW,N,TH,Z,W,BETA,NPTSQ,QW) FVAL(K) = ARGF - ARGW ELSE IF (KPOLE.GT.0) THEN IF ((IX.EQ.1).AND.(BETA(N-1).EQ.D1)) THEN K = KPOLE + 1 I = N-1 WI = WMAP(Z(I),I,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) FVAL(K) = CDABS(WI) - CDABS(W(I)) ELSE IF ((IX.EQ.1).AND.(BETA(N-1).NE.D1)) THEN CONTINUE ELSE IF ((IX.EQ.0).AND.(BETA(N-1).NE.-D1)) THEN IF (I.EQ.N-1) THEN K = KPOLE I = N ELSE IF (BETA(N-1).EQ.D1) THEN K = IVMOD I = N-1 ELSE K = KPOLE + 1 I = N-1 ENDIF CALL RELARG (I,ARGF,ARGW,N,TH,Z,W,BETA,NPTSQ,QW) FVAL(K) = ARGF - ARGW ENDIF ELSE IF (BETA(N-1).EQ.D1) THEN K = N-2 DW = CDABS(W(N-1)) - CDABS(W(N-2)) IF (DABS(DW).GT.TOL) THEN I = N-1 WI = WMAP (Z(I),I,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) FVAL(K) = CDABS(WI) - CDABS(W(I)) ELSE I = N CALL RELARG (I,ARGF,ARGW,N,TH,Z,W,BETA,NPTSQ,QW) FVAL(K) = ARGF - ARGW ENDIF ENDIF C C CALCULATE CONSTAINT EQUATION C TSUM = D0 DO 160 I = 1,N-1 TSUM = TSUM + BETA(I)*TH(I) 160 CONTINUE C TSUM = DMOD(TSUM,PI) IF (DABS(TSUM).GT.PI/D2) TSUM = TSUM - DSIGN(PI,TSUM) FVAL(N-1) = TSUM C RETURN END C********************************************************** C* RELARG GEARLIKE SUBROUTINE RELARG *** C********************************************************** SUBROUTINE RELARG (I,ARGF,ARGW,N,TH,Z,W,BETA,NPTSQ,QW) C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION TH(0:N),Z(0:N),W(0:N),BETA(0:N),QW(1) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ PI = D4*DATAN(D1) C ARGW = DIMAG(CDLOG(W(I))) - DIMAG(CDLOG(W(I-1))) ARGW = DMOD(ARGW,D2*PI) IF (ARGW.LE.TOL) ARGW = ARGW + D2*PI IF (I.EQ.1) THEN IM1 = N ELSE IM1 = I-1 ENDIF DTH = TH(I) - TH(I-1) Z1 = ZQUAD(C0,0,Z(IM1),IM1,N,Z,BETA,NPTSQ,QW) Z2 = ZQUAD(C0,0,Z(I),I,N,Z,BETA,NPTSQ,QW) ARGF = DIMAG(Z2-Z1) + DTH IF (ARGF.LT.D0) ARGW = ARGW - D2*PI C RETURN END C********************************************************** C* YZTRAN GEARLIKE SUBROUTINE YZTRAN *** C********************************************************** SUBROUTINE YZTRAN (N,Y,Z) C C TRANSFORMS Y(K) TO Z(K). SEE COMMENT IN SUBROUTINE GLSOLV. C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Y(N),Z(0:N) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ PI = D4*DATAN(D1) C DTH = D1 THSUM = DTH DO 100 K = 1,N-1 DTH = DTH / DEXP(Y(K)) THSUM = THSUM + DTH 100 CONTINUE C DTH = D2 * PI / THSUM THSUM = DTH Z(1) = CDEXP(CI*THSUM) DO 200 K = 2,N-1 DTH = DTH / DEXP(Y(K-1)) THSUM = THSUM + DTH Z(K) = CDEXP(CI*THSUM) 200 CONTINUE C RETURN END C********************************************************** C* THETA GEARLIKE SUBROUTINE THETA *** C********************************************************** SUBROUTINE THETA (N,Z,TH) C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),TH(0:N) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ PI = D4*DATAN(D1) C TH(0) = D0 TH(N) = D2*PI C DO 100 K = 1,N-1 TH(K) = DIMAG(CDLOG(Z(K))) IF (TH(K).LE.D0) TH(K) = TH(K) + D2*PI 100 CONTINUE C RETURN END C********************************************************** C* WMAP GEARLIKE PRIMARY SUBROUTINE WMAP *** C********************************************************** FUNCTION WMAP (ZZ,KZZ,Z0,W0,K0,N,C,Z,BETA,NPTSQ,QW) C C COMPUTES THE 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 C AND 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,BETA,NPTSQ,QW AS IN GLSOLV (INPUT) C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),BETA(0:N),QW(1) DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ C ZTEMP = ZQUAD (Z0,K0,ZZ,KZZ,N,Z,BETA,NPTSQ,QW) C C PROTECT AGAINST OVERFLOW NEAR SINGULARITIES C IF (DREAL(ZTEMP).GT.75.0D0) THEN ZTEMP = 75.0D0 * ZTEMP / DREAL(ZTEMP) ENDIF C IF (Z0.EQ.C0) THEN WMAP = C * ZZ * CDEXP(ZTEMP) ELSE WMAP = (W0/Z0) * ZZ * CDEXP(ZTEMP) ENDIF C RETURN END C********************************************************** C* WMAPP GEARLIKE PRIMARY SUBROUTINE WMAPP *** C********************************************************** FUNCTION WMAPP (ZZ,Z0,W0,K0,N,C,Z,BETA,NPTSQ,QW) C C COMPUTES MAP DERIVATIVE DW(ZZ)/DZZ C C CALLING SEQUENCE C C ZZ POINT IN THE DISK AT WHICH DW/DZ IS DESIRED (INPUT) C C Z0 NEARBY POINT IN THE DISK AT WHICH W(Z0) IS KNOWN C AND 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,BETA,NPTSQ,QW AS IN GLSOLV (INPUT) C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),BETA(0:N),QW(1) DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ C IF (ZZ.EQ.C0) THEN WMAPP = C ELSE ZP = ZPROD (ZZ,K0,N,Z,BETA) WZ = WMAP (ZZ,0,Z0,W0,K0,N,C,Z,BETA,NPTSQ,QW) WMAPP = ZP * (WZ/ZZ) ENDIF C RETURN END C********************************************************** C* ZMAP GEARLIKE PRIMARY SUBROUTINE ZMAP *** C********************************************************** FUNCTION ZMAP (WW,IGUESS,ZGUESS,Z0,W0,K0,EPS,IERR, & N,C,Z,BETA,NPTSQ,QW) C C COMPTUES THE INVERSE MAP Z(WW) C C CALLING SEQUENCE C C WW POINT IN GEARLIKE DOMAIN AT WHICH Z(WW) IS C DESIRED (INPUT) C C IGUESS CONTROL PARAMETER WHICH DETERMINES WHETHER C AN INITIAL GUESS HAS BEEN SUPPLIED C IGUESS == 0 ----------> NO INITIAL GUESS C IGUESS >< 0 ----------> INITIAL GUESS SUPPLIED C C ZGUESS INITIAL GUESS FOR Z(WW). IGNORED IF IGUESS = 0. C SHOULD NOT BE A PREVERTEX Z(K). (INPUT) C C Z0 POINT IN THE DISK NEAR Z(WW) AT WHICH W(Z0) IS KNOWN C AND FINITE (INPUT) C C W0 W(Z0) (INPUT) C C K0 K IF Z0 = Z(K) FOR SOME K, OTHERWISE 0 (INPUT) C C TOL DESIRED ACCURACY FOR ANSWER Z(WW) (INPUT) C C IERR INPUT: CONTROL (SUPPRESS ERROR MESSAGES) C IERR == 0 --------> PRINT MESSAGES C IERR >< 0 --------> SUPPRESS MESSAGES C OUTPUT: ERROR CODES C IERR == 0 -------> SUCCESSFUL RUN C IERR >< 0 -------> FAILURE TO CONVERGE C C N,C,Z,BETA,NPTSQ,QW AS IN GLSOLV (INPUT) C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) COMMON /PARAM2/ Z2(0:20),C2,CDWDT,W02,Z02,BETA2(0:20), & QW2(100),K02,N2,NPTSQ2 DIMENSION Z(0:N),BETA(0:N),QW(1) DIMENSION SCR(142),ISCR(5) EXTERNAL ZODE LOGICAL ODECAL DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ C ODECAL = .FALSE. ZI = ZGUESS C C THE VALUE OF IGUESS DETERMINES WHETHER ODE WILL BE CALLED C IF NO INITIAL GUESS IS GIVEN (IGUESS = 0), THEN ODE IS CALLED C ALSO, IF AN INITIAL GUESS IS GIVEN, BUT CONVERGENCE (AT THE NEWTON C STEP) FAILS THEN IGUESS IS SET EQUAL TO 0 AND PROGRAM CONTROL IS C TRANSFERED BACK TO THE BRANCH POINT 100 SO THAT ODE WILL BE CALLED C 100 IF (IGUESS.EQ.0) THEN C C CALL ODE TO GET INITIAL GUESS C C COPY INPUT DATA VALUES TO /PARAM1/ FOR ZODE C (THIS IS NECESSARY BECUASE ODE REQUIRES A FIXED CALLING C SEQUENCE IN SUBROUTINE ZODE.) C N2 = N NPTSQ2 = NPTSQ C2 = C Z02 = Z0 W02 = W0 K02 = K0 C DO 125 K = 0,N Z2(K) = Z(K) BETA2(K) = BETA(K) 125 CONTINUE C DO 150 K = 1,NPTSQ*7 QW2(K) = QW(K) 150 CONTINUE C CDWDT = WW C TSTRT = 0.0D0 TEND = 1.0D0 ZI = C0 IFLAG = -1 RELERR = 0.0D0 ABSERR = 1.0D-4 CDWDT = WW C CALL ODE (ZODE,2,ZI,TSTRT,TEND,RELERR,ABSERR, & IFLAG,SCR,ISCR) IF ((IFLAG.NE.2).AND.(IERR.EQ.0)) THEN C WRITE (6,20) IFLAG WRITE (20,20) IFLAG ENDIF ODECAL = .TRUE. C ENDIF C C NEWTON STEP C DO 200 K = 1,10 WZ = WMAP (ZI,0,Z0,W0,K0,N,C,Z,BETA,NPTSQ,QW) WZP = WMAPP (ZI,Z0,W0,K0,N,C,Z,BETA,NPTSQ,QW) ZI = ZI + (WW - WZ)/ WZP C C CONTROL NEWTON STEP SIZE C IF (CDABS(ZI).GT.1.1D0) ZI = 0.9D0*ZI/CDABS(ZI) C C TEST FOR CONVERGENCE C IF (CDABS(WW-WZ).LT.EPS) GOTO 300 C 200 CONTINUE C C CONVERGENCE HAS NOT OCCURRED. IF ODE WAS NOT CALLED, START OVER C IF (.NOT.ODECAL) THEN IGUESS = 0 GOTO 100 ENDIF C C WRITE ERROR NOTICE C IF (IERR.EQ.0) THEN C WRITE (6,21) WRITE (20,21) ENDIF C 300 ZMAP = ZI C 20 FORMAT (/4X,'Nonstandard return form ODE: IFLAG = ',I2/) 21 FORMAT (/4X,'Possible error in ZMAP: No convergence after'/ & 4x,'10 iterations. New guess ??'/) C RETURN END C********************************************************** C* ZODE GEARLIKE SUBROUTINE ZODE *** C********************************************************** SUBROUTINE ZODE (T,ZZ,ZDZDT) C C COMPUTES DERIVATIVE ZDZDT C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) COMMON /PARAM2/ Z(0:20),C,CDWDT,W0,Z0,BETA(0:20), & QW(100),K0,N,NPTSQ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ C WZP = WMAPP (ZZ,Z0,W0,K0,N,C,Z,BETA,NPTSQ,QW) ZDZDT = CDWDT/WZP C RETURN END C********************************************************** C* ZQUAD GEARLIKE SUBROUTINE ZQUAD *** C********************************************************** FUNCTION ZQUAD (ZA,KA,ZB,KB,N,Z,BETA,NPTSQ,QW) 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 REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),BETA(0:N),QW(1) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ C C IF ((CDABS(ZA).GT.1.1D0).OR.(CDABS(ZB).GT.1.1D0)) WRITE(6,120) IF ((CDABS(ZA).GT.1.1D0).OR.(CDABS(ZB).GT.1.1D0)) WRITE(20,120) C ZMID = (ZA + ZB)/D2 ZQUAD = ZQUAD1 (ZA,ZMID,KA,N,Z,BETA,NPTSQ,QW) & - ZQUAD1 (ZB,ZMID,KB,N,Z,BETA,NPTSQ,QW) C 120 FORMAT (/' **** WARNING IN ZQUAD: Z OUTSIDE THE DISK') C RETURN END C********************************************************** C* ZQUAD1 GEARLIKE SUBROUTINE ZQUAD1 *** C********************************************************** FUNCTION ZQUAD1 (ZA,ZB,KA,N,Z,BETA,NPTSQ,QW) 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 REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),BETA(0:N),QW(1) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ C DATA RESPRM /1.D0/ C C CHECK FOR ZERO-LENGTH INTEGRAND C IF (CDABS(ZB-ZA).EQ.D0) THEN ZQUAD1 = C0 RETURN ENDIF C C STEP 1: ONE-SIDED GAUSS-JACOBI QUADRATURE FOR LEFT ENDPOINT: C R = DMIN1(D1,DIST(ZA,KA,N,Z)*RESPRM/CDABS(ZB-ZA)) ZAA = ZA + R * (ZB - ZA) ZQUAD1 = ZQSUM (ZA,ZAA,KA,N,Z,BETA,NPTSQ,QW) C C STEP 2: ADJOIN INTERVALS OF PURE GAUSSIAN QUADRATURE IF NECESSARY: C 100 IF (R.EQ.D1) RETURN R = DMIN1(D1,DIST(ZAA,0,N,Z)*RESPRM/CDABS(ZB-ZAA)) ZBB = ZAA + R * (ZB - ZAA) ZQUAD1 = ZQUAD1 + ZQSUM (ZAA,ZBB,0,N,Z,BETA,NPTSQ,QW) ZAA = ZBB GOTO 100 C END C********************************************************** C* DIST GEARLIKE SUBROUTINE DIST *** C********************************************************** FUNCTION DIST (ZZ,KS,N,Z) C C DETERMINES THE DISTANCE FROM ZZ TO THE NEAREST SIGNULARITY Z(K) C OTHER THAN Z(KS). C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N) C DIST = 99.D0 DO 100 K = 1,N IF (K.EQ.KS) GOTO 100 DIST = DMIN1(DIST,CDABS(ZZ-Z(K))) 100 CONTINUE C RETURN END C********************************************************** C* ZQSUM GEARLIKE SUBROUTINE ZQSUM *** C********************************************************** FUNCTION ZQSUM (ZA,ZB,KA,N,Z,BETA,NPTSQ,QW) C C COMPUTES THE COMPLEX LINE INTEGRAL OF ZPROD FROM ZA TO ZB BY C APPLYING A ONE-SIDE GAUSS-JACOBI FORMULA WITH POSSIBLE C SINGULARITY AT ZA. C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),BETA(0:N),QW(1) DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ C R = CDABS(ZB) C ZS = C0 ZH = (ZB-ZA)/D2 ZC = (ZB+ZA)/D2 C IF (KA.EQ.0) GOTO 100 IF (DABS(BETA(KA)).EQ.DH) GOTO 300 C 100 DO 200 I = 1,NPTSQ ZARG = ZC + ZH*QW(I) ZS = ZS + QW(NPTSQ+I) * (ZPROD(ZARG,0,N,Z,BETA)-C1) / ZARG 200 CONTINUE ZQSUM = ZS*ZH RETURN C 300 IOFFST = NPTSQ*2 IF (BETA(KA).GT.D0) IOFFST = IOFFST*2 DO 400 I = 1,NPTSQ ZARG = ZC + ZH*QW(IOFFST+I) ZS = ZS + QW(IOFFST+NPTSQ+I) * ZPROD(ZARG,KA,N,Z,BETA) / ZARG 400 CONTINUE ZQSUM = ZS*ZH FAC = ((D1-R)/D2)**BETA(KA) ZQSUM = ZQSUM*FAC - DLOG(R) C RETURN END C********************************************************** C* ZPROD GEARLIKE SUBROUTINE ZPROD *** C********************************************************** FUNCTION ZPROD (ZZ,KS,N,Z,BETA) C C COMPUTES THE INTEGRAND C C N C PROD (1-ZZ/Z(K))**BETA(K) K >< KS C K=1 C C NOTE: IN PRACTICE THIS IS THE INNERMOST SUBROUTINE C IN GLSOLV CALCULATIONS. THE COMPLEX LOG C CALCULATIONS BELOW MAY ACCOUNT FOR AS MUCH AS C HALF OF THE TOTAL EXECUTION TIME C IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) DIMENSION Z(0:N),BETA(0:N) DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ C ZSUM = C0 DO 100 K = 1,N IF (K.NE.KS) THEN ZTMP = C1 - ZZ/Z(K) ZSUM = ZSUM + BETA(K)*CDLOG(ZTMP) ENDIF 100 CONTINUE ZPROD = CDEXP(ZSUM) C RETURN END //GO.SYSIN DD gldrv.f echo gl.f 1>&2 sed >gl.f <<'//GO.SYSIN DD gl.f' 's/^-//' -C********************************************************** -C* GEAR-LIKE SAMPLE CALLING PROGRAM GEAR-LIKE *** -C********************************************************** - - IMPLICIT REAL*8 (A,B,D-H,O-V,X,Y), COMPLEX*16 (C,W,Z) - DIMENSION Z(0:20),W(0:20),BETA(0:20),QW(100),TH(0:20) - DATA D0,DH,D1,D2,D4 /0.D0,0.5D0,1.D0,2.D0,4.D0/ - DATA C0,C1,CI /(0.D0,0.D0),(1.D0,0.D0),(0.D0,1.D0)/ - CHARACTER*80 CMMNT - CHARACTER*40 FN,TITLE - CHARACTER*8 PLTCLR,CIRCLR,RADCLR - PI = D4*DATAN(D1) -C - 10 WRITE (6,30) - READ (5,20) FN - IF (FN(1:1).EQ.' ') GOTO 10 -C - IPER = INDEX (FN,'.') - IF (IPER.GT.0) FN = FN(:IPER-1) -C - OPEN (UNIT = 30, FILE = FN//'.TXT', STATUS = 'OLD') -C - READ (30,*) CMMNT - READ (30,*) TITLE - READ (30,*) CMMNT - READ (30,*) ITITLE - READ (30,*) CMMNT - READ (30,*) IFOOT -C - OPEN (UNIT = 20, FILE = FN//'.OUT', STATUS = 'NEW') - WRITE(20,120) TITLE -C -C PARAMETER WC ADDED FOR COMPATIBILITY OF CALLING SEQUENCES WITH -C THE SCHWARZ-CHRISTOFEL DRIVER -C - READ (30,*) CMMNT - READ (30,*) X, Y - WC = DCMPLX(X,Y) -C -C INITIALIZE NUMBER OF VERTICES AND VERTEX VALUES AND POWERS -C - READ (30,*) CMMNT - READ (30,*) N -C - READ (30,*) CMMNT - DO 50 I = 1,N - READ (30,*) BETA(I), RAD, ARG - W(I) = RAD * CDEXP(CI*ARG*PI) - 50 CONTINUE -C -C COMPUTE NODES AND WEIGHTS FOR PARAMETER PROBLEM: -C PASS N & BETA FOR COMPATIBILITY OF CALLING SEQUENCE WITH -C SCHWARZ-CHRISTOFEL DRIVER -C - READ (30,*) CMMNT - READ (30,*) NPTSQ - CALL QINIT (N,BETA,NPTSQ,QW) -C -C SOLVE PARAMETER PROBLEM: -C - IPRINT = 0 -C - READ (30,*) CMMNT - READ (30,*) IGUESS - IF (IGUESS.NE.0) THEN - READ (30,*) CMMNT - DO 100 I = 1,N - READ (30,*) ARGZI - Z(I) = CDEXP(CI*ARGZI*PI) - 100 CONTINUE - ELSE - READ (30,*) CMMNT - ENDIF -C - TOL = 10.0D0**(-(NPTSQ+1)) - CALL GLSOLV (IPRINT,IGUESS,TOL,ERREST, - & N,C,Z,WC,W,BETA,NPTSQ,QW) -C -C COMPUTE IMPAGE BOUNDARY VALUES BETWEEN VERTICES W(K) -C - TTOL = TOL*10.0D0 - READ (30,*) CMMNT - READ (30,*) IEDGE - IF ((ERREST.LT.TTOL).AND.(IEDGE.EQ.1)) THEN -C - CALL THETA (N,Z,TH) -C - DO 300 K = 1,N - WRITE(6,320) K-1,K - WRITE(20,320) K-1,K - DO 300 J = 1,9,2 - T = TH(K-1) + DFLOAT(J)*(TH(K)-TH(K-1))/1.D1 - ZT = CDEXP(CI*T) - WZT = WMAP (ZT,0,C0,WC,0,N,C,Z,BETA,NPTSQ,QW) - WRITE(6,321) T/PI,WZT,CDABS(WZT), - & DIMAG(CDLOG(WZT))/PI - WRITE(20,321) T/PI,WZT,CDABS(WZT), - & DIMAG(CDLOG(WZT))/PI - 300 CONTINUE - ENDIF -C -C CALL PLOTTING SUBROUTINE FOR PLOTTER (FORWARD MAP) -C - READ (30,*) CMMNT - READ (30,*) IPLOTF - READ (30,*) CMMNT - READ (30,*) PLTCLR - READ (30,*) CMMNT - READ (30,*) PAPER - READ (30,*) CMMNT - READ (30,*) SCAL - READ (30,*) CMMNT - READ (30,*) DMESH - READ (30,*) CMMNT - READ (30,*) INC - READ (30,*) CMMNT - READ (30,*) CIRCLR - READ (30,*) CMMNT - READ (30,*) INR - READ (30,*) CMMNT - READ (30,*) RADCLR -C - IF ((ERREST.LT.TTOL).AND.(IPLOTF.EQ.1)) THEN - CALL PLTTER (PAPER,SCAL,DMESH,INC,INR, - & FN,TITLE,ITITLE,IFOOT, - & PLTCLR,CIRCLR,RADCLR, - & N,C,Z,WC,W,BETA,NPTSQ,QW) - ENDIF -C -C CALL PLOTTING SUBROUTINE FOR PLOTTER (REVERSE MAP) -C - READ (30,*) CMMNT - READ (30,*) IPLOTR - READ (30,*) CMMNT - READ (30,*) PLTCLR - READ (30,*) CMMNT - READ (30,*) PAPER - READ (30,*) CMMNT - READ (30,*) DMESH - READ (30,*) CMMNT - READ (30,*) INC - READ (30,*) CMMNT - READ (30,*) CIRCLR -C - IF ((ERREST.LT.TTOL).AND.(IPLOTR.EQ.1)) THEN - CALL PLTTR2 (PAPER,DMESH,INC, - & FN,TITLE,ITITLE,IFOOT, - & PLTCLR,CIRCLR,TTOL, - & N,C,Z,WC,W,BETA,NPTSQ,QW) - ENDIF -C - CLOSE (UNIT = 20, STATUS = 'KEEP') - CLOSE (UNIT = 30, STATUS = 'KEEP') -C - 20 FORMAT (A40) - 30 FORMAT (//4X,'Enter parameter filename: ') - 120 FORMAT (//10X,'PLOT TITLE = ',A40//) - 320 FORMAT (/10X,'TH(',I2,') TO TH(',I2,')'/) - 321 FORMAT (3X,'TH/PI, W(Z), POLAR(W(Z)): ', - & F14.8,2(5X,'(',F14.8,',',F14.8,')')) -C - STOP - END - - - //GO.SYSIN DD gl.f echo gl.txt 1>&2 sed >gl.txt <<'//GO.SYSIN DD gl.txt' 's/^-//' -From ttacs1.ttu.edu!kkent Tue Dec 12 21:19:47 1989 -Date: 1 Dec 89 16:43:00 CST -From: "Pearce, Kent" -To: "ehg" - -'THE FOLLOWING LINE SETS THE VALUE OF THE TITLE HEADING ' -'FIGURE ' -'THE FOLLOWING LINE SETS THE VALUE OF ITITLE ' -1 -'THE FOLLOWING LINE SETS THE VALUE OF IFOOT ' -1 -'THE FOLLOWING LINE SETS THE VALUE OF WC ' -0.0D0 0.0D0 -'THE FOLLOWING LINE SETS THE VALUE OF N ' -6 -'THE FOLLOWING N LINES SET BETA (-EXTERIOR), RADIUS & ARGUMENT/PI AT W(I) ' -1.0 1.0 .25 --.5 1.5 .25 --.5 1.5 -.35 -1.0 .5 -.35 --.5 2.0 -.35 --.5 2.0 .25 -'THE FOLLOWING LINE SETS THE VALUE OF NPTSQ ' -8 -'THE FOLLOWING LINE SETS THE VALUE OF IGUESS ' -0 -'IF IGUESS >< 0, THEN THE NEXT N LINES SET THE VALUES OF [ARG Z(I)]/PI ' -'THE FOLLOWING LINE SETS THE VALUE OF IEDGE ' -0 -'THE FOLLOWING LINE SETS THE VALUE OF IPLOTF ' -1 -'THE FOLLOWING LINE SETS THE VALUE FOR THE PLOT COLOR' -'RED' -'THE FOLLOWING LINE SETS THE VALUE OF PAPER ' -4.0 -'THE FOLLOWING LINE SETS THE VALUE OF SCAL ' -2.0 -'THE FOLLOWING LINE SETS THE VALUE OF DMESH ' -0.01 -'THE FOLLOWING LINE SETS THE VALUE OF INC (PLOT CIRCLES) ' -1 -'THE FOLLOWING LINE SETS THE VALUE FOR THE CIRCLES COLOR ' -'GREEN' -'THE FOLLOWING LINE SETS THE VALUE OF INR (PLOT RADIALS) ' -1 -'THE FOLLOWING LINE SETS THE VALUE FOR THE RADIALS COLOR ' -'CYAN' -'THE FOLLOWING LINE SETS THE VALUE OF IPLOTR ' -0 -'THE FOLLOWING LINE SETS THE VALUE FOR THE PLOT COLOR' -'RED' -'THE FOLLOWING LINE SETS THE VALUE OF PAPER ' -3.0 -'THE FOLLOWING LINE SETS THE VALUE OF DMESH ' -0.01 -'THE FOLLOWING LINE SETS THE VALUE OF INC (PLOT CIRCLES) ' -1 -'THE FOLLOWING LINE SETS THE VALUE FOR THE CIRCLES COLOR ' -'GREEN' - - //GO.SYSIN DD gl.txt