Circular arc polygons, numerical conformal mappings, and moduli of quadrilaterals

We study numerical conformal mappings of planar Jordan domains with boundaries consisting of finitely many circular arcs and compute the moduli of quadrilaterals for these domains. Experimental error estimates are provided and, when possible, comparison to exact values or other methods are given. The main ingredients of the computation are boundary integral equations combined with the fast multipole method.


Introduction
Many applications of complex analysis make use of conformal mappings or harmonic functions in an essential way. Examples include modeling of airfoils, analysis of turbulence, design of dams, and mathematical theory of electricity [1]. A first example is a mathematical model of a physical condenser, which is a planar domain G together with a compact set E ⊂ G [39]. The pair (G, E) is called a condenser and its capacity is defined as [14] cap(G, E) = inf u∈A G |∇u| 2 dm, where A is the family of all harmonic functions with u(x) ≥ 1 for all x ∈ E and u(x) → 0 when x → ∂G. In concrete applications the sets E and ∂G have a simple geometry, both have a finite number of components, each component being a piecewise smooth curve. In this case, it is known that the infimum is attained by a harmonic function. A second example is to solve the following Dirichlet-Neumann boundary value problem for the Laplace equation in a planar Jordan domain G with a piecewise smooth boundary ∂G = ∪ 4 k=1 ∂G k where the sets ∂G k occur in positive order-the sets ∂G k and the domain G define a quadrilateral. This problem is 0, on G, u = 1, on ∂G 1 , u = 0, on ∂G 3 , ∂u/∂n = 0, on ∂G 2 , ∂u/∂n = 0, on ∂G 4 .
If u is a solution function to this problem, then the expression Q |∇u| 2 dxdy defines the modulus of the quadrilateral. Sometimes, conformal mappings can be applied to simplify the geometry and, in fact, to solve this Dirichlet-Neumann problem. By the Riemann mapping File: circa20210723.tex, printed: 2021- 11-12, 7.52 theorem we know that a Jordan domain G is conformally equivalent to the unit disk, but this theorem does not give a clue for finding the conformal map. In the case of polygonal domains, widely used numerical methods based on the Schwarz-Christoffel transformation have been developed by D. Gaier [15], L.N. Trefethen and T.A. Driscoll [13], and N. Papamichael and N. Stylianopoulos [37]. Jordan domains where the boundary is a union of finitely many circular arcs, have been studied by P. Brown and M. Porter [9], [40], by U. Bauer and W. Lauf [5], and, in particular, by D. Crowdy [11]. See also [6,22,41]. We give new numerical methods for this same case and present our results in the form of numerical tables, graphics, and analysis of algorithm performance. The method is based on boundary integral equations (R. Kress, [26]) as developed and implemented by M.M.S. Nasser in a series of papers during the past two decades, see e.g. [23], [29]- [32]. The method uses the fast multipole method implementation from [17] for the speed-up of solving linear equations. In a recent series of papers ( [32]- [34]), the method was applied for the capacity computation of planar condensers and for the study of isoperimetric problems for capacity. In particular, we will make use of the very recent results in [32].
The structure of this paper is as follows. In Section 2, we give the basic facts about the boundary integral method and its efficient implementation. Section 3 describes the method of Kress [24] for the treatment of non-smooth boundary points in numerical integration and a refinement principle in trapezoid rule. Section 4 contains applications of the presented method to several circular arc polygonal quadrilaterals. Several numerical examples for the computation of the exterior modulus of quadrilaterals are presented in Section 5. Finally, in Section 6, the proposed method is applied to gear domains.

Preliminary notions
The capacity of a condenser defined in the introduction can be defined in many equivalent ways as shown in [14], [16]. First, the family A may be replaced by several other families by [16,Lemma 5.21,p. 161]. Furthermore, cap (G, E) = M(∆(E, ∂G; G)), (2.1) where ∆(E, ∂G; G) is the family of all curves joining E with the boundary ∂G in the domain G and M stands for the modulus of a curve family [16,Thm 5.23,p. 164]. For the basic facts about capacities and moduli, the reader is referred to [14,16,20].

