Distribution of temperature in a single lens due to absorption of light and heat conduction: an adaptive solver

We develop an algorithm for the solution of the stationary heat-equation in a single lens due to absorption of light, heat-conduction and transfer of the heat to the environment while we assume rotational symmetry for the whole situation. The proceeding is based on an easy to implement ﬁnite difference scheme, which is best suited for rectangular areas. Therefore, we have to transform the heat equation and the boundary conditions from the original domain, i.e. the surface of section of the lens by the aid of tensor methods to a rectangle. So the algorithm generates a grid, which adopts automatically to the actual shape of the lens. In this sense, we characterize the method as adaptive. In the examples, we investigate the effect of a high-transmission glass on the distribution of temperature and further demonstrate the adjustment to a realistic lens shape with a strong deviation from a spherical surface in form of a kink near the edge. We compare the results with a simple model for the distribution of temperature and show the strong dependency of the results on the transmission of the materials.


INTRODUCTION
Depending on the application, raising the power of light transferred through optical systems may have very different advantages, reaching from the reduction of machining times in laser-material processing to a more intense experience of movies due to modern digital cinema projectors [21].
Thus, the search for optical glasses for high power light flux applications has a history of at least two decades, see e.g.[6] and the literature cited there, but due to new powerful sources of light their importance is increasing, so that new high-power glasses are developped and offered [21].This high power flux, however, does not only raise the temperature in the optical elements homogeneously due to the absorption of light, but is also the source for an inhomogeneous distribution of temperature, which in turn influences considerable the optical quality of the system.Measurements of such effects are challenging but of great interest and are therefore an issue in latest research [17].
Here we want to calculate the distribution of temperature in a single lens due to the absorption of light, heat conduction and the thermal interaction with the environment.The relevant scientific basics in the sense of a phenomenological description have been known in principle for a long time [42], even if some of the historical explanations are antiquated today.The development of this formalism marks a glorious era in early classical theoretical and experimental physics, beginning, very simplified speaking, with works of I. Newton and the legendary Théorie analytique de la chaleur of J.B. Fourier, which was called the 'bible of the mathematical physicist' by A. Sommerfeld [56], and culminating finally in the discovery of the conservation of energy.In an earlier work [41], we discussed a few of the basic topics of the theory very briefly and derived approximations for a thin single lens.
The task to be solved here is given by Poisson's equation, which is one of the most fundamental equations in mathematical physics.Hence, for the solution of this class of problems, there is no absence of numerical schemes (see, e.g.[7], [8], [20], [23], [29], [30], [32], [37], [48], [62], or the literature on electrostatic field calculation).In the following, we will choose a way which is taylored for the shape of lenses, may they be aspheres or not.
Generally, we will neglect the transport of heat in the interior of the lens by radiation which may get important for higher differences of temperature.First we will set up the stage and state the problem in mathematical form, i.e., we write down the differential equation and the boundary conditions.In the second step, we will transform the coordinates in such a way that the new domain is a rectangle.For the description of this transformation we use a few elements of tensor calculus which will be briefly explained.Then the problem will be discretized, leading to a system of linear equations which may be solved by different methods.Here, one of the iterative strategies is suggested to avoid the need of additional linear algebra subprograms.
We will also look at a few examples, and compare the results with a simple approximation.
FIG. 1 The domain where the problem has to be solved.Boundary no. 1 is the optical axis, boundary no. 2 and 4 are the upper ("t: top") and lower ("b: bottom") optical surfaces, respectively, and boundary 3 is the edge of the lens.D is the center thickness and R is half the diameter of the lens.n t is is the outer unit normal of the upper (t:"top") surface.For later convenience, we will sometimes write x 1 for ρ and x 2 for z.
Finally, we discuss hints for further generalizations.

