# 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" <kkent@ttacs1.ttu.edu>
-To: "ehg" <ehg@research.att.com>
-
-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" <kkent@ttacs1.ttu.edu>
-To: "ehg" <ehg@research.att.com>
-
-'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