Quadrilaterals.
A Jordan domain in the complex plane C is a domain with boundary homeomorphic to the unit circle. A quadrilateral is a Jordan domain D together with four distinguished points z 1 , z 2 , z 3 , z 4 ∈ ∂D which define a positive orientation of the boundary. In other words, if we traverse the boundary, then the points occur in the order of indices and the domain D is on the left hand side. The quadrilateral is denoted by (D; z 1 , z 2 , z 3 , z 4 ). The modulus of the quadrilateral is a unique positive number h such that D can be conformally mapped by some conformal map f onto the rectangle with vertices 0, 1, 1 + ih, ih such that The modulus is denoted mod(D; z 1 , z 2 , z 3 , z 4 ). The following basic formula is often used: mod(D; z 1 , z 2 , z 3 , z 4 ) mod(D; z 2 , z 3 , z 4 , z 1 ) = 1 .
A simple example of a quadrilateral is the case when the domain D is a rectangle with sides a and b and the points z 1 , z 2 , z 3 , z 4 are its vertices. Depending on the labeling of the vertices, the modulus is either a/b or b/a. An alternative equivalent definition is based on the Dirichlet-Neumann problem mentioned in the introduction L.V. Ahlfors [2,Thm 4.5,p. 63].
2.4. Quadrilateral modulus and curve families. The modulus of a quadrilateral (D; z 1 , z 2 , z 3 , z 4 ) is connected with the modulus of the family of all curves in D, joining the opposite boundary arcs (z 2 , z 3 ) and (z 4 , z 1 ), in a very simple way, as follows 2.6. Grötzsch ring and elliptic integrals. A ring domain, or briefly a ring, in the plane is a domain whose complement has exactly two components. Such a domain can be conformally mapped onto {z ∈ C : r < |z| < 1}. The number log(1/r) is called the modulus of the ring. One can also consider a ring R with complementary components E, F ⊂ C as a condenser (C \ E, F ) and define the capacity of the ring as The formula (2.1) gives now the connection between the modulus of a ring and the modulus of the family of curves joining its boundary components. The so called Grötzsch ring G r = B 2 \ [0, r], 0 < r < 1, is frequent in the study of capacities. In the case n = 2, r ∈ (0, 1), the following explicit formulas hold for the capacity of this ring [20, (7.18), p. 122], where K(r) and K (r) are the elliptic integrals of the first kind The elliptic integrals E(r) and E (r) of the second kind are 2.10. History of numerical conformal mapping. The Schwarz-Christoffel method described in the introduction has a long history which goes back to the nineteenth centure, see [7]. During the past fifty years, the development of computational methods has revolutionized the applications of numerical conformal mapping, see [15,13,37,27].

Boundary integral method
3.1. The integral equation. We assume that the boundary Γ = ∂G is a smooth Jordan curve parametrized by a 2π-periodic function η : [0, 2π] → Γ which is twice continuously differentiable and satisfies η (t) = 0 for all t ∈ [0, 2π] (piecewise smooth boundaries will be considered in the next subsection). The boundary Γ is oriented such that G is to the left of Γ, i.e., Γ is oriented counterclockwise for bounded G and clockwise for unbounded G. We denote by H the space of all Hölder continuous real-valued functions on the boundary Γ.
where α is a given auxiliary point in the domain G. The generalized Neumann kernel N (s, t) is defined by [43] The kernel N (s, t) is continuous on [0, 2π] × [0, 2π] and hence, the integral operator N defined on H by is compact. The integral equation involves also the kernel which is singular and has the representation Here the kernel M 1 is continuous on The integral operator M defined on H by is singular, but is bounded on H [43].   A MATLAB function fbie for solving the integral equation (3.10) in O(n log n) operations, where n is the number of nodes in the interval [0, 2π], is presented in [30]. In the function fbie, the integral equation is discretized using the Nyström method with the trapezoidal rule and then using the GMRES method to solve the obtained linear system. The matrix-vector product in the GMRES method is computed using the MATLAB function zfmm2dpart in the toolbox FMMLIB2D [17]. Let et, etp, A, gam, and rho, be n × 1 discretization vectors of the functions η(t), η (t), A(t), γ(t), and ρ(t), respectively. Then, the n × 1 vectors rho and h are computed by Theoretically, all elements of the vector h are equal to a constant h in (3.11). Numerically, we will approximate the constant h by the arithmetic mean of the elements of the vector h. For the other parameters in fbie, we choose iprec = 5, gmrestol = 0.5 × 10 −14 , restart = [ ], and maxit = 100. This means that the tolerances of the methods are 0.5 × 10 −15 for FMM and 0.5 × 10 −14 for GMRES. Moreover, the GMRES is used without restart, and the maximum number of GMRES iterations is 100. The auxiliary points α in (3.2) need to be chosen sufficiently far away from the boundary Γ. For some domains (see Examples 4.10 and 4.12 below), we need to choose α carefully to ensure the convergence of the method.

