Design & Optimization of Large Cylindrical Radomes with Subcell and Non-Orthogonal FDTD Meshes Combined with Genetic Algorithms

: The word radome is a contraction of radar and dome . The function of radomes is to protect antennas from atmospheric agents. Radomes are closed structures that protect the antennas from environmental factors such as wind, rain, ice, sand, and ultraviolet rays, among others. The radomes are passive structures that introduce return losses, and whose proper design would relax the requirement of complex front-end elements such as ampliﬁers. The radome consists mostly in a thin dielectric curved shape cover and sometimes needs to be tuned using metal inserts to cancel the capacitive performance of the dielectric. Radomes are in the near ﬁeld region of the antennas and a full wave analysis of the antenna with the radome is the best approach to analyze its performance. A major numerical problem is the full wave modeling of a large radome-antenna-array system, as optimization of the radome parameters minimize return losses. In the present work, the ﬁnite difference time domain (FDTD) combined with a genetic algorithm is used to ﬁnd the optimal radome for a large radome-antenna-array system. FDTD uses general curvilinear coordinates and sub-cell features as a thin dielectric slab approach and a thin wire approach. Both approximations are generally required if a problem of practical electrical size is to be solved using a manageable number of cells and time steps in FDTD inside a repetitive optimization loop. These approaches are used in the full wave analysis of a large array of crossed dipoles covered with a thin and cylindrical dielectric radome. The radome dielectric has a thickness of ~ λ /10 at its central operating frequency. To reduce return loss a thin helical wire is introduced in the radome, whose diameter is ~0.0017 λ and the spacing between each turn is ~0.3 λ . The genetic algorithm was implemented to ﬁnd the best parameters to minimize return losses. The inclusion of a helical wire reduces return losses by ~10 dB, however some minor changes of radiation pattern could distort the performance of the whole radome-array-antenna system. A further analysis shows that desired speciﬁcations of the system are preserved.


Introduction
This work presents the design procedure of a large radome, to cover an array with a large number of antennas, up to 48, whose height is of the order of 24 wavelengths and with small elements for adaptation of the order of 0.002 wavelengths. We present the technique we consider most suitable, finite difference time domain (FDTD), with sub-cell The main issue in FDTD calculations is the correct introduction of boundary conditions. FDTD solves the Maxwell's curl equations and boundary conditions must describe the structure under study, otherwise the numerical results are not correct. Because typical FDTD technique uses a rectangular cell, curved boundaries are modeled with a staircase shape [25] approaching the true curved boundary. When curved boundaries have a small radius a large mesh is necessary. The use of curvilinear coordinates circumvents the staircase issue [29], however, the number of field components is duplicated with the introduction of covariant and contravariant components [33,34]. Furthermore, new matrices are needed to obtain contravariant components from covariant ones. The metric tensor associated to each curvilinear cell is the matrix that transforms covariant into contravariant components. The main advantage of curvilinear cells is the easy way to apply the boundary conditions, which are introduced in a simple natural manner [29].
The choice between a simple rectangular cell with simple updating FDTD equations and a sophisticated FDTD with a large number of matrices and memory allocation would depend on the shape and geometry of the boundary conditions of the structure under study. There is a trade-off between memory consumption versus expected accuracy, and this is strongly related to geometrical complexity of the modelled antenna. For example, planar antennas are always efficiently managed using rectangular cells [24], and also quasi planar antennas [25].
In this work, the FDTD in curvilinear coordinates is revisited, [16,34], in order to show that it is a powerful tool in the full wave analysis of the radome-antenna structure. The curvilinear FDTD permits the establishment of boundary conditions in arbitrarily shaped surfaces, in an easy and natural way. Moreover, we introduce non-orthogonal sub-cell features to include thin dielectric slabs and thin wires. These sub cell features, which are conveniently inserted in the FDTD algorithm, increase the efficiency of the FDTD implementations because they considerable reduce the number of cells required to give a proper resolution to the radome thickness and metal inserts used for its tuning.
With the proposed technique we face the three-dimensional electromagnetic modeling of a large radome-antenna-array system. The technique is employed in the analysis of an array of 20 to 48 crossed dipoles [35] separated λ/2, covered with a thin cylindrical dielectric radome. The thickness of the radome is negligible in terms of wavelength, the thickness below λ/10, and dimensions of the overall structure, whose height is around 24λ, 11.28 m at the central frequency of the working band. Although the proposed technique is powerful in the analysis of radome-antenna structures, there is still a major problem: the challenge to find the optimal metal insert that maximizes the overall transmission for a given frequency band. In our case, it is proposed to introduce a helical metal wire of diameter~0.002λ in the dielectric cover to cancel the capacitive behavior of the dielectric.
The FDTD code running on a Graphical Processing Unit (GPU) is fast enough to be combined with an optimization algorithm [36], and we use a genetic algorithm (GA) as optimization algorithm. GAs are optimization algorithms inspired by natural selection processes observed in nature [37]. GAs are a powerful tool when searching for an optimal or semi-optimal solution in complex optimization problems [38][39][40], typically used in phased array design [41], or circuits [42]. A large number of comercial software packages exist also as open sources, however is difficult to find reliable FDTD or finite element method (FEM) software with this numerical power: Curvilinear mesh, sub-cell features, and optimization procedures for design purpose. This is the case, for example, of Ansys using Finite Element Method (FEM), [43,44]. There is no optimization procedure, and it works frequency by frequency, in FDTD we obtain results in the whole band using the Fast Fourier Transform. In Ansys the design optimization is by hand, we have to modify the boundaries gently until we obtain the desired response, and finding a proper design is not easy [43]. The programs used are our own, programmed by us since almost the beginning of the FDTD technique in the last 30 years, in Fortran, C, and Matlab, these have been updated and improved over time [18]. This manuscript is organized as follows. Section 2 revisits the curvilinear FDTD technique. The objective of this section is to present the curvilinear FDTD technique in a clear and pedagogical way, providing the information of managed matrices and memory allocation. Section 3 presents the sub-cell dielectric features that are introduced in the original curvilinear FDTD technique. The sub-cell correction of static field in the proximity of a thin wire is presented in Section 4.
Section 5 presents the GA optimization algorithm to obtain the best parameters for minimizing return losses and, finally, results, analysis and discussion are presented in Section 6 for a Radome-Array of Crossed Dipoles system (RACD) working at UHF frequencies.