MATHEMATICAL STATEMENT OF THE PROBLEM
The theoretical basics and sign conventions as well as the notation of the thermal quantities can be found in more detail in [34], [41].
The original problem is a 3D one, but due to symmetry around the z-axis, which we assume here, it is sufficient to consider the surface of section of the lens shown in Figure 1.
In the inner of the volume, the distribution of temperature is governed by the stationary heat equation [34]: where the last equation is in cylindrical coordinates, having dropped the additional term − 1 ∂φ 2 T(r) due to our assumption of rotational symmetry.
Here, λ is the heat conductivity in W K•m , and the source term ν(r) in W m 3 is the spatial density of the heat produced by absorption of light.
Due to the theory of partial differential equations, the statement of the problem is only complete with boundary conditions, which we will consider in the sequence given by the numbering of the boundaries in Figure 1: 1.On the optical axis, we have: One way of comprehension bases on Gauss' theorem.Consider a cylindrical volume-element in the inner of the lens around the optical axis with radius ρ and height dz and integrate over the surface resp.volume of this cylinder: For small ρ, the left hand side of this equation is where our assumption of rotational symmetry around the optical axis causes the actual form of the last term.Now, ν shall be a smooth function and therefore for small ρ, the right hand side is Since the right hand side of this equation vanishes like O(ρ 2 ) for ρ → 0, the LHS does so, too.This gives 2. If there is no transfer of heat to the environment through the optical surfaces, the current density of heat and thus the gradient of temperature has to be parallel to the tangent vector of the surface.In other words, the dot product between the unit normal vector n t (which is shown schematically in Figure 1) directing outwards of the surface and the gradient of temperature vanishes.
As a generalization, one can take the exchange of heat with the environment into account.To do this, we first need the definition of the current density of heat [34]: The component of the current density in normal direction is therefore given by: This definition of the current density implies a flow of heat outside the volume if j n is positive.If there are no sources of heat at the surface, j n has to be continuous across the boundary [50].
A well established and probably the most simple way to model the exchange of heat with the ambient surroundings is by an empirical law, the so-called Newton's law of cooling (see, e.g.[33], [44], [50], [59], [60], and many others); it is also used for optical components e.g.[4], [24], [61].
According to e.g.[47], the combined action of conduction, convection and radiation can be described by where α N is the Newtonian heat-transfer coefficient W K•m 2 and T s the temperature of the ambient surroundings.
Finally, we find at the optical surfaces: with α N = 0, if the transfer of heat through the optical surfaces is surpressed.
In the case of the upper surface no. 2, we have with the outer unit normal n t ("t" : top) Here, I ρ and I z are the orthonormal cylindrical base vectors in ρ and z-direction, respectively.
3. For simplicity we will choose although more complicated boundary conditions could be treated as well.Notice, however, that the solution will only be fixed uniquely, if not all the boundary conditions are given by derivatives of the temperature; otherwise one could add any constant value to the distribution of temperature in the lens, i.e. to one solution, to get a new one.
4. Here, with the outer unit normal n b ("b": bottom), we have the same conditions as in item no. 2, i.e., The finite difference method which we will use here is best suited for rectangular domains, and in a former calculation of the author the treatment of boundary points and even the creation of an accurate grid caused problems for some lens shapes.Thus, it seemed likely to transform the surface of section of the lens to a rectangle.The price to pay is the transformation of the differential equation and the boundary conditions, which will be done in the next sections.

TRANSFORMATION OF THE GEOMETRY
For convenience, we will sometimes write x 1 for ρ and x 2 for z, and will use upper indices here to be consistent with the tensor notation later.They should not be confused with powers; the latter will be marked by brackets like [.] 2 throughout the whole article.
We considered the surface of section of the lens in cylindrical (ρ|z) resp.pseudo-cartesian (x 1 |x 2 ) coordinates.Now we want to transform this domain into the logical rectangle Mathematically this could be interpreted as a mapping (the "alibi" oder "active" point of view), but we will take this transformation as a description of the original points by the aid of new coordinates (the "alias" oder "passive" point of view) [26].So here the points do not change their location (that may eventually be senseful e.g. for the treating of thermal expansion, depending on the way of description), but they change their name, i.e. coordinates.Such a change of the coordinate system is common practice in physics and engineering.The temperature T(q(x)) at a fixed point in space is the same, regardless of the way of denomination of the coordinates.Therefore, scalar quantities like the temperature need not be transformed.The derivatives, however, obviously do depend on the chosen coordinate system which will be treated by using the cain rule of multi-varibale calculus.The transformation could be done directly from 3D cartesian coordinates, but it seems to be simpler to choose cylindrical coordinates as a starting point.

