A simple algorithm for multiplereflection-refraction of rays from smooth surface

P.Bruscaglioni, A.Ismaelli. Department of Physics. University of Florence, Italy

1.Introduction.

This document describes a simple algorithm for tracing trajectories of rays (e.m. or optical radiation) incident and reflected on a smooth surface, for which geometrical optics can be used for reflection and refraction.

The surface is determined by a series of points whose co-ordinates are with reference to a system X, Y, Z..

The points x, y , are at regular intervals with distance

xm-1 -. xm = Hx, and yn-1 - yn- = Hy

In the rest of the document the case Hx = Hy = H will be considered,. However, an extension to the general case can be easily made.

Moreover, H will be considered as unity for the distances, and will be taken as H = 1.

The ray is assumed to start at a position P0 (x0, y0, z0) with director cosinuses Sx, Sy, Szwith respect to the axes X, Y, Z., with Sz < 0

z0 is assumed to be equal to the maximum value of the z co-ordinate of the surface.

Initially the point P0 is found to be placed in a square a,b,c,d on a plane x y, as is shown by Fig.1.

The vertex of the square have the co-ordinates

1 )

a : xa = m-1, ya = n, za ( they stand for (m-1)H, nH, (za)H

b : xb = m, yb = n, zb

c : xc = m-1, yc = n-1, zc

d : xd = m, yd = n-1, zd

Given Po (x0,y0,z0), one has :

2) m = Integer (x0/H) + 1 , n = Integer (y0/H) + 1

As shown in Fig.1, for the sequel, the squares will be divided in two parts : I and II. by the diagonal c b. The two parts of the square are assumed to be plane. Thus the whole surface is represented as an ensemble of plane triangles.

It can be noted that the division of the square in two parts by the diagonal a d would result in a different approximation of the real surface. That is two slightly different surfaces could be obtained when H is small with respect to some characteristic distance where the surface changes substantially the z co-ordinate.

2 Determination of the point of incidence.

From the point P0 the ray is made to move by a distance L along its direction. to the point where it may incide on the plane of the triangle I of the square of Fig.1.

Three cases can occur.

A The point T of incidence cannot be found. The ray trajectory is parallel to the plane,

Or L is negative : the trajectory forms an angle smaller than  / 2. with the normal to the plane directed toward the external medium.

B L is positive, but the point T is found, but it is out of the triangle.

C. The point T is found on the triangle

Let the director cosines of the normal n to the surface( directed toward the external medium) be U1, U2, U3

In the cases B, C, the distance L is

3) L = ( U1 (xa-x0) + U2 ( ya-y0) + U3 (za-z0) ) / (U1 Sx + U2 Sy +U3 Sz)

For the case A Eq.1 shows L = . The ray trajectory is perpendicular to the surface normal.

If the cases A, B are found, the exam is repeated for the triangle II.

We can again have three cases A, B, C, with the distance L is given by Eq.3 with the director cosinuses U1,U2, U3 of the plane b c d

Let us assume that the case C occur, for the triangle I or the triangle II

Let the point of intersection be with the co-ordinates XT, YT, ZT

4)

XT = x0 + L Sx
YT = y0 + L Sy

ZT = z0 + L Sz

The cosinus directors of the normal to the surfaces I and II, equal to

For the triangle I

5a)

U1 = (za-zb) / R U2 = (zc - za) / R U3 = 1 / R

With R = (1 + (za-zb)2 + (zc - za)2 ) 1/2

And for the triangle II.

5b)

U1 = (zc - zd) / R U2 = (zd - zb) / R U3 = 1 / R

and R = = (1 + (zc-zd)2 + (zd - zb)2 ) 1/2

3 One of the triangles of Fig. 1 is interceped.

Let us assume that the triangle I is intercepted at Point T.

a) There is reflection and the new trajectory start at point T whose co-ordinates become the new xo, yo, zo, with new cosinuses Sx, Sy, Sz, where now Sz can be positive.

The trajectory is continued to decide whether the triangle II can be intercepted, following the procedure of Sect.2. Due to the different orientation of the triangles, the consecutive encounters with the two are possible.

Analogously, if the triangle I, is not intercepted, and the triangle II is, the procedure a) is followed for the couple II  I.

4. The square of Fig.1 is not intercepted.

For both the triangles I and II, case C of Sect.2 does not occur.

The trajectory is then moved forward in its direction, and it comes over a next square, which is determined by the trajectory direction as follows.

Fig.2 is used where the central square A is that of FIG.1.

From the square A the trajectory can proceed to be over the squares B,C,D,E

With no reflection on the square A, the co-ordinates of the point P0 where the trajectory starts are still are xo, y0, z0.. The cosinus directors of the trajectory with respect to the X Y Z axes are still Sx, Sy, Sz.

