LASER BEAM SHAPING AND TRANSFORMATION USING FREEFORM SURFACE
DESIGN
by
Xiangchao Zhu
A thesis submitted to the Graduate Faculty of
Auburn University
in partial fulfillment of the
requirements for the Degree of
Master of Science
Auburn, Alabama
May 3, 2014
Copyright 2014 by Xiangchao Zhu
Approved by
Hulya Kirkici, Chair, Professor of Electrical and Computer Engineering
Stanley Reeves, Professor of Electrical and Computer Engineering
Michael Baginski, Associate Professor of Electrical and Computer Engineering
ii
Abstract
Lasers are used in many areas including material processing (lithography, semiconductor
manufacturing), welding, cutting, drilling, medical procedures (eye surgery and cosmetic skin
treatments), and optical data processing. Laser diodes are the most developed lasers in recent
years, due to their multiple features like low cost, compactness, electronic compatibility, broad
range of wavelengths, and high pulse prepetition frequency values. On the other hand, the
irradiance distribution of a laser-diode is elliptical Gaussian, not circular. In order to be widely
used commercially, the beams from the laser diodes must be shaped to a desired shape through
optical elements by the redistribution of energy profile into other beam irradiance profiles [3].
This work focuses on redistribution of Gaussian profile laser beam irradiance to other
forms propagating through an optical media. In this work, single lens system is designed by
getting surface parameters which are the coordinates of the surface points through solving the
first order partial differential equations based on Snell?s law and energy conservation principle.
Meantime, a lens array is applied to the multi-laser source with the Gaussian distribution. The
simulation results show that this method can be very well used for beam shaping. For the two
lens system, the surface parameters are solved by the same method used in one freeform lens
design. Meanwhile, another new method based on the separated variable mapping is also used in
acquiring the surface parameters. This separated variable mapping means that the surface
parameters of the source plane and target plane can be numerically specified based on ?source-
to-target? [4] or ?target-to-source? [5,6]. And the mapping is calculated based on energy
iii
conservation principle, Snell?s law, and constancy of optical path length. According to
comparison of the simulation results from two different methods, the method based on solving
partial differential equations is more effective and efficient to transform a collimated Gaussian
laser beam into a rectangular ?flat-top? one, and the uniformity is almost 90%. Also, a lens array
is applied to the multi-laser source with Gaussian distribution. The simulation results again show
that the method based on solving partial differential equations is better than the other one.
iv
Acknowledgments
The author would like to thank his advisor Dr. Hulya Kirkici, who has been giving so
much patience and guidance during the past two years that allows the author to finally finish this
project. Without her support this project could not come to the current level of performance.
The author would like to thank the other two committee members, Dr. Baginski who taught him
a lot of knowledge of electromagnetic fields and Dr. Reeves who guided him to gain the
knowledge of using MATLAB to do the simulation during the classes. The author also would
like to express appreciation to his research group members, thanks for all the helps that they have
given to support him. Finally, the author wishes to thank his parents, Zhian Zhu and Hongying
Qian, for their consistent encouragement and care.
v
Table of Contents
Abstract ............................................................................................................................................ ii
Acknowledgments........................................................................................................................... iv
Table of Contents .............................................................................................................................v
List of Tables ................................................................................................................................. vii
List of Illustrations ........................................................................................................................ viii
List of Abbreviations .......................................................................................................................x
CHAPTER 1 ................................................................................................................................... 1
BACKGROUND AND RESEARCH OBJECTIVES .................................................................... 1
1.1 Laser diode .................................................................................................................... 1
1.2 Non-imaging optics....................................................................................................... 3
1.3 Research Objectives ...................................................................................................... 4
CHAPTER 2 ................................................................................................................................... 7
DEVELOPED TECHNIQUES AND TOOLS ............................................................................... 7
2.1 Runge-kutta method ...................................................................................................... 7
2.2 First-order quasi-linear partial differential equations (PDEs) ...................................... 9
2.3 Optical path length (OPL)........................................................................................... 10
2.4 Fermat's principle ........................................................................................................ 11
2.5 Vector Snell?s Law...................................................................................................... 12
2.6 Non-uniform rational basis spline............................................................................... 14
vi
2.6.1 NURBS CURVE.............................................................................................. 14
2.6.2 NURBS SURFACE ......................................................................................... 15
2.7 Top-hat beam .............................................................................................................. 17
2.8 Monte-Carlo method ................................................................................................... 21
CHAPTER 3 ................................................................................................................................. 23
METHODOLOGY........................................................................................................................ 23
3.1 One freeform lens........................................................................................................ 23
3.2 Two freeform lens ....................................................................................................... 26
3.2.1 Based on Partial Differential Equations (PDEs) .............................................. 26
3.2.2 Based on separated variables mapping ............................................................ 28
CHAPTER 4 ................................................................................................................................. 32
DESIGN EXAMPLES AND RESULTS ...................................................................................... 32
4.1 One freeform lens........................................................................................................ 32
4.1.1 Incident beam is circular Gaussian beam......................................................... 33
4.1.2 Incident beam is elliptical Gaussian beam ....................................................... 35
4.2 Two freeform lens ................................................................................................................... 38
4.2.1 Incident beam is circular Gaussian beam......................................................... 38
4.2.2 Incident beam is elliptical Gaussian beam ....................................................... 42
CHAPTER 5 ................................................................................................................................. 51
CONCLUSIONS........................................................................................................................... 51
REFERENCES ............................................................................................................................. 52
Appendix ....................................................................................................................................... 58
vii
List of Tables
Table 4.1 ..................................................................................................................................... 30
Table 4.2 ..................................................................................................................................... 31
Table 4.3 ..................................................................................................................................... 33
Table 4.4 ..................................................................................................................................... 36
Table 4.5 ..................................................................................................................................... 40
viii
List of Illustrations
Figure 1.1 ...................................................................................................................................... 2
Figure 1.2 ...................................................................................................................................... 4
Figure 2.1 .................................................................................................................................... 11
Figure 2.2 .................................................................................................................................... 14
Figure 2.3 .................................................................................................................................... 15
Figure 2.4 .................................................................................................................................... 17
Figure 2.5 .................................................................................................................................... 18
Figure 3.1 .................................................................................................................................... 22
Figure 3.2 .................................................................................................................................... 24
Figure 3.3 .................................................................................................................................... 27
Figure 4.1 .................................................................................................................................... 32
Figure 4.2 .................................................................................................................................... 32
Figure 4.3 .................................................................................................................................... 33
Figure 4.4 .................................................................................................................................... 34
Figure 4.5 .................................................................................................................................... 35
Figure 4.6 .................................................................................................................................... 35
Figure 4.7 .................................................................................................................................... 38
Figure 4.8 .................................................................................................................................... 40
Figure 4.9 .................................................................................................................................... 41
ix
Figure 4.10 .................................................................................................................................. 43
Figure 4.11 .................................................................................................................................. 44
Figure 4.12 .................................................................................................................................. 45
Figure 4.13 .................................................................................................................................. 46
Figure 4.14 .................................................................................................................................. 46
Figure 4.15 .................................................................................................................................. 47
Figure 4.16 .................................................................................................................................. 47
Figure 4.17 .................................................................................................................................. 48
x
List of Abbreviations
VCSEL Vertical Cavity Surface Emitting Laser diode
ODE Ordinary Differential Equation
PDE Partial Differential Equation
OPL Optical Path Length
NURBS Non-uniform rational basis spline
1
CHAPTER 1
BACKGROUND AND RESEARCH OBJECTIVES
1.1 Laser diode
A laser diode is electrically a P-I-N diode. The active region of the laser diode is in the
intrinsic (I) region, and the carriers, electrons and holes, are pumped into it from the N and P
regions respectively. While initial diode laser research was conducted on simple P-N diodes, all
modern lasers use the double-heterostructure implementation, where the carriers and the photons
are confined in order to maximize their chances for recombination and light generation. Unlike a
regular diode used in electronics, the goal for a laser diode is that all carriers recombine in the (I)
region, and produce light. Thus, laser diodes are fabricated using direct bandgap semiconductors.
The laser diode epitaxial structure is grown using one of the crystal growth techniques, usually
starting from an N doped substrate, and growing the (I) doped active layer, followed by the P
doped cladding, and a contact layer. The active layer most often consists of quantum wells,
which provide lower threshold current and higher efficiency [1]. A semiconductor laser diode is
shown schematically in Figure 1.1, where the gain medium in the middle is surrounded by the
guiding layers. The refractive index of the guiding layers is somewhat greater than that of the
surrounding regions (substrate and cladding) so that the light can be confined to a relatively
narrow region by total internal reflection. The two ends of the diode in the cross direction (x-y
plane) are cleaved to form perfectly smooth, parallel edges, forming a Fabry-Perot resonator, and
thus an optical cavity is formed. When the electrical current is injected through the positive
2
electrode and collected at the base-plate on the bottom side of the junction (ground electrode),
the original stable state of thermal equilibrium of the junction will be disturbed, and electrons
can absorb energy either from heat caused by the collision among the massive ly increased
amount of electrons or from the photons caused by the spontaneous emission which is necessary
to initiate the laser oscillation. As a result, there are more electrons in higher energy states than
in lower energy states, thus population inversion is achieved in the gain medium within the
optical cavity.
When an electron is excited from a lower to a higher energy level, it will not stay that
way forever; instead it will go back to the ground state and emit photons. These photons will
travel along the waveguide and be reflected several times from the end faces before they are
emitted out. As the light wave passes through the cavity, it is amplified by the stimulated
emission which generates another photon of the same frequency, travelling in the same direction,
with the same polarization and phase as the first photon. This is how laser diode works and
generates laser light.
Figure 1.1: A semiconductor laser diode with the p-n junction structure.
3
1.2 Non-imaging optics
Non-imaging optics (also called anidolic optics) [2] is the branch of optics concerned
with the optimal transfer of light radiation between a source and a target. Unlike traditional
imaging optics, the techniques involved do not attempt to form an image of the source; instead
they intend to form an optimized optical system for optical radiative transfer from a source to a
target is desired.
The two design examples that non-imaging optics solves better than imaging optics are
solar energy concentration and illumination. Especially, the illumination part is the core concept
in this paper. The illumination part is controlling the distribution of light, typically so it is
"evenly" spread over some areas and completely blocked from other areas. Typical variables to
be optimized at the target include the total radiant flux, the angular distribution of optical
radiation, and the spatial distribution of optical radiation. These variables on the target side of the
optical system often must be optimized while simultaneously considering the collection
efficiency of the optical system at the source.
Examples of non-imaging optical devices include optical light guides, non-imaging
reflectors, non-imaging lenses or a combination of these devices. Common applications of non-
imaging optics include many areas of illumination engineering (lighting). Examples of modern
implementations of non-imaging optical designs include automotive headlamps, LCD backlights,
illuminated instrument panels, fiber optic illumination devices, and LED Lights.
4
1.3 Research Objectives
The laser diode example shown above is the basic edge-emitting semiconductor laser,
only one of three kinds of laser diodes. The other two includes double hetero-structure laser
diode and vertical cavity surface emitting laser diode (VCSEL), as shown in Figure 1.2. Though
different from each other in the structure, these laser diodes work on the same principle,
functioning as an optical oscillator by stimulating a chain reaction of photon emission inside a
tiny chamber. Thus they produce the same beam characteristics and can share the same beam
shaping methods. Throughout this thesis, edge-emitting laser diodes will be referred to as an
example by the author to demonstrate how to design and analyze the beam shaping systems
involved.
Figure
1.2: Three different structures of laser diodes [7]
As mentioned earlier, collimation, uniform intensity and pre-determined beam shape are
the three main features in the regions of laser beam shaping applications. There are many papers
presenting either to collimate the beam, to make it circular or rectangular, or to bring the
irradiance profile to a uniform irradiance. And most of them use the traditional Galilean or
Keplerian beam shaping system. These two kinds of refractive beam shaping systems are
constructed by dual aspheric lens [54-55] based on energy conservation principle, constancy of
optical path length, and Snell?s law. Specifically, the two different off-axis mirror system and
5
two different lens system are developed to circularize, collimate, and expand the elliptical beam
of edge-emitting semiconductor laser based LIDAR systems [3].
Traditional beam shaping optics design techniques usually utilize rotational or
translational symmetries, where the calculated cross section curve is swept around or along its
symmetry axis to generate the 3-D geometry shape. However, the need for beam shaping optics
that can distribute light beam in a non-rotational or non-translational manner has strongly
increased. This leads to the concept of freeform optical surfaces, which can provide much more
controlling freedom. The rapidly advancing manufacture technologies of freeform surfaces
stimulate the development of freeform beam shaping elements [8]. The freeform surfaces design
method has advantages with high degrees of design freedom that can be used to generate simpler
designs with better performance and produce some special illuminations [9]. There are many
kinds of freeform surfaces design methods in nonimaging optics: multiparameter optimization
[10-15], the ?Tailored Freeform Surface Method? [16-24], and the simultaneous multiple surface
(SMS) method [25-27].
Inspired by these works, this thesis will design optical systems to transform a collimated
elliptical or circular Gaussian laser beam into a uniform rectangular one with/without changing
the wavefront. And also this thesis will realize the large area uniform illumination by using 9
lens units. The equations specific to this system are derived analytically based on Snell?s Law
and Energy Conservation. Then Matlab is used to solve these equations. Matlab codes are
written to calculate and determine the parameters of the optical components. After getting the
discrete data of the freeform surface by calculation, Rhinoceros are applied for the 3D molding.
And then the Tracepro are applied to simulate the illumination system. The author expects that
this beam-shaping system can be used and also easily modified, and the analysis methods
6
developed in this thesis would be a valuable resource for the scientific community interested in
beam shaping systems.
7
CHAPTER 2
DEVELOPED TECHNIQUES AND TOOLS
2.1 Runge-kutta method
In mathematics, an ordinary differential equation (ODE) is an equation containing a
function of one independent variable and its derivatives. The term "ordinary" is used in contrast
with the term partial differential equation which may be with respect to more than one
independent variable.
Generally, there are two ways to solve the ODE, one is Euler method, and the other is
Runge-kutta method. However, the latter is more widely used in solving ODE. There are several
reasons that Euler?s method is not recommended for practical use, among them, (i) the method is
not very accurate when compared to other methods run at the equivalent step size, and (ii)
neither is it very stable.
In numerical analysis, the Runge-Kutta methods are an important family of implicit and
explicit iterative methods, which are used in temporal discretization for the approximation of
solutions of ordinary differential equations. These techniques were developed around 1900 by
the German mathematicians C. Runge and M. W. Kutta.
The formula for the fourth order Runge-Kutta method (RK4) is given below. Consider
the problem [28]:
0
( , )()y f t yyt ????? ?
?
(2.1)
8
Define h to be the time step size and ti = t0 + ih. Then the following formula:
0
1
1
2
2
3
43
1 1 2 3 4
( , )
( , )
22
( , )
22
( , )
1
( 2 2 )
6
ii
ii
ii
ii
ii
w
k hf t w
kh
k hf t w
kh
k hf t w
k hf t h w k
w w k k k k
?
?
??
?
?
?
?
? ? ??
?
?
? ? ??
?
? ? ??
?
? ? ? ? ? ?
?
(2.2)
Computes an approximate solution, that is wi ? y(ti).
9
2.2 First-order quasi-linear partial differential equations (PDEs)
If a PDE can be expressed as follows, then it can be called the quasi-linear PDE:
( , ) ( , )( , , ( , ) ) ( , , ( , ) ) ( , , ( , ) )u x y u x yx y u x y b x u u x y c x y u x yxy? ?????? (2.3)
The general solution to this kind of PDE is Finite Difference Method (FDM) and Finite
Element Method (FEM). However, the Runge-kutta method can be used to solve the first-order
quasi-linear PDE when the boundary conditions are certain. Especially, the fourth order Runge-
kutta method is more widely applied to this kind of quasi-linear PDEs.
10
2.3 Optical path length (OPL)
In optics, optical path length (OPL) or optical distance is the product of the geometric
length of the path that light follows through the system, multiplied by the index of refraction of
the medium through which it propagates. A difference in optical path length between two paths
is often called the optical path difference (OPD). Optical path length is important because it
determines the phase of the light and governs interference and diffraction of light as it propagates.
Optical path difference results in phase shift between two previously coherent sources
when passed through different mediums. For example a wave passed through glass will appea r to
travel a greater distance than an identical wave in air. This is because the light beam in the glass
will have travelled a greater number of wavelengths due to the higher refractive index of the
glass.
The OPD can be calculated from the following equation:
1 1 2 2OPD d n d n?? (2.4)
where d1 and d2 are the distances of the ray passing through medium 1 or 2, n1 is the greater
refractive index (e.g., glass) and n2 is the smaller refractive index (e.g., air).
In a medium of constant refractive index, n, the OPL for a path of physical length d is just
OPL nd? (2.5)
If the refractive index varies along the path, the OPL is given by
? ?LOPL n l dl?? (2.6)
where n(l) is the local refractive index as a function of distance, l, along the path L. Fermat's
principle states that the path that light takes between two points is the path that has the minimum
optical path length [29].
11
2.4 Fermat's principle
In optics, Fermat's principle or the principle of least time is the principle that the path
taken between two points by a ray of light is the path that can be traversed in the least time. This
principle is sometimes taken as the definition of a ray of light [30]. However, this version of the
principle is not general; a more modern statement of the principle is that rays of light traverse the
path of stationary optical length with respect to variations of the path [31]. In other words, a ray
of light prefers the path such that there are other paths, arbitrarily nearby on either side, along
which the ray would take almost exactly the same time to traverse [32].
Fermat's principle can be used to describe the properties of light rays reflected off mirrors,
refracted through different media, or undergoing total internal reflection. It follows
mathematically from Huygens ' Principle (at the limit of small wavelength). Fermat's text
Analyse des r?fractions exploits the technique of adequality to derive Snell's law of refraction
[33] and the law of reflection.
12
2.5 Vector Snell?s Law
In Figure 2.1 we have an interface between two materials with different indexes of
refraction n1 and n2. n1 is the index of refraction of the material you come from, and n2 of the
material you go to (in case of refraction).
The direction vector of the incident ray (= incoming ray) is i , and we assume this vector
is normalized. The direction vectors of the reflected and transmitted rays are r and t . These
vectors will be normalized as well. We also have the normal vector N , orthogonal to the
interface and points towards the first material n1.
The direction vectors of the rays can be split in components orthogonal and parallel to the
interface. We call these the normal component v? and the tangent component ||v of a vector v .
Figure 2.1: Situation for refraction and reflection [34]
The calculation of the refracted ray starts with Snell?s law [35] which tells that the
products of the index of refraction and sines of the angles must be equal:
12sin sinitnn??? (2.7)
13
Then we should try to find a formula for t . The first thing is to split it up into tangent
and normal parts:
||t t t??? (2.8)
According to the formula (2.7), the norms of the tangent parts happen to be equal to sines.
Hence:
1|| ||
2ti
??? (2.9)
Since ||t and ||i are in the same direction, this becomes:
11|| ||
22 c o s it i i N
?? ???? ? ??? (2.10)
The other one can be calculated by using the Pythagoras theorem:
2||- 1-t t N? ? (2.11)
Then substitute (2.10) and (2.11) in (2.8) to get the refracted direction vector:
211
||22 c o s 1innt i t N???? ? ? ?????
(2.12)
It is easy to know that the vector ||t equals to sin?t. Finally we get:
211
22
c o s 1 s initnnt i N????? ? ? ????? (2.13)
By Snell?s law, (2.13) can be expressed:
222 2 2
11
22
s in = s in ( 1 c o s )t i inn? ? ?? ? ? ???? ? ? ?? ? ? ? (2.14)
The last two equations will be used to calculate the refracted direction vector.
14
2.6 Non-uniform rational basis spline
Non-uniform rational basis spline (NURBS) is a mathematical model commonly used in
computer graphics for generating and representing curves and surfaces. It offers great flexibility
and precision for handling both analytic (surfaces defined by common mathematical formulae)
and modeled shapes.
NURBS are commonly used in computer-aided design (CAD), computer-aided
manufacturing (CAM), and computer-aided engineering (CAE) and are part of numerous
industry wide standards, such as IGES, STEP, ACIS, and PHIGS. NURBS tools which are found
in various 3D modelling and animation software packages.
They can be efficiently handled by computer programs and yet allow for easy human
interaction. NURBS surfaces are functions of two parameters mapping to a surface in three-
dimensional space. The shape of the surface is determined by control points. NURBS surfaces
can represent simple geometrical shapes in a compact form. T-splines and subdivision surfaces
are more suitable for complex organic shapes because they reduce the number of control points
twofold in comparison with the NURBS surfaces.
2.6.1 NURBS CURVE
In general, editing NURBS curves and surfaces is highly intuitive and predictable.
Control points are always either connected directly to the curve/surface or act as if they were
connected by a rubber band. Depending on the type of user interface, editing can be realized via
an element?s control points, which are most obvious and common for B?zier curves, or via
higher level tools such as spline modeling or hierarchical editing [36]. A B?zier curve is a
parametric curve frequently used in computer graphics and related fields.
The general form of a NURBS curve is [36-38]:
15
, *0 ,
0,
0
()( ) ( )
()
n n
i k i ii k i i
n i
i k ii
w N u Cr u N u C
w N u? ????
? ?? (2.15)
Here,
,* ,
,0
()()
()i k iki n i k ii
w N uNu
w N u?? ?
(2.16)
U={u0,?, um} is a nondecreasing sequence of real numbers, i.e., ui?ui+1, i=0,?, m-1. The ui are
called knots, and U is the knot vector. The ith B-spline basis function of p-degree (order p+1),
denoted by Ni,p(u), is defined as
1,0 1() 0 iii if u u uNu o th e r w is e ????? ??
1, , 1 1 , 1
11
( ) ( ) ( )ipii p i p i p
i p i i p i
uuuuN u N u N uu u u u??? ? ?
? ? ? ?
?????? (2.17)
In this, n is the number of control points Ci and wi are the corresponding weights. The
denominator is a normalizing factor that evaluates to one if all weights are one. The Figure 2.2
shows the NURBS curve, the control points are intersections of the lines.
Figure 2.2: NURBS curve [36]
2.6.2 NURBS SURFACE
A B-spline surface is obtained by taking a bidirectional net of control points, two knot
vectors, and the products of the univariate B-spline functions.
16
, , , ,00
, , ,00
( ) ( )( ) , 0 , 1
( ) ( )
nm
i j p i q j i jij
nm
i j p i q jij
w N u N u Cr u u v
w N u N u
??
??
? ? ????? (2.18)
With
11
11
11
11
{ 0 , .. ., 0 , , .. ., , 1 , .. ., 1 }
{ 0 , .. ., 0 , , .. ., , 1 , .. ., 1 }
p r p
pp
q s q
qq
U u u
V u u
? ? ?
??
? ? ?
??
?
?
(2.19)
U has r+1 knots, and V has s+1.The Figure 2.3 shows the NURBS surface.
Figure 2.3: NURBS surface
17
2.7 Top-hat beam
In optics, a top-hat (or top-hat) beam such as a laser beam or electron beam has a near-
uniform fluence (energy density) within a circular disk. It is typically formed by diffractive
optical elements from a Gaussian beam. Top-hat beams are often used in industry, for example
for laser drilling of holes in printed circuit boards. They are also used in very high power laser
systems, which use chains of optical amplifiers to produce an intense beam. Top-hat beams are
named for their resemblance to the shape of a top hat [39].
The irradiance distribution of the top-hat beam is:
m a x, 2
m a x
( / )()o u t T H c ir c R RIR R?? (2.20)
where Rmax labels the radius of the transverse size of the top-hat beam, and circ(R/Rmax) is equal
to unity for 0?R?Rmax and to zero otherwise [40].
Normally, the super-Gaussian, flattened-Gaussian, Fermi-Dirac, and super-Lorentzian are
included in the top-hat distributions.
A general flattened output irradiance profile can be written as a product of its
normalization constant, I0(?,R0), which is determined by requiring the input beam profile to be
normalized according to energy conservation, and its functional dependence f(?, R/R0); which
depends on the shape parameter ? and on the ratio of the radial coordinate, R; to the beam width,
R0, of each distribution
0 0 0 0( , ) ( , ) ( , )I R R I R f R R? ? ??? (2.21)
The width parameter, R0, defines a length scale over which the profile decreases to some
significant value for the specific profile distribution, such as, e-2 of its axial value. The shape
parameter ? species the shape of the profile distribution, such as, the power of the radial
coordinate of a super-Gaussian profile.
18
Explicitly, the Fermi-Dirac function is defined by:
1
, 0 ,( , ) ( , ) 1 e x p 1o u t F D F D F D F D FD
RI R R I R R? ? ? ?????????? ? ???????
??????
(2.22)
where RFD corresponds to R0 in Eq. (2.21) and is the radius at which the output irradiance falls to
half of its axial value. ? is a dimensionless shape parameter for the Fermi-Dirac function. When
? increases, the output profile approaches a more uniform distribution where there is a
continuous variation of irradiance from zero to its maximum value on axis.
The super-Lorentzian function is defined by
0,, ( , )( , ) 1 ( )S L S Lo u t S L S L M
SL
I M RI M R R RR? ???
??
(2.23)
where RSL is the radius at which the output irradiance falls to half of its axial value, and M is a
dimensionless shape parameter for the super-Lorentzian function. When M increases, the output
profile approaches a more uniform distribution where there is a continuous variation of
irradiance from zero to its maximum value of I0,SL (M,RSL).
Figure 2.4: (a) Fermi-Dirac and (b) super-Lorentzian output irradiance for shape parameters ?
and M. All profiles are normalized such that the total energy contained within each irradiance
19
distribution over all space is equal to unity. The cross-sectional distance is expressed in the
normalized units of the beam width.
The super-Gaussian function is defined by
, 0 ,( , ) ( , ) e x p 2
p
o u t S G S G S G S G SG
RI p R R I p R R???????? ??
????
(2.24)
where RSG is the radius at which the output irradiance falls to e-2 of its axial value, and p is a
dimensionless shape parameter for the super-Gaussian function. When p increases, the output
profile approaches a more uniform distribution where there is a continuous variation of
irradiance from zero to its maximum value of I0,SL (M,RSG).
The flattened-Gaussian function is defined by
? ? ? ? 22
, 0 , ,0
1( , ) ( , ) e x p 2 ( 1 ) ( )
!!
nm
N FG
o u t F G F G F G F G F G mn
N R RI N R R I N R N R R
nm
?
?
???????? ? ?
?? ?(2.25)
where RFG is the radius at which the output irradiance falls to ?(N + 1,N + 1)2 / ?(N + 1)2 of its
axial value. N is a dimensionless shape parameter for the flattened-Gaussian function. When p
increases, the output profile approaches a more uniform distribution where there is a continuous
variation of irradiance from zero to its maximum value of I0,FG (N,RFG).
20
Figure 2.5: (a) super-Gaussian and (b) flattened-Gaussian output irradiance for shape parameters
p and N. All profiles are normalized such that the total energy contained within each irradiance
distribution over all space is equal to unity.
21
2.8 Monte-Carlo method
Monte Carlo methods (or Monte Carlo experiments) are a broad class of computational
algorithms that rely on repeated random sampling to obtain numerical results; typically one runs
simulations many times over in order to obtain the distribution of an unknown probabilistic
entity. The name comes from the resemblance of the technique to the act of playing and
recording your results in a real gambling casino. They are often used in physical and
mathematical problems and are most useful when it is difficult or impossible to obtain a closed-
form expression, or infeasible to apply a deterministic algorithm. Monte Carlo methods are
mainly used in three distinct problems classes: optimization, numerical integration and
generation of draws from a probability distribution [41].
In physics-related problems, Monte Carlo methods are quite useful for simulating
systems with many coupled degrees of freedom, such as fluids, disordered materials, strongly
coupled solids, and cellular structures (see cellular Potts model). Other examples include
modeling phenomena with significant uncertainty in inputs, such as the calculation of risk in
business; and, in math, evaluation of multidimensional definite integrals with complicated
boundary conditions. In application to space and oil exploration problems, Monte Carlo-based
predictions of failures, cost overruns and schedule overruns are routinely better than human
intuition or alternative "soft" methods [42].
The modern version of the Monte Carlo method was invented in the late 1940s by
Stanislaw Ulam, while he was working on nuclear weapons projects at the Los Alamos National
Laboratory. It was named by Nicholas Metropolis, after the Monte Carlo Casino, where Ulam's
uncle often gambled [43]. Immediately after Ulam's breakthrough, John von Neumann
22
understood its importance and programmed the ENIAC computer to carry out Monte Carlo
calculations.
23
CHAPTER 3
METHODOLOGY
3.1 One freeform lens
Based on the Snell?s law and the energy conservation, we can deduce from the
characteristics of the source and desired illumination, to get a set of first-order partial differential
equations [44-48]. These two scientific laws are the basic principle of the freeform lens design.
The Snell?s law tells how the rays are refracted and reflected, which means the redistribution of
the luminous energy. With regard to the energy conservation, the relationship between input
energy and output one can be defined, and, furthermore, the input-to-output irradiance mapping
can be obtained based on this law using the variable separation method [49]. Then the collimated
light freeform lens design method will be built in the Cartesian coordinate system.
Assume I is the vector of the collimated incident light at point p, while O is the vector of
the refracted light by the freeform lens from point p to point t, and the normal vector at point p of
the freeform lens is N. The Cartesian coordinate system is built as shown in Figure 3.1, z is the
optical axis.
24
Figure 3.1: Schematic representation of the one freeform lens
The coordinates on p and t are (x,y,z(x,y)) and (tx,ty,tz), respectively. N is the normal unit
vector to the freeform surface. I and O are the unit vectors of the incident and transmitted beams,
respectively. Using the concept of differential calculus of multivariate, they can be expressed as:
22
1 ( , , 1 )1 xy
xy
N z zzz? ? ??? (3.1)
where, zx and zy are the first-order partial derivatives of z with respect to x and y, respectively.
2 2 2
1 ( , , )( ) ( ) ( ) x y z
x y z
O t x t y t zt x t y t z? ? ? ?? ? ? ? ? (3.2)
After using Snell?s law at point p, we have the equations as follows:
22 2 ( )o i o i o in O n I N n n n n O I? ? ? ? ? (3.3)
Here, no and ni indicate the refractive index of the incidence and emergence medium. Further, we
can get:
22 2 ( )o x i x x o i o in O n I N n n n n O I? ? ? ? ? (3.4)
25
22 2 ( )o y i y y o i o in O n I N n n n n O I? ? ? ? ? (3.5)
22 2 ( )o z i z z o i o in O n I N n n n n O I? ? ? ? ? (3.6)
Here, Ix, Iy, Iz are the three components of I. Ox, Oy, Oz are the three components of O. Nx, Ny, Nz
are the three components of N. Moreover, I=(0,0,1), and tz is a constant, which is defined as the
distance between source plane A and target plane C. Then we can get:
2 2 20
()( ) ( ) ( ) ( )oxx
o z i x y z
n x tz n t z n t x t y t z?? ? ? ? ? ? ? ? (3.7)
2 2 20
()
( ) ( ) ( ) ( )oyy o z i x y z
n y tz
n t z n t x t y t z
??
? ? ? ? ? ? ?
(3.8)
Equations (3.7) and (3.8) are the first-order partial differential equations about zx and zy, which
can be solve by using the Runge-kutta method.
According to light transmission energy conservation condition, the output of source is
equal to the flux incident in the target plane.
( , ) ( , )i n o u t x y x yI x y d x d y I t t d t d t? (3.9)
Therefore, the relationship between x, y and tx, ty can be acquired from equation (3.9).
26
3.2 Two freeform lens
3.2.1 Based on Partial Differential Equations (PDEs)
According to the previous section, we know how to design freeform lens based on the
Snell?s law and energy conservation principle. And the two freeform lens design is similar to the
one freeform lens design. One more principle is needed, which is the constancy of optical path
length.
The two freeform refractive surface optical system is shown schematically in Figure 3.2.
Figure 3.2: Schematic representation of the two freeform lens
Assume I is the unit vector of the collimated incident light at point p, R is the unit vector
of the ray through p and q, O is the unit vector of the refracted light by the second freeform lens
from point q to point t, and O is parallel to the optical axis. n1, n0 and n2 are set as the refractive
indices of the mediums on the right side of the first freeform surface, between the two freeform
surfaces and on the left side of the second freeform surface, respectively. Similar to the one
freeform lens design, I, R, O, N and the relationships among them can be expressed as:
(0,0,1)I ? (3.10)
(0,0,1)O? (3.11)
27
2 2 2
1 ( , , )( ) ( ) ( ) q p q p q p
q p q p q p
R x x y y z zx x y y z z? ? ? ?? ? ? ? ? (3.12)
22
1 ( , , 1 )1p p x p y
p x p y
N z zzz? ? ??? (3.13)
220 1 0 1 0 12 ( )pn R n I N n n n n R I? ? ? ? ? (3.14)
Where, xp, yp, zp and xq, yq, zq are the Cartesian coordinates of points p and q on the first and
second freeform surfaces, respectively. And Np is the normal unit vector at point p of the first
freeform lens. Then substitute (3.10), (3.11), (3.12) and (3.13) in (3.14) to get the partial
differential equations:
0
2 2 201
()
( ) ( ) ( ) ( )qppx q p q p q p q p
n x xz
n z z n x x y y z z
???
? ? ? ? ? ? ?
(3.15)
0
2 2 201
()
( ) ( ) ( ) ( )qppy q p q p q p q p
n y yz
n z z n x x y y z z
???
? ? ? ? ? ? ?
(3.16)
The relationship between xp, yp, zp and xq, yq, zq can be calculated based on the constancy of
optical path length.
? ? ? ? ? ? 12 2 2 21 0 2 ( ) .p q p q p q p t qn z n x x y y z z n z z C o n s t??? ? ? ? ? ? ? ? ????? (3.17)
Because I and O are parallel to the optical axis, then xp=xs, yp=ys, xq=x t, yq=yt. zs and zt are
constants. On the other hand, the relationship between xs, ys and xt, yt can be calculated based on
the energy conservation.
( , ) ( , )i n s s s s o u t t t t tI x y d x d y I x y d x d y? (3.18)
28
3.2.2 Based on separated variables mapping
This method is to produce desired irradiance distribution whilst forming a prescribed
wavefront from a given input beam. To simplify the computation, the input-to-output irradiance
mapping is firstly obtained based on Energy conservation using the separated variables method
[50]. Thus, both the input and output irradiance distributions must fulfill the restriction that they
can be factorized in two orthogonal transverse coordinates, which is still suitable for many laser
beam shaping applications. Then, the two freeform optical surfaces are generated simultaneously
and point by point corresponding to the input and output ray sequences defined by the first step.
The method is fast and can be easy understood by optical engineers since it avoid solving the
underlying Monge?Amp?re equation [49].
Consider a two-freeform-refractive-surface optical system shown schematically in Figure
3.3 (one or both of the two optical surfaces could also be reflective). The input and output beam
are supposed as propagating toward the positive z direction and their wavefronts can be
represented as position vectors S = (xs,ys,zs) and T = (xt,yt,zt), respectively. The points on the two
freeform surfaces are denoted by P = (xp,yp,zp) and Q = (xq,yq,zq), respectively. n1, n0 and n2 are
set as the refractive indices of the mediums on the right side of the first freeform surface,
between the two freeform surfaces and on the left side of the second freeform surface,
respectively.
29
Figure 3.3: Geometrical construction of the double freeform optical surfaces for achieving a
specified input-output ray mapping [49].
Let Iin and Iout denote the prescribed irradiance distributions over the planes perpendicular
to the z-axis near the input and output wavefronts, respectively. Their relationship can be
described by Energy conservation written as Eq. (3.19):
( , ) ( , )i n s s s s o u t t t t tI x y d x d y I x y d x d y? (3.19)
Assume that Iin and Iout can be separated into a product of two one-dimensional irradiance
distributions: Iin(xs,ys) = Iin,x(xs) Iin,y(ys) and Iout(xt,yt) = Iout,x(xt) Iout,y(yt). Thus, both (xs, ys) and (xt,
yt) can be numerically specified based on ?source-to-target? [50] or ?target-to-source? [51,52]
variable separation mapping strategies if one of them is predefined. Take ?source-to-target
variable separation mapping strategy for example, if (xs, ys) are equidistantly divided into n?m
rectangular grids: xs=xs,j and ys=ys,i, i=1,2,? n, j=1,2,? m, then xt,j and yt,i can be calculated by
satisfying Eq. (3.20) and Eq. (3.21), respectively.
, , , ,
,1 ,1 ,1 ,1, , , ,( ) ( ) ( ) ( )
s j s n t j t n
s s t t
x y x yi n x s s i n y s s o u t x t t o u t y t t
x y x yI x d x I y d y I x d x I y d y?? ? ? ?
(3.20)
30
, , , ,
,1 ,1 ,1 ,1, , , ,( ) ( ) ( ) ( )
s m s i t m t i
s s t t
x y x yi n x s s i n y s s o u t x t t o u t y t t
x y x yI x d x I y d y I x d x I y d y?? ? ? ?
(3.21)
Then, the input and output unit ray vectors Ini,j and Outi,j can be obtained as the unit normal
vectors at Si,j(xs,j,ys,i,zs,i,j) and Ti,j(xs,j,ys,i,zs,i,j), respectively, as shown in Eq. (3.22) and Eq. (3.23):
, , ,, , ,,, ,,
22
, ,
,
, , 1 1
s j s i x ys j s i x ys j s i s j s i
s s s sij
xy xy
z z z z
x y x y
?? ? ? ? ?? ? ? ?? ? ? ???? ? ? ? ?? ? ? ?? ? ? ?
? ? ? ?? ? ? ?? ? ? ???In
(3.22)
, , ,, , ,,, ,,
22
, ,
,
, , 1 1
t j t i x yt j t i x yt j t i t j t i
t t t tij
xy xy
z z z z
x y x y
?? ? ? ? ?? ? ? ?? ? ? ???? ? ? ? ?? ? ? ?? ? ? ?
? ? ? ?? ? ? ?? ? ? ???O u t
(3.23)
Thus, the input ray sequence is defined both by Ini,j and Si,j, and the output ray sequence is
defined both by Outi,j and Ti,j, i = 1,2,? n, j = 1,2,? m.
The next step is to calculate the data points Pi, j and Qi, j on the two freeform optical
surfaces which are desired to transform the input ray sequence into the output ray sequence.
From two starting points P1,1 and Q1,1 e.g. the central points of the first and second freeform
surfaces, the normal vector N1,1 at P1,1 is firstly calculated using the vector form of Snell?s law
so that the ray emitted from P1,1 can reach Q1,1, as shown in Eq. (3.24):
? ? 221 ,1 0 1 ,1 1 1 ,1 0 1 0 1 1 ,1 1 ,12 ( )n n n n n n? ? ? ? ?N R I n R I n (3.24)
wherein R1,1 denotes the unit vector of the ray through P1,1 and Q1,1. P2,1 can be obtained by
intersecting its corresponding input ray to the tangent plane of P1,1 [53], which can be
formulated as Eq. (3.25) and Eq. (3.26):
? ?2 ,1 1,1 1,1 0? ? ?P P N (3.25)
? ?2 ,1 2 ,1 2 ,1 2 ,1 2 ,1??P S P S = I n (3.26)
Then, Q2,1 is obtained by equaling the OPL of the ray through P2,1 and itself to the OPL of the
ray through P1,1 and Q1,1, as shown in Eq. (3.27) and Eq. (3.28):
31
1 2 , 1 2 , 1 0 2 , 1 2 , 1 2 2 , 1 2 , 1 1 1 , 1 1 , 1 0 1 , 1 1 , 1 2 1 , 1 1 ,1, , , , , ,n n n n n n? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?S P P Q Q T S P P Q Q T
(3.27)
? ?
2 ,1 2 ,1 2 ,1 2 ,1 2 ,1? ? ?T Q T Q O u t
(3.28)
wherein [X,Y] represents the distance between two arbitrary points X and Y. Then, the required
normal vector N2,1 at P2,1 can be computed so that the ray through P2,1 can be refracted to Q2,1.
The entire starting curves on the two surfaces can be generated by repeating the above process.
All the points of the second curve on the first freeform surface can be obtained on the tangent
planes of the starting curve?s points, written as Eq. (3.29) and Eq. (3.30):
? ?,2 ,1 ,1 0i i i? ? ?P P N (3.29)
? ?, 2 , 2 , 2 , 2 , 2i i i i i??P S P S = I n (3.30)
Then, points of the second curve on the second freeform surface can be computed by equa ling
their OPLs to that of the ray through P1,1 and Q1,1, expressed as Eq. (3.31) and Eq. (3.32):
1 , 2 , 2 0 , 2 , 2 2 , 2 , 2 1 , 1 , 1 0 , 1 , 1 2 , 1 , 1, , , , , ,i i i i i i i i i i i in n n n n n? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ?S P P Q Q T S P P Q Q T
(3.31)
? ?
, 2 , 2 , 2 , 2 , 2i i i i i? ? ?T Q T Q O u t
(3.32)
All the normal vectors at the points of the second curve on the first freeform surface can then be
obtained. We can repeat this process to generate all the surface points and interpolate them to
get the entire freeform surfaces [38].
32
CHAPTER 4
DESIGN EXAMPLES AND RESULTS
This chapter provides design examples to the systems discussed in chapter 3. The
parameters of the design examples are calculated with Matlab codes developed. The examples
will be discussed in two main parts: one freeform lens, and two freeform lens. The Matlab and
Tracepro results are presented and discussed separately.
4.1 One freeform lens
Parameters for the design are listed in Table 4.1, which are the values of the basic known
parameters of our design example.
Second, we use the equations derived in chapter 3 to obtain the parameters of the
freeform surface, which are realized in the software Matlab. In detail, the parameters are the
three-dimensional Cartesian Coordinates of the data points on the surface. Then the surface can
be reconstructed in the software Rhinoceros based on the idea of NURBS which is stated in
chapter 2.
Table 4.1 Parameters of Laser shaping system
Parameter Value
Wavelength 870 nm
The glass used BK7
The refractive index at the wavelength 1.509493
33
According to above design steps, four different schemes for the freeform lens with
nonrotational symmetry are constructed.
4.1.1 Incident beam is circular Gaussian beam
The other parameters for this beam shaping system is shown in Table 4.2
W 4mm
L 4mm
Incident Beam radius 2mm
Incident beam waist w0 (e-2 of beam profile) 1mm
Dx 6mm
Dy 6mm
z0 500mm
W and L are the width and length of the freeform lens, respectively. Dx and Dy are the
length respectively in x and y direction of the output beam in the target plane. z0 is the distance
between the first surface of the freeform lens and the target plane.
In the design, the freeform surface is reconstructed with 201?201 points based on
NURBS. Ray tracing is implemented based on the Monte-Carlo method. About 1,000,000 rays
are traced in this beam shaping system, the irradiance distribution of the output beam in the
target plane is shown in Figure 4.1.
Table 4.2 Parameters of one freeform lens system
34
Figure 4.1: Irradiance distribution on the target plane.
From the ray tracing results, the ratio of the power on the desired region to the total
power is about 90.9%.
When the source consists of 3 by 3 laser diode arrays, each lens array is composed of 9
lens, which is used to shape the corresponding laser diode. Figure 4.2 is the layout of the system.
Figure 4.2: Layout of the lens array.
About 10,000 rays are traced in this beam shaping system, the irradiance distribution of
the output beam in the target plane is shown in Figure 4.3.
35
Figure 4.3: Irradiance distribution of the lens array on the target plane.
From the ray tracing results, the ratio of the power on the desired region to the total
power is about 91.52%. According to this result, large area illumination can be realized by using
the idea of lens array.
4.1.2 Incident beam is elliptical Gaussian beam
The other parameters for this beam shaping system is shown in Table 4.3
W 4mm
L 4mm
Incident Beam radius in x direction 2mm
Incident Beam radius in y direction 1.6mm
Incident beam perpendicular waist w0x 1mm
Incident beam parallel waist w0y 0.8mm
Dx 4mm
Dy 4mm
z0 500mm
Table 4.3 Parameters of one freeform lens system
36
Also, W and L are the width and length of the freeform lens, respectively. Dx and Dy are
the length respectively in x and y direction of the output beam in the target plane. z0 is the
distance between the first surface of the freeform lens and the target plane.
In the design, the freeform surface is reconstructed with 201?201 points based on
NURBS. Ray tracing is implemented based on the Monte-Carlo method. About 1,000,000 rays
are traced in this beam shaping system, the irradiance distribution of the output beam in the
target plane is shown in Figure 4.4.
Figure 4.4: Irradiance distribution of the freeform lens on the target plane.
From the ray tracing results, the ratio of the power on the desired region to the total
power is about 91.51%.
When the source consists of 3 by 3 laser diode arrays, each lens array is composed of 9
lens, which is used to shape the corresponding laser diode. Still, about 2,000,000 rays are traced
in this beam shaping system, the irradiance distributions of the incident and output beam in the
input and target plane are shown in Figure 4.5 and Figure 4.6, respectively.
37
Figure 4.5: Irradiance distribution of the freeform lens on the input plane.
Figure 4.6: Irradiance distribution of the freeform lens on the target plane.
From the ray tracing results, the ratio of the power on the desired region to the total
power is about 92.33%.
38
4.2 Two freeform lens
As mentioned in chapter 3, there are two different methods used in designing two
freeform lens. Same parameters are used both in these two methods.
4.2.1 Incident beam is circular Gaussian beam
Parameters for the design are listed in Table 4.4, which are the values of the basic known
parameters of our design example.
Table 4.4 Parameters of two freeform lens system
Incident beam radius 20mm
Incident beam waist (e-2 of beam profile) 10mm
Desired output beam dimension 30mm?30mm
Distance between input and output plane 310mm
In the design, the freeform surface is reconstructed with 101?101 points based on
NURBS. Ray tracing is implemented based on the Monte-Carlo method. About 1,000,000 rays
are traced in this beam shaping system.
If the surface parameters are calculated by PDEs, the irradiance distribution of the output
beam on four receivers placed at 0mm, 200mm, 500mm and 1000mm away from the flat emitted
surface of the secondary lens are shown in Figure 4.7.
39
40
Figure 4.7: In the PDEs method , irradiance distribution on the receivers are placed at (a) 0mm,
(b) 200mm, (c) 500mm, (d) 1000mm away from the emitted surface of the second lens.
If the surface parameters are calculated based on variable separation mapping, then the
irradiance distribution of the output beam on four receivers placed at 0mm, 200mm, 500mm and
1000mm away from the flat emitted surface of the secondary lens are shown in Figure 4.8.
41
42
Figure 4.8: In the variable separation mapping method, irradiance distribution on the receivers
are placed at (a) 0mm, (b) 200mm, (c) 500mm, (d) 1000mm away from the emitted surface of
the second lens.
From the ray tracing results, in the PDEs method, the ratio of the power on the desired
region to the total power is about 90% on four receivers. However, in the variable separation
mapping method, the ratio is only about 60%. According to these, we can find that the previous
method is more effective to design the two freeform lens to realize uniform illumination.
Moreover, the computation is much faster in the PDEs method.
4.2.2 Incident beam is elliptical Gaussian beam
Parameters for the design are listed in Table 4.5, which are the values of the basic known
parameters of our design example.
Table 4.5 Parameters of two freeform lens system
Incident beam radius in x direction 26mm
Incident beam radius in y direction 8mm
Incident beam waist w0x 12.7mm
Incident beam waist w0y 3.5mm
Desired output beam dimension 60mm?60mm
Distance between input and output plane 310mm
In the design, the freeform surface is reconstructed with 101?101 points based on
NURBS. Ray tracing is implemented based on the Monte-Carlo method. About 1,000,000 rays
are traced in this beam shaping system.
43
The irradiance distribution of the incident and output beam in the source is shown in
Figure 4.9.
Figure 4.9: Irradiance distribution of the freeform lens on the input plane.
If the surface parameters are calculated by PDEs, the irradiance distribution of the output
beam on four receivers placed at 0mm, 200mm, 500mm and 1000mm away from the flat emitted
surface of the secondary lens are shown in Figure 4.10.
44
45
Figure 4.10: In the PDEs method, irradiance distribution on the receivers are placed at (a) 0mm,
(b) 200mm, (c) 500mm, (d) 1000mm away from the emitted surface of the second lens.
If the surface parameters are calculated based on variable separation mapping, then the
irradiance distribution of the output beam on four receivers placed at 0mm, 200mm, 500mm and
1000mm away from the flat emitted surface of the secondary lens are shown in Figure 4.11.
46
Figure 4.11: In the variable separation mapping method, irradiance distribution on the receivers
are placed at (a) 0mm, (b) 200mm, (c) 500mm, (d) 1000mm away from the emitted surface of
the second lens.
From the ray tracing results, in the Partial Differential Equations (PDEs) method, the
ratio of the power on the desired region to the total power is about 91% on four receivers.
However, in the variable separation mapping method, the ratio is only about 82%. According to
these, also we can find that the previous method is more effective to design the two freeform lens
to realize uniform illumination. Moreover, the computation is much faster in the Partial
Differential Equations (PDEs) method.
47
Above all, the Partial Differential Equations (PDEs) method is more effective and
efficient to design freeform lens. And the reconstruction points wouldn?t be a big problem in
designing. However, more data points should be used in the variable separation mapping method
if good illumination result is desired, for example, 400?400 points are normally used in
reconstructing the freeform surfaces.
Then Partial Differential Equations (PDEs) method can be used for 3 by 3 laser diode
arrays, each lens array is composed of 9 lens, which is used to shape the corresponding laser
diode. And the incident beam of circular and elliptical Gaussian shape will be reshaped
according to this system.
For circular Gaussian beam, the layout of the system is shown in Figure 4.12, and the
irradiance distributions of the incident and output beam in the source and target plane are shown
in Figure 4.13 and Figure 4.14, respectively, the target plane is placed at 1000mm away from the
emitted surface of the second lens.
Figure 4.12: The layout of the system with 3-by-3 laser diode arrays, the distribution of incident
beam is circular Gaussian.
48
Figure 4.13: Irradiance distribution of the freeform lens on the input plane.
Figure 4.14: Irradiance distribution of the freeform lens on the target plane 1000mm away from
the emitted surface of the second lens.
For elliptical Gaussian beam, the layout of the system is shown in Figure 4.15, and the
irradiance distributions of the incident and output beam in the source and target plane are shown
in Figure 4.16 and Figure 4.17, respectively, the target plane is also placed at 1000mm away
from the emitted surface of the second lens.
49
Figure 4.15: The layout of the system with 3-by-3 laser diode arrays, the distribution of incident
beam is elliptical Gaussian
Figure 4.16: Irradiance distribution of the freeform lens on the input plane
50
Figure 4.17: Irradiance distribution of the freeform lens on the target plane 1000mm away from
the emitted surface of the second lens.
From the stimulated irradiance distributions on the receivers placed at 1000mm away
from the emitted surface of the second lens, which are shown in Figure 4.14 and 4.17, the
uniformities are 90% and 94%, respectively. In consequence, we can know that the PDEs
method is very effective to design the lens array which will reshape irradiance of the laser diode
arrays.
51
CHAPTER 5
CONCLUSIONS
To summarize, the irradiance distribution of the laser diode was mainly concentrated in
this work. And the freeform surfaces design method was proposed for transforming this
irradiance distribution into uniform one. One freeform lens and two freeform lens were applied
to this beam shaping system.
For one freeform lens design, the first order partial differential equations were solved to
obtain the parameters of the freeform surface. The simulation results show that this one freeform
lens system can effectively transform a collimated Gaussian laser beam into a rectangular ?flat-
top? one. Furthermore, the lens array was used for the large area illumination, and the simulation
results are very good, which means this lens design method is very advantageous. With the help
of the lenslet technology, we can simplify the multi-sources shaping problem to the single laser
diode shaping problem.
Meantime, for the two freeform lens system, there are two different methods which can
be used to calculate the surfaces parameters. One is based on solving partial differential
equations, the other is on the basis of variable separation mapping. According to the simulation
results, the previous one, the PDEs based method, is more effectively and efficiently in
transforming a collimated Gaussian laser beam into a ?flat-top? rectangular one with a long
depth of field and small divergence.
52
REFERENCES
[1] Larry Coldren, Scott Corzine, Milan Mashanovitch. ?Diode Lasers and Photonic
Integrated Circuits? (Second Ed.), John Wiley and Sons (2012).
[2] Julio Chaves, Introduction to Nonimaging Optics, CRC Press, 2008.
[3] Mert Serkan, ?Laser beam shaping optical system design methods and their application
in edge-emitting semiconductor laser-based LIDAR systems,? doctoral dissertation.
[4] L. Wang, K. Y. Qian, and Y. Luo, ?Discontinuous free-form lens design for prescribed
irradiance,? Appl. Opt. 46(18), 3716?3723 (2007).
[5] Y. Han, X. Zhang, Z. Feng, K. Qian, H. Li, Y. Luo, X. Li, G. Huang, and B. Zhu,
?Variable-separation three dimensional freeform nonimaging optical system design based on
target-to-source mapping and micro belt surface construction, ? Sciencepaper Online 1?9(2010).
http://www.paper.edu.cn/en/paper.php?serial_number=201002-443.
[6] Y. Luo, Z. Feng, Y. Han, and H. Li, ?Design of compact and smooth free-form optical
system with uniform illuminance for LED source,? Opt. Express 18(9), 9055?9063 (2010).
[7] Computer Desktop Encyclopedia@ 2000 The computer language Co. Inc.
[8] Zexin Feng, Lei Huang, Mali Gong,* and Guofan Jin. ?Beam shaping system design
using double freeform optical surfaces,? Optics Express, Vol. 21, No. 12 / June 2013.
[9] Rengmao Wu, Zhenrong Zheng,* Haifeng Li, and Xu Liu. ?Constructing optical
freeform surfaces using unit tangent vectors of feature data points,? J. Opt. Soc. Am. A / Vol. 28,
No. 9 / September 2011.
53
[10] W. J. Cassarly and M. J. Hayford, ?Illumination optimization: the revolution has begun,?
Proc. SPIE 4832, 258?269 (2002).
[11] T. L. R. Davenport, T. A. Hougha, and W. J. Cassarlya, ?Optimization for illumination
systems: the next level of design,? Proc. SPIE 5456, 81?90 (2004).
[12] T. R. Davenport, ?3D NURBS representation of surfaces for illumination,? Proc. SPIE
4832, 293?301 (2002).
[13] J. Bortz, N. Shatz, and M. Keuper, ?Optimal design of a nonimaging TIR doublet lens for
an illumination system using an LED source,? Proc. SPIE 5529, 8?16 (2004).
[14] M. G. Turner and K. J. Garcia, ?Optimization using rational B?zier control points and
weighting factors,? Proc. SPIE 7061, 70610H (2008).
[15] W. Z. Zhang, Q. X. Liu, H. F. Gao, and F. H. Yu, ?Free-form reflector optimization for
general lighting,? Opt. Eng. 49, 063003 (2010).
[16] Y. Ding, X. Liu, Z. R. Zheng, and P. F. Gu, ?Freeform LED lens for uniform
illumination,? Opt. Express 16, 12958?12966 (2008).
[17] R. Winston and H. Ries, ?Nonimaging reflectors as functionals of the desired
irradiance,? J. Opt. Soc. Am. A 10, 1902?1908 (1993).
[18] A. Rabl and J. M. Gordon, ?Reflector design for illumination with extended sources: the
basic solutions,? Appl. Opt. 33, 6012?6021 (1994).
[19] H. R. Ries and R. Winston, ?Tailored edge-ray reflectors for illumination,? J. Opt. Soc.
Am. A 11, 1260?1264 (1994).
[20] P. T. Ong, J. M. Gordon, and A. Rabl, ?Tailoring lighting reflectors to prescribed
illuminance distributions: compact partialinvolute designs,? Appl. Opt. 34, 7877?7887 (1995).
54
[21] W. Tai and R. Schwarte, ?Design of an aspherical lens to generate a homogenous
irradiance for three-dimensional sensors with a light-emitting diode source,? Appl. Opt. 39,
5801?5805 (2000).
[22] H. Ries and J. Muschaweck, ?Tailored freeform optical surfaces,? J. Opt. Soc. Am. A 19,
590?595 (2002).
[23] V. Oliker, ?Optical design of freeform two-mirror beam-shaping systems,? J. Opt. Soc.
Am. A 24, 3741?3752 (2007).
[24] Z. R. Zheng, X. Hao, and X. Liu, ?Freeform surface lens for LED uniform illumination,?
Appl. Opt. 48, 6627?6634 (2009).
[25] P. Ben?tez, J. C. Mi?ano, J. Blen, R. Mohedano, J. Chaves, O. Dross, M. Hern?ndez, J. L.
Alvarez, and W. Falicoff, ?SMS design method in 3D geometry: examples and applications,?
Proc. SPIE 5185, 18?29 (2004).
[26] F. Mu?oz, P. Ben?tez, O. Dross, J. C. Mi?ano, and B. Paikyn, ?Simultaneous multiple
surface design of compact air-gap collimators for light-emitting diodes?,Opt. Eng. 43, 1522?
1530 (2004).
[27] O. Dross, R. Mohedano, P. Ben?tez, J. C. Mi?ano, J. Chaves, J. Blen, M. Hern?ndez, and
F. Mu?oz, ?Review of SMS design methods and real world applications,? Proc. SPIE 5529, 35?
47 (2004).
[28] Hazewinkel, Michiel, ed. (2001), "Runge-Kutta method", Encyclopedia of Mathematics,
Springer, ISBN 978-1-55608-010-4
[29] http://en.wikipedia.org/wiki/Optical_path_length
[30] Arthur Schuster, An Introduction to the Theory of Optics, London: Edward Arnold, 1904
[31] Ghatak, Ajoy (2009), Optics (4th ed.), ISBN 0-07-338048-2
55
[32] Feynman, Richard, The Feynman Lectures on Physics, Vol. 1, pp. 26?7
[33] Katz, Mikhail G.; Schaps, David; Shnider, Steve (2013), "Almost Equal: The Method of
Adequality from Diophantus to Fermat and Beyond", Perspectives on Science 21 (3): 7750,
arXiv:1210.7750, Bibcode:2012arXiv1210.7750K
[34] B de Greve, ?Reflections and Refractions in Ray tracing.?
http://graphics.stanford.edu/courses/cs148-10-summer/docs/2006--degreve--
reflection_refraction.pdf
[35] Eric Weisstein?s World of Physics, Snell?s Law, http://scienceworld.wolfram.com
[36] http://en.wikipedia.org/wiki/Non-uniform_rational_B-spline
[37] http://web.cs.wpi.edu/~matt/courses/cs563/talks/nurbs.html
[38] Piegl, L and Tiller, W. ?The NURBS book.? New York: Springer, 1997: 188-212.
[39] http://en.wikipedia.org/wiki/Tophat_beam
[40] David L. Shealya and John A. Hoffnagleb, ?Beam shaping profiles and propagation.?
SPIE Conf. Laser Beam Shaping VI, Proc. 5876-13, 2005.
[41] http://en.wikipedia.org/wiki/Monte_Carlo_method
[42] Anderson, H.L. (1986). "Metropolis, Monte Carlo and the MANIAC". Los Alamos
Science 14: 96?108.
[43] Baeurle, Stephan A. (2009). "Multiscale modeling of polymer materials using field-
theoretic methodologies: A survey about recent developments". Journal of Mathematical
Chemistry 46 (2): 363?426. doi:10.1007/s10910-008-9467-3.
[44] H. Ries and J. Muschaweck, ?Tailoring freeform lenses for illumination,? Proc. SPIE
6338, 633808 (2001).
56
[45] H. Ries and J. Muschaweck, ?Tailored freeform optical surfaces,? J. Opt. Soc. Am. A 19,
590-595 (2002).
[46] Y. Ding and P. F. Gu, ?The Freeform Reflector for Uniform Illumination,? Acta Optica
Sinica 27, 540-544 (2007).
[47] Y. Ding, X. Liu, H. F. Li, and P. F. Gu, ?The design of the freeform reflector for uniform
illumination,? in Proceedings of Asia Display 2007, Volume 1. (Shanghai, China, 2007), pp.
735-738.
[48] Y. Ding, X. Liu, Z. R. Zheng, and P. F. Gu, ?Freeform LED lens for uniform
illumination,? Opt. Express 16(17), 12958?12966 (2008).
[49] Zexin Feng, Lei Huang, Mali Gong, and Guofan Jin, ?Beam shaping system design using
double freeform optical surfaces,? Opt. Express Vol. 21, No. 12 (2013).
[50] L. Wang, K. Y. Qian, and Y. Luo, ?Discontinuous free-form lens design for prescribed
irradiance,? Appl. Opt. 46(18), 3716?3723 (2007).
[51] Y. Han, X. Zhang, Z. Feng, K. Qian, H. Li, Y. Luo, X. Li, G. Huang, and B. Zhu,
?Variable-separation three dimensional freeform nonimaging optical system design based on
target-to-source mapping and micro belt surface construction,? Sciencepaper Online 1?9(2010).
[52] Y. Luo, Z. Feng, Y. Han, and H. Li, ?Design of compact and smooth free-form optical
system with uniform illuminance for LED source,? Opt. Express 18(9), 9055?9063 (2010).
[53] W. B. Elmer, The optical design of reflectors, 2nd ed. (Wiley, New York, 1980).
[54] F. M. Dickey, S. C. Holswade, and D. L. Shealy, eds., Laser Beam Shaping Applications
(CRC Press, 2005).
57
[55] Haotong Ma, Zejin Liu, Pengzhi Jiang, Xiaojun Xu, and Shaojun Du, ?Improvement of
Galilean refractive beam shaping system for accurately generating near-diffraction-limited
flattop beam with arbitrary beam size?, Opt. Express Vol. 19, No. 14.
58
Appendix
1. Two freeform lens design based on partial differential equations
clear all,clc
format long
global z0 nO nI W L N
e1=0;
t1=cputime;
W = 20;
L = 20;
Lx = L;
Ly = W;
w0x = 5;
w0y = 5;
nO = 1;
nI = 1.5095; % BK7 -- 870nm
n1 = nI ;
n2 = nI ;
N = 200;
Dx = 40; % [mm] the size of the output beam in the target plane
Dy = 40;
I0 = 1; % the intensity for the input beam
fzx = @(xp,yp,xq,yq,zp,zq) nO*(xp-xq)/(nO*(zq-zp)-n1*sqrt( (xq-xp)^2 + (yq-yp)^2 + (zq-
zp)^2 )) ;
fzy = @(xp,yp,xq,yq,zp,zq) nO*(yp-yq)/(nO*(zq-zp)-n1*sqrt( (xq-xp)^2 + (yq-yp)^2 + (zq-
zp)^2 )) ;
fun_x_t=@(x,y) exp(-2*(x./w0x).^2).*exp(-2*(y./w0y).^2);
% fplus = @(x,y,tx,ty,z) fzx + fzy;
% % X = zeros(1,N+1);
% % Y = zeros(1,N+1);
% % Z = zeros(N+1,N+1);
% % % % h1 = Lx/N ;
% % % % h2 = Ly/N ;
% % % %
% % % % X = -Lx/2 : h1 : Lx/2 ;
% % % % Y = -Ly/2 : h2 : Ly/2 ;
dx = L/2/N;
h1 = 2*dx;
59
h2 = 2*dx;
for i=1:N+1
for j=1:N+1
%Space axis
X(j)=2*((j-1)*dx-L/4);
Y(i)=2*((i-1)*dx-L/4);
end
end
E0 = I0 / Dx / Dy * quad2d(fun_x_t,-Lx/2,Lx/2,-Ly/2,Ly/2); % normalization coefficient ---
for the output beam
% tx = I0 / Dy / E0 * quad2d(fun_x_t,-Lx/2,X,-Ly/2,Ly/2);
% ty = I0 / Dx / E0 * quad2d(fun_x_t,-Lx/2,Lx/2,-Ly/2,Y);
for j = 1:N+1
tx(j) = I0 / Dy / E0 * quad2d(fun_x_t,-Lx/2,X(j),-Ly/2,Ly/2);
end
for i=1:N+1
ty(i) = I0 / Dx / E0 * quad2d(fun_x_t,-Lx/2,Lx/2,-Ly/2,Y(i));
end
tx = tx - Dx/2;
ty = ty - Dy/2;
%% P & Q & S & T
x_s = X ;
y_s = Y ;
x_p = X ;
y_p = Y ;
x_q = tx ;
y_q = ty ;
x_t = tx ;
y_t = ty ;
%% Boundary condition
z_p(1,1) = 8.5 ;
z_p(1,N+1) = 8.5 ;
z_p(N+1,1) = 8.5 ;
z_p(N+1,N+1) = 8.5 ;
z_q(1,1) = 295.5 ;
z_q(1,N+1) = 295.5 ;
z_q(N+1,1) = 295.5 ;
z_q(N+1,N+1) = 295.5 ;
z_s = zeros(N+1,N+1) ;
z_t = 300 * ones(N+1,N+1) ;
% Constancy of the Optical Path Length
Const = n1 * (z_p(1,1) - z_s(1,1)) + nO * sqrt( (x_q(1)-x_p(1))^2 + (y_q(1)-y_p(1))^2 +
(z_q(1,1)-z_p(1,1))^2 ) + n2 * (z_t(1,1)-z_q(1,1)) ;
%% z_p(x,-Ly/2) & z_q(x,-Dy/2)
i=1;
for j = 2:N
60
f1 = feval(fzx,x_p(j-1),y_p(i),x_q(j-1),y_q(i),z_p(i,j-1),z_q(i,j-1));
% step 2
xq2 = I0 / Dy / E0 * quad2d(fun_x_t,-Lx/2,x_p(j-1)+h1/2,-Ly/2,Ly/2) - Dx/2;
zp2 = z_p(i,j-1)+ h1/2*f1 ;
ff2 = @(zq2) ([n1 * (zp2- z_s(i,j-1)) + nO * sqrt( (xq2-x_p(j-1)-h1/2)^2 + (y_q(i)-
y_p(i))^2 + (zq2-zp2)^2 ) + n2 * (z_t(i,j-1)-zq2) - Const]) ;
zq2 = fsolve(ff2,[z_p(i,j-1)]);
f2 = feval(fzx,x_p(j-1)+h1/2,y_p(i),xq2,y_q(i),zp2, zq2);
% step 3
zp3 = z_p(i,j-1)+ h1/2*f2 ;
ff3 = @(zq3) ([n1 * (zp3- z_s(i,j-1)) + nO * sqrt( (xq2-x_p(j-1)-h1/2)^2 + (y_q(i)-
y_p(i))^2 + (zq3-zp3)^2 ) + n2 * (z_t(i,j-1)-zq3) - Const]) ;
zq3 = fsolve(ff3,[z_p(i,j-1)]);
f3 = feval(fzx,x_p(j-1)+h1/2,y_p(i),xq2,y_q(i),zp3, zq3);
% step 4
xq3 = I0 / Dy / E0 * quad2d(fun_x_t,-Lx/2,x_p(j-1)+h1,-Ly/2,Ly/2) - Dx/2;
zp4 = z_p(i,j-1)+ h1*f3 ;
ff4 = @(zq4) ([n1 * (zp4- z_s(i,j-1)) + nO * sqrt( (xq3-x_p(j-1)-h1)^2 + (y_q(i)-y_p(i))^2
+ (zq4-zp4)^2 ) + n2 * (z_t(i,j-1)-zq4) - Const]) ;
zq4 = fsolve(ff4,[x_p(i,j-1)]);
f4 = feval(fzx,x_p(j-1)+h1,y_p(i),xq3,y_q(i),zp4, zq4);
z_p(i,j) = z_p(i,j-1) + h1*(f1+2*f2+2*f3+f4)/6;
ff5 = @(zq5) ([n1 * (zp4- z_s(i,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (y_q(i)-y_p(i))^2 +
(zq5-z_p(i,j))^2 ) + n2 * (z_t(i,j)-zq5) - Const]) ;
z_q(i,j) = fsolve(ff5,[z_p(i,j)]);
end
%% z_p(x,Ly/2) & z_q(x,Dy/2)
i=N+1;
for j = 2:N
f1 = feval(fzx,x_p(j-1),y_p(i),x_q(j-1),y_q(i),z_p(i,j-1),z_q(i,j-1));
% step 2
xq2 = I0 / Dy / E0 * quad2d(fun_x_t,-Lx/2,x_p(j-1)+h1/2,-Ly/2,Ly/2) - Dx/2;
zp2 = z_p(i,j-1)+ h1/2*f1 ;
ff2 = @(zq2) ([n1 * (zp2- z_s(i,j-1)) + nO * sqrt( (xq2-x_p(j-1)-h1/2)^2 + (y_q(i)-
y_p(i))^2 + (zq2-zp2)^2 ) + n2 * (z_t(i,j-1)-zq2) - Const]) ;
zq2 = fsolve(ff2,[z_p(i,j-1)]);
f2 = feval(fzx,x_p(j-1)+h1/2,y_p(i),xq2,y_q(i),zp2, zq2);
% step 3
zp3 = z_p(i,j-1)+ h1/2*f2 ;
ff3 = @(zq3) ([n1 * (zp3- z_s(i,j-1)) + nO * sqrt( (xq2-x_p(j-1)-h1/2)^2 + (y_q(i)-
y_p(i))^2 + (zq3-zp3)^2 ) + n2 * (z_t(i,j-1)-zq3) - Const]) ;
zq3 = fsolve(ff3,[z_p(i,j-1)]);
f3 = feval(fzx,x_p(j-1)+h1/2,y_p(i),xq2,y_q(i),zp3, zq3);
% step 4
xq3 = I0 / Dy / E0 * quad2d(fun_x_t,-Lx/2,x_p(j-1)+h1,-Ly/2,Ly/2) - Dx/2;
zp4 = z_p(i,j-1)+ h1*f3 ;
61
ff4 = @(zq4) ([n1 * (zp4- z_s(i,j-1)) + nO * sqrt( (xq3-x_p(j-1)-h1)^2 + (y_q(i)-y_p(i))^2
+ (zq4-zp4)^2 ) + n2 * (z_t(i,j-1)-zq4) - Const]) ;
zq4 = fsolve(ff4,[z_p(i,j-1)]);
f4 = feval(fzx,x_p(j-1)+h1,y_p(i),xq3,y_q(i),zp4, zq4);
z_p(i,j) = z_p(i,j-1) + h1*(f1+2*f2+2*f3+f4)/6;
ff5 = @(zq5) ([n1 * (z_p(i,j)- z_s(i,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (y_q(i)-y_p(i))^2
+ (zq5-z_p(i,j))^2 ) + n2 * (z_t(i,j)-zq5) - Const]) ;
z_q(i,j) = fsolve(ff5,[z_p(i,j)]);
end
%% z_p(-Lx/2,y) & z_q(-Dx/2,y)
j = 1 ;
for i = 2:N
g1 = feval(fzy,x_p(j),y_p(i-1),x_q(j),y_q(i-1),z_p(i-1,j),z_q(i-1,j));
% step 2
yq2 = I0 / Dx / E0 * quad2d(fun_x_t,-Lx/2,Lx/2,-Ly/2,y_p(i-1)+h2/2) - Dy/2;
zp2 = z_p(i-1,j)+ h2/2*g1 ;
gg2 = @(zq2) ([n1 * (zp2 - z_s(i-1,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (yq2-y_p(i-1)-
h2/2)^2 + (zq2-zp2)^2 ) + n2 * (z_t(i-1,j)-zq2) - Const]) ;
zq2 = fsolve(gg2,[z_p(i-1,j)]);
g2 = feval(fzy,x_p(j),y_p(i-1)+h2/2,x_q(j),yq2,zp2,zq2);
% step 3
zp3 = z_p(i-1,j)+ h2/2*g2 ;
gg3 = @(zq3) ([n1 * (zp3 - z_s(i-1,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (yq2-y_p(i-1)-
h2/2)^2 + (zq3-zp3)^2 ) + n2 * (z_t(i-1,j)-zq3) - Const]) ;
zq3 = fsolve(gg3,[z_p(i-1,j)]);
g3 = feval(fzy,x_p(j),y_p(i-1)+h2/2,x_q(j),yq2,zp3,zq3);
% step 4
yq3 = I0 / Dx / E0 * quad2d(fun_x_t,-Lx/2,Lx/2,-Ly/2,y_p(i-1)+h2) - Dy/2;
zp4 = z_p(i-1,j)+ h2*g3 ;
gg4 = @(zq4) ([n1 * (zp4 - z_s(i-1,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (yq3-y_p(i-1)-h2)^2
+ (zq4-zp4)^2 ) + n2 * (z_t(i-1,j)-zq4) - Const]) ;
zq4 = fsolve(gg4,[z_p(i-1,j)]);
g4 = feval(fzy,x_p(j),y_p(i-1)+h2,x_q(j),yq3,zp4,zq4);
z_p(i,j) = z_p(i-1,j) + h2*(g1+2*g2+2*g3+g4)/6;
gg5 = @(zq5) ([n1 * (z_p(i,j)- z_s(i,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (y_q(i)-y_p(i))^2 +
(zq5-z_p(i,j))^2 ) + n2 * (z_t(i-1,j)-zq5) - Const]) ;
z_q(i,j) = fsolve(gg5,[z_p(i,j)]);
end
%% other points
for j = 2:N+1
for i = 2:N
g1 = feval(fzy,x_p(j),y_p(i-1),x_q(j),y_q(i-1),z_p(i-1,j),z_q(i-1,j));
% step 2
yq2 = I0 / Dx / E0 * quad2d(fun_x_t,-Lx/2,Lx/2,-Ly/2,y_p(i-1)+h2/2) - Dy/2;
zp2 = z_p(i-1,j)+ h2/2*g1 ;
62
gg2 = @(zq2) ([n1 * (zp2 - z_s(i-1,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (yq2-y_p(i-1)-
h2/2)^2 + (zq2-zp2)^2 ) + n2 * (z_t(i-1,j)-zq2) - Const]) ;
zq2 = fsolve(gg2,[z_p(i-1,j)]);
g2 = feval(fzy,x_p(j),y_p(i-1)+h2/2,x_q(j),yq2,zp2,zq2);
% step 3
zp3 = z_p(i-1,j)+ h2/2*g2 ;
gg3 = @(zq3) ([n1 * (zp3 - z_s(i-1,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (yq2-y_p(i-1)-
h2/2)^2 + (zq3-zp3)^2 ) + n2 * (z_t(i-1,j)-zq3) - Const]) ;
zq3 = fsolve(gg3,[z_p(i-1,j)]);
g3 = feval(fzy,x_p(j),y_p(i-1)+h2/2,x_q(j),yq2,zp3,zq3);
% step 4
yq3 = I0 / Dx / E0 * quad2d(fun_x_t,-Lx/2,Lx/2,-Ly/2,y_p(i-1)+h2) - Dy/2;
zp4 = z_p(i-1,j)+ h2*g3 ;
gg4 = @(zq4) ([n1 * (zp4 - z_s(i-1,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (yq3-y_p(i-1)-h2)^2
+ (zq4-zp4)^2 ) + n2 * (z_t(i-1,j)-zq4) - Const]) ;
zq4 = fsolve(gg4,[z_p(i-1,j)]);
g4 = feval(fzy,x_p(j),y_p(i-1)+h2,x_q(j),yq3,zp4,zq4);
z_p(i,j) = z_p(i-1,j) + h2*(g1+2*g2+2*g3+g4)/6;
gg5 = @(zq5) ([n1 * (z_p(i,j)- z_s(i,j)) + nO * sqrt( (x_q(j)-x_p(j))^2 + (y_q(i)-y_p(i))^2 +
(zq5-z_p(i,j))^2 ) + n2 * (z_t(i-1,j)-zq5) - Const]) ;
z_q(i,j) = fsolve(gg5,[z_p(i,j)]);
end
end
%% change the sign of y-axis
X1 = X; % unchangable
Y1 = -Y;
tx1 = tx ;
ty1 = -ty ;
z_p = flipud(z_p);
z_q = flipud(z_q);
%% rhino
for i=1:N+1
for j=1:N+1
x_S(i,j) = X1(j);
x_P(i,j) = X1(j);
x_Q(i,j) = tx1(j);
x_T(i,j) = tx1(j);
y_S(i,j) = Y1(i);
y_P(i,j) = Y1(i);
y_Q(i,j) = ty1(i);
y_T(i,j) = ty1(i);
end
end
z_S = z_s ;
z_P = z_p ;
63
z_Q = z_q ;
z_T = z_t ;
ii=N+1;
jj=ii*ii;
S_rhino=zeros(jj,3);
for j=1:ii
S_rhino((j-1)*ii+1:j*ii,1)=x_S(:,j);
S_rhino((j-1)*ii+1:j*ii,2)=y_S(:,j);
S_rhino((j-1)*ii+1:j*ii,3)=z_S(:,j);
end
P_rhino=zeros(jj,3);
for j=1:ii
P_rhino((j-1)*ii+1:j*ii,1)=x_P(:,j);
P_rhino((j-1)*ii+1:j*ii,2)=y_P(:,j);
P_rhino((j-1)*ii+1:j*ii,3)=z_P(:,j);
end
Q_rhino=zeros(jj,3);
for j=1:ii
Q_rhino((j-1)*ii+1:j*ii,1)=x_Q(:,j);
Q_rhino((j-1)*ii+1:j*ii,2)=y_Q(:,j);
Q_rhino((j-1)*ii+1:j*ii,3)=z_Q(:,j);
end
T_rhino=zeros(jj,3);
for j=1:ii
T_rhino((j-1)*ii+1:j*ii,1)=x_T(:,j);
T_rhino((j-1)*ii+1:j*ii,2)=y_T(:,j);
T_rhino((j-1)*ii+1:j*ii,3)=z_T(:,j);
end
e1=cputime-t1+e1;
figure
pcolor(x_P,y_P,z_P);
figure
pcolor(x_Q,y_Q,z_Q);
2. Two freeform lens design based on separated variables mapping
clear all,clc
format long
% OE_2013_Beam shaping system design using double freeform optical surfaces
w0=10; x0=15; % distance from the respective axes at which the intensity has fallen to half
its axial value
y0=15;
betax=40;
betay=40;
fun_in=@(x,y) exp(-2*(x./w0).^2).*exp(-2*(y./w0).^2);
Q_in=quad2d(fun_in,-20,20,-20,20);
fun_out=@(x,y) ((1+exp(betax*(abs(x)./x0-1))).^-1).*((1+exp(betay*(abs(y)./y0-1))).^-1);
64
Q_out=quad2d(fun_out,-15,15,-15,15);
AxAy=Q_in./Q_out; % [xs,ys]=meshgrid(-10:20/400:10,-10:20/400:10);
% xs=linspace(-10,10,401);
% ys=linspace(-10,10,401);
% xt=zeros(1,401);
% yt=zeros(1,401);
% xt(1)=-20;
% yt(1)=-5;
% % % % % % % % % %
N=100;
% % % % % % % % % %
L=20;
% dx : step size
dx=L/N;
for i=1:N+1
for j=1:N+1
%Space axis
xs(j)=2*((j-1)*dx-L/2);
ys(i)=2*((i-1)*dx-L/2);
% Gau(n,m)=exp(-(x(m)^2+y(n)^2)/(wo^2));
% FDrect(n,m)=((1+exp(betax*(abs(x(m))/x0-1)))^-1)*((1+exp(betay*(abs(y(n))/y0-1)))^-
1);
% r=sqrt(x(m)^2+y(n)^2);
% FDGau(n,m)=(1+exp(beta*(r/R0-1)))^(-1);
end
end
xtt(1)=-15;
xtt(N+1)=15;
ytt(1)=-15;
ytt(N+1)=15;
% for i=1:size(rrr,2)-1
% rrr(rrr==0)=eps;
% RRR(RRR==0)=eps;
% kk1=((2/pi/g0/(w0)^2)*(rrr(i)/RRR(i))*(1+exp(beta*(RRR(i)/R0-1)))*exp(-
2*(rrr(i))^2/w0^2));
% tt=RRR(i)+kk1*h;
% tt(tt==0)=eps;
% kk2=((2/pi/g0/(w0)^2)*(rrr(i+1)/tt)*(1+exp(beta*(tt/R0-1)))*exp(-
2*(rrr(i+1))^2/w0^2));
% RRR(i+1)=RRR(i)+h*(kk1+kk2)/2;
% end
% h=0.0001;
% for m=1:size(xs,2)-1
% kk1=A3*((1+exp(betax*(abs(xt(m))./x0-1))).^-1).*((exp(-2*(xs(m)./w0).^2)).^-1);
% tt=xt(m)+kk1*h;
65
% kk2=A3*((1+exp(betax*(abs(tt)./x0-1))).^-1).*((exp(-2*(xs(m+1)./w0).^2)).^-1);
% xt(m+1)=xt(m)+h*(kk1+kk2)/2;
% end
% A3=A1*A2/AxAy/A4; %
% fun3=@(xx,t) (1+exp(betax*(abs(xx)./x0-1))).^-1;
% xx0=xt(1);
% dxx=0.1;
% yy0=0;
% yy_end=0;
% e1=0;
% while yy_end