The transformation equations of the coordinates
We demand for the following properties of the transformation: 1. Points with x 1 = 0 (i.e. the optical axis) shall be mapped to points with q 1 = 0: 2. A point from the upper surface shall be mapped on a point q = q 1 L 2 : 3. Points with x 1 = R shall be mapped to points with 4. For a point of the lower surface, it shall be: A possible transformation of the coordinates is given by: with and where H t (...), H b (...), and H(...) are considered as functions of x 1 .The inverse transformation is therefore The choice of the values of L 1;2 will be considered later.
As one easily can verify, this transformation maps the points in Figure 2 P → P', Q → Q', and R → R', S → S', respectively, and points on the left, upper, right, and lower boundary are mapped to their corresponding counterparts.

Elements of tensor calculus
Unlike a change to cylindrical or spherical coordinates, our transformation generates a mesh which is not locally orthogonal.Thus, the adequate formalism is that of tensor analysis.Therefore, we want to fix the notation and briefly provide the basics of tensor calculus used here.Mainly the notation bases on [2].That book we found very useful for our purposes here and recommend it for further reading in case of need.
We will make use of the summation convention, i.e., we will sum up over two indices if they are symbolized by the same letter in a term, one as a sub-and one as a superscript.
In this section, it is simpler to consider the general threedimensional case because the distinction between coordinate lines resp.coordinate surfaces is clearer.So here, our indices take the values 1, 2, 3, while later it will be possible to restrict our calculations to the two-dimensional case.
For the ordinary orthonormal cartesian base vectors we will write I k or I m , where the indices like k and m take the values 1, 2, 3.The orthonormality of the cartesian base vectors gives for the dot product where We now insert new coordinates q l by an invertible transformation of the cartesian coordinates x n : the points in space are now functions of the new coordinates r = r(q 1 , q 2 , q 3 ) =   x 1 (q 1 , q 2 , q 3 ) x 2 (q 1 , q 2 , q 3 ) x 3 (q 1 , q 2 , q 3 )   = x 1 (q 1 , q 2 , q 3 ) I 1 + x 2 (q 1 , q 2 , q 3 ) I 2 + x 3 (q 1 , q 2 , q 3 ) I 3 = x m (q 1 , q 2 , q 3 ) I m .
We get the coordinate line belonging to q 1 if the values of the coordinates q 2 and q 3 are fixed, while q 1 is varied.So, we find the tangent vector to the coordinate line by Three such tangent vectors may define a local base, the socalled covariant base (l = 1, 2, 3): where we have defined A coordinate surface for q 1 is given by the (in terms of x m : implicit) function and the vector is the gradient to this surface and therefore orthogonal to it in the point (x 1 , x 2 , x 3 ).By the aid of three such vectors, we define the contra-variant base (k = 1, 2, 3): with From the chain rule we get the following Thus, the dot product between the co-and contra-variant base vectors is: So the bases G l and G k are said to be reciprocal to each other, while the cartesian base is the reciprocal to itself, see Eq. ( 21).Hence, for cartesian systems, a distinction between co-and contravariant indices is neither necessary nor common practice in most circumstances.
Again, from the chain rule, we have If we multiply Eq. ( 30) with α n k , we therefore find or By renaming of indices, Eq. ( 26) can be written as Now, we are able to write down the laws for transformation of the components of a vector V given in cartesian coordinates: i.e., and Finally, we want to find the formula for the partial derivative of a scalar function ϕ(x m ).In cartesian coordinates, the gradient of ϕ is given by: where we have used the comma-notation for the partial derivative.
In the new coordinates, the gradient is according to the chain rule given by: Hence, it is

Differentials of the transformation
Now we will turn back to our effective two-dimensional problem; the indices will take the values 1, 2.
The differentials of the transformation, esp. the β's will become important later, so we calculate the quantities defined in Eq. ( 27) and Eq. ( 31): For the α l m 's we find The β m k 's are given by As one can see easily, our transformation comes true with the requirements Eq. ( 32) and (34), which are eight single equations in total.
Later, we will need the derivatives of the β's, too: (47) From the point of view of tensor calculus, there is a relation between the α m i 's or β m i 's and the metric tensor or the derivatives of the differentials and the Christoffel symbols.Moreover, the Laplacian ∆ in the transformed system may be computed from the metric tensor, see e.g.[26].But here we will do the transformation of the Laplacian by straightforward calculation and such a deep entering into mathematical theory is not necessary, since we do not need e.g.derivatives of a vector [2].

