Direct design of aspherical lenses for extended non-Lambertian sources in three-dimensional rotational geometry

: Illumination design used to redistribute the spatial energy distribution of light source is a key technique in lighting applications. However, there is still no effective illumination design method for extended sources, especially for extended non-Lambertian sources. What we present here is to our knowledge the first direct method for extended non-Lambertian sources in three-dimensional (3D) rotational geometry. In this method, both meridional rays and skew rays of the extended source are taken into account to tailor the lens profile in the meridional plane. A set of edge rays and interior rays emitted from the extended source which will take a given direction after the refraction of the aspherical lens are found by the Snell’s law, and the output intensity at this direction is then calculated to be the integral of the luminance function of the outgoing rays at this direction. This direct method is effective for both extended non-Lambertian sources and extended Lambertian sources in 3D rotational symmetry, and can directly find a solution to the prescribed design problem without cumbersome iterative illuminance compensation. Two examples are presented to demonstrate the effectiveness of the proposed method in terms of performance and capacity for tackling complex designs.


Introduction
One main purpose of illumination design is to redistribute the spatial energy distribution of the light source to produce a prescribed distribution by a means of some elaborately designed optical surfaces [1]. Generally, the illumination design algorithms can be divided into two groups: zero-étendue algorithms [2][3][4][5][6][7][8][9] and algorithms for extended sources [10][11][12][13]. For a zero-étendue algorithm, the light source is assumed as an ideal source (a point source or a parallel beam), in which there is only one single ray passing through each point on the optical surface, as shown in Fig. 1(a). This simplistic assumption allows us to easily achieve a predefined deflection of the optical surface for each incident ray [6,9]. Since the étendue of an actual light source usually cannot be zero in practical applications, these zero-étendue algorithms are invalid in a compact design where the influence of the size or the angular extent of the extended source on the performance of an illumination system cannot be ignored. In this case, we need to resort to the algorithms for extended sources. For an extended source algorithm, the size or the angular extent of an extended source is considered in the design. Consequently, there will be more than one ray passing through each point on the optical surface, as shown in Fig. 1(b), which makes the illumination design for extended sources quite different, and more difficult than that for ideal sources. Illumination design for extended sources in three-dimensional (3D) geometry can be classified into three subcategories: translational design, rotational design and freeform design, as shown in Fig. 2. The translational case is a fundamental design, and the other two cases are more common in practical applications and more difficult to achieve. Most of the existing direct methods developed for extended sources are only effective in translational cases for extended Lambertian sources whose luminance is a constant [10][11][12][13]. That means those methods cannot be applied to either rotational cases or freeform cases, and are not valid for extended non-Lambertian sources whose luminance distribution is usually a function of position and direction. The simultaneous multiple surface (SMS) method is a powerful tool in illumination design for extended sources [14]. The difficulty of this direct method lies in establishing a specific relationship between the input and output wavefronts, and currently the SMS method is only effective for extended Lambertian sources. The only direct method of illumination design for extended non-Lambertian sources in 3D translational symmetry is developed by Wu in our previous work [15]. However, there is still no effective direct method that can be applied to the rotational cases or the freeform cases for extended non-Lambertian sources. Some illuminance compensation approaches [16][17][18][19] can be applied to rotational cases or freeform cases for extended non-Lambertian sources under certain conditions. The key to these compensation methods is iteratively using illuminance compensation to improve the performance of a design created by a zero-étendue method during each iteration. These methods may not be the best choice due to the cumbersome iterative ray-tracing process, and sometimes energy waste caused by total internal reflection cannot be avoided due to the optical surfaces generated by a zero-étendue method. Besides, these methods are less effective in those designs where the ratio of the distance between the lens and the source to the size of the source (i.e. H/d) is less than 3 [20], as shown in Fig. 2(c). Illumination design for extended sources in 3D geometry. (a) In a translational design, the lens is created by translating the lens profile along the y-axis. The performance on each meridian plane is the same. (b) In a rotational design, the lens is created by applying the rotation to the lens profile. (c) In a freeform design, the lens has no rotational symmetry. H represents the distance between the optical surface to the light source, and d represents the size of the extended light source.
To resolve the challenge currently faced by illumination design for extended sources, in this paper we will for the first time develop a direct method of illumination design for extended non-Lambertian sources in 3D rotational geometry. As the rest of this paper will show, both meridional rays and skew rays of an extended source can be well controlled by tailoring the lens profile only in a meridional plane of an aspherical lens, and the solution to a given design problem can be directly found by the proposed method without cumbersome iterative illuminance compensation. Besides, this direct method is effective for both extended non-Lambertian sources and extended Lambertian sources in 3D rotational geometry, and has the performance and capacity to easily tackle complex designs. Usually, a prescribed illuminance design can be converted into a prescribed intensity design as long as the influence of the lens size on the optical performance can be ignored. Thus, we only address the prescribed intensity design in this paper. The rest of this paper is organized as follows. In Section 2, at first we will introduce two principles that are used to develop the direct method for extended non-Lambertian sources, and then we will show how to tailor the lens profile in a meridional plane of an aspherical lens to control both meridional rays and skew rays of the extended source by using these two principles. In Section 3, design examples will be presented to verify this direct method, and elaborate analyses and comparisons of these designs will also be made in this section before we conclude our work in Section 4.

