As you might already guess by the name of this website, the spherical wave horn inspired my work a lot. If we assume expanding spherical wave fronts in round horns and want to stretch the round profile to an ellipse then we must inevitably think about ellipsoidal surface areas and of course the associated mathematics. My own learning phase was not easy either, until I found myself able to mathematically master the challenge. I will describe the results of my work in this post and try to give as many details as possible and describe as much math as necessary. The whole stuff is quite complicated and therefore we will simply start with the two-dimensional part and then gradually come to the ellipsoidal surfaces.

Regarding the transformation of the round mouth basic form into any elliptical variant, the mathematical algorith is based on the general Lamé curves which are defined by the following equation:
\tag{1}\left|\frac{x}{a}\right|^n + \left|\frac{y}{b}\right|^n = 1
Piet Hein introduced the name Superellipse for these curves with exponents above 2 and especially the exponent 2.5 was favored by him. References can be found here: Ref. 1, Ref. 2, Ref. 3, Ref. 4.
It is much easier to use polar coordinates to calculate the desired points on these curves. For the general ellipse using polar coordinates we get:
\tag{2} x = a \cdot \cos{\theta}
\tag{3} y = b \cdot \sin{\theta}
\tag{4} \theta \in 0 \rightarrow 2\pi
For the transformation of the general elliptic equations into the Lamé form, we can recall and transform the well-known relationship:
\tag{5} \cos{\theta}^2 + \sin{\theta}^2 = 1
\tag{6} \left(\left|\cos{\theta}\right|^\frac{2}{n}\right)^n + \left( \left|\sin{\theta}\right|^\frac{2}{n}\right)^n = 1
This results in the sought-after Cartesian coordinates of the general Lamé form:
\tag{7} x = sign\left(\cos{\theta}\right) \cdot a \cdot \left|\cos{\theta}\right|^\frac{2}{n}
\tag{8} y = sign\left(\sin{\theta}\right) \cdot b \cdot \left|\sin{\theta}\right|^\frac{2}{n}
During the calculation of different variations for the elliptical shapes, simply setting n = 2.0 as input parameter gives the general ellipse profile. As a rule of thumb, only values n \geq 2.0 \leq 8 really make sense as values larger than 8 could results in ugly artifacts for the point cloud. When n > 2.0 is selected, the general ellipse slowly changes to a rectangular shape with rounded corners for increasing n. Using this set of equations many different shapes are possible.
Pure elliptical (n = 2.0):

Super Ellipse (n = 2.5), the exponent value favored by Piet Hein:

And a Superellipse n = 8.0 with nearly rectangular shape:

While looking at this series of images, it should be clear that the corresponding surface areas increase with increasing exponent n with Superellipse main axes unchanged. The method used to correct this behavior is simply to find another set of shorter main axis while preserving the provided stretch factor. This can be done by an iterative procedure that minimizes an error function. More about this procedure will be described later in this post.

But how does the transition from round to elliptical or to rectangular shape is implemented in the program code without creating discontinuities in the profile along the horn axis? The transition should eventually be smooth and without abrupt profile changes. Because of these requirements the round throat surface area is always declared as starting point of my algorithm and if the horn axis is called z, then by definition we have z=0 at this point. The radius of the throat surface r_t is an input parameter and must be provided to the program. Then it goes on by determining the radius r_m (z) which is valid at a given point on the horn axis according to the respective horn-specific formulas. The associated circular flat surface area results to:
\tag{9} S_{round} (z) = \pi \cdot r_m^2 (z)
If we want to convert the circle into an ellipse, this can be done as a first approximation by holding the flat surface areas the same. The surface area of the general ellipse is defined as:
\tag{10} S_{ellipse} (z) = \pi \cdot a(z) \cdot b(z)
For a practical implementation of the algorithm in my programs, stretching means that the major axis a gets longer and the minor axis b gets shorter. The corresponding stretch factor is denoted here as f_s(z). My decision was to not stretch the major axis but to compress the minor axis by being divided by f_s(z):
\tag{11} b(z) = \frac{r_m(z)}{f_s(z)}
\tag{12} a(z) = \frac{r_m^2(z)}{b(z)}
Once the two main axes of the ellipse have been calculated, the above formulas (7) and (8) are used to determine all points on the ellipse we are looking for our point cloud. Strictly speaking, this only applies to the general ellipse with n = 2.0. Therefore, in the program code, the variables just described are used only as start parameters as we intend to yield the surface area of an triaxial ellipsoid cap which is per definition assumed to be equal to the spherical surface cap area derived from the horn formulas (for more details see below).