Transformation of the differential equation
Clearly, the functional dependency on the variables is a different one in T(q) resp.T(x), and thus it might be formally more correct to write for the transformation (see, e.g.[11]) in what follows.But here, we want to emphasize the fact that at one and the same point in space, the value of the temperature is the same and thus take the same symbol T for both sets of coordinates, as customary in physics and engineering.
We are starting with the chain rule: For a function F q k (x m ) it is therefore: which we will use in the next calculation.
The second derivatives can be calculated by: In the last formula (55), we are especially interested in the case i = m for the further calculation where Now we are ready to transform the differential equation (1): or, by replacing ρ → x 1 and z → x 2 : Hence, we find the transformed differential equation: Without use of the summation convention, the transformed differential equation reads: This could be further simplified, since β 1 2 = 0 and

Transformation of the boundary conditions
Now we have to transform the boundary conditions.We will start with the dot product of the outer unit normal vector and the gradient of temperature: The transformation laws for the vectors are given by Eq. (40) and Eq.(44): where according to Eq. ( 41) the components of the normal in the q-system are: i.e., Thus, from Eq. ( 10), we get at the optical surfaces: A charming property of tensor calculus and the dot product is that one can change directly from Eq. ( 10) to Eq. (66).
First, we consider the boundary no.2: The vector of location of the upper surface is given by: Thus we find a secant resp.a tangent vector of the upper surface: In cartesian coordinates, the upper unit normal (which directs to the positive zresp.x 2 -direction for ρ → 0 i.e. it is an outer normal of the lens at this boundary) of the surface is therefore: Hence the cartesian components of the normalized upper normal are: , and Incidentally, we note that as defined in Eq. (26), is collinear to the tangent vector (68) of the upper surface q 2 = L 2 , while as defined in Eq. (30), is collinear to the the unit normal vector if one takes x 2 = H t (x 1 ) for the upper surface.
From Eq. (65), we have Since the argument x of the β's has to be evaluated at the fixed upper surface x 2 = H t (x 1 ), it will later be sufficient to consider N 1;2 t (x 1 ) as functions of x 1 alone.
Together with Eq. (66), we find the boundary condition at the upper optical surface: We proceed with boundary 4: The lower normal (which directs to the negative zresp.x 2direction for ρ → 0) of the surface is in cartesian coordinates: The cartesian components of the lower normal are: , and Hence the boundary condition at the lower surface is given by: with Finally, we treat the boundary no.1: In this case we have: The cartesian normal is given by The cartesian components of the left normal are: Here, we use the normal pointing inside the volume, but this will make no difference in the following.
According to Eq. ( 53) resp.the chain rule, we have: since β 2 1 (x) vanishes for an ordinary lens at the optical axis due to ∂H/ ∂ρ ρ→0 −→ 0. This finally gives: In case one wishes to investigate an axicon, one has further to treat the last but one equation in (82).If the original problem is given by a planparallel plate, the surface of section is already a rectangle.As one can verify easily, in that case, the transformed problem (i.e. the differential equation with boundary conditions) reduces to the original one as it should.

DISCRETIZATION AND METHOD OF SOLUTION
Now we will discretize the rectangular domain.We choose M 1 + 1 points along the q 1 -axis, and M 2 + 1 points along q 2 .The distance between the grid points shall be δ in both directions.This results in the dimensions L 1;2 of the rectangle, see Figure 3: We introduced the scale factors L 1;2 in the transformation Eq. ( 18) to have a quadratical grid in the q-space for numerical convenience, while in the original geometry, we can choose the subdivision independently in ρ− resp.z−direction by taking appropriate values of M 1;2 .
The grid points are: Hence, the grid points in the original (x 1 |x 2 )-coordinate system are according to Eq. ( 20):

Discretization of the differential equation
The inner points are given by: while for the boundaries it is either I = 0 (first boundary), J = M 2 (second boundary), I = M 1 (third boundary), or J = 0 (fourth boundary).
Since there is no penury of numerical differentiation formulas, there are a number of schemes for discretization of Eq. ( 61).
The most simple way would probably be the usage of approximations like the following ones given e.g. in [28]: As we remark in passing: in this book and in further publications (e.g., [3], [15] or [26]), there are even formulas and methods for domains which are bound by curved surfaces like in our original task.But here we want to overcome problems caused by individual lens shapes by our transformation to a rectangle.Now, for the discretization, we use a method discussed in [15] which is a little more advanced; it seems to be a good compromise between flexibility, stability, accuracy, pace and ease of implementation, while there are still even more sophisticated algorithms.
In [3] we find the following forward difference formula: With Eq. ( 83) this gives the discrete version of the boundary condition at the boundary no.1: In the q 1 -direction, we use the discrete approximation of the first derivative given in [1]: Therefore, we need a left and a right neighbour point, i.e., this formula may only be used for 0 < I < M 1 .