Principles and design methodology
The proposed direct method of illumination design for extended non-Lambertian sources is built upon the following two important principles. Assume an infinitesimal light beam, 1, emitted from the extended light source is refracted by the optical surface at point P, as shown in Fig. 3. We suppose the optical system is loss-less. Then, the luminance of beam 1 will be conserved in this loss-less system. That means the luminance of beam 1, L s , will be equal to that of its outgoing beam, 2, after the refraction of the optical surface. Thus, the luminance of each outgoing ray can be obtained as long as its incident ray is known. Fig. 3. The luminance of an incident beam is conserved in a loss-less system.
Given an optical surface, a set of incident rays (both edge rays and interior rays) emitted from the extended non-Lambertian source can be found, which take the same direction after the refraction of the optical surface, as shown in Fig. 4(a). Assume O = (sinβ i ,0,cosβ i ) is the unit vector of the outgoing rays. Here β i is the angle formed between the outgoing ray and the z-axis (the optical axis) in a meridional plane of the lens. Extending these outgoing rays yields a set of intersection points in the xy-plane, as shown in Fig. 4(b). Then, we can obtain a domain Ω i formed by these intersection points in the xy-plane. Since both positions and directions of the incident rays are known and the luminance is conserved in this loss-less system, a luminance distribution defined by these outgoing rays on the domain Ω i can be obtained. Suppose this luminance distribution can be represented mathematically by L i = L i (x,y). Then, the output intensity formed by these outgoing rays at the direction of the unit vector O can be calculated by To address the challenge of illumination design for extended non-Lambertian sources in 3D rotational symmetry, next we will show how to tailor the lens profile only in a meridional plane of an aspherical lens by using these two principles. Figure 5 gives the flow chart of the proposed direct method. There are two design phases in this method, as shown in Fig. 5. First, we need to find an initial curve of the lens profile and then calculate the rest of the lens profile during the second phase. The key steps of both these design phases are finding a set of edge rays and interior rays emitted from the extended source which exit the aspherical lens toward a given direction, and calculating the output intensity at this direction according to the luminance distribution of those outgoing rays at this direction. More details about the design process of this proposed direct method will be given below. We assume the extended non-Lambertian source is a disk light source with a diameter of d 1 and an angular range of emission between 0 and φ max . We further assume the optical axis (i.e. the z-axis) intersects with the source at the center of the disk and the optical surface is rotationally symmetric around the z-axis. Assume the extended source is surrounded by a medium of refractive index n, and let L(r,φ) denote the luminance function of the source. Here, r is the radius of the point where the incident ray emanates from the source, and φ is the incident angle formed between the incident ray and the z-axis. For an extended source, the relationship between the luminance of the source and the power, Φ 1 , emitted from the source can be expressed as where E denotes the étendue of the extended light source in 3D geometry. Integrating Eq. (2) yields the total power of the extended source, which is given by ( ) Let I t (β) represent the prescribed output intensity distribution with an angular range between 0 and β max . We suppose that the total flux of the outgoing beam is equal to that of the incident beam in a loss-less system. That means Φ 1 should also satisfy ( ) Although both the meridional rays and the skew rays of the extended source need to be well controlled in a 3D rotational case, we can still perform the given design only in a meridional plane of an aspherical lens. Figure 6 shows a meridional plane (the xz-plane) of the aspherical lens. Here, S 1 S 2 represents the extended source in the meridional plane. Due to the nature of an extended source, an initial curve is required before we start the design. In Fig. 6(a), C 1 C 2 represents the initial curve and C 2 is a mirror point of C 1 about the z-axis. One requirement for the initial curve is that the refracted rays of the two edge rays, S 1 C 1 and S 2 C 2 , by the initial curve should take the same resulting direction angle β = 0°. A common way to find an initial curve to meet this requirement is to employ polynomial fitting techniques, such as a 2nd-order polynomial [1]. However, the polynomial fitting techniques may not be able to generate a stable design [21]. In our previous work, an alternative technique was presented by Wu [21]. This technique predefines a linear relationship between the x-coordinates of those points on the initial curve and the x-coordinates of those points on the light source. Assume that (x N , h) are the coordinates of C 1 in the xz-plane. According to this technique, we define a set of equally spaced data points along the x-axis on C 1 C 2 and S 1 S 2 , respectively, and let the incident rays Q i P i (i = 0,1,…,N) exit the initial curve toward the direction β = 0°, as shown in Fig. 6(a). As a result, a set of outgoing rays are obtained which are parallel to the z-axis between the two edge rays, 1 and 2. Since the incident rays of these outgoing rays are predefined, the luminance of these outgoing rays is known. By applying the rotation to C 1 P N/2 (here, P N/2 is the vertex of the aspherical surface.), we can obtain an initial surface patch. According to the technique for designing the initial curve, it is not difficult to find all incident (skew and meridional) rays which will take the resulting direction β = 0° after refraction of the initial surface patch, because those incident skew rays can be easily found just by rotating those incident meridional rays Q i P i about the z-axis. After that, we can get a domain Ω 0 formed by all the outgoing rays at the direction β = 0° in the xy-plane, and also the luminance distribution, L 0 = L 0 (x,y), defined by these outgoing rays on the domain Ω 0 . Then, the output intensity, I(0), formed by these outgoing rays at the direction β = 0° can be calculated by Eq. (1), or alternatively the output intensity at this direction can be approximately calculated by where, x i is the x-coordinate of point P i on the lens profile, and f 0 (x i ) represents the luminance of the incident ray Q i P i . Then, we only need to choose an appropriate value of x N for point C 1 such that the output intensity, I(0), at the direction β = 0° equals the prescribed intensity, I t (0). Here, we use the bisection method to find the value of x N . When the initial curve is ready, the following procedure is taken to calculate the rest of the lens profile. Take the calculation of point P N+i as an example (Assume the previous point P N+i-1 is already obtained.). The first step is to trace an edge ray, S 2 P i and calculate the direction angle β i of its outgoing ray, the ray 1, as shown in Fig. 6(b). The second step is to trace a ray from S 1 to calculate the new point P N+i of the lens profile. As shown in Fig. 6(b), the ray from S 1 satisfies two conditions: (1) its intersection with the tangent line of its previous point, P N+i-1 , is at the point P N+i , and (2) its refracted ray exits the lens toward the direction β = β i . The normal vector at point P N+i can be calculated by application of the Snell's law. Then, by applying the rotation to P N+i P N/2 , we can obtain a surface patch. The next step is to find a set of incident rays (edge rays and interior rays of the extended source) which will take the same resulting angle β = β i after the refraction by this surface patch. Finding these edge rays is a non-trivial process in a 3D rotational geometry. As shown in Fig. 6(b), there are only two edge rays, S 2 P i and S 1 P N+i , in the meridional plane that exit the surface patch toward the direction β = β i . However, there exist the skew edge rays which exit the surface patch toward the direction β = β i . Due to the symmetry of the 3D rotational design, we only need to find the skew edge rays in the first and second quadrants. Figure 7(a) shows one situation where i<N/2. In this figure, point P N-i is a mirror point of P i about the z-axis on the lens profile. We calculate a set of discrete data points M j (j = 1,…,N1) between points P N-i and P N+i by using linear interpolation where M 1 and M N1 correspond to P N-i and P N+i , respectively. Assume (x, 0, z) are the coordinates of point M j and (N x , 0, N z ) is the unit normal vector of the surface patch at point M j . We rotate the point M j about the z-axis by the angle θ j , and get a new point R j on the surface patch, as shown in Fig.  7(a). The coordinates of point R j are obtained as (xcosθ j , xsinθ j , z), and the unit normal vector of the surface patch at point R j are denoted as N R = (N x cosθ j , N x sinθ j , N z ). A ray traced inversely at direction β = β i passes through the surface patch at point R j and intersects the xy-plane at point S Rj . Let I denote the unit vector of this incident ray, and we have I = (-sinβ i , 0, -cosβ i ). According to the Snell's law, we can get the unit vector, O = (O x , O y , O z ), of the outgoing ray, which is calculated by where the parameter T 1 is given by Then, we can obtain the coordinates of point S Rj , which are defined by cos , sin Changing θ j from 0 to π rad in Eq. (8) produces a trajectory of point S Rj in the xy-plane, as shown in Fig. 8(a). In this figure, the j-th red solid curve represents the trajectory of point S Rj . Similarly, if we rotate point P N-i (M 1 ) about the z-axis by the angle θ 1 from 0 rad to π rad, we can also get a trajectory in the xy-plane which is represented by the 1st red solid curve in Fig. 8(a). From this figure we can find a trajectory that will have no intersection point with the boundary of the extended source when the x-coordinate of point M j is included in the set [0, x N-i ) [Here, x N-i is the x-coordinate of point P N-i (M 1 ).]. That is to say, all the points located between points P N/2 and P N-i will not be able to generate an intersection point with the boundary of the extended source in the xy-plane. Also, from this figure we see there is only one intersection point between each trajectory and the boundary of the light source for each point M j . These intersection points are just what we want. To find these intersection points, we let the coordinates of point S Rj satisfy Eq. (9).  Equation (9) is a highly nonlinear equation for which numerically finding solutions might be very difficult. Figure 8(a) suggests that the solution to Eq. (9) is unique. Thus, we can use the bisection method to find the root. From Fig. 8(a) we can see P N-i (M 1 ) will produce an intersection point with the boundary of the extended source when P N-i is rotated about the z-axis by the angle of π rad, and P N+i (M N1 ) will produce an intersection point when P N+i is rotated about the z-axis by the angle of 0 rad. After we get the angle of rotation, θ j , for each point M j (j = 1,…,N1), we can obtain a set of the skew edge rays which will exit the surface patch toward direction β = β i .
Next, we need to calculate a set of the interior rays that will exit the surface patch also toward the direction β = β i , by using those edge rays obtained above. At each point M j , we define a set of data points between points P N/2 and M j by using linear interpolation. Then, we rotate these points about the z-axis by the angle θ j which is just the angle of rotation obtained above, and get a set of new points Q j,k (k = 1,…,N2) on the surface patch, as shown in Fig. 8(b). Obviously, these new points will lie on the curve P N/2 R j on the surface patch. By applying Eqs. (6)-(8) to these new points Q j,k , we obtain a set of intersection points A j,k (k = 1,…,N2) on the source in the xy-plane for the rotation angle θ j , as shown in Fig. 8(b). Repeating the numerical calculation, we obtain a net of boundary intersection points and interior intersection points on the source, as shown in Fig. 8(c). The green solid square represents the interior intersection points, and the blue solid circle represents the boundary intersection points. According to these intersection points, we can calculate the incident rays traced from the source, which will take the direction β = β i after the refraction of the surface patch. Extending these outgoing rays and calculating the intersection points of these outgoing rays in the xy-plane, we can obtain a domain Ω i formed by these outgoing rays in the xy-plane, as shown in Fig. 8(d). Since the incident rays are known and the luminance is conserved in this loss-less system, we can get a luminance distribution defined by these outgoing rays on the domain Ω i . Then, the output intensity produced by these outgoing rays at the direction β = β i can be calculated by Eq. (1). Sometimes, direct calculation of Eq. (1) might not be that easy. Here, we give an alternative way to calculate the output intensity. Let L j,k denote the luminance of the (j-th,k-th) outgoing ray. For a triangle formed by three neighboring intersection points on domain Ω i , the output intensity produced by this triangle at the direction β = β i is given by where (Area) j,1 represents the area of the triangle, as shown in Fig. 9(a). For a quadrangle formed by four neighboring intersection points on domain Ω i , the output intensity produced by this quadrangle will be Then, the output intensity formed by all the outgoing rays at the direction β = β i can be calculated by ( ) Increasing the values of N1 and N2 will undoubtedly reduce the deviation between the results of Eqs. (1) and (12), however, the computation speed should also be taken into account. Thus, we need to define appropriate values of N1 and N2 to make a trade-off between the calculation accuracy and the computation speed. What we need to do next is only to adjust the position of the point P N+i on the tangent line of the previous point P N+i-1 so that the output intensity I(β i ) equals the prescribed intensity, I t (β i ). Similarly, this calculation can be finished by the bisection method. The design process presented above is for the situation where i<N/2, and it will be a little different when i≥N/2. When x i ≥0 (x i is the x-coordinate of point P i .), we calculate a set of discrete data points M j (j = 1,…,N1) between points P i (M 1 ) and P N+i (M N1 ) by using linear interpolation, as shown in Fig. 7(b). Similarly, we rotate the point M j about the z-axis by the angle θ j and get a trajectory of point S Rj in the xy-plane, as shown in Fig. 8(e). From 8(e) we see the angle of rotation for point P i (M 1 ) equals 0 rad, and that for point P N+i (M N1 ) equals 0 rad too, which is quite different from that shown in Fig. 8(a). After we get the angle of rotation, θ j , for each point M j (j = 1,…,N1), we can obtain a set of skew edge rays which exit the surface patch toward the direction β = β i . At each point M j , the set [0, θ j ] is uniformly discretized into a set of discrete angles {θ j,k | θ j,k = (k-1) θ j /(N2-1), k = 1,…,N2}. Then, rotating the point M j about the z-axis by the angles θ j,k yields a set of new points Q j,k (k = 1,…,N2) on the surface patch, as shown in Fig. 8(f). By applying Eqs. (5)-(8) to these new points Q j,k , similarly we obtain a set of intersection points A j,k (k = 1,…,N2) on the source, as shown in Fig. 8(f). Further, we can obtain a net of boundary intersection points and interior intersection points on the source, as shown in Fig. 8(g). With these intersection points, we could find the incident rays which take the direction β = β i after the refraction of the surface patch and their corresponding outgoing rays. Extending these outgoing rays yields a domain Ω i formed by these outgoing rays in the xy-plane, as shown in Fig. 8(h). According to Eqs. (10)-(12), we can calculate the output intensity produced by these outgoing rays at the direction β = β i .
When the coordinates of point P N+i are obtained, the normal vector of the lens profile at point P N+i can be calculated by application of the Snell's law, and then we perform the whole calculation described above until the direction angle of the incident ray from S 1 satisfies φ >φ max . In Fig. 10, M denotes the right end point of the lens profile. Assume the direction angle of the incident ray S 1 M equals φ max , and that of its outgoing ray, the ray 6, equals β m . We calculate the incident ray S 2 P i1 with a direction angle φ m , which exits the lens also toward direction β = β m . Since only one surface is used here, we have β m <β max with no more degrees of freedom to control those rays emitted from S 2 between φ m <φ≤φ max [1]. It means the incident rays from S 2 between φ m <φ≤φ max cannot be well controlled, and consequently we have β m <β max (here, we define β m as the maximum effective angle we can obtain.). Then, an arbitrary ray emitted from S 2 with a direction angle between φ m <φ≤φ max will take a resulting direction angle between β m <β≤β t (β t is the direction angle of the ray 7. Usually, β t >β max ). When the lens profile is ready, the 3D model of the aspherical lens is generated by applying the rotation to the lens profile. Fig. 10. The rays from S 2 between φ m <φ≤φ max cannot be well controlled.