The stretch factor itself must change steadily and smoothly along the horn axis. It should be clear that at the horn throat f_s(z_0) = 1 and that f_ {s}(z) at a distinct point on the horn axis should reach the value of the program input parameter f_{s, inp}. It generally makes sense here that the maximum stretch factor is reached exactly at the horn mouth. But for the spherical wave horn, however, the stretch factor can also be set as a function of the progressing horn wall radius, and the maximum value is then achieved in the rollback region of the horn at the widest extension of the profile. For the ease of use within the program code, I personally found it useful to rewrite the formula of the stretch factor in a more practical form with respect to programming efficiency. Assuming that a stretch factor of f_{s,inp} = 1.5 was defined as an input parameter, then one can also write:
\tag{13} f_{s,inp} = 1 + 0.5
\tag{14} {f_{s,inp}^{'}} = 1 – f_{s,inp}
The stretch factor at a certain point of the horn axis then becomes:
\tag{15} f_{s}(z) = 1 + {f_{s,inp}^{'}} \cdot f_{func}(z)
But where is the benefit for the program workflow by using this cumbersome spelling? Now, by using f_ {func}(z), we are able to apply any desired transition function to influence the stretch factor, if this function is normalized between 0 and 1 within distinct boundaries along the horn axis and therefore fulfills all of our previously defined constraints. I’m a big fan of using continuous mathematical functions in terms of the stretch factor transition. This provides smooth progression and avoids abrupt changes for the horn profile. Finally, that was one of the major design goals. I would like to give now a few examples of how some of the useful transition functions look like. This selection is of course a rather subjective assessment and many other functions are also possible which leaves room for future fine tuning. The simplest case arises by using a linear function:
\tag{16a} f_{s, lin}(z) = \frac{z}{z_{mouth}}; \quad z \leq z_{mouth}
\tag{16b} f_{s, lin}(z) = \frac{z_{max}-z}{ z_{max}-z_{mouth}}; \quad z > z_{mouth}
Furthermore, the already mentioned dependence on the respective horn radius r_m at the current z value where r_t is the throat radius and r_k the construction vector of the spherical wave horn:
\tag{17} f_{s, r_m}(z) = \frac{r_m(z)-r_t}{r_k-r_t}
However, two trigonometric functions have also proven themselves, for example:
\tag{18a} f_{s, sin}(z) = \sin\left(\frac{z}{z_{mouth}}\cdot \frac{\pi}{2}\right); \quad z \leq z_{mouth}
\tag{18b} f_{s, sin}(z) = \sin\left(\frac{z_{max}-z}{ z_{max}- z_{mouth}} \cdot \frac{\pi}{2}\right); \quad z > z_{mouth}
\tag{19a} f_{s, tanh}(z) = \tanh\left(\frac{z}{z_{mouth}} \cdot \pi\right); \quad z \leq z_{mouth}
\tag{19b} f_{s, tanh}(z) = \tanh\left(\frac{z_{max}-z}{ z_{max}- z_{mouth}} \cdot \pi\right); \quad z > z_{mouth}
In addition, the squares of these trigonometric functions shown are also well working alternatives. The overall behavior of the proposed stretch functions can be seen in the following picture:

The functions are derived from the spherical wave horn formulas, so the normalization takes place at maximum roll-back of the profile here. What we see in the diagram may seem a bit confusing because as coordinate for the abscissa the top of the wave front was selected. Of course, the resulting horn profile has the well-known rollback. What we can recognize very quickly is the linear function as straight black line with the distinct knee. I have only incorporated the linear function as reference as it divides the curves in above and below the linear behavior. For all curves the minimum value is always 0.0 at mouth and 1.0 at maximum roll-back which is exactly what was intended. Now, one or the other reader might wonder why the curves at the maximum z-value go back to zero? This is easy to explain because the profile of the spherical wave horn closes somewhere behind the throat and therefore the function must turn back again to zero. With help of some initial BEM simulations an assessment was done and my favorites are currently the trigonometric functions, especially the dashed curves. These functions stretch the horn profile above the linear behavior, so a more pronounced widening of the major axis occurs at the throat region which leads to a higher preference for horizontal radiation and that’s exactly what was intended. The following example was taken from an 400 Hz cut-off spherical wave horn calculation where a stretch factor of 1.5 and the pure \sin stretch function (16a, 16b) was applied. The next two figures show, first the resulting flat horn profile:
and second the same profile with stereographic projection applied and a coherent radius:

It should be obvious that the blue curves represent the major axes and the yellow curves the minor axes of the ellipse. The black vertical dashed line marks the mouth plane of the round horn. The circle with center at 0.0 is nothing else than the projection sphere with coherent radius exactly at the termination condition of the horn profile which was the maximum horn length in this case.

If it were just a question of planar wave fronts, the process would be described as rather simple. Now, however, the Kugelwellenhorn (spherical wave horn) is based on the assumption of curved wave fronts, which expand exponentially. This implies that already at the throat the wave front is assumed to be slightly curved. The initial wave front surface area S_ {0,round} must therefore be calculated as a dome surface or spelling different as surface area of a spherical cap:
\tag{20} S_{0,round} = 2 \pi \cdot r_{k} \cdot h_0
Where h_ {0} is the height of the dome and r_ {k} is the construction vector of the horn. If it were to stick with the pure round shape, then the math would remain quite simple as just shown. For an elliptical shape, however, the surface of the wave front becomes an triaxial ellipsoid! The calculation of it’s surface area, however, no regular formulas have been found in mathematics so far. For a quite long period I thought about if it is worth to spent more time into this issue until I had familiarized myself with the method of numerical integration. So, if we do not have any analytical formulas to calculate the surface area, it can be estimated by a numerical integration using elliptical coordinates. See here for a detailed description of the equations (Ref. 5). The following figure gives a more insight view of the coordinates and parameters used for the numerical integration process:

In grey color we see the cut surface through the ellipsoid. Everything above this surface has to be numerically integrated with help of the following formula. Within the program code flow a double loop is doing this job:
\tag{21} S_{ecap} = \int_{\theta=0}^{\theta_{max}} \int_{\phi=0}^{2\pi} \left|\sin{\theta}\right| \cdot \\ \sqrt{r_c^2 \cdot \sin{\theta}^2 \cdot \left( r_b^2 \cdot \cos{\phi}^2 + r_a^2 \cdot \sin{\phi}^2 \right) + r_a^2 \cdot r_b^2 \cdot \cos{\theta}^2} \cdot \partial\phi  \, \partial\theta
\tag{22}r_c(z) = r_{swh}(z)
\tag{23}r_b(z) = \frac{r_{swh}(z)}{f_s(z)}
\tag{24}r_a(z) = \frac{r_{swh}^2(z)}{r_b(z)}
\tag{25}\theta_{max} = \arccos\left(\frac{1-h}{r_c(z)}\right)
The challenge here is, first, to have sufficient computer power available and, on second, the definition of suitable control parameters, in particular for the incremental step size of the angles during numerical integration. The best way is of course, if there is a reference available to which the result can be adjusted with respect to the desired numerical accuracy. This reference is not surprisingly easy to find, namely the spherical surface area for which we have an analytical expression (9). So, by selecting a round horn profile and by calculating the spherical surface cap analytically and then with help of the numerical integration procedure and by comparing both results the step sizes for angular increments can be adjusted to ensure the desired accuracy.

Up to this point, the transformation of the spherical wave fronts into ellipsoidal wave fronts has been seamlessly implemented in the sense of the original spherical wave horn theory. Although, the method of transforming a horn from round into an elliptical or rectangular shape is not new, i.e. Ref. 6, however, the desired flat rectangular mouth surface area along the horn axis was set equal to it’s corresponding flat round surface area. As far as I know, never before have ellipsoidal wave fronts been laid on the basis of a horn theory for the spherical wave horn.

Last bot not least I will describe shortly the overall workflow of my program applying a stretch factor and an optionally additional stereographic projection. Everything always starts with the calculation of a round profile following the given horn formula. Besides the numerical step size along the horn axis the most important input parameters that need to be provided to the program are the throat radius r_{t}, the horn cut-off frequency, the stretch factor f_{s,inp} and a projection factor f_{stp} which scales the projection base radius. Following the horn formulas by using polar coordinates we only need to calculate the radius r_m(z) to the wall at point z in this first step. With polar coordinates it is very easy to adjust how many points should be calculated around the circle simply by dividing 2\pi through the desired number of steps. Once we have determined the radius r_m(z) and intend to stretch the profile we do not need to calculate points on the circle. With stretching applied the next step is to calculate the ellipse minor and major axis but only as a first approximation by assuming the flat surface areas are equal. This is done with eq. (11) and (12). Now, we have a first set of start parameters for the ellipse minor and major axis. Wait please! Why do you say a first set of start parameters? Well, not to disappoint the reader but if we use these start parameters a_{start}, b_{start} and numerically integrate the surface area of the ellipsoid cap (21) the result differs from the spherical cap. We could also say it a little bit clearer that the result is wrong, not much wrong, but it is wrong if we want to follow the given horn theory consequently. The only solution is to do further number crunching and implement an iterative process that minimizes the error function we have to define:
\tag{26}F_{err} = S_{ellipse}(z)-S_{round}(z)
What we need now is a method to minimize this error function. A well know algorithm to accomplish this is the so-called Newton method (Ref. 7). I will describe it in more detail as it is used inside the program also for other tasks and has proven to be very valuable and efficient and generally converges in a few steps to the desired result. Before we loop in the iterative process the first error function value needs to be calculated with eq. (26). During the whole process the minor axis b(z) is held fix to the value calculated with eq. (11). All what we have to do is to find the fitting major axis a(z) that minimizes the error function. Following the Newton method the next optimized step is defined as
\tag{27}a_{n+1}(z) = a_n(z)-\frac{F_{err}}{F_{err}^{'}}
\tag{28}{F_{err}^{'}} = \frac{\partial}{\partial a}F_{err}
Eq. (28) is another mathematical nightmare! It is impossible to solve it analytically, at least for me. But assuming that we have enough computer power we simply solve this numerically by an approximation, also called as infinitesimal calculation. If we find a reasonable small enough step size s and set the desired point of the first derivative exactly in the mid of this interval then the following equation provides a reasonable good approximation of the first derivative:
\tag{29}{F_{err}^{'}} \approx \frac{F_{err}(a+\frac{s}{2})-F_{err}(a-\frac{s}{2})}{s}
The next major axis value a_{n+1}(z) can now be calculated withe eq. (27) and a new error function that we hope has become smaller than a pre-defined threshold. If the error function is still above the threshold we need another complete run through the procedure. It has proven that this method converges generally fast in a few steps but it is always a wise decision to implement a max loop counter because if the Newton method diverges in some rare cases then we end up in an infinite loop. Once the error function is below the defined threshold we have found the best fitting a_{n}(z) value that delivers the surface area of the ellipsoid cap we have searched for.

If additionally a stereographic projection has been requested, then we leave a little bit the path of the only truth as another Newton method on top of the already existing Newton method would be too complicated to implement and is no longer in proportion to the effort needed to implement this. I have done some test calculations with a special subroutine to numerically integrate the resulting curved mouth plane after the projection applied (with the center point gathered from the round mouth) and to my surprise the difference to the flat ellipse surface that we received as mouth surface without any projection applied was very small. What we know about the effect of the stereographic projection is that the overall horn profile encounters a slight compression as already described in a previous post. But this is very similar to just select a slightly smaller cut-off frequency, So, I do not worry much about that issue. The overall procedure described in this post will be used explicitly for the Kugelwellenhorn (spherical wave horn) and additionally for several other horns as control reference.

What I hope to conclude that I have not lost too many readers on the way to this point. If you are not very interested in mathematics and the background presented here, you can also wait for the next posts when the spreadsheets are presented and then simply use them to create your own horn profiles. So, stay tuned because besides the well know Kugelwellenhorn also some completely new horn profiles will be presented in forthcoming posts!