SUBROUTINE H{2 or 3 or A } LC(
complex wavenumber,
point (p
and the unit vector up, if necessary),
geometry of the panel
(vertices which define panel),
quadrature rule (weights and abscissae for the standard element),
validation and control parameters,
discrete Helmholtz integral operators (solution),
work space ).
View files H2LC.FOR, H3LC.FOR and H3ALC.FOR [33] in order to observe the subroutines.
The logical variable LVALID enables the user to switch on ( .TRUE.) or off ( .FALSE.) the validation of the input parameters. Setting LVALID= .TRUE. puts a series of tests on the input data into effect. The use of this facility will have a computational cost but could be useful when initially constructing a program that calls the subroutines or to facilitate the diagnosis of an error that appears in the calling program or in the subroutine.
The logical variables LLK, LMK, LMKT and LNK allows the user to choose the specific operators of Lk, Mk, Mkt or Nk that are required. The LLK, LMK, LMKT and LNK that have .TRUE. values results in the values of the corresponding discrete operator to appear in DISLK, DISMK, DISMKT and DISNK respectively, on exit from the subroutine.
Following from the techniques of defining the surfaces in Chapter 2, in H2LC and H3ALC two points QA and QB, each one-dimensional real arrays in Fortran with two components, define the panel. In H2LC, QA and QB are the edges of the straight line panel. In H3ALC, QA and QB define the outer rims of the truncated conical shell panel, the two points lying on the edges of the generator of the panel. In H3LC, three points QA, QB and QC, each one-dimensional real arrays in Fortran with three components, are the vertices that define the planar triangular panel.
The normal to the panel is computed within the subroutines. In H2LC and H3ALC, the normal is defined to be [ QA(2)- QB(2), QB(1)- QA(1)] normalised. Hence the normal is to the right on the line QA- QB. Note that in H3ALC the normal is the normal on the generator of the panel. In H3LC the normal is defined to be [QB-QA] × [QC-QA] (normalised), where × denotes the vector cross product. Hence the normal is presumed to be outward (to the observer) when QA- QB- QC are arranged anti-clockwise around the triangle.
The integrals that define Lk and Nk require special treatment when the point p lies on the panel. The logical variable LPONEL is used to inform the subroutine as to whether the point p lies on the panel ( .TRUE.) or not ( .FALSE.). In the case LPONEL= .TRUE., methods based on the expressions for the integral operators in equations (39) and (40) are used.
For the computation of the discrete Mkt and Nk operators the vector up needs to be defined. In the subroutines it is denoted VECP and it must represent a unit vector. Results for a non-unit vector up can be obtained by scaling the vector up before calling the subroutine (so that it is a unit vector) and suitably scaling the results DISMKT and DISNK afterwards. If the point p ( P) lies on the panel then up ( VECP) must be the unit normal to the panel, as defined above.
When the point p does not lie on the panel, the integrand is continuously differentiable and hence standard quadrature rules are satisfactory. However, when p lies on the panel, the integrands are continuous (after the subtractions described by equations (39) and (40)), but they are generally not differentiable at the point p. If the input quadrature rule does not take this property into account then it could lead to an inaccurate evaluation of the discrete forms.
In a practical computation, one of the subroutines will be called many times. The subroutines allow different quadrature rules to be input for each call. For example when p lies on the panel a special quadrature rule is advised, as discussed above. Also if one panel is relatively large or relatively close to the point p then a relatively more accurate quadrature rule is generally desirable.
In subroutines H2LC and H3LC the parameter MAXNQ sets the limit on the number of quadrature points in the quadrature rule defined in the generator and in the angular direction. The values of these parameters must not be changed between calls of the subroutines. The actual number of quadrature points, denoted NQ, must be less than or equal to MAXNQ but can be altered between calls. For the H3ALC subroutine the number of quadrature points consists of NGQ along the generator and NTQ in the angular direction. The values of NGQ and NTQ may be changed between calls of H3ALC, but they must always be less than MAXNGQ and MAXNTQ respectively. Hence the total number of quadrature points used in H3ALC is NGQ× NTQ.
In subroutine H2LC, the parameters AQ and WQ store the abscissae and the weights of the quadrature rule. The parameters AGQ, WGQ and ATQ, WTQ store similar quadrature rules for the generator and angular directions in H3ALC. Note that in H3ALC the panel is defined as the result of rotating the generator QA- QB about the z axis. In practice, the size of the panel in the angular direction will vary greatly between calls of H3ALC. Since it is good practice to distribute the quadrature points fairly evenly over the panel, the number of quadrature points should be used in the angular direction should be roughly proportional to the `radius' of the panel.
Generally the most efficient quadrature rules to use in subroutines H2LC and H3ALC are standard Gauss-Legendre rules. The weights and abscissae for such rules are listed in Stroud and Secrest [79] and can also be generated using NAG routine D01BBF, for example. When the point p lies on the panel then an efficient quadrature rule is obtainable by dividing the domain at p and using suitably scaled Gauss-Legendre rules on the separate regions. Note that such division is only necessary in the generator direction in H3ALC.
In subroutine H3LC, the parameters XQ, YQ, WQ store the x coordinates, the y coordinates and the weights of the quadrature rule over the standard triangle. Quadrature rules of the Gauss-Legendre type are listed in Laursen and Gellert [56]. In the case when p lies on the panel then a suitable quadrature rule can be obtained by dividing the triangle into three regions, by connecting the point p to the vertices of the panel, and map suitably scaled standard quadrature rules onto the three separate triangular regions.
If LVALID is .TRUE. then the parameters to the subroutines EGEOM and EQRULE set the maximum absolute tolerance allowed in the input geometrical and quadrature rule data. For example if the components of QA and QB are all within EGEOM of each other (the points coincide) or the sum of the weights in the quadrature rules is greater than NQ× EQRULE (allowing for an error of EQRULE in each component) then an error message will be output to the error file. In such cases the subroutine will be aborted and the parameter LFAIL will register .TRUE..