If reflection has occurred, xo, y0, z0. become the coordinates of T, and Sx, Sy, Sz become the director cosinuses of the reflected ray.

A is the square where reflection occurs. The co-ordinates of its corners are those indicated by Eq.1

From x0, y0.z0 the trajectory is made to proceed till a point over one of the squares B,C,D,E is reached.

The decision for which of them is taken as follows, based on simple geometrical consideration.

Case 1, Sx, Sy > 0.

Let us put

Ly = (n- y0) / Sy L'y = (n+1 -y0) / Sy

Lx = (m- -x0) / Sx L'x = (m+1 - x0) / Sx

If Lx > Ly Then the distance to be run is L:

L = 1/2 (Ly + Min (L'x, Lx))

And the trajectory is made to proceed to a point of square B

If Lx< Ly, Then :

Lp = 1/2 ( Lx + Min (L'x,Ly)

And the trajectory is made to proceed to a point of square C

Case Sx > 0, Sy < 0.

Ly = (n-1 - y0) / Sy L'y = (n-2-y0) / Sy

Lx = (m -x0) / Sx L'x = (m+1 - xo) / Sx

If Lx > Ly Then the distance to be run is L:

L= 1/2 (Ly + Min (L'x, Lx))

And the trajectory is made to proceed to a point of square C

If Lx< Ly, Then :

L = 1/2 ( Lx + Min (L'x,Ly) / 2

And the trajectory is made to proceed to a point of square E

The two cases above, together with the other two cases of the signs of sSx, and Sy, can be summarized by the following formulas :

By defining :

Lx = (m -x + (sign ( Sx) -1) / 2) / Sx

L'x =Lx + 1 / Abs(Sx)

Ly = (n - y + (sign(Sy) -1)/2)( / Sy

L'y = Ly + 1 / Abs(Sy)

Distance L to be run :

If Lx< Ly L = Lx + Min (L'x, Ly)) /2

If Lx > Ly L = (Ly + Min (L'y, Lx) /2.

The co-ordinates of the final point are

x' = x0 + L Sx, y' = y0 + L Sy, z' = z0+ L Sz

Then, one of the squares B,C,D,E of Fig 2 is reached.and it is determined by the new coordinates of the corners, which can be again named as a,b,c,d:

a = Integer (x'), Integer (y' +1)

b = Integer (x'+1), Integer (y' +1)

c =Integer (x'), Integer (y')

d =Integer (x' +1) Integer (y')

The square is divided in two trangles and the intersection point is searched for.

If no intersection is found. The trajectory is continued to the next squares , repeating the procedure of Sect.2, 3,4

The point obtained by the procedure of Sect 4. belongs certainly to a square next to that of reflection. A of Fig.1. Let this square be named Ad and the end point after the L distance be Pd.

The direction and polarization of the reflected ray are given by geometrical optics

One has to examine if and where the ray actually encounters the square Ad.

For this, as was shown in Sect. 2 the square is divided in the two triangles.

Since the point Pd may be beyond one of triangle planes, a negative distance can be accepted for the encounter point.

If the encounter point is not found the ray is made to proceed to a next square, starting from Pd and searching for a next square.

Reflection at the incident point

The direction vector U, normal to the plane, for the triangles I and II are given as follows, with reference to Fig.1 and positions 1

Triangle I

UI = ( k - (zb-za) i - (za-zc) j ) / ( 1 + (zb-za)2- (za-zc) 2)1/2

Triangle II

UII = ( k - (zc-zd) i - (zd-zb) j ) / ( 1 + (zc-zd)2- (zd-zb) 2)1/2

The z components are positive (toward the external medium)

Given the direction vector Si of the incident radiation :

Si = Sx x + Sy y + Sz z,

one has for the reflection direction Sr :

Sr = Si - 2 (Si . U)U,

with U according to which triangle is encountered.

Polarization of the reflected ray.

To apply the Fresnel coefficients, electric vector of incident radiation must be divided in the two components parallel and normal to the reflection plane (plane of Si and Sr)

Let the incident electric vector Ei have the components Ex, Ey, Ez., and the components of the normal vector un to the surface at the reflection point be U1, U2, U3

The direction vector urn normal to the reflection plane is

urn = si un / si  un

The component of E normal to the reflection plane is

Ein = Eiurn and the component on the plane is

Eip = Ei - Ein.

To Ein the reflection coefficient Rte has to be applied to obtain the component of the reflected E normal to the plane.

To Eip the Rpe coefficient has to be applied.

Whre Rte and Rpe are the Fresnel coefficients.

The components of the reflected field Er are referred to the X Y Z system, and Er has no component along the propagation direction.

One has to take into acount that Rte and Rpe are generally complex numbers.

FIG.1


FIG.2