Design examples
In this section, two challenging designs will be given to demonstrate the effectiveness of the proposed direct method in illumination design for extended non-Lambertian sources in 3D rotational symmetry. In the first design, the luminance distribution of the extended non-Lambertian light source is given in Eq. (13). This equation shows clearly that the luminance distribution of the light source is a function of both position and direction. The prescribed output intensity distribution is a piecewise function, which is given in Eq. (14). Here, the value of K 1 can be calculated by applying the energy conservation, as shown in Eqs. (3) and (4). The other design parameters are given in Table 1. Figure 11(a) shows the luminance distribution of the light source at the direction φ = 0°, which is a Gaussian distribution with a beam waist of 1mm. According to the proposed direct method, we calculate the luminance distribution of the outgoing rays at the direction β = 0° in the meridional plane, as shown in Fig.  11(b). Figure 11(c) gives the actual intensity distribution generated by the proposed method, which is represented by the solid red line.  The fractional RMS is employed here to quantify the difference between the actual intensity and the prescribed one, which is defined by where, Num denotes the number of the sample points, I tj1 is the target intensity of the j1-th point and I aj1 is the actual intensity of the j1-th point. A smaller value of RMS represents less difference (i.e. a better agreement) between the actual intensity and the prescribed one. Due to the limitation of one single surface, inevitably there will be a region of abrupt intensity change near β C which denotes the maximum effective angle obtained from ray-tracing, as shown in Fig. 11(c). Since this region cannot be eliminated, a design can be considered acceptable as long as good agreement is achieved for the region within [0,β C ]. Let Num = 400. According to the actual intensity distribution, we have β C = 53.09° and RMS = 0.0027. Obviously, good agreement is achieved between the actual intensity and the target. The lens model and the lens profile are both given in Fig. 11(d). Let H denote the z-coordinate of the vertex of the lens, and we have H/d = 2.69 in this design. Such a compact system might be a huge challenge for most illuminance compensation methods which use zero-étendue methods in the design, because those methods are usually less effective when the ratio H/d is less than 3 [18,20]. Besides, we observe that all the incident rays have passed through the aspherical surface. That means total internal reflection does not take place on the aspherical surface in this compact design. According to the actual intensity distribution, we know that the outgoing rays within [0,β C ] account for 81.19% of the total energy, and those within [0,β max ] account for 98.12% of the total energy. Obviously, the energy efficiency within the required region is pretty high. When the luminance distribution of the extended non-Lambertian source is both direction-invariant and space-invariant, the extended non-Lambertian source will degenerate into an extended Lambertian source. For an extended Lambertian source, the design process presented in Section 2 can be simplified (We only need to find the boundary intersection points.) because the luminance of a Lambertian source is a constant. In this second design example, except that an extended Lambertian source is utilized, the target intensity distribution and the design parameters remain identical to the first design example. Following the simplified design process, we obtain a new lens design with both the lens profile and 3D model shown in Fig. 12(a). From this figure we can see the lens is quite compact and we have H/d = 2.62. In this compact design, we also observe all the incident rays have passed through the aspherical surface without total internal reflection taking place on the aspherical surface. Figure 12(b) plotted the actual intensity distribution in solid red line while the target intensity in dashed blue line. From this figure we have the fact that β C = 53.87° and RMS = 0.0016. Obviously, the actual design achieved very good agreement between the actual intensity and the prescribed one. According to the actual intensity distribution, the outgoing rays within [0,β C ] account for 81.15% of the total energy, and those within [0,β max ] account for 97.10% of the total energy. To make a comparison with the first design example, Fig. 12(c) plots the intensity distribution of the extended Lambertian source and the intensity distribution of the extended non-Lambertian source. From this figure we observe that the two intensity distributions are quite different. Figure 12(d) shows the sag difference between the lens profile designed for the non-Lambertian source and that for the Lambertian source on the region [-2.3674, 2.3674]. From this figure we see the two different intensity distributions yield obvious difference between the two lens profiles. In our previous work [21], we developed a feedback method for extended Lambertian sources in 3D rotational geometry. By this method, a prescribed 3D intensity design is first converted into a two-dimensional intensity design for the extended Lambertian source, and then a feedback strategy is employed to improve the performance of the aspherical lens in 3D rotational geometry. This method also has the performance and capacity to easily tackle complex designs. For the purpose of comparison, here, we will employ this feedback method to achieve the same prescribed intensity design given in Eq. (14) while maintaining the same light source and design parameters as those utilized in the second design example. Figure 13(a) plotted the actual intensity distribution obtained by the feedback method in solid red line while the target intensity in dashed blue line. This figure shows that good agreement is achieved between the actual intensity and the target with the fact that β C = 53.82° and RMS = 0.0026. In this design, we once again observe that all the incident rays have passed through the aspherical surface. Overall the outgoing rays within [0,β C ] account for 82.88% of the total energy, and those within [0,β max ] account for 96.97% of the total energy. Figure 13(b) shows the sag difference between the lens profile generated by the feedback method and that generated by the direct method for the Lambertian source on the region [-2.3674, 2.3674]. Less than 8 μm sag difference is achieved across the lens aperture. By making comparisons of the optical performance and the lens profiles of the two designs, we can see the two designs generated by the proposed direct method and the feedback method are almost the same. This indicates that the proposed direct method and the feedback method are both effective for extended Lambertian sources in 3D rotational symmetry.

Conclusion
In conclusion, we present the first direct method that works robustly for extended non-lambertian light sources in 3D rotational geometry. This direct method is built upon two important and basic principles, and uses a novel strategy to simultaneously control both meridional rays and skew rays of the extended source in the meridional plane. The proposed direct method is effective for both extended non-Lambertian sources and extended Lambertian sources. For a non-Lambertian source, a set of edge rays and interior rays which take the given direction after refraction of the aspherical surface need to be found. For a Lambertian source, however, only a set of edge rays are required. The output intensity at a given direction is equal to the integral of the luminance function of all the outgoing rays at this direction. A quite good solution to the given design problem can be directly found by this direct method without cumbersome iterative illuminance compensation. This direct method allows all the incident rays to pass through the aspherical surface without being totally reflected internally on the surface, which results in high energy efficiency, as demonstrated by the design examples. This work is important from both mathematical and practical standpoints. The current study is focused on 3D rotational designs, and it would be interesting to generalize this direct method to freeform designs in future work.