Kress method.
In this paper, we shall assume that Γ is a piecewise smooth Jordan curve with a finite number of corner points such that each of these corner points is not a cusp. We assume that the tangent vector of the boundary has only the first kind discontinuity at each corner point where the left tangent vector at each corner point is considered as the tangent vector at this point. For such boundaries Γ, the integral operator with the generalized Neumann kernel (3.10) is not compact, but this operator can be written as a sum of a compact operator and bounded non-compact operator with norm less than one in suitable function spaces [33]. Hence, we can apply the Fredholm theory to the integral equation with the generalized Neumann kernel although the operator is not compact [24].
Using the method described above to solve the integral equation (3.10) when Γ is a piecewise smooth Jordan curve yields only poor convergence since the solution of the integral equation (3.10) has a singularity in its first derivative in the vicinity of the corner points [24,33]. To achieve a satisfactory accuracy, we first remove the discontinuity of the derivatives of the solution of the integral equation at the corner points by using a suitable substitution [24,25]. Then, the transformed equation can be solved using the above method.
Following Kress [24,25], we define a bijective function w : The function w is strictly monotonically increasing and infinitely differentiable function, and the integer p ≥ 2 is the grading parameter. In our numerical experiments below, we choose p = 3.
3.14. Parametrizing the boundaries of polygonal domains. In this paper, we assume that the boundary Γ is a polygon consisting of a finite number of finite segments, circular arcs, or both segments and circular arcs. The boundary Γ will be parametrized as described in the preceding subsection. We discretize the interval [0, 2π] by n equidistant nodes, Then the parametrization of the boundary Γ is discretized by η(s 1 ), η(s s ), . . . , η(s n ) where s k = δ(t k ), k = 1, 2, . . . , n. Similarly, the derivative of the parametrization of the boundary is discretized by η (s k )s k where s k = δ (t k ), k = 1, 2, . . . , n. A MATLAB function plgsegcirarcp.m for computing such a parametrization can be downloaded from https://github.com/mmsnasser/ circa. To use this MATLAB function, assume that Γ is a polygon with m finite vertices v 1 , v 2 , . . . , v m with v m+1 = v 1 . For the parametrization η(t) of the boundary Γ, it follows from (3.13) and (3.15) that For k = 1, 2, . . . , m, we assume that the center of the portion of the boundary between v k and v k+1 is c k if the portion is a circular arc and c k = ∞ if the portion is a segment. Further, if the portion of the boundary between v k and v k+1 is a circular arc, then we introduce an indicator d k for it with value 1 if the arc is positively oriented with respect to the pertaining circle center c k and −1 otherwise. For segment portions of the boundary this indicator is d k = 0. We assume here that n is an integer multiple of m so that each side of the polygon will be discretized by n/m points. Define the vectors Then, discretizations of the parametrization of the boundary and its derivative can be computed using the MATLAB function plgsegcirarcp.m by calling 3.17. Conformal mapping onto the unit disk. In this subsection, we review a numerical method for the computation of the conformal mapping w = Φ(z) from a polygonal domain G onto the unit disk D = Φ(G) [35,31]. For a bounded domain G, let γ = − log |η(t) − α|, let ρ be the unique solution of the integral equation (3.10), and let the constant h be given by (3.11). Then, the mapping function Φ with normalization can be written for z ∈ G ∪ Γ as The boundary values of the conformal mapping are now given by is a one-to-one increasing function on [0, 2π] with S(2π) − S(0) = 2π. The function S(t) is known as the boundary corresponding function [21, p. 380]. Differentiation yields When G is an unbounded domain, we assume γ = log |η(t) − z 1 | where z 1 is a given point in the exterior domain of G, i.e., z 1 is in the bounded domain enclosed by Γ. We assume also that ρ is the unique solution of the integral equation (3.10) and the constant h is given by (3.11). Then, the mapping function Φ with the normalization where c = lim z→∞ (zΦ(z)) > 0 and the function f (z) is analytic in G with the boundary values Hence, the boundary values of the conformal mapping are given by The mapping function Φ maps the boundary Γ onto the unit circle ∂D = Φ(Γ). Although the orientation of Γ is clockwise (G is on the left of Γ), the orientation of ∂D will be counterclockwise. Then |Φ(η(t))| = 1 and hence Φ(η(t)) can be written as in (3.20) where is a one-to-one increasing function on [0, 2π] with S(2π) − S(0) = 2π. By differentiation, we obtain For both bounded and unbounded G, the function ρ(t) is computed by solving the integral equation (3.10) and its derivative ρ (t) is computed by approximating ρ(t) by a trigonometric interpolating polynomial and then differentiating the interpolating polynomial. This polynomial can be computed using FFT [42]. Then the boundary corresponding function S(t) is computed through (3.21) for bounded G and by (3.23) for unbounded G. The boundary values of the mapping function Φ can be computed through (3.20). Thus, the function is a parametrization of the unit circle. The values of the mapping function w = Φ(z) for z ∈ D as well as the values of the inverse mapping function z = Φ −1 (w) for w ∈ D can be computed using the Cauchy integral formula. For the direct mapping, we have For the inverse mapping, we have The mapping function Φ maps the vertices v k , k = 1, 2, . . . , m, of the circular arc polygon domain D onto points w k = Φ(v k ), k = 1, 2, . . . , m, on the unit circle. In literature, the points w k , k = 1, 2, . . . , m are know as the preimages of the vertices v k , k = 1, 2, . . . , m. In view of (3.16), it follows from (3.24) that the preimages are given by w k = ζ(t (k−1)n/m+1 ) = e iS(t (k−1)n/m+1 ) , k = 1, 2, . . . , m.
The implementation of the method presented in this subsection is given in the following MATLAB function mapdisk.m. All computer codes of the computations presented in this paper can be found in the link https://github.com/mmsnasser/circa.