The Curvilinear FDTD
The curvilinear finite difference time domain technique, [12,13], is deployed in a general curvilinear mesh with non-orthogonal cells. The mesh is defined with coordinates (u 1 , u 2 , u 3 ), see Figure 1, specially developed for the structures under study, i.e., the boundary conditions are imposed along u 1 , u 2 , or u 3 lines or u i , u j , planes (i, j = 1,2,3). The cylindrical or spherical coordinates are special cases of the general curvilinear coordinates [12]. The curvilinear mesh can be related to an orthogonal Cartesian mesh through a mapping Q that establishes a relationship between the nodes (u 1 , u 2 , u 3 ) and a Cartesian coordinate system (x, y, z), that is (u 1 , u 2 , u 3 ) = (u 1 (x, y, x), u 2 (x, y, z), u 3 (x, y, z)). The discretization of the curvilinear mesh is similar to the one of orthogonal mesh, by setting the vertices of each cell along the (u 1 , u 2 , u 3 ) directions. The difference with orthogonal FDTD is the cell: in a curvilinear mesh the cell is a hexahedron, a volume enclosed in six faces with twelve edges and eight vertices. is not easy [43]. The programs used are our own, programmed by us since almost the beginning of the FDTD technique in the last 30 years, in Fortran, C, and Matlab, these have been updated and improved over time [18]. This manuscript is organized as follows. Section 2 revisits the curvilinear FDTD technique. The objective of this section is to present the curvilinear FDTD technique in a clear and pedagogical way, providing the information of managed matrices and memory allocation. Section 3 presents the sub-cell dielectric features that are introduced in the original curvilinear FDTD technique. The sub-cell correction of static field in the proximity of a thin wire is presented in Section 4.
Section 5 presents the GA optimization algorithm to obtain the best parameters for minimizing return losses and, finally, results, analysis and discussion are presented in Section 6 for a Radome-Array of Crossed Dipoles system (RACD) working at UHF frequencies.