Solution of the discretized problem
Now we have for the boundary points Eqs. ( 108), ( 113), ( 115), (118), while for inner points Eq. ( 106) holds.This is a system of linear equations for the unknown T(I; J)'s.Since every T(I; J) only depends directly on values of the temperature of a few points which are spatially near, this system is said to be sparse, and a direct solution would at least demand for a lot of memory and the usage of subroutines for doing numerical linear algebra.To avoid this, it is customary to solve such systems by iterative methods; here, one can make use of the relaxation method (SOR), which may even be parallelized.For further details, see e.g.[3], [15], or [45].Faster schemes based on numerical linear algebra routines can be found in e.g.[8].

An approximation of the solution
In the examples, we want to compare the results of our our algorithm with an approximative method, taylored for thin lenses, which will be discussed now.

Model for the generation of heat due to absorption
As we are interested in the heat produced in the lens by the absorption of light, we have to speak about optical radiation issues.For the purposes of this article, i.e. the description of a numerical algorithm for the distribution of temperature, a simplified discussion is sufficient.More realistic simulations would demand for more sophisticated methods for the determination of the light field.In case of need, readers are referred to e.g.[43] for the basics, as well as e.g.[9], [13], [19], [27], [36], [39], [40], [49], [54], [58], or the pertinent literature given in [41] for a selection of methods or examples, whereas the problems dealt with here are are still issues in active academical and industrial research [5], [14].It is assumed, that for the simulation of an optical system, this problem may be solved by ray-tracing or beam-propagation software.
Light transports energy, and so the basic quantity is the flux Φ [W] of a light source or a beam.For monochromatic light of frequency f 1 s , one can write where Ṅ 1 s is the number of photons per second and h = 6.6 • 10 −34 Ws 2 is Planck's constant.The irradiance I W m 2 is the flux per surface of section A of the beam and is given by [46]: If absorption or scattering in a medium takes place, the number of photons lossed is proportional to the acute number of photons if the events are statistically independent, i.e., if the intensity is not too high (for the description of further effects, including solarization, see e.g.[6] or [38]).So, it is dΦ ∝ Φ, and thus we find the law of Lambert-Beer: where ds [m] is an element of the arc length along the path of the photons, while α 1 m is a parameter which describes the weakening of the beam.Since we will not consider the details of propagation, we will make a rather crude model for the distribution of irradiance in the medium here, where the light is assumed to propagate in positive z-direction: I(ρ, 0) may be given by a Gaussian or a rectangular distribution for example, while our algorithm is able to treat any distribution of the sources of heat ν(ρ, z) in principle.
The loss of power in the beam during propagation in the lens is where the first term is caused by scattering and the second one by absorption.The absorbed power is assumed to change into heat, thus the source-density in the heat equation is given by: where the approximation in the last step is valid for which is met normally very well for the situations in lens systems.For the implentation of our algorithm, though, we do not need that approximation.
Due to the lack of data, we will neglect the effect of scattering in our examples (α sc = 0) and will determine the values of α = α abs from the transmission of the materials, while for e.g.fused silica there are data available and scattering might be included [53].

An effective heat equation
In [41], we described an approximation called 'the modified heat equation' for the radial distribution of temperature in a thin lens which we will use here for comparison.The incorporation of the transfer of heat to the environment was discussed for the simple case of a lens with one flat optical surface, but this can easily be generalized, compare [18]: The bar over the quantities shall symbolize that they belong to the approximation, and in some sense are mean values of the corresponding 3D quantities along the z-direction.So the source density is given by: Technically, the main advantage of Eq. ( 128) is that it is an ordinary differential equation for T(ρ), which is much easier to solve numerically than the partial differential Eq. ( 1).The meaning of each of its terms can be traced back to the elementary physical effects.It seems to be typical for this approximation that it under-estimates the temperature of the lens a little bit due to the neglection of the curvature of T in the zdirection.
For the comparison with the approximation, we calculated the mean values of the temperatures in z-direction of the 3Dsolution for every value of I.