Modulus of quadrilaterals.
Consider the quadrilateral (D; z 1 , z 2 , z 3 , z 4 ) where D is a bounded simply connected polygonal domain and z 1 , z 2 , z 3 , z 4 are four distinguished points on ∂D with counterclockwise orientation. The modulus mod(D; z 1 , z 2 , z 3 , z 4 ) can be computed in two steps. In the first step, we map the domain D using the conformal mapping w = Φ(z) described in § 3.17 onto the unit disk D. The boundary ∂D is then mapped onto the unit circle. Here, the points z 1 , z 2 , z 3 , z 4 need not to be vertices of the polygon. Assume that z k = η(t k ) wheret k ∈ [0, 2π], k = 1, 2, 3, 4. Then the four points z 1 , z 2 , z 3 , z 4 will be mapped onto four points w 1 , w 2 , w 3 , w 4 on the unit circle where, by (3.20), w k = e iS(t k ) , k = 1, 2, 3, 4. By the conformal invariance of the modulus, we have mod(D; z 1 , z 2 , z 3 , z 4 ) = mod(D; w 1 , w 2 , w 3 , w 4 ).
In the second step, the modulus mod(D; w 1 , w 2 , w 3 , w 4 ) will be computed using the exact formula [37, (2.6.1)] where the absolute (cross) ratio |w 1 , w 2 , w 3 , w 4 | is defined by [20, p. 33] (3.27)  As in the bounded case, the unbounded domain D − can be mapped using the conformal mapping w = Φ(z) described in § 3.17 onto the unit disk D so that the four points z 1 , z 2 , z 3 , z 4 on ∂D (in clockwise orientation) are mapped onto four points w 1 , w 2 , w 3 , w 4 on the unit circle ∂D (in counterclockwise orientation). Assume that z k = η(t k ) wheret k ∈ [0, 2π], then (3.20) implies that w k = e iS(t k ) , k = 1, 2, 3, 4. By the conformal invariance of the modulus, the exterior modulus of the quadrilateral (D; z 1 , z 2 , z 3 , z 4 ) is equal to mod(D; w 1 , w 2 , w 3 , w 4 ) which can be computed using the exact formula (3.26).