The Curvilinear FDTD
The curvilinear finite difference time domain technique, [12,13], is deployed in a general curvilinear mesh with non-orthogonal cells. The mesh is defined with coordinates (u 1 , u 2 , u 3 ), see Figure 1, specially developed for the structures under study, i.e., the boundary conditions are imposed along u 1 , u 2 , or u 3 lines or u i , u j , planes (i, j = 1,2,3). The cylindrical or spherical coordinates are special cases of the general curvilinear coordinates [12]. The curvilinear mesh can be related to an orthogonal Cartesian mesh through a mapping Q that establishes a relationship between the nodes (u 1 , u 2 , u 3 ) and a Cartesian coordinate system (x, y, z), that is (u 1 , u 2 , u 3 ) = (u 1 (x, y, x), u 2 (x, y, z), u 3 (x, y, z)). The discretization of the curvilinear mesh is similar to the one of orthogonal mesh, by setting the vertices of each cell along the (u 1 , u 2 , u 3 ) directions. The difference with orthogonal FDTD is the cell: in a curvilinear mesh the cell is a hexahedron, a volume enclosed in six faces with twelve edges and eight vertices.   The Faraday's law circulates the electric field components along the edges of a face to calculate the time derivative of the magnetic field. The calculated magnetic field component is perpendicular to the surface enclosed by the circulation path. The electric field components along the edges are covariant components, and the magnetic field components are contravariant components. Contravariant components are perpendicular to the faces of the hexahedron, therefore are not aligned along the cell edges in general curvilinear meshes.
The Ampère's law circulates the magnetic field components along the edges of each face, these are the covariant magnetic field components, used to calculate the time derivative of the electric field orthogonal to circulation plane, in doing this, the contravariant electric field components are obtained.
The implementation of Curvilinear FDTD differs from rectangular FDTD in the use of covariant and contravariant components, inherent to a hexahedron cell. In rectangular FDTD, covariant and contravariant components are coincident because at every cell one of the edges is perpendicular to the surface defined by the other two edges. This is not the case when general hexahedron cells are used [29]. Twice the field components are necessary in curvilinear FDTD when compared with rectangular FDTD. However, with a curvilinear mesh the boundary conditions are set very accurately using a coarse mesh; and a coarse mesh is a smaller mesh that consumes lower computer resources. Thus, if the expected curvilinear mesh is smaller than half the equivalent rectangular mesh, the use of non-orthogonal cells is justified.
In curvilinear FDTD the covariant fields are needed to calculate the contravariant fields, and extra computer resources in terms of calculation and memory allocation are necessary in its implementation: Extra memory in curvilinear FDTD: -Covariant and contravariant fields at each cell node -Metric tensor g ij , is a 3ˆ3 symmetrical matrix defined each cell to carry out contravariant to covariant transformation -Calculation effort in terms of CPU/GPU in curvilinear FDTD -Contravariant to covariant transformation using the metric tensor g ij Contravariant fields must be down-converted to covariant ones and this requires a matrix operation [33] using the metric tensor of the mesh, g ij . These matrices are position dependent in a general curvilinear mesh and have to be calculated at each cell. The memory requirements are increased in comparison with a rectangular scheme. The choice between both techniques should analyze a trade-off accuracy + computer + complexity. There is an extra element to make the choice, the complexity in doing a curvilinear mesh to fit the device under study. These are the key factors to decide whether to use a very fine and large rectangular mesh, or a coarse and small, but complex, conformal mesh.
The mesh points are given in the mesh file. Indexing them (i, j, k) as a numbering of adjacent cells along the lines of the mesh, the (i, j, k) indexing follows the numbering of a Cartesian mesh, and the mesh points are defined as a mapping between the Cartesian index and the mesh points: ri " 1 . . . n 1 , j " 1 . . . n 2 , k " 1 . . . n 3 s Þ Ñ " Ñ r pu 1 pi, j, kq, u 2 pi, j, kq, u 3 pi, j, kqq The covariant basis vectors are obtained at each node of the FDTD mesh along the cell edges, or mesh lines (see Figure 1), And the contravariant vectors are normal to the FDTD cell faces.
Electronics 2021, 10, 2263 6 of 26 The metric tensor is: g is cell volume: The covariant and contravariant formulation of Maxwells's curls equations arises natural from original Ampère and Faraday-Henry Laws in integral form [16] in general curvilinear coordinates. The contravariant components of the electric field are obtained by projecting the curl on the Ñ A i , contravariant vectors: The two remaining components have similar expressions, and are obtained from Equation (6) by permutation of indices.
The contravariant (E i , H i ) and covariant (E j , H j ) field components are obtained by projecting the corresponding vectors on the corresponding basis vectors: Following Equation (7), is obtained the following expression involving the contravariant electric field components E 1 , and the H 2 , H 3 covariant components: Similar expressions are derived for the other field components E 2 , E 3 by index permutation.
From Faraday's law, the following expression is derived that relates H 1 with E 2 , E 3 .
Furthermore, similar expressions are derived for the other field components H 2 , H 3 , by index permutation.
Equations (10) and (12) are discretized using the centered difference scheme of FDTD, with the standard notation: The basis vectors are obtained, which lately are used to calculate g ij : Discretization of Equations (10) and (12) is as follows for E 1 , H 1 , and similarly for E 2´3 , H 2´3 by index rotation: To stagger the Ampere and Faraday laws in the numerical FDTD scheme, the conversion from contravariant to covariant components is required, and this conversion is carried out by means of the metric tensor g ij [33], where Einstein notational convention of summation over repeated indexes is followed: The FDTD-GCC implementation is summarized in the following steps: - An important parameter to keep the stability of the time marching FDTD is the choice of time step ∆t. Some derivations for the Courant condition in curvilinear coordinates have been addressed in our previous work [17,29]: min`ag mm pi, j, kq2 c , m " 1, 2, 3; pi, j, kq " p1 : n1, 1 : n2, 1 : n3q

Sub-Cell Modeling of a Thin Dielectric Slab with Losses
We assume a thin dielectric slab located along the surface u 2 -u 3 , (u 1 = constant). On this, the equations for the field components crossing the slab are modified. The model introduces a new and additional field node located at the position of the thin slab, adapted to the special cell [16]. The slab has thickness d in the u 1 direction (see Figure 2).

Sub-Cell Modeling of a Thin Dielectric Slab with Losses
We assume a thin dielectric slab located along the surface u2-u3, (u1 = constant). On this, the equations for the field components crossing the slab are modified. The model introduces a new and additional field node located at the position of the thin slab, adapted to the special cell [16]. The slab has thickness d in the u1 direction (see Figure 2). The thin material slab has conductivity σ, dielectric constant ε, and permeability μ = μ0. The above FDTD procedure is used for all the cells with exception of the special cells where the thin slab is located. In these special cells, the electric field normal to the slab, E 1 , is split in sub-components E 1 (i*, j, k) and E 1 (i0, j, k), to take in account the discontinuity in the dielectric. The i0, accounts for the line u1 = constant, where the thin slab is located, and i* is the index of the new field node located at the position of the thin slab. The other field components do not need to be split because they are tangential to the slab and are continuous functions. The update equations for E 1 (i0, j, k) involve H3 and H2 that are outside the dielectric. The E 1 (i*, j, k), located in the special nodes of the dielectric are updated using the same magnetic field components as E 1 (i0, j, k), because the H3 and H2 components are continuous across the dielectric. However, the closed path of Ampere's law used in the updating of E 1 (i*, j, k) is inside the dielectric. Therefore, when calculating E 1 (i*, j, k), dielectric constant changes and losses are considered, The updated equations for the electric fields tangential to the slab are obtained by considering that only a portion d of the enclosed path is inside the material, the factor f is defined as: E 2 , and E 3 , are obtained using the modified conductivity and permittivity σm, εm: The thin material slab has conductivity σ, dielectric constant ε, and permeability µ = µ 0 . The above FDTD procedure is used for all the cells with exception of the special cells where the thin slab is located. In these special cells, the electric field normal to the slab, E 1 , is split in sub-components E 1 (i*, j, k) and E 1 (i 0 , j, k), to take in account the discontinuity in the dielectric. The i 0 , accounts for the line u 1 = constant, where the thin slab is located, and i* is the index of the new field node located at the position of the thin slab. The other field components do not need to be split because they are tangential to the slab and are continuous functions. The update equations for E 1 (i 0 , j, k) involve H 3 and H 2 that are outside the dielectric. The E 1 (i*, j, k), located in the special nodes of the dielectric are updated using the same magnetic field components as E 1 (i 0 , j, k), because the H 3 and H 2 components are continuous across the dielectric. However, the closed path of Ampere's law used in the updating of E 1 (i*, j, k) is inside the dielectric. Therefore, when calculating E 1 (i*, j, k), dielectric constant changes and losses are considered, The updated equations for the electric fields tangential to the slab are obtained by considering that only a portion d of the enclosed path is inside the material, the factor f is defined as: E 2 , and E 3 , are obtained using the modified conductivity and permittivity σ m , ε m : Equations for E 2 , and E 3 are: The new field component E 1 (i*, j, k) is located in the cell fraction that contains the dielectric slab, and is introduced in the equations where this lies in the closed surrounding path of Ampere's law. Thus, the equation for H 1 is not modified (normal to the slab). However, H 2 , and H 3 , (tangential to the slab) include the new field component E 1 (i*, j, k), because it lies in the closed surrounding path that contains the slab.
The covariant fields are calculated with the corresponding metric tensor, whereas the special covariant E 1 (i*, j, k), is calculated using the metric tensor of E 1 (i 0 , j, k). This approach is consistent with our assumption of a thin material slab having a width much smaller than the cell dimensions.

Thin Wire Sub-Cell Model
Typical FDTD implementations model wires of finite radius as infinitesimal thin wires. Depending on the structure under study the infinitesimal modeling would give good results or not. Traditionally, the establishment of the hard boundary of zero tangential fields has a slight error of discretization that may lead to small numerical inaccuracies which mainly affect resonant structures, giving rise to small displacements in frequency with respect to experimental results. These small deviations appear in the FDTD analysis of wire structures, or narrow apertures in conductors. Moreover, resistive thin wires have known impedance that depends on the finite radius and conductivity of the wire, which is neglected in case of infinitesimal assumption. One way to take into account the finite radius of the wires, and improve the numerical FDTD results, is the incorporation of the known behavior of the electromagnetic field in the close vicinity of conductors. This is especially important in the treatment of thin wires [45][46][47][48][49], as well as in various matching structures to improve the behavior of the radomes. The electric field is perpendicular near a conducting surface, and the magnetic field rotates around the conductors, both fields having a dominant behavior 1/r, where r is the distance to conductor. This is a quasi-static behavior, and its introduction into Maxwell's equations allows the modeling of small conductive structures within the cells.
Let's assume a conductive thin wire of radius r 0 along the u 3 -axis, being r 0 much smaller than the cell dimensions. The thin wire is assumed to be in the u 3 -axis, where the E 3 (i = i 0 , j = j 0 , k) nodes are located ( Figure 3); therefore the typical boundary condition in the wire is E 3 (i 0 , j 0 , k) = 0. To take into account the finite dimension of the radius r 0 inside the curvilinear cell, being r 0 << a g 11 pi, j, kq, we assume a quasi-static behavior of the electric and magnetic fields in the proximity of the wire. Near the conductor tangential electric field and normal magnetic field approach 0, and normal electric field and circulating magnetic field approach static field behavior 1/r, where r is the distance to the center of the wire. Therefore electromagnetic fields can be written in terms of known field values and 1/r, dependence: The substitution of the above field dependence in the Maxwell's curl equations produces the following expressions, Equation (24) leads to the update equation for H 2 :  The substitution of the above field dependence in the Maxwell's curl equations produces the following expressions, BH 2 pi,j,kq n Bt¨?
Equation (24) leads to the update equation for H 2 : A similar expression is derived for the H 3 contravariant component. The covariant components are calculated from the contravariant components using the same Equation (14), with the metric tensor g ij .

Finding the Optimal Radome Design for a Large Antenna array
Full-wave analysis of large antennas covered by radomes is confronted with the coexistence of large structures with very small elements of reduced dimensions, and these small elements are critical because they are decisive in their operation. The design of a radome for an array of crossed dipoles in the frequency band [470 MHz, 806 MHz] is proposed. The radome is a cylindrical dielectric thin shell, made on fiber glass, dielectric constant ε r = 4.25, whose thickness is a fraction of the wavelength,~λ/10 in the frequency band, i.e., around 10-20 mm, which is chosen for mechanical and manufacturing reasons. The cylinder has 1600 mm diameter to comfortably contain the λ/2 crossed dipoles. To improve the transmission through the radome an inductive element is introduced, a thin metal wire, whose diameter is also a negligible fraction of the wavelength, which is around 0.002λ. The use of the technique presented in the previous sections seems perfectly indicated for the analysis of this type of structure, whereas the overall antenna-radome structure has a height of~24λ, and a diameter~λ.
The presented technique is powerful and useful in the analysis of the antenna-radome structure, now it remains to address another important problem; to find the optimal radome structure, within the imposed requirements, that will minimize the reflection, or return loss in the given operating band. In our case the objective is to introduce a helical silver wire, conductivity 6.14ˆ10 7 Ω´1¨m´1, having diameter around 0.002λ, that is~1 mm diameter. The helical wire will be wrapped in the external surface of the dielectric shell to cancel the capacitive behavior of the dielectric. The inductive behavior of the helical wire is expected to cancel the capacitive behavior of the dielectric shell. The goal is finding the optimum design that maximizes transmission, or conversely, to minimize return losses. Inside the radome, along the cylinder axis, there is an array of 20 to 48 crossed dipoles spaced λ/2. The crossed dipoles are transversal to the cylinder axis. Each pair of crossed dipoles is mounted horizontally attached to a common mast. The two dipoles are fed 90 degrees out of phase to generate an omnidirectional pattern. The array structure is expected to radiate broadside with a high gain a horizontal polarized field in the horizontal plane of the crossed dipoles. This antenna itself is expected to provide full coverage of the UHF commercial broadcast bands.

Radome Design Strategy Using GA
Three parameters must be obtained in the optimal design, which will be chosen from among those that give the best radiating behavior of the radome-antenna system. These parameters are: -Dielectric shell thickness, which will be chosen in the range 15-50 mm, for mechanical and manufacturing reasons. -Diameter of the silver metallic wire that will constitute the helix that envelops the radome. -Spacing between turns of the helical wire enveloping the radome.
We use a genetic algorithm (GA) as the optimization algorithm. GAs are optimization algorithms inspired by natural selection processes observed in nature. GAs have been successfully employed in a variety of optimization problems in electromagnetics [37], new algorithms GA are being developed [38], for complex optimization problems [39], in antennas are typically used for array design and optimization [41]. Other techniques could be implemented, but this one works reasonably well and we have tested it in other optimization problems [42].
The GA follows the following scheme: (1) Random production of an original population whose number of individuals is N.
(2) Produce the next generation by crossovers and mutations among the individuals.
(3) Produce the new population of N individuals from the generation of (2) by selection.
(4) Produce the next population by repeating steps (2) and (3) until obtaining the individual that satisfies the optimal conditions defined.
In the GA, the individuals of the next population are selected according to the value of the defined fitness function, which is calculated for each individual. In GA, the individuals whose values of the fitness function are higher than the others are kept, and individuals with low fitness are discarded. The value of the fitness function determines the ending condition of the GA search process [41]. If the values of the fitness function do not continue to improve, the algorithm is terminated and the individual with the highest fitness is presented as the optimal solution. Another ending condition of GA is the maximum number of generations, which is also an input. For example, a maximum of 500 generations in the populations can be set, thus another termination condition of the algorithm is when it reaches 500 generations. In that case the algorithm ends and the individual with the highest value of the objective function of the last generation is presented as the optimal solution.
In GA, crossover and mutations are carried out by means of binary 0-1 operations, therefore the coding of individual parameters is necessary expressing the individual parameters by binary strings of 0 s and 1 s. In this optimization procedure each Radome-Array of Crossed Dipoles system (RACD) is an individual having three parameters: Thickness, diameter, and spacing. These three design parameters are coded by three 8-bit string of 0 s and 1 s. That is, each parameter is a real number and is converted to binary format to set the three chromosomes of the individuals: The GA population has N individuals as above, each one with three the chromosomes, and each chromosome with eight genes. The objective is to minimize reflection within the operating bandwidth [f 1 , f 2 ] of the RACD. We look for a good transmission at the center frequency, which is f c = (f 1 + f 2 )/2. Thus, the fitness function is defined as the average transmission coefficient of the average frequency response in the band with the response at central frequency. In this fitness we give equal importance to the average frequency response over the entire band and the frequency response at the center frequency.
The objective to minimize will be the average reflection coefficient within the whole band, also giving importance to the center frequency.
The mesh for this optimization procedure is always 42ˆ36ˆ500 in u 1 -u 2 -u 3 respectively. The FDTD mesh is terminated with Liao Absorbing Boundary Conditions (ABC) [50] because of its simplicity and fast calculation in a curvilinear FDTD mesh. Although PML could be used [51] this are avoided in this modeling to reduce the time spent in the large number of numerical calculations involved in the GA. The reflection coefficient S 11 , return loss, is calculated by using the total-field reflected-field formulation with a synthetic excitation using the near field [21], covering the [f 1 , f 2 ] frequency band. The RACD system is designed specifically for UHF broadcast, and the manufacturing specifications are that the radome should minimally alter the original cross-dipole array pattern, especially outside the broadside direction.

FDTD-GA Optimization Scheme
The FDTD program is introduced inside the loop of the GA to calculate the fitness function. The optimization scheme is as follows: Compares individuals in the population to find the one with the best fitness F(S 11 ). This is compared against previous optimum value, and the best one between them is selected. Finds the individual with the worst fitness and this is substituted by the individual with the best fitness. With this process, the individuals with the worst fitness are eliminated from the population and replaced by those with the best fitness. The best individual is obtained by comparison when fitness converges or by the end of generations.
The individuals with the worst fitness are eliminated from the population and replaced by those with the best fitness. The program ends when the best individual is obtained or at the end of the number of generations. In the latter case the optimum value is the best of the final population, taking into account that this population has undergone natural selection and is the one with the highest fitness. Therefore, at the end of the number of generations we could say that the optimum is obtained from the best individual of the best population.
The FDTD-GA optimization scheme was implemented using Matlab TM R2021a, doing parallelization to take the capabilities of the Graphical Processing Unit (GPU), a NVIDIA GTX-970 GPU card. In doing this the FDTD-GA optimization scheme was able to give a useful design in terms of a few hours, [36].

Return Loss
The cylindrical dielectric cover is made on fiberglass, ε r = 4.25. Diameter is 1600 mm. The conductivity of the helical wire is σ = 6.14ˆ10 7 Ω´1¨m´1, and the other parameter to be found are thickness of dielectric cover, diameter of the wire, and spacing between turns of the helix, which are bounded in the following intervals given in meters: [0.01, 0.05]; [0.10, 0.50], and [0.1ˆ10´3, 1.0ˆ10´3] respectively. The radome completely covers the antenna array, with an extra length of half a wavelength below the first antenna of the array and half a wavelength above the last antenna of the array (see Figure 4).  The FDTD mesh is 42 × 36 × 500 in the directions (u1, u2, u3), which account for the (r, ϕ, z). The overall structure, including the helical wire, has not any symmetry because of the helix. The program runs for 6000Δt. Electromagnetic fields inside and outside of the radome are sampled, stored and processed. The excitation is a half sine function with 90° phase shift between both dipoles. The transmission is calculated as the ratio of the field with radome and without radome, and the return losses accordingly, by subtraction of incident and total field inside the radome. Some preliminary calculations are carried out, and Figures 5 and 6 show the magnitude of the E-field with the radome at a given time instant. Return losses are shown in Figures 7 and 8   The FDTD mesh is 42ˆ36ˆ500 in the directions (u 1 , u 2 , u 3 ), which account for the (r, φ, z). The overall structure, including the helical wire, has not any symmetry because of the helix. The program runs for 6000∆t. Electromagnetic fields inside and outside of the radome are sampled, stored and processed. The excitation is a half sine function with 90˝phase shift between both dipoles. The transmission is calculated as the ratio of the field with radome and without radome, and the return losses accordingly, by subtraction of incident and total field inside the radome. Some preliminary calculations are carried out, and Figures 5 and 6 show the magnitude of the E-field with the radome at a given time instant. Return losses are shown in Figures 7 and 8  radome are sampled, stored and processed. The excitation is a half sine function with 90° phase shift between both dipoles. The transmission is calculated as the ratio of the field with radome and without radome, and the return losses accordingly, by subtraction of incident and total field inside the radome. Some preliminary calculations are carried out, and Figures 5 and 6 show the magnitude of the E-field with the radome at a given time instant. Return losses are shown in Figures 7 and 8 for the frequency band under study [470 MHz, 806 MHz] for different values of the spacing. The return losses for the radome without helical wire are in the interval [−10.6, −6.5] dB, the wrapping helix with spacing s = 0.22 m improves return losses by more than 10 dB around 600 MHz, and 3-15 dB in the operating band. This shows that the procedure is adequate, and that it is now feasible to search for an optimal solution.      The program runs comfortably and fast on a pc with a NVIDIA GeForce GTX 970 GPU card, i7-8700 CPU at 3.2 GHz, and 32 Gb of RAM, windows 10, 64 bits.
The optimization scheme is set as follows. The GA is set to run 500 generations working with a population of 100 individuals: Individuali=1..100 = (Thickness, diameter, spacing)i=1..100i = (00101111, 01010101, 11110000) The probability of mutation is set pm = 0.8, based in previous work and trials [42]; and the fitness function is set to the return losses as in (26) evaluated from the tangential fields. Figure 9 shows the fitness versus generation. The fitness function as Equation (29)   The program runs comfortably and fast on a pc with a NVIDIA GeForce GTX 970 GPU card, i7-8700 CPU at 3.2 GHz, and 32 Gb of RAM, windows 10, 64 bits.
The optimization scheme is set as follows. The GA is set to run 500 generations working with a population of 100 individuals: 100i = (00101111, 01010101, 11110000) The probability of mutation is set pm = 0.8, based in previous work and trials [42]; and the fitness function is set to the return losses as in (26) evaluated from the tangential fields. Figure 9 shows the fitness versus generation. The fitness function as Equation (29)  However if we re-define the fitness as the best behavior at the central frequency: It achieves fitness = −26.9 dB, see Figure 11, and obtained design parameters are: -Thickness = 0.0100 m The GA can be modified introducing changes in fitness, to analyze differences, for instance to take in account only the average transmission in frequency band [470 MHz, 806 MHz] with the fitness: The GA achieves fitness =´21.80 dB, see Figure  However if we re-define the fitness as the best behavior at the central frequency: It achieves fitness = −26.9 dB, see Figure 11, and obtained design parameters are: -Thickness = 0.0100 m However if we re-define the fitness as the best behavior at the central frequency: It achieves fitness =´26.9 dB, see Figure 11, and obtained design parameters are: Convergence is observed to the same optimal values in any chosen fitness function. Results obtained with three different objective functions are shown in Table 1, and optimal parameters have negligible differences. Inside each population at each generation there is a large variability in fitness performance, as is shown in Figure 12, however there is a deep minimum for the fitness function when convergence is achieved. In the population of Figure 12 fitness oscillates between −5 dB and −27 dB. Obviously we look for the deepest minimum that is −27 dB. Convergence is observed to the same optimal values in any chosen fitness function. Results obtained with three different objective functions are shown in Table 1, and optimal parameters have negligible differences. Inside each population at each generation there is a large variability in fitness performance, as is shown in Figure 12, however there is a deep minimum for the fitness function when convergence is achieved. In the population of Figure 12 fitness oscillates betweeń 5 dB and´27 dB. Obviously we look for the deepest minimum that is´27 dB. We observe a tendency of the optimization procedure to minimize the thickness of the dielectric material and the wire diameter. The spacing between turns of the helix is adapted to these parameters maintaining a value close to 0.15 m.
If the boundaries are changed to increase the plastic cover thickness, for example in the interval [0.025, 0.05] m; the optimization chooses the smallest thickness. Therefore the limiting factor in the optimization is the dielectric thickness. Convergences with several bounds which are presented in Table 2 are shown in Figure 13. We observe a tendency of the optimization procedure to minimize the thickness of the dielectric material and the wire diameter. The spacing between turns of the helix is adapted to these parameters maintaining a value close to 0.15 m.
If the boundaries are changed to increase the plastic cover thickness, for example in the interval [0.025, 0.05] m; the optimization chooses the smallest thickness. Therefore the limiting factor in the optimization is the dielectric thickness. Convergences with several bounds which are presented in Table 2 are shown in Figure 13.    Figure 13. Fitness versus generation. Fitness defined in (29). Bounds in GA search of Table 2, cases (1)-(2)-(3)-(4)-(5)-(6)of Table 2.

Radiation Pattern
The RACD system is designed specifically for UHF broadcast, and the manufa specifications are that the radome should minimally alter the original cross-dipol pattern, especially outside the broadside direction. There is concern about the pos of some degradation of the radiation pattern due to the presence of the helica Figure 13. Fitness versus generation. Fitness defined in (29). Bounds in GA search of Table 2, cases (1)-(2)-(3)-(4)-(5)-(6)-(7) of Table 2.

Radiation Pattern
The RACD system is designed specifically for UHF broadcast, and the manufacturing specifications are that the radome should minimally alter the original cross-dipole array pattern, especially outside the broadside direction. There is concern about the possibility of some degradation of the radiation pattern due to the presence of the helical wire. Although the helical wire improves the transmission, it works as a passive parasitic element which induced currents may influence the final radiation pattern of the whole RACD. The helical wire would operate as a monofilar helix antenna, and helix antenna may have two radiating operating modes: Axial mode and normal mode. Usually helix antenna operates in axial mode, in which it radiates along its axis circular polarization. This is the usual operating mode of the helix antenna for which it is designed. In normal mode the radiation is normal to the axis. From the original work of Kraus [52,53], we can guess that our design is correct. The normalized circumference of the helix is C λ = 2π r/λ = 10.7, and the optimal normalized spacing s λ = s/λ is between 0.2-0.6, therefore is out of the axial mode region, very far on top of the line C λ = 2 ' (s λ + 1). The pitch angle = atan(s λ /C λ ) is between 1˝-3˝, therefore the absence of an axial mode is almost guaranteed. The helix works similarly as an array of circular loops, which are expected to radiate horizontal polarization in the direction normal to the axis. However we make some numerical calculations to confirm our hypothesis and analyze any undesired effect.
To analyze any far field pattern distortion we calculate the radiation pattern (RP) of the array with and without radome to analyze differences. Induced currents in the helical wire work as a helical antenna, thus re-radiating and may be distorting the far field of the array of crossed dipoles. The current in the helical wire can be calculated from the incident near field from the FDTD simulation. The tangential electric field induces superficial skin currents which depend on the wire diameter, conductivity and frequency. The near field is FDTD calculated at each point of the helical wire tangential to the wire, and the near to far field calculation [53] is carried out at 638 MHz at the center of 470-806 MHz band. Figure 14 shows the amplitude of induced current versus path along the discrete points of the helical wire using the optimum parameters of Table 2-(1). The current along the helix oscillates~10% and the number of oscillations is the nearest integer to the number of antennas of the RACD divided by three. The radiation of a RACD of 48 antennas is presented in Figure 15. The RACD is expected to be used in terrestrial broadcast radiating broadside, at 90˝, horizontal polarized electric field (E-field), and parasitic radiation from currents induced in the helix may be a problem. In Figure 15 is shown the RP of the array without helix, the RP of helix induced currents and total RP. The shown RP s are H-plane, i.e., the vertical plane that contains the axis of the RACD that will be vertically mounted. The RP's start from 0˝in the end-fire direction (upside direction), end at 180˝in the ground direction. The RP of the array is symmetrical, and the RP of the helix is not, as array geometry is symmetrical in the horizontal plane, and the helix is not. Although the current has a symmetrical distribution along the wire, the position of each small, discrete current, along the wire is not symmetrical. The helix in the radome covers the whole array starting at λ/2 below the first antenna and ending λ/2 on top of the last antenna. The helix RP has three lobes, the "main lobe" at 90˝, but also has two lobes 1.3 dB higher than the main lobe at 35˝and 145˝having a width of 8˝from maxima. The helix has no axial mode [48], as we expected, but has three radiating lobes. The results of Figure 15 are normalized to the maximum radiated power, and although the helix radiates two undesired lobes, the power of them is 7.8 dB below the power radiated by the array, and its contribution to the total RP is low, therefore RP distortion is negligible. The total RP (array + helix) is shown in blue color in Figure 15, the total RP overlaps with the desired RP from the array (in green color) having a negligible deviation from the desired RP at 35˝and 145˝.
expected, but has three radiating lobes. The results of Figure 15 are normalized to the maximum radiated power, and although the helix radiates two undesired lobes, the power of them is 7.8 dB below the power radiated by the array, and its contribution to the total RP is low, therefore RP distortion is negligible. The total RP (array+ helix) is shown in blue color in Figure 15, the total RP overlaps with the desired RP from the array (in green color) having a negligible deviation from the desired RP at 35° and 145°.   Finally, the RACD is intended to operate horizontally polarized E-field at 90°, the field originally radiated by the array of crossed dipoles at 90°. In Figure  presented the RP of horizontal polarized field (Eφ), and polarized field Eθ wi without radome. The helix has a clear influence on polarization. This is due to the v component of the current, due to the vertical projection of the helical wire. The does not radiate pure horizontal polarization, however at ±7°, from the maximu vertical polarization is 10 dB below the horizontal polarization, at 90° is 30 dB dow Finally, the RACD is intended to operate horizontally polarized E-field at 90˝, that is, the field originally radiated by the array of crossed dipoles at 90˝. In Figure 16 are presented the RP of horizontal polarized field (E ϕ ), and polarized field E θ with and without radome. The helix has a clear influence on polarization. This is due to the vertical component of the current, due to the vertical projection of the helical wire. The RACD does not radiate pure horizontal polarization, however at˘7˝, from the maximum the vertical polarization is 10 dB below the horizontal polarization, at 90˝is 30 dB down, for the intended purpose of the RACD the improvement in return losses that provides the helical wire outperforms the losses in power by undesired cross polarization.
The results for a RACD of 20 antennas are also presented to provide a better graphical view of other aspects of the RACD performance, because the number of lobes increases with the number of antennas, making visual inspection of the plots difficult. We simulate the RACD with 20 antennas to provide a clear view of the influence of the helical wire in the radiation pattern, also to clearly see the main lobes and nulls for interpretation. The radome starts at λ/2 below the first antenna and ends λ/2 on top of the last antenna for any number of antennas, therefore depending on the number of antennas we have a different height and volume of the whole structure. Results are shown in Figures 17 and 18. The number of antennas changes the RP's. A smaller number of antennas lead to lower directivity. However, where it is noticeable is especially in the RP of the helix, where four secondary lobes arise in a pattern without symmetry. However, the levels of these lobes although relatively high within the helix RP have a very low power with respect to what the dipoles radiate. The effect of the helix is negligible, on the order of 20 dB below the dipole radiation, which is why in Figure 16 the green line is superimposed on the blue line, and blue line is not visible. The RACD of 20 antennas does not radiate horizontal polarization, with exception of 90˝, as in previous case, however at˘7˝, from the maximum the vertical polarization is 10 dB below the horizontal polarization, at 90˝is below 30 dB.
Electronics 2021, 10, x FOR PEER REVIEW 2 the intended purpose of the RACD the improvement in return losses that provid helical wire outperforms the losses in power by undesired cross polarization. The results for a RACD of 20 antennas are also presented to provide a better gra view of other aspects of the RACD performance, because the number of lobes inc with the number of antennas, making visual inspection of the plots difficult. We sim the RACD with 20 antennas to provide a clear view of the influence of the helical w the radiation pattern, also to clearly see the main lobes and nulls for interpretatio radome starts at λ/2 below the first antenna and ends λ/2 on top of the last anten any number of antennas, therefore depending on the number of antennas we h different height and volume of the whole structure. Results are shown in Figures 1  18. The number of antennas changes the RP's. A smaller number of antennas lead to directivity. However, where it is noticeable is especially in the RP of the helix, wher secondary lobes arise in a pattern without symmetry. However, the levels of these although relatively high within the helix RP have a very low power with respect to the dipoles radiate. The effect of the helix is negligible, on the order of 20 dB belo dipole radiation, which is why in Figure 16 the green line is superimposed on th line, and blue line is not visible. The RACD of 20 antennas does not radiate hori polarization, with exception of 90°, as in previous case, however at ±7°, fro maximum the vertical polarization is 10 dB below the horizontal polarization, at below 30 dB.    If we analyze the Table 2, depending on the bounds, the GA may converge t unacceptable return losses. This is the case of Table 2-(2), return losses perform worse. Th results of introducing the helix worsen the performance of the radome, its introductio If we analyze the Table 2, depending on the bounds, the GA may converge to unacceptable return losses. This is the case of Table 2-(2), return losses perform worse. The results of introducing the helix worsen the performance of the radome, its introduction does not have any advantage with respect to the only-dielectric radome. This bad performance also appears in the RP's and polarization, shown in Figures 19 and 20. The helix has also a stronger contribution to the total RP (see Figure 19), changing the null distribution, but also increasing 10 dB the radiated power to 0˝and 180˝. Polarization is shown in Figure 20, the results in circles plot the polarization without helix, and the solid lines with helix. The results for horizontal polarization E 2 = E ϕ , overlap in the interval 60-120˝, with negligible differences out of this. The E θ polarization has slightly more differences, in the interval 50-70˝that increase to almost 5 dB near 0˝and 180˝, to ground and zenith of the RACD. does not have any advantage with respect to the only-dielectric radome. This bad performance also appears in the RP's and polarization, shown in Figures 19 and 20. The helix has also a stronger contribution to the total RP (see Figure 19), changing the null distribution, but also increasing 10 dB the radiated power to 0° and 180°. Polarization is shown in Figure 20, the results in circles plot the polarization without helix, and the solid lines with helix. The results for horizontal polarization E2 = Eφ, overlap in the interval 60°-120°, with negligible differences out of this. The Eθ polarization has slightly more differences, in the interval 50°-70° that increase to almost 5 dB near 0° and 180°, to ground and zenith of the RACD.  The helix increases E3 = Eθ polarization, this is due to the current induced alo axial direction in the helix. However, in every case, in a beam of 14 degrees wide th polarization is 10 dB down, also, in a very narrow beam of a few degrees in bro direction this level is below 30 dB, i.e., at broadside direction the horizontal polar is kept.
Obviously, the helix has influence on the radiation pattern of RACD, althou influence is negligible for near broadside radiation, and is outweighed by the ben better matching of around 10 dB in the whole frequency band. However, a careful a is necessary to avoid the appearing of undesired side lobes that may influenc performance.

Conclusions
In this work we address a complex electromagnetic problem, the desig optimization of a large radome to protect a large array of crossed dipoles. The numerical technique using general curvilinear coordinates seems appropriated combined with a search optimization algorithm, as it is the genetic algorithm. We u cell features inside the curvilinear FDTD to give proper modeling of thin dielectr thin wires. Numerical FDTD runs inside GA are feasible using parallelization in G helical wire is used to give a better matching to the dielectric cover. The pro structure to cover the array of crossed dipoles is technically feasible by imp transmission around 10 dB, however, this introduces low levels of cross polarizatio undesired side lobes, depending on the parameters. Further analysis shows that intended use of the RACD system, which is to provide high gain at broadside di with horizontal polarization, these are minor effects. The benefit of introducing outweighs the detriments of cross-polarization and side lobes by doing a rigorous a and striving for an optimized design.  The helix increases E 3 = E θ polarization, this is due to the current induced along the axial direction in the helix. However, in every case, in a beam of 14 degrees wide the cross polarization is 10 dB down, also, in a very narrow beam of a few degrees in broadside direction this level is below 30 dB, i.e., at broadside direction the horizontal polarization is kept.
Obviously, the helix has influence on the radiation pattern of RACD, although this influence is negligible for near broadside radiation, and is outweighed by the benefit of better matching of around 10 dB in the whole frequency band. However, a careful analysis is necessary to avoid the appearing of undesired side lobes that may influence total performance.

Conclusions
In this work we address a complex electromagnetic problem, the design and optimization of a large radome to protect a large array of crossed dipoles. The FDTD numerical technique using general curvilinear coordinates seems appropriated when combined with a search optimization algorithm, as it is the genetic algorithm. We use sub-cell features inside the curvilinear FDTD to give proper modeling of thin dielectrics and thin wires. Numerical FDTD runs inside GA are feasible using parallelization in GPU. A helical wire is used to give a better matching to the dielectric cover. The proposed structure to cover the array of crossed dipoles is technically feasible by improving transmission around 10 dB, however, this introduces low levels of cross polarization, also undesired side lobes, depending on the parameters. Further analysis shows that for the intended use of the RACD system, which is to provide high gain at broadside direction with horizontal polarization, these are minor effects. The benefit of introducing a helix outweighs the detriments of cross-polarization and side lobes by doing a rigorous analysis and striving for an optimized design.