Examples
Now, we will look at a few examples and we will compare the results of the method deveolped here with the approximation.
Since our primary concern here is the demonstration of the numerical algorithm, we will show the raw data of calculation in the 3D plots.
6.1 Comparison of a lens from N-BK7 resp.

N-BK7HT
As mentioned in the introduction, there are new glasses on the market with enhanced transmission.Now we want to compare the effects on the temperature for one such glass, N-BK7HT.
We consider a biconvex lens with the following parameters: Upper radius of curvature: R t = −300 mm, lower radius of curvature: R b = 2500 mm, mean thickness of the lens: D = 15 mm, half of the diameter: R = 80 mm.The heat conductivity is given by λ = 1.114W m•K for both, N-BK7 and N-BK7HT [51], [52].The illumination is given by a Gaussian distribution Eq. ( 123) with I 0 = 1.0 • 10 5 W m 2 and b = 50 mm.The temperatures at the edge of the lens and that of the ambient surroundings are assumed to be 20 • C; for the Newtonian heat-transfer-coefficient we guess α N = 10 W m 2 •K .
The results are shown in Figure 4 and Figure 5.
From these calculations we see that the rise of temperature in the center of the lens relative to its edge is proportional to the value of the exponent in the Lambert-Beer-law, α, under otherwise similar circumstances.This agrees with our earlier approximation [41], although the latter was without heat exchange to the environment.
For N-BK7, we used τ i = 0.997, as was given by Schott.So, the actual value of τ i can be assumed to be between 0.9965 and 0.9974, and the actual value of α should be between 0.26 1 m and 0.35 1  m .Therefore, we have an uncertainty of 30% in total in α and thus in the the value of the temperature in the middle of the lens only due to that harmless looking rounding of the transmission data, because it is not τ i but 1 − τ i which causes the effect of heating.

Realistic lens form
Now we want to give an example for a realistic lens form with a strong deviation from a spherical surface, in our case a kink near the edge of the lens.In Figure 6 the form of our example is shown, while other typical lens forms could be treated as well.We see, that our algorithm is even able to treat such a form of a lens with kink.However, a problem may arise if one of the grid points lies directly on the kink, because e.g. the derivative in Eq. ( 50) is no more defined in these cases.Here, one has slightly to change the value of M 1 to shift the grid points.Visually, the curves of the approximation and the average of the numerical solution are almost identical in this case.

Lack of experimental data
We extracted the density of heat sources ν from transmission data.But in reality, there may be scattering in the volume [53] or scattering resp.absorption in the coatings.Surface effects may cause damage and are considered e.g. in DIN EN ISO 21254, [16], [22], [25], or [57].Reliable, systematical and standardized data of absorption in the coatings seem to be rare although there are a number of works dealing with methods of measurement.One estimate from literature is an absorption of 0.12% of the power of the beam for 1064 nm per surface of a coated quartz lens [35]; an other is the assumption, that the absorption in the coating approximately equals the absorption in 10 mm bulk of quartz glass [12].This indicates that the contributions of the coatings to the absorption and thus for the heating of the lens may be of great importance.
In the examples, we pointed out that the rise of temperature in the center of the lens depends sensitively on the value of transmission given for a material and a wavelength.This is not a problem of our algorithm, which is able to calculate the distribution of temperature from the density of heat sources, but is due to the knowledge of that density.
The validity of the Lambert-Beer law is limited to linear effects, i.e. it will fail for very high intensities.Those points were investigated and data published for some materials (mainly those of interest in microlithography), but a more systematically base of data for optical glasses seems to be missing.The same is true for data like the coefficient in Newton's law or similar boundary conditions.Therefore we neglected the effects here or tried estimated values.This problem, however, does not concern our method of solution of the mathematical problem, at least not in principle: the results of any software will suffer by barely known input parameters.
Our aim here was to present a way which allows the calculation of the distribution of temperature.Perhaps, it may be the base for more advanced schemes.In any case, the development of the algorithm made it clearer, at which places information from assumptions or models occur in the calculation.Thus, comparisons with measured data are highly desirable.Hopefully, examples and methods like the one given here may be a help for a better experimental determination of such quantities.