Half-disk. Let
Then (A; z 1 , z 2 , z 3 , z 4 ) where z 1 = r, z 2 = s, z 3 = e iσ , z 4 = e iβ , is a quadrilateral and its modulus can be computed by means of elementary conformal mappings and it is The approximate numerical values of the modulus mod(A; z 1 , z 2 , z 3 , z 4 ) are computed by the above described method with n = 2 13 and the exact values are presented in Table 1. The relative error in the approximated values are also presented in Table 1  , The above method with n = 2 13 is used to compute approximate values of the modulus mod(T ; z 1 , z 2 , z 3 , z 4 ) for several values of L ∈ (1, 5]. The relative error in the computed values is presented in Figure 3. The values of the function µ and its inverse are computed as described in [36]. Figure 3 presents also the relative error in the approximate values obtained using the SC Toolbox [12]. As we see from the figures, for L > 2.5, the relative errors of the two methods are almost identical. For small L, the results obtained with SC Toolbox are better than the results obtained by the proposed method.  where P + is the polygon with vertices i, 0, b, and a + i (see Figure 4 (left)). If b = a + 1, then by (4.5), the exact value of the modulus is given by For b = a + 1, the proposed method is used with n = 2 13 to compute approximate values of mod(P + ; i, 0, z 3 , z 4 ) for several values of a. The relative error in the computed values is presented in Figure 4 (right). For other values of b, let the real function u(a, b) be defined for u(a, b) = 2mod(P + ; i, 0, b, a + i).
The values of the function u(a, b) are computed for several values of a, b using the above proposed method and using the SC Toolbox. The contour lines of the values of the function u for the proposed method are presented in Figure 5 (left). In Figure 5 (right), we present the absolute value of the difference between the values of the function u obtained by the above method and by the SC Toolbox.

