Efficient Crosspolar Optimization of Shaped-Beam Dual-Polarized Reflectarrays Using Full-Wave Analysis for the Antenna Element Characterization

Efficient Crosspolar Optimization of Shaped-Beam Dual-Polarized Reflectarrays Using Full- Wave Analysis for the Antenna Element Characterization Daniel R. Prado, Manuel Arrebola, Marcos R. Pino, Rafael Florencio, Rafael R. Boix, José A. Encinar, and Fernando Las-Heras Abstract—A method for the optimization of the crosspolar component of dual-polarized reflectarrays using full-wave analysis at the element level is described and demonstrated. The reflectarray full-wave analysis is based on local periodicity (LP) and integrated within the optimization process in order to accurately characterize the crosspolar far field. The proposed method is based on the generalized Intersection Approach framework using the Levenberg-Marquardt Algorithm as backward projector, and the employed full-wave analysis is based on the Method of Moments assuming Local Periodicity (MoM-LP). Several strategies to accelerate the computations are exploited, such as the parallelization of all the algorithm building blocks. To minimize the impact of MoM-LP in the optimization process, a strategy to reduce the number of MoM-LP calls is described, further accelerating the algorithm. Moreover, the convergence is improved by working with the squared field amplitude, alleviating the trap problem of local optimizers. This method allows to optimize the crosspolar component in the whole visible region or only in the coverage zone to facilitate the convergence, reducing computing time and memory usage. Two test cases are provided to validate the technique, one with an isoflux pattern for global Earth coverage and another with European coverage for direct broadcast satellite application. Index Terms—Crosspolar optimization, direct broadcast satellite (DBS), full-wave analysis, generalized Intersection Approach, isoflux, Levenberg-Marquardt algorithm (LMA), reflectarray, shaped beam antenna. I. INTRODUCTION S INCE the first reflectarray antenna prototype based on waveguide elements introduced in 1963 [1] and its popularization in the 1980s due to the development of low-profile printed antennas [2], reflectarrays have demonstrated to be a feasible alternative to parabolic antennas, allowing to design high gain antennas with similar efficiency. Printed reflectarrays have also demonstrated their capabilities to produce contoured beams, for instance, for direct broadcast satellite (DBS) applications [3], where a single-layer reflectarray was used. However, single-layer reflectarrays present bandwidth limitations, although they can be overcome by other types of broadband elements, such as several layers of stacked rectangular patches [4], several parallel dipoles [5], [6], or other types [7]. Stacked rectangular patches were used to design reflectarrays with enhanced bandwidth for DBS applications [8], [9] with very stringent requirements. The most common technique for the design of shapedbeam reflectarrays consists in a phase-only synthesis (POS), which can be carried out with several algorithms, such as Intersection Approach [10], [11], Levenberg-Marquardt [12], or particle swarm optimization [13], among others. Then, the geometry of each cell is adjusted in order to fit the required phase shift previously synthesized. This technique has proved to be very effective to design reflectarrays that fulfill very stringent requirements [8], [9], but since it is a POS, it has no control over the crosspolar far field during the synthesis process. Some techniques for lowering the crosspolarization of reflectarrays have been described, which rely on positioning the elements in certain configurations [14] or by adjusting the elements dimensions forcing a null in the amplitude and a change of sign in the phase of the crosspolar reflection coefficients [15]. However, these techniques present limitations and cannot be used to synthesize shaped-beam reflectarrays. In [16], an extension of the Intersection Approach algorithm of [11] was developed to include crosspolar requirements in the synthesis while keeping the efficiency of the algorithm by using the fast Fourier transform (FFT) in both projectors. The result was a distribution of reflection coefficient matrices. The refiectarray geometry should be adjusted in such a way that the element behavior matched the required synthesized matrix, which can be difficult to accomplish. In [17], a similar approach to [16] was used, adjusting the refiectarray element geometry to match the reflection coefficient matrix by applying certain geometrical transformations to the unit cell described in [7]. Another technique for the crosspolar optimization of refiectarrays has been explained in [18]-[20] where the Intersection Approach is used along with the Broyden-Fletcher- Goldfarb-Shanno (BFGS) optimization algorithm. However, only single-polarized refiectarrays are considered with very few elements (225 in total), since the algorithm was slow. Finally, in [21]-[23], a general optimization technique for printed refiectarrays is presented and used, which employs a gradient minimax algorithm for the direct optimization of the refiectarray geometry. It allows to synthesize both copolar and crosspolar components using several degrees of freedom. Nevertheless, it uses a lookup table of scattering parameters to speed up computations during the optimization process, which needs to be generated for each substrate and unit cell. However, as the number of optimizing variables is increased, the time required to generate the table rapidly increases as well. Furthermore, to create the table some variables, which are not needed during the optimization process (i.e., they are not optimized) are considered, such as the pair of incident angles, since for a given refiectarray geometry, varying the incident angle can modify the reflection coefficients [2]. Also, because interpolation is used to calculate the scattering parameters from the lookup table, accuracy decreases. This fact is important, since the crosspolar behavior is highly nonlinear and its absolute values are very low with regard to those of the copolar component. In this paper, a new alternative technique for the crosspolar optimization of dual-polarized refiectarrays is presented. The refiectarray is modeled using a full-wave analysis based on local periodicity (LP) and the spectral domain method of moments LP (MoM-LP) [24]-[26] to analyze the antenna element. The full-wave analysis is integrated within the optimization process in order to accurately characterize the crosspolar far field, which is being optimized. The technique is based on the generalized Intersection Approach (generalized IA) framework of [27] and it uses the Levenberg-Marquardt Algorithm (LMA) as its backward projector. Several strategies are employed that allow the use of MoM-LP directly in the optimization process within acceptable computing time. These strategies comprise the parallelization of all the algorithm building blocks in order to scale well with the number of available processors, taking advantage of modern CPU units and workstations; and also strategies to minimize the number of MoM-LP calls. In order to improve the algorithm convergence, the squared field amplitude is used, which allows to obtain good results in less iterations in comparison when working with the complex field, as has been recently demonstrated in [28] where both approaches are compared. The algorithm can handle several thousands of optimizing variables, and they can be managed along with the number of far-field points to control computing time and memory usage. Two test cases Fig. 1. Geometry of a single-offset printed refiectarray. are provided (an isoflux pattern for global Earth coverage and a European coverage for DBS application), which show the capabilities of the developed technique, maintaining the copolar pattern shape within specifications while lowering the crosspolar levels for both polarizations at the same time. This paper is organized as follows. Section II provides a brief overview of the refiectarray analysis and the computation of the far fields. Section III details the algorithm employed in the crosspolar optimization along with some techniques to improve its performance. Section IV is devoted to the two test cases to validate the presented algorithm. Section V discusses the results obtained in Section IV. Finally, Section VI contains the conclusions. II. ANALYSIS OF THE REFLECTARRAY A. Field in the Aperture The scheme of an offset printed refiectarray is shown in Fig. 1. The refiectarray is illuminated by a feed, which generates an incident electric field, Emc, on the surface of the refiectarray. This tangential field can be expressed for both polarizations in compact notation as ÍI/F, . EÍZ Oc,y)x + E^Yy (x,y)y (1) where the superscripts indicate the feed polarization and the subscripts the component of the field with regard to the refiectarray coordinate system (e.g., EYX would be the x-component of the projected field over the refiectarray surface, in the refiectarray coordinate system, when the feed radiates in Y polarization, in the feed coordinate system). Note that each polarization has, in general, two components of incident field, the desired and the crosspolar. The crosspolar component of the incident field appears, since the feed is not ideal and also due to the projection of the field radiated by the feed onto the refiectarray surface. In a similar way, the reflected field on the refiectarray surface can be written as EX/Y ^ref (x,y) E*iYAx, y)x + Ef'Y(x, y) y. (2) Jref,x Jref,y The relation between the incident and reflected fields at each element (m,n) of the reflectarray is given by a matrix of reflection coefficients that characterize that element EX/Y where \Xm j yn ) R" Rmn-EfniY(xm,yn) ( ,mn ,mn \ rxx rxy 1 ,mn ,mn I ryx lJyy / (3) (4) and (xm, yn) are the coordinates of the (m, n)th element. The components of Rm" are complex numbers and fully characterize the behavior of the cell. This method of analysis considers three sources of crosspolarization: the feed, the projection of the incident field onto the surface of the reflectarray, and the elements through the use of Rm". The reflection coefficient matrix for each element is computed using MoM assuming LP (i.e., it embeds each element in an infinite array comprised of the same element) [24]-[26] and it approximately considers the losses due to the substrate, mutual coupling between elements, the field reradiated by the element (metallizations), and reflected by the substrate and ground plane. Out of the three crosspolarization sources, the Rm" matrix is the most important for the crosspolar optimization, since both the feed and incident field projection contributions remain fixed during the optimization, while the Rm" matrix changes for each element at each iteration of the algorithm. Hence, a correct characterization of Rm" using a full-wave technique based on LP, such as MoM-LP, during the optimization process is required for a correct prediction of the crosspolar far field. The LP assumption can be improved by an extended LP, which also considers the real surrounding elements [29], [30], although the extra computational burden introduced by this technique might not be suitable for optimization purposes [29]. In any case, in this paper, MoM-LP [24]-[26] will be used to compute Rm". This technique [25] has already been demonstrated by full-wave simulations of the whole antenna [31] and measurements [32], showing good agreement both in copolar and crosspolar patterns. On the other hand, the tangential reflected magnetic field is computed at each reflectarray element assuming a locally incident plane wave coming from the feed, using the following expression: ¿X/Y H. ref \Xm, yn) kref x ^ref \^m> 3W OJ/UQ (5) with m = 2K fo and the components of km{ kx = -ko sin 9mn cos y imn •\mn , —ko sin 9""' sin < +k0cos9mn mn mn (6) where (9mn,<pmn) is the incident angle of the plane wave for the (m, n)th element. km{ corresponds to the reflected propagative wave in the specular reflection direction, in the absence of grating lobes. In order to compute the cross product in (5), all three components of the electric field are needed, but (3) only provides two, namely, Ex and Ey. The Ez -component can be obtained through the plane wave relation £ref • -Eref = 0. Solving (7), it immediately follows that: ~kmi,x Elefx ~ %ef,)i Emfy Eieti kiei,z (7) (8) which allows to solve (5), obtaining the reflected tangential magnetic field on the reflectarray surface. B. Far-Field Patterns Once the tangential electric field has been obtained using (3), the radiated far fields can be efficiently computed using the FFT algorithm [2]. It has been demonstrated [32], [33] that the use of the First Principle of Equivalence allows a correct characterization of both the copolar and crosspolar components, in contrast to the Second Principle, which might not give accurate predictions for the crosspolar far field. The radiated far field by a planar aperture using the First Principle of Equivalence is [34] jk0e-Jkor 7X/Y EX/Y Anr nX/Y Px COS <p ,X/Y sin< - rjcosO ( Q jk0e-Jkor X/Y Anr (c sin i XIY » - Qy cos#>! X / Y cos ip + Qy sin cos 0 (l ,X/Y sin i 3z/y COS tp (9) where n = fioco is the vacuum impedance and Px, Py, Qx, and Qy are the spectrum functions, which can be calculated as the Fourier transforms of the tangential electric and magnetic fields in the aperture as PT{U,V) = (J E^xl(x,y)e^+°yUxdy x/y G*//(«,»)= HIxJYx/y(x,y)e^ux+vyUxdy (10) with u = sin 9 cos q>, v = sin 9 sin q> and S the reflectarray surface, which is already discretized into a regular grid of periodicity a x b. Please note that with the compact notation in (10), there are two spectrum functions per polarization, making a total of eight considering both fields. Considering the tangential fields constant within each element (m,n), the spectrum functions take the form P*'J(u, v) = K Y.Y. EfJiYx/y(x™> yn)' m n QxJ¡Yy{u, v) = K • Y, X HHí,x/y(x>n> y») ,jko(uxm+vyn) Jko(uxm+vy„) (ID (12) being K a constant . . /k0ua\ . /k0vb\ K = ab sine I I sine I I v 2 ; i 2 ; that accounts for the amplitude element pattern (when the array factor approximation is used, K = 1 [35]). The double sum in (11) can be written as a 2-D inverse discrete Fourier transform (IDFT2), which can be efficiently computed by a 2-D Inverse Fast Fourier Transform [2]. In that case, the spectrum functions are calculated as pX'/iu, o) = N-K. IDVT2[E^Yx/y(xm, yn)] Q*¡* (u,v) = N-K- IDFT2[7/rgx7/); (xm, y„)] (13) where N is the total number of reflectarray elements. Once the far field has been obtained in spherical coordinates using (9), the copolar (£Cp) and crosspolar (Exp) radiation patterns given by Ludwig's third definition [36] are obtained for both polarizations. The gain (Gq, and Gxp) is estimated by computing the total power radiated by the feed [2] (input power), which is usually modeled with a cos* 6 function [37]. The input power remains constant through the synthesis. This method of analysis considers the spillover and illumination efficiencies as well as substrate losses [2]. The far-field pattern specifications can be given in form of templates with top and bottom bounds [2], [10], [12], 7max(M, v) and Tmin(u, v), which can be defined in the whole UV grid or in the area of interest. These templates should be provided for copolar and crosspolar patterns for both polarizations and are usually given in gain [2]. In this regard, the phase variation in the far field is not of concern, even within the main beam, since the developed technique is for broadcast applications only. The control of the radiated farfield phase could be of interest in some synthesis processes, since the spatial dispersion produces different delays in receptors situated at different points within the shaped region. However, this fact does not affect broadcast applications, and thus, far-field phase variation over the coverage zone is not critical. III. CROSSPOLAR SYNTHESIS BASED ON IA-LMA In this section, a description of the generalized IA algorithm using the LMA as backward projector is provided. It directly optimizes the geometry of the elements using MoM-LP as analysis tool. For this reason, a special effort is made in order to accelerate the algorithm, implementing several strategies in the LMA. They are aimed to reduce the impact of MoM-LP operations through minimizing its use and parallelizing all the building blocks of the algorithm. The generalized IA considers two sets, the set of the radiation patterns that can be obtained with the reflectarray optics (set 1Z) and the set of radiation patterns that comply with the specifications (set M). At each iteration i of the algorithm, the following operation is performed: £¿+i = B[F(Ei)] (14) where E is the radiated far field by the antenna. T is the forward projector, which projects the far field radiated by the antenna (belonging to the 1Z set) onto the set of fields that comply with the specifications (M set). Conversely, B is the backward projector, which projects the far field that fulfills the specifications onto the set of far fields that can be radiated by the reflectarray. The goal of (14) is to find a radiation pattern that belongs simultaneously to both sets or, if that is not possible, to find a radiation pattern whose distance to the set of patterns that fulfill the requirements is minimal. Also, the generalized IA is a local optimizer and as such the starting point is of the utmost importance, since it will determine the correct converge of the algorithm [27]. The potential local minima reached by local optimizers are commonly known as traps. In the generalized IA framework, there are two main sources of traps [27]: the number of degrees of freedom, which in this case will be the number of optimizing variables; and the nonconvexity of the involved sets. The latter problem can be addressed by working with the squared field amplitude instead of the field itself. The gain is also suitable, since it is proportional to the squared field amplitude Gcc|£|2 (15) and the specifications are usually given in gain. The proportionality factor between the gain and squared field amplitude remains constant through the synthesis, so it does not change the degree of nonlinearity of the optimization problem. Also, as stated before, the reflectarray is analyzed, considering approximately the losses by means of a full-wave analysis of the unit cell and also the illumination efficiency and spillover [2]. Using the gain or squared field amplitude improves the convergence of the algorithm, since it convexities the set M of radiation patterns that comply with the specifications (this set is denoted as Yc in [27], where Yc is defined as convex). However, this leads to a redefinition of the backward projector in [11] and [16], since the FFT cannot be used anymore to recover the field at the aperture. This causes the algorithm to lose the efficiency derived from the FFT being the most time-consuming operation. Now, to recover the field at the aperture, an optimization algorithm is needed, which in the case of the present work is the LMA. For the first problem, the solution is to reduce the number of variables in the first steps of the synthesis process. This has been implemented in the reflectarray POS case with a fictitious reduction of variables modifying the illumination taper of the feed [11]. The idea behind the fictitious reduction of variables is that an array with a very low illumination at the edges is equivalent to using a smaller array, since the outer array elements barely contribute to the far fields, thus effectively reducing the number of variables, although, in reality, the number of variables is the same (hence the fictitious reduction) [11]. Nevertheless, this technique is not used here, since the reflectarray geometry is optimized with the LMA, and the reduction of the number of variables is not fictitious anymore, since one can choose which variables will be optimized and which will not. In order to alleviate the notation, the compact notation for both polarizations will be avoided and generic gains and templates will be used when possible. A. Forward Projection The forward projector is divided in two steps, the first one being optional. In the first step, the templates are normalized to the value of the gain in a given direction. Then, the fields are trimmed using those normalized templates, obtaining a far field, which complies with the specifications, but cannot, in general, be radiated by a passive refiectarray. The requirements of the copolar radiation pattern are usually given in gain, with a maximum and minimum value. If these templates remain unaltered (i.e., the first step is skipped), the synthesis is carried out infixed gain. This method can be useful to further refine a previous synthesis. Alternatively, they can be normalized to the gain in a given (uo, »o) direction, which belongs to the maximum gain region. It is then said that the synthesis is carried out in float gain [11]. This method is more useful when a synthesis is started from scratch (e.g., when doing POS), since it adapts automatically the templates to the maximum gain that can be obtained with the refiectarray. In this case, a normalizing constant is defined as GCp(uo, VQ) Cn (16) where Tav is the average value of the maximum and minimum copolar templates in the (UQ, DQ) direction 'cpjinax (Uo, Vo) + Tcpmin(Uo, VQ) (17) L C„ is defined for both polarizations independently to adapt the template to the current gain of each polarization (compact notation for both polarization was avoided to alleviate the notation, but Cf and CYn are obtained using the gains and templates from both polarizations). Then, it is applied to both copolar and crosspolar templates as T"p(u,v) = Tcp(u,u) • Cn T"p(u,v) = Txp(u,v) -Cn-(18) The normalization in (18) is applied to the minimum and maximum templates of both polarizations. This normalization process sets the middle point of the templates at (uo, »o) at the same level of the gain in that direction. Once the templates have been normalized, both components of the radiated field should accomplish the following condition in both polarizations: rcp,min(M>u) ^ Gcp(u,v) < T^maK(u,v) Txv.mm(u>v) ^ Gxp(u,v) < r" m¡DÍ(u,v). (19) This condition can be fulfilled by the Pr operator, defined for all (U,D) as follows for a generic gain G: PAG) rmax(M>u)> Tmax(u'v) < G(u,v) Tmm(U>v)> G(U,V) < T^{n(u,v) G(U,D), otherwise. The result of this operation is the gain G' G'{u,v) = Pr(G(u,v)). (20) (21) For the sake of simplicity, a generic G gain component and normalized template T" were used in (20) and (21), but they should be applied to the copolar and crosspolar components of both polarizations with their respective normalized templates. Also, the forward projector only works with the gain (which is proportional to the squared field amplitude), and the far-field phase does not play any role, being able to take any value. In contrast to [11] and [16], where the far-field phase remained unchanged, now there is no such constraint, improving the convergence of the algorithm. B. Backward Projection Since the forward projection works with the squared field amplitude (gain) instead of the field itself, the FFT to recover the field at the aperture cannot be used, since for the use of the FFT, the complex field (amplitude and phase) is needed to perform the inverse operation of (13), as in [10], [11], and [16]. In any case, since the goal is to directly optimize the geometry of the refiectarray to improve the crosspolar far field, an optimization algorithm will be used. Now, the B operator is defined as the minimization of the distance of an element m e M to the set 1Z dist(m, 11). (22) For the minimization of the distance in (22), at this stage, we have the element m as the trimmed gain G'(u,v) and the current refiectarray geometry, which generates a gain pattern that belongs to the 1Z set, G(u, v). As a distance definition, the Euclidean norm for square-integrable functions can be used, which can be easily implemented by the weighted Euclidean metric d = dist2(G'(ti, u), G(u, u)) w(u,v)(G'(u, D) - G(u,v))2du dv (23) where w(u,v) is a weighting function and Q is the area belonging to the visible region (u2 + v2 < 1) where the far fields are optimized. This area can be the whole visible region or located around a more limited zone, and can be used as a way of controlling memory and computing time resources allocated in the synthesis (i.e., a bigger Q area means more allocated computer memory and slower computations). Since the radiated fields are already discretized in the UV grid, the integral in (23) can be approximated by a sum for the UV points, which lie in Q d = ^ w(u, v)(G'(u, v) - G(u,v))2A u AD u,v = ^J_C{u, D)(G'(U, D) - G(u, D))]2 where C(u, D) = TJU)(U,D)AUAD (24) (25) and AM and Av are the steps in the discretized UV grid in u and D, respectively. On the other hand, the LMA minimizes cost functions of the form [12] M F(x) = YJ(n(x))2 (26) ¿=i where r(x) is known as residual, which is discretized in M points. Comparing (24) with (26), it is clear that r(u, D) = C(u, D)(G'(U, D) - G(u, v)) (27) so the LMA naturally minimizes the distance between the two gain patterns. The generic gains G and G' represent the copolar and crosspolar components for both polarizations, which are all optimized at the same time. Each far-field component is discretized into M points, hence having a total of AM points in the UV grid. Also, since the relative value of the crosspolar component is significantly lower than the copolar gain, in order to properly reduce its starting value, it might be convenient to scale the crosspolar residual by means of the weighting function contained in C(u,v). Finally, it is not necessary to attain a minimum (generally local) of (22) [27], only to decrease the distance at each iteration of the generalized IA, so very few iterations of the LMA are needed. C. Minimization Algorithm for the Geometry Optimization As stated in Section III-B, due to the similarity between (24) and (26), the LMA is a natural choice to minimize the distance in (22), and in fact will be the minimizing algorithm used here. However, other gradient-based algorithms are also suitable, such as steepest descent [38], Gauss-Newton [39], or self-scaled BFGS [40] by setting up an appropriate cost function [28]. The strategies devised to accelerate the computations of the LMA described below are also suitable for those algorithms, and in particular, also for the self-scaled BFGS, which has been used before in the context of reflectarray optimization [18]-[20], [28]. One advantage of the BFGS over the LMA is the fact that the BFGS does not perform a matrix multiplication (only matrix-vector multiplications) and it can also directly work with the inverse of the Hessian, so it does not need to solve an equation system. These two points allow the BFGS to save time with those operations. However, it presents the disadvantage that it must perform further cost function evaluations using the full-wave analysis tool in a line search to choose a proper step length [39], which wastes the time saved by the matrix multiplication and equation system solver. Moreover, the BFGS updating formula is considered the most effective quasi-Newton formula [39] when working in a standalone mode (i.e., not as part of a more general framework as in the present case). For the first iteration, the identity matrix is usually employed, which makes the first step equivalent to that of the steepest descent [39]. However, since only a few iterations of the gradient-based algorithm are necessary per iteration of the generalized IA [27], the Hessian approximation is continuously restarted and the convergence improvement of the BFGS due to the updating formula is lost. This is not a problem for the LMA, which uses always the same Hessian approximation and whose convergence is controlled by a real positive number. In any case, in this context is the generalized IA the algorithm controlling convergence [27] and the choice between the LMA and BFGS may have little impact in this fact. The LMA described in [12] is used here. It includes a number of strategies in order to improve the computing times of the algorithm. In short, the computation of the Jacobian matrix is parallelized, computing one column per available thread and minimizing the error of the finite difference by an appropriate choosing of its step (which depends on the machine epsilon [39]) to evaluate derivatives; the matrix-vector and matrix multiplications are performed by the Intel Math Kernel Library (MKL) library [41], which takes advantage of highly optimized, fully parallelized algorithms, and low-level hardware operations in order to improve their performance and computing time; a Cholesky factorization-based solver is used, which is the fastest exact solver for this type of problems (also from the MKL library). Also, the result of a big matrix multiplication involving the Jacobian is symmetric, so only a triangular part of it is computed, with a further reduction of computation time. However, in [12], the LMA was used to perform POS, and since now, the crosspolar far field is optimized by directly analyzing the reflectarray geometry with MoM-LP, more efficient strategies are included in the cost function and the Jacobian matrix evaluation (which are the building blocks calling the MoM-LP routine). 1) Cost Function: For the POS case, the most timeconsuming operation when computing the far fields from the tangential field of the reflectarray is the FFT However, now that the reflectarray is analyzed with MoM-LP, the evaluation of the tangential field is slower, since an MoM-LP call is required per reflectarray element. However, since the evaluation of a reflectarray element is independent of the rest (since the analysis assumes LP), this operation can be easily parallelized using OpenMP [42]. This way, the evaluation of the tangential field was sped up, and can take advantage of modern CPU units, scaling well with the number of available processors. This parallelization of the cost function routine only applies when it is called the first time at the beginning of each LMA iteration, and not from the Jacobian evaluation routine, where another kind of parallelization is used. 2) Jacobian Evaluation: Despite the improvement in the Jacobian evaluation described in [12], it is still quite slow when using MoM-LP in the analysis of the reflectarray. In order to reduce the number of calls to the MoM-LP routine, first, a one-sided-difference is used in the evaluation of each column of the Jacobian, instead of the central difference. This reduces the calls to the cost function (and hence the calls to the MoM-LP routine) by half. Still, when calling the cost function in order to evaluate the lateral difference, all elements of the reflectarray are processed with MoM-LP. Assuming that all reflectarray elements are optimized by the LMA, the Jacobian matrix has sN columns (where N is the number of reflectarray elements and s is the number of optimizing variables per reflectarray element), and considering MoM-LP as the dominant operation, the overall time cost of evaluating the full matrix would be O (OMOM-LP -S-N-N) (28) with OMOM-LP = O {kxZ2 + k2Z3) (29) where Z is the number of unknowns in MoM-LP, k\Z2 is the term for filling the MoM-LP matrix, and ^Z3 for the inversion. In the present case, five basis functions per dipole are used (see Section IV for a description of the unit cell), thus having a 40 x 40 MoM-LP matrix, which is very fast to invert. However, in this case, the bottleneck is the matrix filling, which is much slower than the inversion [25]. If L is the number of available processors, since the evaluation of the Jacobian matrix is massively parallelizable due to the fact that each column can be independently computed [12], computing times derived from (28) can be divided up to L times, considering a 100% parallelization efficiency, although in reality gains will be smaller. The time cost in (28) is not suitable for practical implementation. If analyzing one element with MoM-LP takes an average time of 0.1 s and the reflectarray has N = 900 elements with s = 6 optimization variables per element, the time to evaluate the Jacobian matrix would be of the order of 135 h (5.6 days) for L = 1 and 13.5 h for L = 10. If the reflectarray has roughly 7000 elements as in [9], the time would escalate to 340 and 34 days for L = 1 and L = 10, respectively. The computing times shown as described earlier are only a coarse approximation of the time cost to evaluate the Jacobian matrix, and illustrate the necessity of improving (28). In light of these results, it is obvious that the times involved in the Jacobian matrix evaluation are of no practical use. However, they can be greatly reduced by further minimizing the number of MoM-LP calls. Each column of the Jacobian is obtained by a single call to the cost function, which returns a vector with all the components of the column. Furthermore, each column is calculated by deriving the cost function with respect to just one variable, resulting in just one element being modified. Hence, there is no need to recompute Rmn for all elements using MoM-LP, but just one element per Jacobian column, reusing the rest of the reflection coefficients from the first call to the cost function at the beginning of each iteration of the LMA. This can be easily implemented by branching the cost function routine, detecting when it is used to evaluate the Jacobian matrix, and then computing only one element with MoM-LP. With this improvement, the new time cost of evaluating the Jacobian matrix is ^(OMOM-LP -S • N). (30) Using the same examples as before, the new improved computing times would be 9 and 70min for L = 1, which represent an important reduction from the previous 5.6 and 340 days, respectively. With these two optimizations in the computation of the cost function and the Jacobian matrix, apart from others introduced in [12], the LMA is able to directly optimize the geometry of large reflectarrays using MoM-LP for the crosspolar optimization in a reasonable amount of time, since the algorithm will scale well with the number of available processors, taking advantage of farm servers. Finally, the implemented algorithm is able to handle thousands of optimizing variables obtaining good results due to its improved convergence properties derived from the use of the squared field amplitude and the elimination of the far-field phase constraint (which now can take any value), provided a suitable starting point. Also, in order to control computing times and memory usage, the number of optimizing variables can be reduced as well as the number of points in the Fig. 2. Employed reflectarray unit cell based on parallel and coplanar dipoles in two layers of metallization. UV grid where the optimization is performed, carrying it out, for instance, only around the coverage area instead of in the whole visible region. IV. VALIDATION In order to validate the algorithm described in Section III-C2, two shaped patterns for satellite applications are optimized in order to lower their crosspolar levels while maintaining the copolar pattern within specifications. The patterns are an isoflux pattern for global Earth coverage [43] and a European coverage for DBS application [44]. In both cases, a POS is done using the IA [11], and then, a design is obtained following [2], using the reflectarray element shown in Fig. 2. The unit cell is comprised of two layers of metallizations with four parallel and coplanar dipoles for each polarization, and is described in detail in [32]. These two designs are used as starting point for the crosspolar optimization. The copolar template is the same as the one used in the POS, and the crosspolar template is defined as a certain constant value below the maximum copolar level, for the whole region where the optimization is carried out. A. Isoflux Pattern The first test case is an isoflux pattern for global Earth coverage. Reference [43] details how to obtain a parametrized isoflux pattern, which in this case has been particularized for P = 20°, a geostationary orbit, side lobe level 19dB below the maximum copolar level, and 0.35 dB of allowable ripple in the coverage area. The reflectarray is circular, with 1020 total elements placed in a rectangular grid with 36 elements along the diameter in both reflectarray axes. The feed horn is modeled as a cos^ 6 function [37] with a g-factor of 14.8, which provides an illumination taper of — 12dB at the reflectarray edges. The feed horn is placed at rf = (40,0, 195) mm with regard to the reflectarray center (see Fig. 1). The working frequency is 30 GHz and the period of the unit cell is 5 mm x 5 mm (0.5/i). The substrate for both layers of the unit cell is the same, 19 18 17 l 1 / .» /F W \ ***» <^^ ***-» * —^s*i • -»••• ^ 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 I \\ - - - X pol. (initial) — X pol. (optimized) - - - Y pol. (initial) — Y pol. (optimized) -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 (a) 19 18 17 <i l f —^N ^ *•... • >... ^*^:-*' -* k f 1 1 1 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 1 1 1 X pol. (initial) — X pol. (optimized) — - - Y pol. (initial) — Y pol. (optimized) * 1 * 1 ll If ^* **^ f -11 1 I : A -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 (b) Fig. 3. Main cuts for the copolar pattern for both polarizations before and after the crosspolar optimization. The copolar gain for Y polarization is increased during the optimization, (a) Cut along u for v = 0. (b) Cut along v for u = 0.34. with a height of 0.787 mm and a complex relative permittivity er =2.33-;3.029- 10"3. The optimization is carried out in the whole visible region, with a resolution of the far fields of 256 x 256 points for the FFT, having a total of 51543 points in the UV grid. All of the refiectarray elements are optimized at the same time, considering all dipoles as independent optimizing variables (s = 8), so the total number of optimizing variables (a) (b) Fig. 4. Cuts along u and v of the crosspolar far field before and after the optimization. The cuts are taken at the point where the crosspolar level is maximum for each pattern, (a) X polarization, (b) Y polarization. is 8160. The optimization is carried out in fixed gain, setting the required gain to 18.35 dB in the center of the coverage zone. The crosspolar template is set 35 dB below the maximum copolar gain template, and the crosspolar component is scaled by a factor of 100 in natural units with the weighting function. The LMA is set to perform three iterations per iteration of the generalized IA. The results obtained after the optimization are compared with the starting point in Figs. 3 and 4. On the one hand, Fig. 3 shows the results for the copolar pattern. There were two goals: to maintain the copolar specification (isofiux shape) while reducing the crosspolar component, and to increase the gain for Y polarization while maintaining the X polarization gain. As it can be seen, the starting point has around 0.4 dB less in gain for Y polarization than the other in the coverage area. However, after the optimization, the gains for both polarizations become similar while mostly preserving the copolar shape within specifications. On the other hand, Fig. 4 shows some cuts in (u, v) of the crosspolar patterns. The cuts correspond to the points were the -u = — sin V cos ip Fig. 5. European coverage, (u, u) are in the antenna coordinate system. crosspolar is maximum, for both polarizations, even though they are not the same for the patterns before and after the optimization. For the X polarization, the initial maximum crosspolar level is —4.81dBi, and is reduced by 8.36dB to —13.17 dBi after the optimization. For the Y polarization, initially there is a maximum level of —4.12 dBi for the crosspolar component and is reduced an amount of 9.05 dB to — 13.17 dBi. This important crosspolar reduction is achieved while keeping the gain of the copolar pattern (and slightly increasing it in the case of Y polarization) and its shape. The tradeoff is a slight increase in the side lobes, as seen in Fig. 3 and a redistribution of the lobes in the crosspolar component, which is now more uniform toward higher levels in the whole UV grid. B. European DBS Coverage The second test case corresponds to a European DBS coverage shaped-beam, as shown in Fig. 5 [44]. The working frequency is 11.85 GHz, and the satellite is placed in a geostationary orbit in position 10° E longitude, 0° latitude. The minimum gain specified in the coverage area is 28 dBi, which has been enlarged to consider typical pointing errors (0.1° in roll, 0.1° in pitch, and 0.5° in yaw). The coverage shown in Fig. 5 is specified in the antenna coordinate system, and before carrying out the synthesis, it needs to be transformed to the reflectarray coordinate system defined in Fig. 1 [2]. The reflectarray is square and formed by 5180 elements in a regular grid of 74 x 70 cells. The feed horn is modeled as a cos^ 0 function [37] with a g-factor of 23, which provides an illumination taper of — 17.9dB at the reflectarray edges. The feed horn is placed at r j = (358, 0, 1070) mm with regard to the reflectarray center. The period of the unit cell is 14 mm x 14 mm, which is 0.55/1 at the working frequency. Now, the substrate for the bottom layer has a height of 2.363 mm and a complex relative permittivity er = 2.55 — y'2.295 • 10-3, while the top layer has a height of 1.524 mm and a complex relative permittivity er =2.17 - j 1.953 • 10~3. The optimization is carried out in a region defined by u e [0.1,0.45] and v e [-0.15,0.15] around the coverage area, with a resolution of the far fields of 512 x 512 points for the FFT, having a total of 4047 points in the UV grid. 0.15 0.1 0.05 -0.05 -0.1 -0.15, ^LT Coverage 28 dBi [ ' I ~ ! u -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 (a) 0.15 0.05 -0.05 -0.15, -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 (b) Fig. 6. Starting European coverage radiation pattern in gain (dBi) for Y polarization after POS. (u,u) are in the reflectarray coordinate system, (a) Copolar. (b) Crosspolar. All of the reflectarray elements are optimized at the same time, considering six dipoles as independent optimizing variables (s = 6) maintaining the unit cell symmetry, thus having a total number of optimizing variables of 31080. The optimization is carried out in fixed gain, setting the required gain to 31.90dB in the direction (uo, vo) = (0.28, 0). The crosspolar template is set 35 dB below the maximum copolar gain template, and the crosspolar component is scaled by a factor of 105 in natural units with the weighting function. The LMA is set to perform three iterations per iteration of the generalized IA. Although the design and crosspolar optimization have been carried out for both linear polarizations at the same time, only the results for Y polarization will be shown graphically, since the radiation patterns are similar for both polarizations before the optimization, and the Y polarization results after the optimization are slightly worse than for the X polarization. However, numerical values will be provided for both polarizations. Fig. 6 shows the radiation pattern obtained for the design after the POS for Y polarization. The gain requirements are