Generalizations
The columns x 1 = const.are equidistant in the present algorithm, but the transformations (18) are not the most general ones.With a more general transformation of the coordinates, this may be overcome to achieve a very high resolution e.g.near the optical axis or the edge of the lens for the description FIG.7 Lens with kink, raw data and comparison with the approximation.For the parameters see Figure 6 and the text.Notice the different length-scales in ρ-resp.z-direction in the 3D plot.
of effects located there.One possible generalization may be with f (0) = 0 and f (R) = R.However, one has to pay attention to the invertibility of the transformation.
A slightly other example may be the choice If necessary, one can construct a non-uniform grid in x 1direction with this.
By a similar treatment of the transformation along the x 2direction, a streching near the optical boundaries may be done for a better possibility for the modelling of surface effects.
As we have seen, absorption in the coatings or at the surfaces is an important issue for realistic calculations.Of course, that could be taken into consideration in the calculation.A more explicit modelling of the coatings may eventually be possible by the aid of a more sophisticated transformation, but requires a better knowledge of the relevant properties of the materials involved.
For simplicity, and for the comparison with our earlier approximation, we took the same boundary condition for all the points at the boundary no.3 (x 1 = R), i.e. at the edge of the lens.In general this is not necessary, of course.
The present state of our algorithm does not take into account a heat conductivity depending on temperature.The perhaps simplest way to overcome this would be the definition of λ = λ(T) and the determination of that value for the grid points with T(I; J) during each step of iteration, while the classical theory demands for a discretization of the equation [34]: see [45].
In case the distribution of the sources of heat is not radially symmetric, separation of the variables may be helpful.Then our presented theory is part of a more general calculation.
For simplicity, we used the same temperature of the ambient surroundings T s on both sides of the lens.This might not be true since, e.g., one surface of the single lens is at an outer side of the objective, while the other is at an inner one.It is, of course, not difficult to generalize this point to different temperatures on both sides.
A method for the calculation of the distribution of temperature due to conduction and convection is discussed in [10] and may be of interest for the treatment of the air spaces between lenses.
Simple models for the change of the shape of the lens due to the distribution of temperature are given in [55] or [18].Such methods are well suited for our present approach.A more accurate way would be the solution of the equations of thermoelasticity [31].

CONCLUSION
We developed an algorithm for the calculation of the distribution of temperature in a single lens due to absorption of light and heat conduction for the case of a stationary situation with symmetry around the optical axis.First, we remarked that it is sufficient to choose a surface of section of the lens as computational domain, then we transformed this area to a rectangle.Therefore, the differential equation and the boundary conditions had to be transformed, too.This transformed problem is well suited for a treatment by finite differences, regardless of the original shape of the lens.Thus, the task was transferred to a system of algebraic equations.
A simple model for the description of the sources of heat due to absorption of light based on the known data of internal transmission was discussed, and we considered a few examples.We learned, that even a small deviation in the transmission data may have a conspicuous effect in the distribution of temperature.An other example demonstrated the ability of the method to adjust to special but realistic lens forms.Finally, we pointed out suggestions for experimental investigations in the future and gave prospects for further generalizations of this theory.

FIG. 2
FIG.2Transformation of the surface of section of the lens to a rectangle.

FIG. 4 N
FIG. 4 N-BK7, raw data and comparison with the approximation.For the parameters see the text.Notice the different length-scales in ρ-resp.z-direction in the 3D plot.Visually, the curves of the approximation and the average of the numerical solution are almost identical in this case.
Dfor ρ ≤ 60 mm const.for60 mm < ρ ≤ 80 mm .(136) R t = 150 mm is the radius of curvature of the upper surface, while D = 5 mm is the thickness of the lens.The material of the lens shall be N-BK7, so we have λ = 1.114W m•K and α N-BK7 = 0.30 1 m as in the example above.Likewise, the temperatures at the edge of the lens (boundary no. 3) and that of the ambient surroundings are assumed to be 20 • C; for the Newtonian heat-transfer-coefficient we guess once more α N = 10 W m 2 •K .The illumination is given by Eq. (124) with I 0 = 1.0 • 10 5 W m 2 and b = 40 mm.The results are shown in Figure 7.

FIG. 5 N
FIG. 5 N-BK7HT, raw data and comparison with the approximation.For the parameters see the text.Notice the different length-scales in ρ-resp.z-direction in the 3D plot.

FIG. 6
FIG.6Lens with kink at the left surface.