4.8.
Trapezoid with a small curvature at the left side. In this example, we consider the trapezoid T with the vertices z 1 = 0, z 2 = 1 + i tan πσ, z 3 = 1 + i(L + tan πβ), z 4 = iL, and with a small curvature at the left side parametrized the angle π (see Figure 6 (left)). This trapezoid has been considered in [3, p. 14] but no numerical results were presented in [3].
In this paper, we use the above method with n = 2 12 to compute approximate values of the modulus mod(T ; z 1 , z 2 , z 3 , z 4 ) for several values of 0.0001 ≤ ≤ 0. 25. The values of the computed modulus is presented in Figure 6 (right) for L = 2, σ = π/8 and β = π/4.
and Ψ(z) = Γ (z)/Γ(z) is the digamma function. Here, we approximate t 0 with   [37, p. 44]). We consider here v 2 and v 4 as vertices. The L-shaped domain with four vertices (in counterclockwise orientation) is a quadrilateral. There are 280 possible choices of such quadrilaterals. The proposed method is used with n = 2 13 to compute the modulus for these 280 quadrilaterals and the reciprocal error based on the identity (2.3), i.e., the reciprocal error in Tables 2 and 3 is defined by (4.11) |1 − mod(L; z 1 , z 2 , z 3 , z 4 ) mod(L; z 2 , z 3 , z 4 , z 1 )| .
The results are presented in Figure 9 (left). Extensive numerical tests [18] related to capacity computation show that the reciprocal error (4.11) agrees with several other error estimates. However, for the current proposed method, the reciprocal error is not significant because the method is based on mapping the domain D and the four points z 1 , z 2 , z 3 , z 4 on its boundary to the unit disk D with four points w 1 , w 2 , w 3 , w 4 on the unit circle. Then mod(D, w 1 , w 2 , w 3 , w 4 ) is computed using the exact formula. Thus the reciprocal error in mod(D, z 1 , z 2 , z 3 , z 4 ) is the same as the reciprocal error in mod(D, w 1 , w 2 , w 3 , w 4 ). Thus, the reciprocal error in our method measures only the error in the numerical computation of the special function "µ" in the exact formula.
For this example, the exact values of the modulus of the L-shaped quadrilateral for several choices of vertices are given in [15]. Table 2 presents these exact values as well as the approximate values obtained using the proposed method and the relative error in the approximate values. In this example, as well as in the next example, the auxiliary point α in (3.2) need to be chosen carefully to ensure the convergence of the method. In our numerical computation we choose α inside the domain L, sufficiently far from the boundary, and close to the arithmetic mean of the vertices z 1 , z 2 , z 3 , z 4 of the quadrilateral (L; z 1 , z 2 , z 3 , z 4 ).   Table 3. The values of the modulus of the L-shaped (circular arc) quadrilateral for several choices of vertices.
4.12. L-shaped domains: circular arc polygonal boundary. Now, we consider a circular arc polygon with the same vertices as in Example 4.10 (see Figure 8 (right)). This polygon is obtained by replacing each side-segments in the polygon in Figure 8 (left) by a circular arc such that the angle between the segment and the tangent to the circular arc is . We consider here = 1/3. The proposed method is used with n = 2 13 to compute the modulus for the 280 possible choices of quadrilaterals and the reciprocal error based on the identity (2.3). The obtained results are presented in Figure 9 (right). The values of the modulus of this L-shaped quadrilateral for several choices of vertices are given in Table 3. Fig. 6]). This polygon consists of 5 straight segments and 2 circular arcs with centers 0.9 − 2i and −2 + 0.8i, respectively (see Figure 10). There are 140 possible choices of four vertices z 1 , z 2 , z 3 , z 4 to get a quadrilateral. The proposed method is used with n = 7 × 2 10 to compute the modulus for these 140 quadrilaterals and the approximate values of the modulus for some of these choices are presented in Table 4.  Table 4. The values of the modulus of the polygon with 5 straight segments and 2 circular arcs.

Examples: exterior modulus of quadrilaterals
In this section, we consider several numerical examples to illustrate the accuracy of the proposed method for computing the exterior modulus of quadrilaterals. In the first example, the exact value of the exterior modulus is known. For the second example, we compare the above proposed method against the methods presented in [19].

Rectangle.
In this example, we consider the rectangle R with the vertices z 1 = 0, z 2 = ib, z 3 = a + ib, z 4 = a with a, b > 0. The exact value of the exterior modulus of the quadrilateral (R; z 1 , z 2 , z 3 , z 4 ) is given by [37, p. 82] In our numerical examples below, we assume that a = 1 and we choose several values of κ such that 0.02787 ≤ κ ≤ 0.7306, then 0.02 < b = 1/ψ(κ) < 10. For these values of b, the proposed method with n = 2 13 is used to compute approximate values of the exterior modulus of the quadrilateral (R; z 1 , z 2 , z 3 , z 4 ). The relative error in the computed values is presented in Figure 11 where the exact values of the exterior modulus is computed by (5.2). Figure 11. The relative error in the computed approximate values of the exterior modulus vs b for the rectangle with the vertices z 1 = 0, z 2 = ib, z 3 = 1 + ib, z 4 = 1.

Method
(P 1 ; z 1 , z 2 , z 3 , z 4 ) (P 2 ; z 1 , z 2 , z 3 , z 4 ) Our Method 0.9923416331 0.9592571729 SC Toolbox [19] 0.9923416323 0.9592571721 AFEM [19] 0.9923500126 0.9593012739 hp-FEM (Interior) [19] 0.9923416332 0.9592571731 hp-FEM (Exterior) [19] 0.9923416332 0.9592572007 6. Conformal mapping onto gear domains 6.1. Gear domains. A gear domain D is a special case of the circular arc polygons described above. It is a starlike simply connected domain containing the origin and bordered by arcs of circles centered at the origin and segments of lines passing through the origin [8,9,10,38]. Here, we assume that D is bounded. The method presented in § 3.17 can be used to compute the conformal mapping w = Φ(z) for the gear domain D onto the unit disk D normalized by (3.18) as well as its inverse z = Φ −1 (w) from D onto D. Assume that D has m vertices v k , k = 1, 2, . . . , m.
Then the method can be used to compute also the preimages w k , k = 1, 2, . . . , m, of these vertices.

6.2.
A gear domain with 6 vertices. As our first example, we consider a gear domain with 6 vertices. This example has been considered in [38, Fig. 4(b)] (although the vertices of this domain are not given explicitly in [38], we approximate these vertices from Fig. 4(b) in [38]). The vertices are given in Table 6. The method presented in § 3.17 is used with n = 3 × 2 11 to compute the conformal mapping from the gear domain D onto the unit disk and its inverse. Figure 13 (left) shows the images of several circles |w| = r, for r = 0.3, 0.45, 0.6, 0.75, 0.9, in the unit disk under the inverse conformal mapping z = Φ −1 (w). The image of the circle |z| = r or part of the circle for r = 0.19, 0.58, 0.88, 1.13, in the domain D under the conformal mapping w = Φ(z) is shown in Figure 13 (right). The square markers on the unit circle are the preimages of the vertices of the gear domain. The approximate values of the preimages are presented in Table 6. 6.3. One-tooth gear domain. Numerical Conformal Mappings onto one-tooth gear domains have been considered in [8,9,10]. Without loss of generality, a one-tooth gear domain is a circular arc polygonal domain with the vertices v 1 = e −iθ , v 2 = βe −iθ , v 3 = βe iθ , and v 4 = e iθ where 0 < θ < π is the gear angle and β > 1 is the gear ratio [9]. The method presented in § 3.17 is used with n = 2 13 to compute the conformal mapping from the gear domain D onto the unit disk and its inverse for β = 1.5 and θ = π/6. Figure 14 (left) shows the images of several circles |w| = r, for r = 0.09, 0.19, . . . , 0.99, under the inverse conformal mapping z = Φ −1 (w). The image of the circle |z| = r or part of the circle for r = 0.09, 0.19, . . . , 1.49, under the conformal mapping w = Φ(z) is shown in Figure 13 (right). The square markers on the unit circle indicate the preimages of the vertices of the gear domain. Figure 14. The gear domain with one-tooth for θ = π/6 and β = 1.5 (left) and the unit disk (right).
For a fixed β and for 0 < θ < π, the modulus of the quadrilateral (D; e −iθ , βe −iθ , βe iθ , e iθ ) has been computed using the proposed method with n = 2 13 . The results for β = 1.1, 1.25, 1.5, 2 are presented in Figure 15. It is clear from Figure 15 that the modulus approaches zero as θ → 0 or θ → π. Further, the results presented in Figure 15 validate numerically the conjuncture in [10, p. 90]. That is, there are exactly two gears corresponding to two different values of θ with the same modulus except for one value of θ ∈ (0, π) where the modulus has its maximum value. These maximums are marked with squares in Figure 15. The location of the maximum value depends on the value of β and it moves towards π as β increases. Moreover, the maximum value of the modulus increases as β decreases toward 1.

A multitooth gear domain.
We consider a multitooth gear domain with 12 vertices as in Table 7. The method presented in § 3.17 is used with n = 3 × 2 12 to compute the conformal mapping from the gear domain D onto the unit disk and its inverse.  0.95294901093885 + 0.30313063611365i Figure 16. The multitooth gear domain (left) and the unit disk (right).