Compur. Bml. Med.
Pergamon
Press
1977. Vol. 7. pp. 35-70.
RAY-TRACING
Prmted
in Great
Britain.
FOR PHYSIOLOGICAL A TUTORIAL
OPTICS-
WAYNE WIITANEN Department of Biology, University of Oregon, Eugene, OR 97403, U.S.A. (Received 24 October 1975)
Abstract-This tutorial article will introduce the physiologist to the basic techniques of ray-tracing: mathematical formulation, and associated computer program. The techniques are illustrated on schematic eyes of the cat, man, the blowfly, and the honeybee. Instructions for use of the program and a listing of the program are included in the article. Physiological optics
FORTRAN
Tutorial
Eyes
Ray-tracing
INTRODUCTION Perhaps one of the more significant contributions that modern digital computation can make to physiological optics is the mechanization of ray-tracing methods. Because of the high degree of precision involved in charting ray paths by this method it seems that such methods should be relevant to physiological optics. In the past, the application of ray-tracing has been limited because of the difficulty in using the method; but, this objection no longer holds. Experience with ray-tracing used in connection with an investigation of the optics of the honey bee’s eye suggests that the method should be of value to those interested in vision and physiological optics. The general plan of the paper is tutorial. It assumes an acquaintance on the part of the reader with elementary vector analysis at a level comparable to that found in any standard introductory text, such as Lass [l]. The mathematical exposition develops in detail a small part of a more elaborate paper on general ray-tracing methods written by Spencer [2]. Unfortunately this paper is not widely known. A basic knowledge of FORTRAN IV format statements will make understanding the card formats described in the sections concerned with using the ray-tracing program much easier. The FORTRAN IV program itself has been completely rewritten from the IBM version (written by D. H. O’Herren of Goodyear Aerospace), with the exception of the iterative procedure for following rays through refracting surfaces (which is based on the work of Spencer mentioned above). There are, then, three goals for this paper: (1) An exposition of ray-tracing methods, (2) instruction in the use of the computer program, (3) to expose visual physiologists to an available technique at a level of detail sufficient for subsequent use. 1. GEOMETRICAL
RAY-TRACING
Fundamental to the understanding of ray-tracing is Snell’s Law, which describes how a ray of light is bent as it passes from one medium into another. Either from experimental work or from the solution of Maxwell’s equations, subject to the proper boundary conditions, Snell’s Law may be found as n, sin 8, = n2 sin 8,, where n, is the index of refraction of the first medium, nz is the index of refraction of the second medium (that into which light is passing), 8, is the angle that the incident ray (in the medium described by n,) makes with the outward normal to the surface separating the two media, and e2 is the angle that the refracted ray makes with the inward normal to the same surface. The incident ray, the refracted ray, and the normal to the refracting surface all lie in the same plane, i.e. are coplanar.
WAYNE WIITANEN
36
A vector formulation of Snell’s Law is desirable for ray tracing. Let u1 be a unit vector representing the incident ray and u2 be a unit vector representing the refracted ray. The components of these vectors are the directional cosines of the rays. Denote the unit vector normal to the surface by n and let it point into the medium into which light is refracted (the inward normal). Because ur, uz, and n all lie in the same plane advantage can be taken of the co-planar relationship between the vectors, and write u2 = ati1 + bn, where a and b are numbers yet to be determined. of this equation with the normal vector, n, giving u2An
To this end take the cross product
= (au1 + bn)/jn = a(u, An) +
b(n/\n)
= a(uI An). Recalling that,lu,l = lull = InI = 1 and that Ip/jq( = IpI.lql sin(p,q), it is clear that the equation above says that sin (u,, n) = a sin (ul, n) or
a sin 8, = sin 19~,
since e1 is the angle between u1 and n, and Q2 the angle between u2 and n. By comparing this last equation with Snell’s Law it follows that a = n1/n2. To determine b another relationship between the vectors ul, u2, and n is needed. One such relationship is that they are all unit vectors (that is, jull = ju2( = InI = 1). Forming the dot product of u2 = au1 + bn with itself we have that 1 = u2*u2 = (au1 + bn)*(aul + bn) = a2uI.u1 + 2abu,*n + b2n.n = a2 + 2abu,.n
+ b’.
This quadratic equation is easily solved for b, b = --au1 .n f [u2(ul *n)2 - u2 + l]“‘. The sign of the radical is yet to be determined. This is most easily done by considering the case where there is no refraction (nl = nz): a = 1, and u1 = u2, Under these conditions b must be equal to zero, implying that b = -au,*n
+ [u2(uI.n)2]“2.
Finally, the vector form of Snell’s Law becomes: u2 =
nlul/n2 - [nluI*n/n2 - (n~(u,*n)2/n~ - n:/n$ + 1)“2]n.
When a > 1 and u,*n is small it is possible to have u2(u,.n)2 - u2 + 1 < 0, corresponding to the condition for total internal reflection. A more symmetrical form of Snell’s Law may be obtained by taking the cross product of n with u2 = (nl/n2)u1 + bn which is n2(u2 A 4
=
nlh
An).
At each refracting surface the incident electromagnetic radiation of a light wave is split into two components: a reflected wave and a transmitted wave. From the Fresnel formulae [3] the amplitudes of the transmitted and reflected waves may be calculated. However, this cakzulation requires that the angIes of incidence and refraction be known. Let 8i designate the angle that the incident ray makes with the normal, and that which
Physiologicaloptics the transmitted we have that
37
ray makes with the normal as et. Since ur represents the incident ray COSBi = u,-n
or 8i = arccos
(II1an).
u2 = au, +
bn.
For the transmitted ray it follows that
Taking the dot product ofboth
sides of the equation with u2 we have 1 = a(ur *u2) + = u(u1~u2) +
&n-u,)
bcos (n,u,),
or 19~= arccos [l - a(uI .u2)]/b. By considering how the energy of the electromagnetic field of incident light energy is distributed between the reflected and transmitted fields for a given angle of polarization it is possible to define the reflectivity and transmissivity of the refracting surface by means of the equations R = RI, cos2 CI+ RL sin2 LX (a = angle of polarization)
where R 11= tan’ (Oi - @,)/tan’ (Bi + 0,) RI = sin2 (ei - @/sin’ (ei + 6,)
and T = T,, cos2 a + Tl sin’ CI (a = angle of polarization)
where T,, = [sin 2&.sin 28,]/[sin2 (Oi + t9,).cos2 (6$ - O,)] TL = sin 281. sin 26,)/sin2 (Oi + 0,).
Natural light is unpolarized and the corresponding reflectivity, R, and transmissivity, T, may be defined by averaging over all angles of polarization, a. Since the average of sin* a or cos’ a is 3 it follows that R = &(R, + RII) T = +(T, + T,,)
and R+T=l. To obtain the corresponding energies from the reflectivity and transmissivity it is necessary to know the energy of the incident light. Call this energy JCi);then we have that J(r) = RJ(i)
where J(*) is the average energy of the reflected ray, and J(t) the average energy of the transmitted wave. From these formulae it is possible to determine the energy of the light rays being transmitted through successive refracting surfaces by setting jCi) equal to Jft) after the refracted ray leaves the refracting surface. The intensity of the reflected and transmitted rays may be calculated from the formulae S”’ = j”‘/cos
(g.
S(i) = JWpos
e. = Rs(i)
I
1
s(t) = J(ycos 8, =
m(i)
cos BJCOS8 t
WAYNE WIITANEN
38
where $I’ and 3”’ are the average respectively. 2. SPECIFICATION TRACING
intensities of the reflected and transmitted
OF REFRACTING RAYS THROUGH
SURFACES THEM
rays,
AND
To specify the rays and the geometry of the refracting surfaces two coordinate systems will be used, both rectangular. The first system will be called a global (or reference) system and it is in this system that the rays originate and their paths are traced. Components of a vector specified in this system, or of a point in the system, will be indicated by having bars over them (this notation and presentation follows that of Spencer [2] so that the mathematical development can be related directly to the program listing in the Appendix). The second coordinate system is a local one that is used to define the shape of the refracting surfaces. Eventually this local surface will have to be positioned within the global one so that the rays may be traced through it. For the moment it is easier to define surface equations in a system that is independent of the global one. Vectors and points in the local coordinate system (local to the surface being defined) will not have bars over them. A given ray is specified by the coordinates (x0, YO,2,) of a point, say P,. through which it passes and the direction cosines (k, I,%) of that ray. Remember that the bars over a letter indicate that the ray is in the global coordinate system. The origin of the global system is denoted by 0. A surface, S, in the global system is defined by the equation F(X, Y Z) = 0. Because the surface is usually defined in a local coordinate
system by the equation
F(X, Y, 2) = 0, it is necessary to relate the local coordinate system to the global one. The steps involved in the ray-tracing problem may, then, be separated as follows: 1. Transform the ray origin (x0, Y,, Z,) and directional cosines (k, 1,&) to their corresponding values in the local coordinate system (X, Y, Z). 2. Find the point of intersection of the ray with the surface, S, as defined by the equation F(X, Y,Z) = 0.
(B)
Fig. 1. Euler angles.
A-rotation about Y-axis (Euler angle CL).B--rotation angle p). C-rotation about Z-axis (Euler angle 7).
about
X-axis
(Euler
Physiological optics
39
3. Find the change in direction of the ray in the local coordinate system following refraction or reflection at the surface, S. 4. Transform the new ray origin (xb, yb, zb) (point of intersection of the incident ray with the refracting -surface) - and the direction cosines (k’, I’, m’) back into the global coordinate system (X, Y, Z). Steps l-4 are repeated as many times as there are refracting surfaces in the optical system. Each of the steps will now be examined in detail. Step 1. Transformation of rays to the local coordinate system of S. The orientation -- of the - local coordinate system (X, Y, Z) is specified with respect to the global one (X, Y, Z) by means of the Euler angles as shown in Fig. 1. The Euler angles specify the rotation of the local system with respect to the global one by means of three successive rotations: CC,p, and y. a is a rotation about the Y-axis, p is a rotation about the X-axis, and y is a rotation about the Z-axis-all in the counter-clockwise sense. Each rotation may be specified by a 3 x 3 matrix as follows
(sin y about Z-axis)
R, = [ i
;i;
-$ygl
(about x-axis)
(about Y-axis)
Rotation through the three angles c(, /I, and y is the product of the three matrices, R,,, = R,R,R,, and in terms of components this is coscl.cosy
R rPa=
+ sincc.sinfi.siny
coscc.siny - sincr.sinp.cosy [
sin cIcosp
-cos/?*cosy cosp~cosy sin p
-sincr.cosy
+ cosa.sinp*siny
-sincl.siny
- cosa.sinfl-cosy cos c(cos p
1 .
Because matrix multiplication is, in general, non-commutative, it is important to note that when defining the orientation of a local coordinate system with respect to the global one by means of the Euler angles that the local system be positioned first with respect to a, then with respect to /?, and finally with respect to y. An important property possessed by each of the matrices R,,R,, and R, is that they are unitary, that is, that the product of the matrix with its transpose yields the identity matrix. This property is easily verified for the three R matrices. The product of two unit matrices is again a unit matrix. All of this means that the inverse of the unitary matrix R is the transpose of that matrix. This property will be used --- just a little later. To translate the ray origin in the (X, Y, Z), or global, system to the local system in which the surface equation has been formulated, it is necessary to know the position of the origin of the local system (X, Y,Z) with respect to the global one (the first card of group 2 does this for the computer program). If the coordinates of 0 (origin
WAYNE WIITANEN
40
of the local system) with respect to the global system are (x0, To, ZO), then the ray origin in the local system will be (X,, Y,, 2,) according to the transformation equation
where (X0, Y,, Z,) is the ray origin in the global system, and the direction cosines become
il [I R
k
1 = R,,,
i
n?
111
Step 2. Point of intersection of the ray withthe surface S. Recall from vector analysis that the vector equation of a line is given by v = p + au. Using this formula, a given ray may be expressed in parametric form as u1 = PO + su2 or, in terms of components, X = X,, + sk Y=
Y,+sl
z = zo + sm.
where s is a parameter of distance along the ray measured from the transformed ray origin PO(XO, YO,Z,). The next problem is to determine the point of intersection of the line representing the incident ray with the refracting surface, S. The scheme employed is to determine the intersection point of the ray with the tangent plane Z = 0, and then to iterate from that point down to the surface S, thereby obtaining the intersection point of the ray with the surface, S. Eventually, an iterative procedure is found allowing determination of the value of the parameter s that defines the point of intersection of the ray with the refracting surface, S. To start out, the parametric equations of the ray are solved for s when Z = 0. This equation is s0 = -Z,/m. By substituting the value so into the parametric equations for the ray, the coordinates of the point of intersection of the ray with the plane Z = 0 are obtained: X1 =X,,+kso Y, = Yo + Is, z,
= 0.
Redefine s so that it is measured from the point (X,, YI, 0), and rewrite the parametric equations for the ray as X = X, + sk Y = Y, + sl Z = sm.
Physiological
41
optics
Next find a value of s such that the resulting values (X, Y,2) from the above equation will satisfy the surface equation F(X, Y, Z) = 0 (describing the refracting surface, S). Solution by purely algebraic means is not practical because of the complexity of the surface equations that may be used. Instead, a Newton-Raphson technique is employed. The Newton-Raphson process is an iterative method for locating the roots of the equation f(-u) = 0 [4]. Suppose that an estimate of a root of the equation is known. Call the estimate a. A correction to the estimate, h, is desired, which when applied to a, will yield a better estimate of the root: ni + 1 = Ui + hi, where hi is the i-th correction to the i-th estimate, ai, of the true root. a. The equation them becomes f(a + h) = 0. To obtain a formula that will give successive correction terms, hi, the last formula is expanded by means of a Taylor’s series about h: f(a
Becausef(a
+ h) -,f(u)
+ hi’(u) + (h2/2)f”(u
0 < 8 < 1.
+ dh),
+ h) = 0 it follows that f(u) + hf’(u) + (h2/2)f“(u
+ #h) A 0.
If h is small the term containing h2 may be neglected, which gives the simple formula f(a) + 1$‘(a) = 0, or h = -f(u)/j-‘(a). Casting this formula into an iterative form for “corrected”
roots gives
Ui+ 1 = Ui + hi = Ui -f(Ui)/f'(Ui),
This formula indicates that the derivative of the equation must exist in order for solutions to exist in the neighborhood of the true root, a. It is reasonable to ask about the error that is introduced into the process by neglecting the h2 term in the complete expansion. This problem may be investigated by subtracting the abbreviated equation from the complete one, which gives h - hi = -h2f”(ui
+ 812)/2f’(ui),
0 < e < 1.
Because h is the true correction, and hi the approximation, it is clear that h - hi is the error in the approximation to h. For convenience the i subscripts will be dropped and F will be written for the maximum value of f”(u) in the neighborhood of an hi. This gives a somewhat simpler appearing formula: E(lz) = 11- hi = -h2F/2f’(u). Clearing of fractions, transposing, and solving the resulting quadratic for h the equation above yields h = l/F[ -f’(u)
+f’(u)(l
+ 2Fhilf’(u))“2].
Expansion of the last term by means of the binomial theorem gives !Z= hi - Fhf/2f’(u)
+ F2hf/2(f’(a))2
+ ...
as the formula describing the error involved in the estimation of h when the Taylor’s series is truncated after the second term. Again. because hi is usually small, the second term (third degree term) in the binomial expansion may be deleted at no great loss of accuracy. This leads to a general formula for the error in the i-th correction to the root: Ei(hi) I (F. h2/2f”(Ui_ 1) 1.
The next question that might be asked is: will the method converge to a solution? Essentially this means that the sequence of errors between the true root and the succes-
WAYNE WIITANEN
42
sive approximations should converge to zero or to some preassigned (convergence tolerance). Let the equation in question be written as
small number
.x =f(x) and let x0 be the first approximation
to the true root, 2. Then write Xl
as the next approximation equations yields
=f(xo)
to the true root. Forming x - x1 = f (x) -f
the difference between the two
(x0)
as the next approximation to the true root, x. By the mean value theorem of the differential calculus, the right-hand side of the above equation can be written as f(x)
- f (x0) = (x - x(J)f ‘(io),
x < 20 < x0
and the entire equation becomes x - x1 = (x - xo)f’(&) and similarly for succeeding approximations. By writing down the n equations by which the n differences x - x, are defined and multiplying them together followed by division by the common factors x - x1, x - x2,. . ., x - x, one gets x -x,
= (x - xo)f’(~o)f’(~~)...f’(2,-l).
If the maximum absolute value off’(x) is less than 1 throughout then Ix - x,1 I Ix - Xoj’F”
the interval (x,x0)
where F is the maximum value of all of thef’(x,): F = Imax(f’(20),f’(il)
,..., f’(i,_l))I
< 1.
The right-hand side of the inequality will approach zero as n becomes large, which means that the error, x - x,, may be made as small as desired by repetition of the iteration process a sufficient number of times. This means that the condition for convergence is that If’(x)1 < 1 in the neighborhood of the desired root. Applying this iterative method to the Newton-Raphson technique it is apparent that the requirement for convergence is If’(%l)I < 1.
From the definition of the Newton-Raphson f ‘(4
process it is clear that
A f (%)f “(4Mf
‘(%))Z,
and the convergence criterion becomes If(%)f
“(a”)/(f’(a”))2
I< 1
or, regrouped, If(4f
“(4) I s (f ‘(%))2.
The exact conditions for convergence and the estimate of the error will depend upon the form of the equation f(x) = 0, which will vary with the type of surface. The applicability of the Newton-Raphson method might be questioned on the grounds that it seems as if the surface equation F(X, Y, Z) = 0 is a function of three variables, while the formulae developed above are in terms of a single variable. This conflict is quite easily resolved by noting that successive values of X, Y, and Z are determined
Physiological
by the parameter ofs:
optics
s. This means that the equation
43
F(X, Y, Z) = 0 is really a function
F(X(s), Y(s), Z(s)) = 0. From this argument it is easy to see that the Newton-Raphson be written as Sj+ 1 = Sj - F(Xj, ~, Zj)/F’(Xj,
iteration formula can
rj, Zj)
where Xj=
X1 + kSj
rj = Ye + 1Sj Zj = WlSj
and the derivative F’( ) is
Because F( ) is a function of s, by way of the implicit functions X(s), Y(s), and Z(s), the derivative becomes dF ds=
i3F
dX aF dY dX’ds+E’ds+dZ’ds
i3F dZ
where -_= dX
k
-_= dY
1 ’
ds
ds
and
’
dZ ds =m.
These derivatives are obtained from the defining equations for Xj, Yj, and Zj as given above. Making the indicated substitutions yields dF
Ek
ds = x
(Xj,Yj,Zj)
+5’
+F,m (Xj,YjZj)
(Xj,Yj,Z,)
The initial estimate for s1 may be taken as s1 =
0
and the iterative process is terminated when ISf - q-11 -cc, where E is a small predetermined convergence tolerance. In general, the choice of a starting point at the Z = 0 plane will insure convergence; however, there are circumstances under which the process will either fail to converge or will converge to an incorrect value. If the derivative F’() becomes zero at any time in the process, then s will blow up and the process will diverge. In addition it is possible that a ray will strike a surface at two or more points, resulting in convergence to an incorrect value (wrong point). Neither of these circumstances is often encountered in practice. Recalling how the line describing the ray path was defined it follows that the total path of the ray from (X,, Y,, Z,) to the surface, S, is p = s0 + So. If p < 0 then the
WAYNE WHTANEN
44
path is a virtual one, extending along the line of the ray in a direction opposite to the direction of propagation. Let the intersection of the ray with the surface be indicated by the coordinates (X,, Y,, Z,), as determined by the iteration procedure. From these values it is easy to obtain the direction cosines for the surface normal at the point of intersection:
The procedure so far is applicable to any surface expressible by the equation F(X, Y, Z) = 0, providing a wide. variety of choices. By requiring that the ray-tracing procedure by implemented on a digital computer, the number of choices for the form of the surface equation becomes somewhat restricted, because symbolic derivative taking facilities (those that are done with paper and pencil by people) are not widely available on digital computers. Most of the refracting surfaces used in optics are obtained by revolving some conic section about a coordinate axis, with the vertex of the surface at the origin. The equation to be used for describing refracting surfaces has been selected to be F(X,Y,Z)=AXZ+BY’+CZ2+DX+EY+FZ-G=O. By judicious selection of coefficients a large number of surfaces may be described. Table 1 lists the surfaces most commonly encountered in practice, and indicates the nature of the coefficients of the general surface equation. Step 3. Refraction of the ray. From the previous discussion it will be recalled that Snell’s Law may be represented conveniently in vector form
where u1 is a unit vector in the direction of the incident ray with components (k, 1,m), u2 is a unit vector in the direction of the refracted ray, and n is a vector normal to the surface, S, at the ray intersection point and has components (K, L, M). n, is the refractive index of the medium in which the ray is incident, and n2 the index of Table 1. Commonly used equations for refracting surfaces. For each of the surface types named the standard equation is given, its center, and the form in which the ray-tracing program requires the coefficients
Surface
UP+
Center
equation
cz*+ Dx+
Plane
ax+Bytcz+o=o
Sphere
x2+ytz2=r2
O
0
I
I
l/b2
(+h, Zk,O)
(fh,+k&nm)
(O,O,O) (+n,?k,x,rm)
Ev+
I
I A
Fz=
G
I 8
C *2/n
-D
tL2h
iZk
(r2-(hz+k~mZ))
I/c2
t2h/a2
+_2k/b2 +2mki
I-(h’/o’+
Ilb2
0
+2h/o’
+2k/b2
(I-h2/a2-k2/b2)
l/b2
0
?2h/a2
+2k/b2
-c
/b’
I/C2f2h/a2
22X/b’
T2mk
.1/b’
0
t2hd
T2k/62
-c
I/b’
d/C2
+2h/a2
?2k/b2
i2m/c
Ellipmd x~/o~+~~/~~+~~/c~= Ellnptic
x2/a2+y2/b Elliptic
I (thI tk
9371
k 2/b%m2/c2
cylinder,Z-am..... =
I
0
paroboloid,Z-axis.. x2/a2+y2/b2-a=0
Hyperboloid
X4o2+yZ/d-r2/C2=I(~h,~,X,+m) Hyperbolic
+cm
of one sheet, Z-axis..
paraboloid,
x2/02-y’/b*-cz=O
I-( h2/8+k2/b2-n?/c’)
X-axis..... (th,tk,im
+cm
‘one, Z-ax,s. h2/a2+ k2/b2-d/
c2
)
Physiological optics
45
refraction of the medium into which the ray is refracted. Recalling that the notation is S and S’ for u1 and u2 it follows that S’ = pS + rn
where P = n&2 and r is an undetermined multiplier. This can be done because the normal to the surface and the two rays are coplanar. The components-of this vector equation are k’ =
pk + I-K
1’ =
pl + TL
m’ =
Squaring and adding the component
pm + TM. equations, above, a quadratic
in r is obtained:
r2 + 2aT+b=O where Q = p(kK + IL + mM)/(K2 + L2 + M2) and b = (p2 - l)/(K’ + L2 + M’). The quadratic equation may be solved by the usual formula. The condition under which it fails to have a real root is b > a2 corresponding to total internal reflection. At a reflecting surface the same equations may be used with p = 1, and the other root of the quadratic. In this case b = 0 and the undetermined multiplier becomes
r = -2a. With this assignment of values the direction cosines (components) become
of the refracted ray
k’ = k - 2aK 1’ = 1 - 2aL m’ = m - 2aM. Step 4. Transformation of the refracted ray back into the global system. The new ray origin coordinates (X, Y, 2) and direction cosines (k’, I’, m’) may be transformed back into the global reference system by using the inverse forms of the transformations used to put the incident ray into the local coordinate system. These are
[
and
9]-R+;;]
46
WAYNE WIITANEN
Because R,,, is a product of unitary matrices it is itself a unitary means that its inverse and transpose are equal. That is, coscc~cosy + sina.sinfl.siny R$
=
cosa.siny
sin a cos /?
- sincc.sinp.cosy cosg~cosy
-cosp*siny -sincr.cosy
matrix and this
+ cosccsinpsiny 3. FOUR
-sina.cosy
RAY-TRACING
sin p COSG!.cosp I
- cosa.sinp.cosy EXAMPLES
To illustrate the power of the ray-tracing method, I have selected as examples the Gullstrand schematic human eye (Duke-Elder [S]), the schematic cat’s eye (Bishop and Vakkur [6]), the compound eye of the fly (Seitz [7]), and the eye of the honey bee (Varela and Wiitanen [S]). Example 1: The schematic human eye (Duke-Elder [S])
In Figs. 2(a) and 2(b) the geometry of the schematic eye, refractive indices placement of apertures, and sample ray paths within the system are shown. The rays traced are Gullstrands
schematic
human eye
I cm
I
/
(a)
(b)
Fig. are the All
2. Gullstrand’s schematic human eye. (ablements of the dioptric apparatus. Two apertures shown: one at the first vertex having radius 5.5 mm (corresponding to radius of the cornea), other at 3.6mm behind the first surface of radius 4.5 mm (corresponding to the pupil). measurements of surface positions and radii of curvature in millimeters. (bEselected ray paths on the same scale as 3(a).
Physiological
optics
47
paraxial with respect to the Z-axis (optical axis of the system) and are spaced at 0.5, 1.0, 2.0 and 4.0mm away from it. The forward shift of the focal point is apparent as the rays become more distant from the Z-axis (less paraxial). The focal point of the system, as computed by ray-tracing using paraxial rays, was found to be at 2.365 cm from the first surface. The focal point of the Gullstrand schematic human eye is generally accepted to be at 2.2875 cm. From this example we see that it is possible to select the origin and initial trajectory of rays incident on the optical system. We also note that each of the refracting surfaces is spherical, is specified by the equation for a sphere (using the radius of curvature measurements), and is positioned in the global system by means of the position of the first surface of each element. In all seven surfaces have been used in the description: the first surface of the cornea, the back surface of the cornea, the front surface of the lens core, the back surface of the lens core, the rear surface of the lens capsule, and finally the surface of the retina. In the next section the manner in which these rays and surfaces are indicated to the ray-tracing program will be explained and the necessary cards for each example will be given. Example 2: The schematic eye of the cat (Bishop and Vakkur [6])
Again a group of four rays are shown in Fig. 3(b). They are parallel to the Z-axis and are spaced at 0.5, 1.0, 2,0 and 4.0 mm from the axis. The geometry of the eye cot\ eye
Retina \
(b)
Fig. 3. Schematic cat’s eye. (a&geometry of the dioptric apparatus. Two apertures: cornea dia. 8.4 mm (at first surface), pupil dia. 7.3 mm (at third surface, 5.20 mm behind first surface). All surface positions and radii of curvature in mm. (btselected ray paths on the same scale as 4(a). C.8.M. 7/1-D
48
WAYNE WIITANEN
Plane
-_-----
Envelope of crystalline -
cone
---_ L-_--
of rhabdom
1 -_---_
I
II I
0
7.4
16.3
146.0 (a)
(b)
Fig. 4. Eye of the honey bee. (atgeometry micrometers. (btselected
of the dioptric apparatus. All measurements ray paths on the same scale as S(a).
in
is shown in Fig. 3(a). Again the refracting surfaces are spherical and are positioned in the global coordinate system by specifying the positions of the various surfaces. The anterior progression of the focal point, as the rays become less and less paraxial, is interesting. Also worth noting is that the focal point lies within the retina itself but anterior to the later of primary receptors. The focal point, as calculated by ray-tracing, was found to be 2.214cm behind the first surface of the eye. The position of the focal point as calculated by Bishop and Vakkur [6] was found to lie 2.2153 cm behind the first surface of the eye. The two results are in good agreement. As for the previous example, the control cards necessary to “set up” this system will be described in detail below. Example 3: The compound eye of the honey bee (Apis mellifera)
Figures 4(a) and (b) illustrate the geometry of the honey bee’s eye and some ray trajectories through the eye, respectively. As the rays become less paraxial we again observe the forward shift of the focal point. The calculated focal point of the bee’s eye was found to be at 108.54 pm below the first surface. This example illustrates the use of a non-spherical refracting surface, namely, the planes at 16.3 and 48.0 pm. Essentially these are spheres of infinite radius of curvature, but are nonetheless easy to represent to the ray-tracing program. Details of computer set-up will be given below. Example 4: The compound eye of the blowfly (Calliohora
erythrocephala)” Figures 5(a) and (b) show the geometry and the ray paths within the optical system, respectively. As the rays become less and less paraxial the focal point progresses anteriorly although not evident from Fig. 6. Ray-tracing locates the focal point as calculated * Seitz
[7],
Physiological
optics
49
Y
73.56 II 77.68 (a)
(b) Fig. 5. Compound eye of the blowfly. (awioptric apparatus geometry. All measurements in micrometers. (bkselected ray paths on the same scale as 6(a).
are
by Seitz [7] using the thick lens formula places the focal point at 73.40 pm below the first surface of the equivalent thick lens. The discrepancy is 13.47 pm. Details of the computer set-up will be given below.
4. USING
THE
RAY-TRACING
PROGRAM
General considerations Digital computer programs usually require some sort of control information preceding the actual program to direct the actions of the computer’s supervisory system. Because this information varies considerably from machine to machine, and even within different supervisory systems for the same machine, it will not be discussed. The ray-tracing program itself is written in pure FORTRAN IV, in the hope that it will be executable on any computer that will accept the language. For this reason the program should have a wider applicability than would one that does the same computations but is written in assembly language. Card groups and computational sequence In the most general terms, the ray-tracing program reads a group of cards describing the surfaces and rays, performs a series of calculations using these data (see Section
50
WIITANEN
WAYNE
Start Return to system. _ Run terminated
IRead input data (Card grol;ps 1,2,3)
No data
Print out a copy _
Loop
Rintout II error mesMge
TeZ+E
3~oZZons
Yes
No
Read fan data
Setup single ray
zIrnL;E * Have all Ypoints been calculated? No 1 SetupXincrements c( IncrementXpoints c Have allXpoints& been cal:ulated 7
I, Calculate direction cosines to put ray into standard (unit vector) format 2. Calculate elements ofR-matrix. 3. Transform ray vector to loco I system. 4. Calcu\at$ incident ray and tongent plane intersectlon point.
I
-
No incidence Error message
a
No convergence
-
Aperture stop Errwmessage
e
Aperture stop
-Error
No refraction message -
Internal reflection
Itemte
to surface Z
Apply aperture to ray path OK 1 Test for internal reflection
OK 1 refracted my vector components. 2. Calculate transmittance, reflectance, etc. 3. Tansform refracted ray back into global system 4. Print out data forth is surface if required.
I. Calculate
t Have all surfaces been used?Yes I I. Calculate intersection 2. Print out intersectionof No
Fig. 6. Flowchart
indicating
of ray with optica I oxis. ruy with optical axis.
+ Single ray?
the computational
Yes
procedure
c
used
Read another
single
by the ray-tracing
raycars
program.
l), and finally writes out the desired answers before it goes on to the next set of data cards. The sequence of ray-tracing computations is presented schematically in Fig. 6. The data cards used by the ray-tracing program fall naturally into three groups: 1. Set-up cards 2. Surface definitions 3. Ray specifications. Each (1) card, only
of these groups will be described in detail. Input card group 1. The first group (set-up cards) consists of three cards: a title card, and an “initial conditions” card. The first, the switches the first three columns of the card and is punched as follows: column
1:
a switches card, uses
O-ray specification cards (group 3 cards) will be for a fan of rays (see description under group 3 cards). l-ray specification cards (group 3 cards) will be for a series of individual rays, rmt a fan.
Physiological optics
51
&rotational symmetry will not be assumed (for the advantages of rotational symmetry see the discussion under group 3 cards). l-rotational symmetry will be assumed. &the PROCES subroutine will not be called and normal output column 3: will be printed. l-the PROCES subroutine will be called instead of producing the normal printed output. Note: a dummy PROCES subroutine is required even if it is not to be called (see the Appendix). FORMAT (311)
column 2:
Technical note regarding the user-provided subroutine PROCES
In order to provide for the possibility of performing calculations with the output from the ray-tracing program rather than producing a printed table, entry to a user-provided subroutine by CALL PROCES (I, J) has been provided at each point where the main (ray-tracing) program would produce printed output. The parameters Z and J define the conditions under which the subroutine PROCES was called: z=
1: 2:
good data ray does not intersect with a refracting surface, takes a virtual path, hits an aperture, or is not refracted 3: end of this study has been reached J = 1,. . .) n: the number of the refracting surface at which the current calculations have been made.
The two parameters, Z and J, alone are insufficient to be of any value for further computations, and so some method for passing along other computational values must be provided. A COMMON statement is set up as follows for this purpose: COMMON
IRAY, IMl, XA, RPARl, YA, RPARZ ZA, RPAR3, XOBAR(35), CKBAR(35), YOBAR(35), CLBAR(35), ZOBAR(35), CMBAR(35), X0(35), YO(35), ZO(35), ALPHA(35) BETA(35), GAMMA(35), A(35), B(35), C(35), D(35), E(35), F(35), G(35), REFR(35). APl(35). AP2(35), AP3(35), AP4(35), AP5(35), AP6(35), NAP(35), NOUT(35), CK, CL, CM, S
where IRAY IMl XA RPARl YA RPAR2 ZA RPAR3 XOBAR
= = = = = = = = =
YOBAR = ZOBAR = CKBAR CLBAR CMBAR X0 YO ZO
= = = = = =
the number of the ray being computed surface number (same as J in the calling sequence) 1st X-coordinate of ray 2nd X-coordinate of ray 1st Y-coordinate of ray 2nd Y-coordinate of ray 1st Z-coordinate of ray 2nd Z-coordinate of ray X-coordinate of intersection (ray with surface number J) Y-coordinate of intersection (ray with surface number J) Z-coordinate of intersection (ray with surface number J) X directional cosine of ray leaving surface J Y directional cosine of ray leaving surface J Z directional cosine of ray leaving surface J X-coordinate of the vertex of surface J Y-coordinate of the vertex of surface J Z-coordinate of the vertex of surface J
WAYNE WIITANEN
52
ALPHA = Euler angle of surface J with respect to the reference system BETA = Euler angle of surface J with respect to the reference system GAMMA = Euler angle of surface J with respect to the reference system A, B, C, D, E, F, G = surface equation coefficients of surface J REFR = index of refraction of surface J-1 APl, AP2, AP3, AP4, AP5, AP6 = aperture parameters for surface J NAP = aperture type for surface number J CK = X directional cosine CL = Y directional cosine CM = 2 directional cosine S = ray distance from tangent plane to surface number J.
By providing the subroutine PROCES with a COMMON statement identical to that shown above it is possible to refer to the data used in the main program. The second card of the first group is a title card. This card will cause the title punched in columns l-72 to appear at the top of each page of printed output. Its format is columns l-72: title information FORMAT
(12‘46)
[A = alphabetic]
The last card gives information about the number of surfaces, the convergence tolerance of the Newton-Raphson iteration method (used by the program to find the point of intersection of the incident ray with the refracting surface-see Appendix 1 for details), and the index of refraction of the medium in which the object lies. The format of this card is columns
l-l 1: the 12-14: the 1532: the FORMAT
initial index of refraction (Fl 1.5) number of surfaces in the system (13) convergence criterion (E15.8) (F11.5, 13, E15.8).
Here is a typical set-up for the first card group: col.
1
5
1 0
1 5
2 0
010 THIS IS A SAMPLE TITLE CARD 1.0 004+0.10000000E-05
(switches) (title) (i.c.).
The switches card says that a fan of rays is to be used as input, that rotational symmetry is to be assumed, and that the PROCES subroutine is nut to be called. The title card provides a heading for each page of printed output from the program. The initial conditions card says that the index of refraction of the object space is 1.0 (free space), that there are 4 refracting surfaces in the optical system, and that 6 places of accuracy are required of the iteration procedure. (2) Input curd group 2. The second group of cards is used to define the shape of refracting surfaces. The surface equation used by the ray-tracing program is AX2+BY2+CZ2+DX+EY+FZ-G=O.
By selecting appropriate values for the coefficients, A-G, any quadratic surface may be described. By referring to Table 1 appropriate combinations of coefficients may be selected to define most refracting surfaces. Each surface is defined in its own local (or independent) coordinate system, without any reference to the overall coordinate system in which the rays are to be traced (the global, or reference, system). In order
Physiological optics
53
to give the local coordinate system an orientation with respect to the reference system (that is, the one in which the rays are being traced), the origin of the local system with respect to the reference system must be supplied, as well as any rotation of the local system with respect to the reference system. To this end card group 2 is subdivided into three sections: (a) Parameters for placing the local coordinate system of the refracting surface within the framework of the reference coordinate system; (b) the surface equation (relative to the local system); (c) definition and placement of apertures in the local coordinate system. These cards must be supplied for each surface in the optical system, and in the order given. To place a surface in the reference system a card is made up as follows (each number is punched in the format F11.5): card 1:
columns
1-11: 12-22: 23-33 : 34-44: 45-55 : 5666:
FORMAT
X-coordinate of the local system origin with respect to the reference system origin Y-coordinate of the local system origin with respect to the reference system origin Z-coordinate of the local system origin with respect to the reference system origin c( Euler angle p Euler angle y Euler angle
(6Fll S).
The CI,0, and y Euler angles define the rotation of the local coordinate system with respect to the reference system, and they are illustrated in Fig. 1. Bear in mind that the order in which the rotations are done is important: first the CIrotation, then the p, and finally the y rotation. In columns 73-80 an identification may be punched. I usually punch in the label n-A where n stands for the surface number associated with this coordinate placement. Each coefficient of the surface equation is punched (in the F11.5 format) as follows card 2:
columns
FORMAT (6F 11.5) card 3: columns
FORMAT
l-11: 12-22 : 23-33: 34-44: 45-55 : 56-66:
A B
l-11: 12-22 :
G (F11.5) refractive index of the medium following this surface. (F11.5)
C D E F
(F11.5) (F11.5) (Fl 1.5) (F11.5) (F11.5) (F11.5)
(2F11.5).
The coefficients are usually punched starting at the left end of each field. It is convenient to label these cards n-B and WC in columns 73-80, respectively. Apertures are of three types: circular-annular, rectangular-trapezoidal, and hyperbolic. The parameters for specifying the type and dimensions of an aperture are punched on the last card of this group. All are in the F11.5 format except those specifying the aperture code and the output code. The parameters are: card 4:
columns
l-11:
X-coordinate of the coordinate of the center of a circular-annular aperture, or the coordinate of the center of a hyperbolic aperture, or the X lowerbound of the base of a rectangular-trapezoidal aperture.
WAYNE WIITANEN
12-22 : 23-33 :
3444 :
45-55 :
56-66:
67:
68:
FORMAT
Y-coordinate of the same as above. Inner radius of a circular-annular aperture, or the length of the semi-major axis of a hyperbolic aperture, or the Y-lower bound for a rectangular-trapezoidal aperture. Outer radius of a circular-annular aperture, or the length of the semi-minor axis of a hyperbolic aperture, or the Y-upper bound for a rectangular-trapezoidal aperture. Y-lower bound for a hyperbolic aperture, or the reciprocal slope of the left-hand side of a rectangulartrapezoidal aperture. Enter 0.0 for a circular-annular aperture. Y-upper bound for a hyperbolic aperture, or the reciprocal slope of the right-hand side of a rectangulartrapezoidal aperture. Enter 0.0 for a circular-annular aperture. Aperture code : 0: circular-annular 1: rectangular-trapezoidal 2 : hyperbolic. Output code : 0: no output at this surface 1: output required at this surface
(6F11.5, 211).
The various apertures are illustrated in Fig. 7. I usually label this card n-D in columns 73-80, where yt, as usual, is the number of this surface. This group of four cards is repeated for each refracting surface being defined. (3) Input curd group 3. The last group of cards specifies the rays to be traced through the system. There are two subcategories: a fan of rays, and individual rays. The first of these requires two cards, and their formats are: card 1:
columns
1-11: 12-22: 23-33 : 3444: 45-55 :
Fig. 7. Aperture
types.
The maximum X-coordinate of the ray origin The Z-coordinate of the ray origin The increment to be used to obtain successive rays The maximum X-coordinate of the second point on the ray The Z-coordinate of the second point on the ray
A-rectangular
or trapezoidal. (shaded area passes
B--circular rays).
or annular.
C-hyperbolic
55
Physiological optics
r
Fig. 8. Isometric view illustrating generation of fans of rays (see text).
The 2 increment to be used to obtain successive second point on a ray. l-11: The maximum Y-coordinate value of the second point card 2: columns on a ray The Y increment to be used to obtain successive points 12-22: on the ray. FORMAT (6F11.5/2F11.5) (the symbol “/” means “next card”). 5666 :
I usually label these cards FAN-l and FAN-2 in columns 73-80. Figure 8 illustrates how the rays are defined and how the fans are generated. Rotational symmetry (as may be specified on the switches card of group 1) means that the entire lens system is rotationally symmetric about the optical axis (Z-axis) of the reference system. Advantage can be taken of this fact to cut down the number of rays to be traced through the system, without sacrificing any information. Should rotational symmetry not be requested then the rays to be traced are generated in the following order: 1: X-coordinate of the first ray point--(X,,) (AX) (-X,,,). 2: X-coordinate of the second ray point+Xk,,) (AX’) (-XL,,.). 3 : Y-coordinate of second ray point+Y,,,) (AY) (- Y,,,). Note: Y-coordinate of the first point is (0) (0) (0). The stepping scheme starts with the Y-coordinate (3), then the second X point (2) and finally the first X point (1). This sequence gives a total of rays. When rotational symmetry is used, the X-coordinate of the second ray point runs from XL,, to 0 by decrements of AX’, and the Y-coordinate of the second ray point runs from Ym, to 0 by decrements of AY, thus reducing the number of rays traced by one-fourth (and consequently cost). A trick to confine the rays to the XY-plane is done by setting Ym, to zero and the Y decrement (AY) to any positive, non-zero number. Figure 9 illustrates the difference between a full fan of rays and those generated by invoking rotational symmetry. Single rays may be traced through the system by specifying the origin and direction of each ray. This specification requires one card per ray in the format: card 1:
columns
l-11: 12-22 : 23-33 :
X-coordinate of the ray origin in the reference system. Y-coordinate of the ray origin in the reference system. Z-coordinate of the ray origin in the reference system.
WAYNEWIITANEN
56
Fig. 9. Full fan of rays. Compare with Fig. 8 to see advantage of rotational symmetry in generating rays
3444: 45-55 : 5666: 67:
68:
69-7 1:
FORMAT
X-coordinate of a second point on the ray, or the X directional cosine of the ray leaving its origin. Y-coordinate of a second point on the ray, or the Y directional cosine of the ray leaving its origin. Z-coordinate of a second point on the ray, or the Z directional cosine of the ray leaving its origin. Input code for ray specification type: 0: ray is specified by two points 1: ray is specified by one point and directional cosines Optical axis computation code: 0: do not compute the intersection point of the ray with the optical axis 1: compute the intersection point of the ray with the optical axis. Ray identification number (right justified). If this number is less than or equal to zero then it indicates that there are no more rays to be processed, including the one on this card.
(6F11.5, 211, 13).
Individually defined rays are illustrated in Fig. 10. After all the rays of a fan have been processed, or after the last single ray has been processed, then the program will return to its beginning and attempt to read in another set of data cards (groups 1, 2, and 3). If there are no more cards the program then returns to the supervisory system. If there are more cards then the program will perform the ray-tracing studies indicated by them. Computational errors and the output format During computation situations may arise that are indicative of errors or unusual situations of which the user should be apprized. (1) If a fan of rays is being used and if the step size is incorrectly specified, usually resulting in a loop in the program,
Physiological optics
Fig. 10. Illustration of three independently
51
specified rays (1, 2 and 3).
the user will be notified of his error by the program, and the program will go on to the next study, if any. (2) It is possible that certain rays will not strike a surface. When this happens the program informs the user that a certain ray, given by number. has had no incidence with the surface given. This arises when there is no convergence with the current surface during application of the Newton-Raphson iteration, or if the convergence is oscillatory making it impossible to get within the required convergence tolerance. (3) Sometimes a ray may be reflected back along its original path. Here the program advises the user that a virtual path for the ray has been encountered, and it goes on to the next ray. (4) Apertures provide obstacles for rays within the optical system, and rays that hit them are lost to subsequent processing. When a ray strikes an aperture, the user is informed and the program goes on to the next ray. (5) Internal reflection is not handled by the program, and whenever an internal reflection is encountered the user is informed that there has been no refraction, and the program continues with the next ray. The printed output from the program has a number of prefatory lines which repeat the surface and ray parameters used for the particular study in a readable format so that the values the program used can be checked. In addition, there is some information concerning the maximum number of rays to be traced and the maximum number of lines to be printed out. Sometimes the actual number of lines printed is less than the maximum because of rays blocked by apertures or lost. The bulk of the printed data is organized in a columnar format. The columns, from left to right (see Fig. 11) identify the number of the ray being traced, the surface at which the ray has arrived and at which this output is being taken (after refraction of the ray), the name of the coordinate point or directional cosine being output, the value of the object point, the second object point or directional cosine (these two numbers, the first and second object points, or directional cosines, define the path over which the original ray was travelling prior to encountering the first refracting surface), the point at which the ray intersects the surface (X, Y, Z-coordinates in the reference system), the directional cosines of the ray as it leaves the surface, and the reflectance/transmittance data. For each ray there are three lines of output: one for the X-coordinate (or k directional cosine), one for the Y-coordinate (or 1 directional cosine), and one for the Z-coordinate (or m directional cosine). By examining the object point and the second point/directional
58
R.1
216
WAYNE
SURF 1
u&J
PI 0
2
*.“0GCc
21t
‘
C.cdtCc C. cccc~” 4. ‘)Jctc E. “Jcr 3 c.cmCC 4.c7crc
216
5
216
4.51021 -“.O”“O”
‘"7.58626 6.28155 -C.CCCCO ll".eoo"C 3.4229" -c.cccI)o 123.02132 -0.79568 C.CCCOO 173.5t.co0 -,.IU3L C.OCCCF L77.68"""
c.cococ G. ‘~,“CO 4rcPxwJ C. O’)CC?
C.CC0PC C.CbCCC 0.CJC‘C c. cord 0 ‘63.99725 ~.“5CCC 0.3F.W” C.“i”“3
*.ooocc
217
3
2‘7
211
L.000~3 “.OICL‘ *.c;C“) 0.03‘“‘) C.CJCCO 4.W)C‘O U.CSOL~ b.
4.
217
c*,cc‘1
n?ccc
“3CCU ~.Cmcc ~.CJL‘C C.CxJ“3 (I. 0u**n (r.
217
2.C[i”“”
“.CJ0”” 10c.c00c0 r.c3coo
c .coccc
1"C."""C0 C.Cmx!C J.C"COC ~“.C*cCo ~."0C"" rJ.CDCCC ‘C".C*CC~ *.m"FC C.‘,CCC ‘"C.CPOCC ~.CCCCO C.CCCCC 1CC.CJCOC
+.0000c -U.CCCCC
‘m.3dl51 3.56770 -C.EC~C‘ ‘"7.4"294 3.36363 -0.C"CCC llC.8CC30 .?*6"24C -C.CCCOC 123.13536 -1.CC403 O.FCbCC 173.56""" -1.3099C C.OCO"‘ 117.e~‘CC
4.0708" 0.00000 C.99799 -0.06881 C.0"""" 0.997b3 4.07665 0.""""" C.99750 4.08318 0.00000 C.99653 -0.08337 ".0f"0" C.99652 -".oa287 C.CCCOC 0.99656
-0.06‘25 0.09000 0.990‘1 -"."5996 0.0"""" 0.99B2" -r.c%157 0.0"""" 0.95810 4.07,+* P.CCOIIC C.99744 4.‘)7160 C.CPCCO 0.997.3
-"."7118 C.""""" 0.99716
L/JlPI/S~R, 0.366.5E_01 0.35517E41 ".36695E-o, “.44lu)E~~
0.*38136-04 c.+C2llE-O* ".11555E-03 ".lC569E-03 ".16629E-03 ".66459E-"3 ".61296E-03 ".6295CE-03 O.L24~SC-05 ".11813E-"5 O.l‘85+E-05 0.88422E-05 ".8366.E-05 0.8395bL-"5
ll~~Tl/Slll 0.96335E 00 0.9336Q 00 0.99lClE 00 0.99995E 00 0.93lUE 00 ".91729E 00 0.99982E 00 ".9++83E CO 0.9~72oE 0" 0.99939E 00 0.9209Y CO ".949,0E 00 ".L00"?E 01 0.946leE 00 ".94969E 00 o.99999E JO 0.966‘8k 00 0.9.995E CO
".366"3E4, 0.35993E41 ".36603C-01 n.~613*E-O4 ".9*32lE4+ C.**580t-c~ C.‘T55.3E-"3 ".16721i-03 ".lb75‘E43 ".6631‘E-03 C.621822-"3 C.b326lE-"3 ".12179E-"5 0.1‘8(99E-OS C.l‘PlPE45 ".81136E-05 0.8.252245 ".8+468E-"5
0.963KlE 0.9160X 0.9539lE 0.9999SE 0.99832E 0.95101 C.PPPOLE ".95.?lY 0.953952 0.999392 C.931olE ".955,3E 0.10000E 0.95269E
00 0" 00 CO JO 00 00 co 00 CO ,0 00 01 00
0.955I.E
CO
0.3bSBlE-CL ".36268E-01 0.36587E-01 ".4b732E41 ".9969*C-01 ".4+8"9E-"9 0.11556143 ".16812E-03 O.l6834E-03 ".bb2SZE-03 0.629222-03 0.6352lE-03
0.963.Lk 0" O.QSSClE DO 0.958&x 0" 0.9999SE 00 0.95639E CO 0.95888E OJ 0.99902E 00
0.99999E 00 ".95268E 00 C.955‘OE 00
lSC.3951~
217
l .C.Jt‘C
213
218
218
22a
cos
“Ill
~.wccc ‘.cOCr u.c”*Cc
21t
WIITANEN
‘
"."3""3 C.C'XIC 4.C'Vh‘ C.OCO‘O L.CNC.0 ,.CJCCC cl.moFC O.CJCCC 4.000"" G.U')CCO C.CJCCC
Fig. 11. Sample
output
3.050”” C. 000”” 10n.00000 3.w000 0 .COOC’) 1w.00000 3.coccc C.WC0C
1cC.ccccc 3.00000 c.00cco 100.0000t
2.99781 -C.CCCCO
100.21314 2.62604 -C.CCCFO 107.20368
2.,,4425 -0.00000 11".00"00 1.70541 -O.oCoLC 123.22454
from ray-tracing
4.05255 C.0U000 0.99062 -0.05156 0.0000" 0.99867 -0.05295 6.00000 C.9986" -0.06625 0.00""" C.99Bl.9
program.
0.95lwE C.95aleE 0.99932E 0.94911E
0.95919E
00 00 00
00 00
SW text for discussion.
cosine data on each of these lines we may determine the trajectory of the ray before it strikes the jirst refracting surface. By looking at the intersection point data of each line we may obtain the coordinates of the point at which the ray intersects with the current surface (number under “SURF” column) indicated. The three directional cosine values for the ray are the directional cosines of the refracted ray as it leaves the surface (and they are, therefore, the directional cosines of the incident ray on the next refracting surface). Although the reflectance, energy of the refracted/reflected wave, and the intensity of the refracted/reflected wave appear to be related to the coordinate axes, lying as they do on the lines so labelled, they are not related to the axes. The data associated with the X-coordinate line corresponds to the reflectance (R) under the column labelled R/J (R)/S(R) while associated with the Y-coordinate line the energy (J(R)) of the reflected wave is to be found, and associated with the Z-coordinate is the intensity (S(R)) of the reflected wave. The data under the column headed T/J(T)/S(T) follows the same pattern but is for the transmissivity (T), energy (J(T)) of the transmitted ray, and the intensity (S(T)) of the transmitted ray. In addition to this information, whenever a ray crosses the optical axis (Z-axis) the point of that crossing is noted on the printout. Extensions to the program and a word of warning Several extensions to the program recommend themselves. One would be to include an option for plotter output so that the dioptric system could be drawn out and individual rays traced through it. Two different versions might be implemented: a sectional view of the system cut by a selected plane, and an isometric view. Another useful extension would be to increase the generality of the refracting surface equation by including, perhaps, exponential and transcendental functions in a general
Physiological
59
optics
form. This modification would have the disadvantage that the number of coefficients required as input data would go up. Considering the somewhat ad hoc method currently being used for handling apertures, a generalized aperture describing system would be worthwhile. To incorporate such a system the shape of an inner and an outer aperture would have to be defined in much the same way as the surface equations. In addition to the two surface equations that would be required for the inner and outer apertures, a card specifying the coordinate position and the Euler angles of the aperture pair would be required. A word of warning might be in order for persons planning to use the program as it stands on a computer having a word length of fewer than 36 bits: the machine being used may not have sufficient digits of accuracy for convergence to be attained. If this is true, then the program should be compiled in double precision, or relaxation of the convergence tolerance be instituted. To illustrate the use of the control and definition cards for the ray-tracing program, I will describe the cards needed for each of the examples presented earlier. Exanrple 1. The schematic human eye. In Table 2(a) the parameters describing the refracting surfaces used in the Gullstrand model for the human eye are listed. All measurements are in millimeters. In Table 2(b) are the computer-readable forms of the same parameters, arranged in the three card groups described above. (Note that the cards marking off the groups in Table 2(b) would never be present in an actual computer run,) Group 1. The switches card shows that a fan of rays has been selected (zero in column one), that rotational symmetry will apply (one in column two), and that the auxillary subprogram PROCES will not be used (zero in column three). The second card of this group is the title, and the last card defines the initial conditions: the refractive index of object space is 1.0; there are seven surfaces to be defined (Group 2 cards), and an accuracy of lop7 is required in the Newton-Raphson iteration. Group 2. Here are the seven quadruplets of cards that describe the seven refracting surfaces. Card A describes the coordinate position of the vertex of the surface and any coordinate rotation of that surface. Cards B and C contain values of the parameters for this surface (coefficients A, B, C, D, E, F, G) and the refractive index of the surface following it. Card D describes the aperture to be used for this surface. Group 3. This group defines the fan of rays to be used in the tracing process. The origin of the rays in one meter (1000 mm) to the left of the first refracting surface’s vertex (defined to be at the origin of the reference coordinate system). All rays are confined to the XZ-plane (done by setting Y = 0.0 on card FAN-2). Points defined at the ray origin will go from 8.0 to - 8.0 (mm) by steps of 0.5 mm while the points at the first surface (destination of the initial ray segment) will go from 4.0-O.Omm by steps of 0.5 mm. Example 2. The schematic eye of the cat. Table 3(a) illustrates the parameters used to describe the schematic cat’s eye, and in Table 3(b) the computer-readable form of the same parameters is given. The description follows that of the first example with Table 2(a). Gullstrand schematic human of the surface. The radius is the radius with respect to the origin of reference center defines the center of the sphere system in which the refracting surface
I 2 3 4 5 6 7
00 O.? 3.6 4 146 6 565 7.2 244
-77
-6X -100 -7911 +57h +60 +120
iw. (0.0. (0.0. (0.0. (0.0. (0.0. (0.0.
eye. The position of each surface is that corresponding to the vertex of curvature, and a plus sign means that the surface is concave system, while a minus sign means that the surface is convex. The defining the refracting surface with respect to the local coordinate is defined. R.I. is the refractive index of the medium following the surface
- 7 7) - 6.8) - IO.01 7.9 + 5.76) + 6.0) + 12.ot
II)
1.376 I.336 1.3X6 I.406 I.386 336 1.336
I
\.1 + \2 + y2 + 22 + y1 + y2 + $ +
.?’ + ?I + ,,1 + I‘. + ;,z + p + ?2 +
1
;z ;z ;z .5: :;: ;i
= = ; = = = =
ii ,.I ’ :: $ $ ).I
1 3
I I
I100 1 I 0 0
1
I
I I 1 I
I I I I
I
1 I I 1
0
0 0 0 0
0
0 0 0 0
-154 -13.6 ~20.0 -15822 +11.52
+120 +24.0
0 0 0 0 0 0 0
WAYNE
WIITANEN
Table 2(b) GROUP I 010 FAN/SYMMETRY/NO PROCES RAY PATH ANALYSIS OF GULLSTRAND 1.0 007+0.10000MX)E-06 GROUP
EYE - WIDE
PUPIL
5112169 HEADER
2..
.O
0.0 1.0 1.376 0.0 0.0 IO 1.336 0.0 0.0 1.0 1.386 0.0 0.0 1.0 I.406 0.0 0.0 IO 1.386 0.0 0.0 1.0 1.336 0.0 0.0 10 1 336 0.0
1.0 0 0 .O 1.0 .O .O .O 1.0 .O .O .O 1.0 .O 0 .O IO .O .O
0 1.0 0 0 0 1.0 0 .O GROUP 8.0 .O
SCHEMATIC
0.0 1.0
0.0 0.0
0.0 0.0
0.0 - 15.4
0.0 0.5 1.0
5.5 0.0 00
0.0 0.0 0.0
0.0 0.0 - 13.6
01
00 36 1.0
5.5 0.0 0.0
0.0 00 0.0
0.0 0.0 - 20.0
01
0.0 4.46 1.0
4.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 - 15.822
01
0.0 6.565 1.0
4.5 0.0 0.0
0.0 00 0.0
0.0 0.0 Il.52
01
0.0 7.2 1.0
4.5 0.0 0.0
0.0 0.0 0.0
0.0 0.0 12.0
01
0.0 24.4 1.0
4.5 0.0 00
0.0 0.0 00
0.0 0.0 24.0
01
0.0
25.0
0.0
00
01
0.5
40
0.0
0.5
IA IB IC ID 2A 2B 2C 2D 3A 39 3c 3D 4A 4B 4C 4D 5A 5B 5c 5D 6A 68 6C 6D 7A 7B 7c 7D
3.. - lrn.0 I.0
FAN-l FAN-2
Table 3(a). Schematic eye of the cat. Same remarks apply to this table as to the table for the Gullstrand schematic human eye SUrfaCe
Position
I , 3 4 5
0.0 0.68 5.20 13.70 21.83
Radius
Center
-8.57 - 7.89 - 7.20 + 8.05 + 12.03
(0.0, (0.0. (0.0. (0.0. (0.0.
+ +
R.I.
8.57) 7.89) 7.20) 8.05) 12.031
Equation
1.376 1.336 I 554 1.336 (1.3361
x1 x2 .x2 x2 .x2
+ + + + +
yz yz y2 y2 y’
+ + + + +
it 2’ z* 27 it
= = = = =
12 ? rl rz rl
A
B
C
D
I
1
I
0 0
1 1 I
I I I I
0 0 0 0
I 1
1 I
1
F
G
-17.14 - 15.78 - 14.40 + 16.10 f24.06
0 0 0 0 0
E
0 0 0 0
Table 3(b) GROUP I I IO SINGLE RAYS/ROTATIONAL SYMMETRY/NO SINGLE RAYS FOR THE SCHEMATIC CATS EYE. 1.0 txl5+0.ltxmmOE-06 GROUP 0 10 .O .O .O 1.0 .O .O .O 1.0
IO
0.0 0.0
0.0 0.0
0.0 - 17.14
0.0 0.68 1.0
84 0.0 0.0
0.0 0.0 0.0
0.0 0.0 - 15.78
01
0.0 5.2 1.0
8.4 0.0 0.0
0.0 0.0 0.0
0.0 0.0 - 14.4
01
0.0 13.7 I.0
7.3 0.0 0.0
0.0 0.0 0.0
0.0 0.0 16.1
01
0.0 21.83 1.0
4.75 0.0 0.0
0.0 0.0 0.0
0.0 0.0 24.06
01
1.0
1.336 0.0
0.0
12.5
0.0
0.0
01
IcmO 1000.0 - 1000.0 1000.0 0.0
00 0.0 0.0 0.0 0.0
0.5 1.0 2.0 4.0 0.0
0.0 0.0 0.0 0.0 0.0
01001 01002 01003 OlGil4 mm
0.0 1.0 1.376 0.0 0.0 I.0 1.336 00 0.0 1.0
1.5544 0.0 0.0 IO I 336 0.0 0.0
.O
.O 1.0
.O 0 0 IO 0
0
0.0
3.. 0.5 1.0
.O
2.0
0
4.0 0.0
.I)
HEADER
2.
.n
GROUP .O 0
PROCESS CALL I9 DECEMBER 1969
IA IB IC ID 2A 29 2c 2D 3A 38 3c 3D 4A 4B 4c 4D 5A 5B 5c 5D
Physiological Table
4(a). The optical
system
of the honey
Position
Radius
Center
R.I.
I 2 3 4 5 6
0.0 7.4 16.3 48.0 48.0
-43.0 -43.0 * 21.0 m
(O,O, - 43.0) (O,O, - 43.0) (0.0s (0.0.21.0) (O,OmJ
1.4896 1.4524 1.4350 1.3114 1.3477 1.3470
146.0
cl2
(0.O.Z)
61
bee (Apis mdlijka).
Surface
)
optics See the remarks
Eauation ‘iJ + ),” + :z = I2 g + y1 + ;z = ?-I .-=0 .x1 + y2 + ? = rz I=0 0
under
Table
2(a)
A
B
c
D
E
F
G
I
I
1
0
0
0
0
0
0 0
0 0
0 0
-86.0 -86.0 1.0 42.0 1.0 1.0
0 0 0 0 0 0
11100 0 0 1 0 0 0 0
I
I
0 0
Table 4(b)
I
GROUP 010 FAN/NO PROCES/SYMMETRY APIS MELLIFERA - COMPLEX DIOPTRICS 1.o 006+0.1tmtxwE-06 GROUP .O 1.0 .O .O .O 1.0 .O .O .O .O
RAYS FOR FOCAL
POINT HEADER
2. 0.0 1.o 1.4896 0.0 0.0 1.0 1.4524 0.0 0.0 0.0 1.4350 0.0 00
.O
.O .O I .o .O
1.0
1.3114 0.0 0.0 0.0 1.3477 0.0 0.0 0.0 1.3470 0.0
.oo .O .O .O .O .O .O .O .O GROUP .Ol .O
- PARAXIAL
0.0 1.0
0.0 0.0
0.0 0.0
0.0 - 86.0
0.0 7.4 1.0
16.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 - 86.0
01
0.0 16.3 0.0
16.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 1.0
01
0.0 48.0 I.0
16.0 0.0 00
0.0 0.0 0.0
0.0 0.0 42.0
01
0.0 48.0 0.0
16.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 1.0
01
0.0 146.0 0.0
14.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 1.0
01
0.0
2.0
0.0
0.0
01
0.0025
0.01
0.0
0.0025
IA LB IC ID 2A 2B 2c 2D 3A 3B 3c 3D 4A 48 4c 4D 5A 58 5C SD 6A 6B 6C 6D
3.. - loo.0 1.0
FAN-I FAN-2
one exception: the switches card (first card of Group 1) indicates that single rays are to be provided to the system (one in column 1). This means that the cards of Group 3 must be in the single ray format. Four rays will be traced (ray numbers 001, . . ., 004) and tracing will be stopped after ray 004 has been processed-the fifth card in Group 3, having a ray number of 000, will cause the tracing to be terminated. For each of the rays to be traced the ray origin will be one meter to the left of the first refracting surface, and the rays will be specified by means of two coordinate points (zero in column 68). The point of intersection of the ray with the optical axis is to be calculated and printed (one in column 69). The point coordinates are grouped as (X1.Y1,Z1), (x;,Y;,z;). Example 3. The compound eye of the honey bee. In Tables 4(a) and (b) are shown the Darameters that describe the bee’s eye (Varela and Wiitanen I2311 - -. and their computerreadable form. Only one new variation has been introduced: surfaces 3, 5, and -6 are planes, and as such have a radius of curvature equal to infinity. Their equations are particularly simple. Otherwise the card groups and descriptions follow those given for Example 1 (above). Table SlKfZKe
I 2 3 4 5 6
.5(a). The optical
system
of Culliphora
erythrocephda.
See Table
Z(a) for remarks
Position
Radius
Center
R.I.
Equation
A
B
C
D
E
F
G
0.0 7.1 I08 23.3 73.56 77.68
-21.16 -21.16
(0.0. - 21.16) (0.0. - 21.16) (0,0/m) (0.0.21.16) lO.O.cc I (0.O.K 1
1.4730 1.4530 1.4150 1.3370 1.3410 1.3490
.g + $ + .-I = 12 .x2 + ).1 + ? = r2 0 \.2 + J’ + -2 = ,i ;=O ; = 0
I I 0 1 0 0
I I 0 1 0 0
I I 0 1 0 0
0 0 0 0 0 0
0 0 0 0 0 0
-42.32 -42.32 1.0 42.32 1.0 1.0
0 0 0 0 0 0
21716 x x.
62
WAYNE WIITANEN Table
GROUP I 010 FAN/NO PROCES/SYMMETRY ANALYSIS OF RAY PATHS IN THE EYE OF CALLIPHORA 1.0 GROUP .O 1.0 .O .O .O 1.0 .O .O .O .O .O .O .O I .o .O .O
(SEITZ
CASE 19)
HEADER
OO6+O.IOOmmE-M 2.. 0.0 IO 1.4730 0.0 0.0 1.0 I .4530 0.0 0.0 00 I.4150 0.0 0.0 1.0 1.3370 0.0 0.0 0.0 1.3410 0.0 0.0 0.0 1.3490 0.0
.O .O
.O .O .O 0 .O .O GROUP 4.0 .O
S(b)
00 1.0
0.0 0.0
0.0 0.0
0.0 - 42.32
0.0 7.1 1.0
II.75 00 0.0
0.0 0.0 0.0
0.0 0.0 -42.32
01
0.0 10.8 0.0
Il.75 0.0 0.0
0.0 0.0 0.0
0.0 0.0
01
0.0 23.3 1.0
II.75 00 0.0
0.0 0.0 0.0
0.0 0.0 42.32
01
0.0 73.56 0.0
10.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0
01
0.0 77.68 0.0
3.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 1.0
01
0.0
2.94
0.0
0.0
01
0.5
4.0
0.0
0.5
1.0
IO
IA IB IC ID 2A 2B 2c 2D 3A 3B 3c 3D 4A 4B 4c 4D SA 5B 5C 5D 6A 68 6C 6D
3.. 1ocQ.o 10
FAN- I FAN-?
Example 4. The compound eye of the blowfly. Tables 5(a) and (b) present the parameters that describe the fly’s eye (Seitz [7]) and their computer-readable form. The card groups follow the description given for Example 1.
5. DISCUSSION In the preceding sections the details involved in using a computer-based ray-tracing technique have been presented along with some elementary examples illustrating its use. The intent was to introduce the worker in physiological optics to the method and to build a degree of confidence in using the program. A byproduct of the choice of illustrative examples was the discrepancy between the focal point of the Gullstrand schematic eye as calculated by the Gaussian lens formula and ray-tracing methods, and that between the focal point of the compound eye of Calliphora as calculated by the Gaussian formulae and ray-tracing methods. Should one wish to pursue the problems of ray-tracing further using the examples given above, some interesting things might be done: the effect of lens accommodation in the two simple eyes could be studied readily, as could the effects of changing the pupil diameter. In the latter, a better placement of the aperture could be effected by introducing a dummy plane (having the same refractive index as the medium in which it is embedded) into the system at the point where the aperture is to be introduced. The interested reader will find further examples of ray-tracing with respect to visual systems and model eyes in Lotmar [9] and Katzper [lo]. 6. SUMMARY In this tutorial paper the author has presented the basic processes of ray-tracing using digital computers. Illustrations of the technique were drawn from the literature about animal visual systems. A complete listing of the program is provided along with indications for its use. REFERENCES 1. H. Lass, Vt~ctor and Trrmr Amdysis. McGraw-Hill, NY (1950). 2. G. H. Spencer. A Germal Ray Tracing Prqram. IBM Corp. Research Paper RC-549 (1961).
Physiological optics
63
3. M. Born and E. Wolf, Principles ofOptics.Pergamon Press, NY (1965). 4. J. B. Scarboroueh. Numerical Mathematical Analvsis. Johns Hopkins Press, Baltimore (1962). 5. Sir W. S. DukerElder. A Textbook of Ophthalmology. Mosby, NY (1941). 6. P. 0. Bishop and G. J. Vakkur, The schematic eye in the cat, Vision Res. 3. 357-381 (1963). 7. G. Seitz. Der Strahlengang im Appositionsauge von Calliphora erythrocephula (Meig.), Z. uergl. Physiol. 59. 205-231 (1968). 8. F. G. Varela and W. Wiitanen, The optics of the compound eye of the honeybee (Apis me[lifera), J. Gem Physiol. 55, 336358 (1970). 9. W. Lotmar. Theoretical eye model with aspherics, J. Opt. Am. Sot. 61, 1522-1529 (1971). 10. M. Katzper, Simulation of optics experiments in Computers in Undergraduate Science Education (Commission on College Physics) (1971). About the Author-WAYNE A. WIITANENstudied mathematics at the Massachusetts Institute of Technology from 1953 to 1957, and received an A.B. in extension studies in 1967, an M.A. in Biology in 1969 and a Ph.D. in Biology in 1972, all from Harvard University. From 1954 to 1968 Dr. Wiitanen worked in the computer field in a variety of capacities: from hardware design to programming language design and implementation. Dr. Wiitanen is currently an Assitant Professor of Biology with a joint appointment in the Department of Computer Sciences at the University of Oregon, where he has been since 1971. He is also program director of the Systems Physiology Training Program in the Department of Biology at the University of Oregon. His research interest is the application of mathematics and computer science to biological problems with particular emphasis on neurobiology. He is a member of the American Association for the Advancement of Science. the Institute of Electrical and Electronics Engineers, and the Society for Computer Simulation.
APPENDIX Listing of FORTRAN program. DISCLAIMER Although the author of the following program has checked it extensively, neither the author nor the publisher makes any warranty, either implied or expressed, for the correct performance of the program. C
C
803 9095 813 C 2 9000 9001 loo0
GENERAL RAY TRACING PROGRAM FOR FORTRAN IV INTEGER SSl COMMON IRAY,IMl,XA,RPAR1,YA,RPAR2,ZA,RPAR3,XOBAR,CKBAR,YOBAR. *CLBAR.ZOBAR,CMBAR,XO,YO,ZO,ALPHA,BETA,GAMMA,A.B,C,D.E,F.G,REFR, *APl.AP2.AP3,AP4,APS.AP6,NAP,NOUT DIMENSION XO(35).YO(35) DIMENSION Z0(35),ALPHA(35) DIMENSION BETA(35),GAMMA(35) DIMENSION A(35),REFR(36) DIMENSION AP1(35),AP2(35) DIMENSION AP3(35),AP4(35) DIMENSION AP5(35),AP6(35) DIMENSION NAP(35),NOUT(35) DIMENSION CKBAR(35),CLBAR(35) DIMENSION CMBAR(35),XOBAR(35) DIMENSION YOBAR(35),ZOBAR(35) DIMENSION B(35),C(35) DIMENSION D(35),E(35) DIMENSION F(35).G(35) DIMENSION TITLE(12) DATA CONV/O.O1745329/,IPG/1/,LINE/O/ DATA IPROS/O/ DEFINE THE TANGENT TAN(A) = SIN(AlICOS(A) READ (5.9095) SSl ,IRS,IPROS FORMAT(311) IPG= 1 LINE=0 KKK=2 READ TITLE CARD READ (5.2) TITLE FORMAT( 12A6) WRITE (6,900O)TITLE,IPG FORMAT(lH1,25X,12A6,25X,5HPAGE 13) WRITE (6.9001) FORMAT(lH0.53X,lOHINPUT DATA//) IF (SSl) 1000,1001,1000 WRITE (6,9002)
64 9002 1001 9003 1002 9021 1013 9022 1020 C 102 C C 9004
C 103 C C 113 C C
C C C
C C
9005
9004
1003 9007
1004 9008
1005 9009
1006 9010 1007 9011
WAYNE WIITANEN FORMAT(43X,31HSWITCH SET FOR SINGLE RAY INPUT ) GO TO 1002 WRITE (6,9003) FORMAT(43X,32HSWITCH SET FOR FAN OF RAYS INPUT ) IF (IRS. NE. 1) GO TO 1013 WRITE (6,9021) FORMAT(43X,27HROTATIONAL SYMMETRY ASSUMED ) IF (IPROS. NE. 1) GO TO 1020 WRITE (6.9022) FORMAT(43X,69HNO OUTPUT OF RAYS WILL BE DONE, BUT SUBROUTINE PROCE *S WILL BE CALLED. ) CONTINUE JB=50 READ HEADER CARD FORMAT (F11.5.13.El5.8) READ (5,102)REFR(l),NOSUR,TOL REFR = INITIAL INDEX OF REFR. NOSUR = NO. OF SURFACES TOL = MIN. LIMIT ON ITERATION INCREMENT WRITE (6,9004) REFR( l),NOSUR,TOL FORMAT(43X.30HINITIAL INDEX OF REFRACTION = F11.5/43X, *21HNUMBER OF SURFACES = 12/43X.24HCONVERGENCE TOLERANCE = El2.5/,‘) NOSUR = NOSUR + 1 NOSR=O DO 203 I=Z,NOSUR READ SURFACE DATA FORMAT (F11.5,F11.5,F11.5,F11.5,F11.5,F11.5,11,11,13) READ (5,lO3)XO(I),YO(I,,Zo(I),ALPHA(I).BETA(I).GAMMA(I) X0, YO, AND ZO ARE LOCAL SYSTEM SURFACE VERTEX COORDINATES. ALPHA, BETA, AND GAMMA ARE EULER ROTATIONAL ANGLES. READ (5,113) A(I),B(I),C(I),D(I),E(I),F(I).G(I).REFR(I) FORMAT(6F11.5) A, B. C, D; E. F, AND G ARE SURFACE COEFFICIENTS REFR = INDEX OF REFRACTION FOLLOWING SURFACE READ (5,103) APl(I),AP2(I,,AP3(I),AP4(I).AP5(I),AP6(I),NAP(I), *NOUT(I) API, AP2, AP3, AP4. AP5, AND AP6 SPECIFY APERTURE DIMENSIONS, SEE INPUT WRITE UP. NAP = 0 FOR CIRCULAR-ANNULAR APERTURE, = 1 FOR RECT-TRAPEZOIDAL APERTURE, = 2 FOR HYPERBOLIC APERTURE. NOUT = 0 FOR NO OUTPUT AT SURFACE, = 1 OTHERWISE. NSR=I-1 WRITE (6,905) NSR,XO(I).YO(I),ZO(I),ALPHA(I),BETA(I),GAMMA(I) FORMAT(lX.RHSURFACE 12.5X.29HLOCAL COORDINATES (X,Y.Z) = (FI 1.5. *lH.FI 1.5,1H.F11.5,1H)/16X,25HEULER ANGLES (A.B.C.) = ( F11.5, *lH,Fll.5.1H,Fl 1.5.IH)) WRITE (6.9006) A(I),B(I),C(I),D(I),E(I),F(I),G(I),REFR(I) FORMAT(I6X.33HSURFACE EQUATION COEFFICIENTS, A= E15.8,3H.B=E15.8, *3H,C=E15.8,3H,D=El5.8/47X,2HE=El5.8,3H,F=E15.8.3H.G=El5.8/ *16X,22HINDEX OF REFRACTION = F11.5) IF (NAP(I)1) 1003,1004,1005 WRITE (6.9007) AP3(I),AP4(I),APl(I),AP2(I),ZO(I) FORMAT(16X,39HAPERTURE IS AN ANNULUS. INNER RADIUS = F11.5, *16H OUTER RADIUS = F11,5/16X,35HCOORDINATES OF APERTURE (X,Y,Z) = *(F11.5,lH,Fl1.5,1H,F11.5.1H)) GO TO 1006 WRITE (6.9008) APl(I),AP3(I),AP2(I),AP4(I).AP5(I).AP6(1) FORMAT(l6X,37HAPERTURE IS A RECTANGLE OR TRAPEZOID./l6X,24H(X,Y)-L *OWER BOUNDS ARE (F11.5.1H,F11.5,1H),2X,24H(X,Y)-UPPER BOUNDS ARE ( *F11.5,1H.F11.5,1H)/l6X,26HRECIPROCAL SLOPE OF LHS = F11,5,7H RHS = *F11.5) GO TO 1006 WRITE (6.9009) APl(I).AP2(I),AP3(I),AP4(I),AP5(I),AP6(1) FORMAT(16X,47HAPERTURE IS A HYPERBOLA, CENTER IS AT (X.Y) = ( *Fll.S,lH,Fl 1.5,1H)/16X.l8HSEMI-MAJOR AXIS = Fl 1.5,2X.l8HSEMI-MINOR *AXIS = F11.5.2X,30H(Y -LOWER,Y -UPPER) BOUNDS ARE (F11.5.1H.Fll.5, *lH)) IF (NOUT(1) .NE. 0 .AND. IPROS .NE. 1) GO TO 1007 WRITE (6.9010) FORMAT(16X,25HNO OUTPUT AT THIS SURFACE //) GO TO 1009 WRITE (6.9011) FORMAT(16X.22HOUTPUT AT THIS SURFACE //) NOSR = NOSR + 1
Physiological optics 1009
203 C 2050 9012 C 205 C C C C C C C C
9013
1008 9014 1010 9015
1011 9016
201 C C C C 116
9017
1014 1015 9019
IF (MOD&S) .NE. 0) GO TO 203 IPG=IPG+ 1 WRITE (6,900O) TITLEJPG WRITE (6,100) CONTINUE READ RAY DATA WRITE (6,9012) FORMAT( lH0,20HRAY DESCRIPTION DATA //) SWITCH 1 ON TO USE SINGLE RAY CARD INPUT IF (SSl) 205,201,205 READ (5,103)XA,YA,ZA,RPARl,RPAR2,RPAR3.NIN.NAXIN.IRAY AN IRAY OF ZERO OR LESS SIGNIFIES END OF SINGLE RAY JOB IF (IRAY .LE. 0) GO TO 803 XA, YA, AND ZA ARE RAY COORDINATES AT FIRST POINT RPARl, RPAR2, RPAR3 ARE RAY COORDINATES AT SECOND POINT OR RAY DIRECTION COSINES AT FIRST POINT NIN = 0 IF RAYS ARE SPEC. BY 2 POINTS NIN = 1 IF RAYS ARE SPEC. BY 1 POINT AND DIR. COSINES, NAXIN = 1 FOR COMPUTATION OF INTERSECTION OF RAY WITH OPTICAL AXIS, = 0 OTHERWISE IRAY = RAY NUMBER. IF (NAXIN .EQ. 1) GO TO 1008 WRITu6.9013) IRAY FORMAT(lOX.1 1HRAY NUMBER 13,5X,56H NO COMPUTATION OF INTERSECTION *OF RAY WITH OPTICAL AXIS ) GO TO 1010 WRITE(6,9014) IRAY FORMAT(lOX,llHRAY NUMBER 13,5X,71 H COMPUTATION OF INTERSECTION OF *RAY WITH OPTICAL AXIS WILL BE PERFORMED ) IF (NIN .EQ. I) GO TO 1011 WRITE (6,901s) XA.YA,ZA,RPARl.RPARZ.RPAR3 FORMAT(24X,38H(X,Y,Z) COORDINATES OF FIRST POINT = (F11.5,lH.Fll.5 *.l H.FI I .5.1H)/24X.39H(X.Y.Z) COORDINATES OF SECOND POINT = (F11.5, *IH.FI 1.5.1H.Fl 1.5.IH)) LINE=0 IPG=IPG+l WRITE (6,900O) TITLEJPG WRITE (6,l) GO TO 206 WRITE(6.9016) XA,YA,ZA,RPARl_RPAR2,RPAR3 FORMAT(24X,38H(X,Y,Z) COORDINATES OF FIRST POINT = (F11.5,lH,F11.5 *,lH,F11.5,lH)/24X,3lH(X.Y,Z) DIRECTIONAL COSINES = (F11.5,lH,Fll.5 *,lH,Fll.5,lH)) LINE=0 IPG=IPG+ 1 WRITE (6,900O) TITLEJPG WRITE (6.1) GO TO 206 READ (5,116)FXUl,FZl,XGAPl,FXU2,FZ2,XGAP2 READ (5,116)FYU2,YGAP2 RAY INPUT DATA. FXUl = MAX.XR, FZl = ZR, FXU2 = MAX. X AT POINT 2 XGAPl = X SPACING AT POINT 1. XGAP2 = X SPACING AT POINT 2 FZ2 = Z COORDINATE AT POINT 2 FYU2 = MAX Y COORDINATE AT 2ND POINT, YGAP2 = Y SPACING AT POINT 2 FORMAT (F11.5,F11.5,F11.5,F11.5,F11.5,F11.5) TEMP = - FXU 1 TEMP2= - FXU2 TEMP3 = - FYU2 WRITE(6,9017) FXUl,XGAPl,TEMP,FZl,FXU2,XGAP2,TEMP2,FZ2,FZl,FYUl. *YGAP2,TEMP3,FZ2 FORMAT(lOX,34HFAN OF RAYS... IST X-COORDINATE = Fll.5,1H(Fll.5,lH) *Fl1.5,5X.7HAT Z = Fll5/25X,l9H2ND X-COORDINATE = F11.5,1H(F11.5, *lH)Ft 1,5,5X,7HAT Z = Fll5/25X,24HlST Y-COORDINATE = 0(0)0,5X,7HAT * Z = Fll5/25X,l9H2ND Y-COORDINATE = Fll.5,1H(Fll.5,lH)Fll.5, *SX,7HAT Z = F11.5) IF (IRS .NE. 1) GO TO 1014 NRAYS=(2*IFIX(FXUl/‘XGAPl)+ l)*(IFIX(FXU2/XGAP2)+ l)*(IFIX(FYU2/ *YGAP2)+ 1) GO TO 1015 NRAYS=(2*IFIX(FXUl/XGAPl)+ 1)*(2*IFIX(FXU2/XGAP2)+1)*(2*IFIX( *FYU2/YGAP2)+ 1) WRITE (6,9019) NRAYS FORMAT(lOXJ7HNUMBER OF RAYS = 15) IF (IPROS .EQ. I) GO TO 1012
65
66
9020
9018
1012
812
8120 8030
814
4 804
8050 805
206 208
207
209
WAYNE WIITANEN NOSR=3*NOSR*NRAYS+?*NRAYS WRITE (6,902O) NOSR FORMAT(lOX,36HMAXIMUM NUMBER OF LINES OF OUTPUT = 16) IF (FXUl .GE. 0.0 .AND. XGAPl .GT. 0.0 .AND. FXU2 .GE. 0.0 .AND. * XGAPZ .GT. 0.0 .AND. FYU2 .GE. 0.0 .AND. YGAP? .GT. 0.0) GO TO * 1012 WRITE (6,9018) FORMAT(69H X OR Y MAXIMUM VALUE(S) AND ITS (THEIR) INCREMENTS ARE *INCONSISTENT. / 64H X AND Y MAXIMA MUST BE GE. 0.0 AND INCREMENTS * MUST BE .GT. 0.0 j ’ GO TO 803 IPG=IPG+l WRITE (6,900O) TITLE,IPG WRITE (6,l) IRAY =0 GAP3=0. RPAR?=FYU2 RPARZ=RPARZ-GAP3 IF (IRS .EQ. 0) GO TO 8120 IF fRPAR2) 8030.814.814 IF (RPARZ’+FYU~) k030,814.814 IF (IPROS .NE. 1) GO TO 803 ICODE = 3 CALL PROCES(ICODE,I) GO TO 803 GAP3 = YGAP? GAP1 =O. XA=FXUl XA=XA-GAPI IF (XA+FXU1)812.804,804 GAP1 = XGAPl GAPZ=O. RPARl = FXU2 RPARl = RPAR 1- GAP2 IF (IRS .EQ. 0) GO TO 8050 IF (RPARl) 4,805,805 IF (RPARl + FXUZ) 4,805,805 GAP2= XGAPZ YA = 0.0 ZA=FZl RPAR3 = FZ2 NIN=O NAXIN = 1 IRAY = IRAY + 1 NIN=NIN+ 1 GO TO (207,208),NIN CKBAR( 1) = RPAR 1 CLBAR(i)=RPARZ CMBAR(l)=RPAR3 GO TO 209 XD=RPARl -XA YD=RPARZ-YA ZD=RPAR3-ZA RALEN=SQRT(XD*XD+YD*YD+ZD*ZD) CKBAR( 1) = XD/RALEN CLBAR(l)=YD/RALEN CMBARll 1= ZD/RALEN XOBAR(i)-XA ’ YOBAR(l)=YA ZOBAR(I)=ZA SI = 1.0 DO 15 I=Z,NOSUR IMl=I-1 SA = SIN(CONV*ALPHA(I)) CA = COS(CONV*ALPHA(I)) SB=SIN(CONV*BETA(I)) CB = COS(CONV*BETA(I)) SG = SIN( CONV*GAMMA(I)) CG = COS(CONV*GAMMA(I)) Xl =XOBAR(IMl)-X0(1) Yl =YOBAR(IMl)-YO(1) Zl =ZOBAR(IMl)-ZO(1) Rll=CA*CG+SA*SB*SG
Physiological optics
311 210
302 211 5
213
8
217 10
1021 104
11 218
219 220 9 301 300 221
1016 105
R12= -CB*SG R13=CA*SB*SG_SA*CG R21 =CA*SG-SA*SB*CG R22= CB*CG R23 = -(SA*SG+ CA*SB*CG) R31= SA*CB R32=SB R33=CA*CB XO=Xl*Rll+YI*R12+ZI*R13 YO=Xl*R21+Yl*R22+Zl*R23 ZO=Xl*R31+YI*R32+Zl*R33 CK=CKBAR(IMl)*Rll +CLBAR(IMl)*R12+CMBAR(IMl)*R13 CL=CKBAR(IMl)*R21 +CLBAR(IMl)*R22+ CMBAR(IMl)*R23 CM = CKBAR(IMl)*R31+ CLBAR(IMl)*R32+ CMBAR(IMl)*R33 J=O IF (CM)210.302.210 X0 = X0 - CK*ZO/CM YO=YO-CL*ZOICM z2=0. GO TO 211 z2=20 s=o. J=J+I X=XO+CK*S Y =YO+CL*S Z=Z2+CM*S FCN=A(I)*X*X+ B(I)*Y*Y +C(I)*Z*Z+ D(T)*X+ E(I)*Y + F(I)*Z -G(I) FX = 2.0*4(1)*X + D(1) FY = 2.0*B(I)*Y + E(1) FZ = 2.O*C(Ij*Z + F(1) DETMT=CK*FX+CL*FY+CM*FZ AMAGN = SQRT(FX*FX+ FY*FY + FZ*FZ) IF (DETMT)218,217,218 IF (FCN) 10.9.10 IF (IPROS .NE. I) GO TO 1021 ICODE= CALL PROCES(ICODE,I) GOT0 II WRITE (6.104) IRAY,IMl,XA,RPARl,YA,RPAR2,ZA,RPAR3 FORMAT(lX,I4,13,13H NO INCIDENCE 2X,8HX OR K 2(Fl1.5,1X)/ *23X,8HY OR L 2(Fl lS,lX)/23X,8HZ OR M 2(Fll,S,lX)//) LINE=LINE+6 IF (LINE .LT. S1) GO TO 11 IPG=IPG+ 1 LINE=0 WRITE (6.9ooO) TITLEIPG WRITE (6.1) CONTINUE IF (SSI) 205.3,205 DELS = - FCNiDETMT DELSZ = DELS*DELS TOL2 = TOL*TOL IF (DELSZ-TOL2)9.9,219 IF (J-JB)220,10,10 S=S+DELS GOT05 IF (CM) 300,301,300 IF (S)221,222,222 IF (S-ZO/CM)221.222,222 IF (IPROS .NE. 1) GO TO 1016 ICODE = 2 CALL PROCES(ICODE.1) GOT0 11 WRITE (6.105) IRAY,IMl,XA,RPARl,YA,RPAR2,ZA,RPAR3 FORMAT(lX.I4.13,13H VIRTUAL PATH 2X,8HX OR K 2(Fll.S,lX)/ *23X,8HY OR L 2(F11,5,1X)/23X,8HZ OR M 2(F11,5,1X)//) LINE = LINE + 6 IF (LINE .LT. 51) GO TO 11 LINE=0 IPG=IPG+ I WRITE (6.9000) TITLE,IPG WRITE (6.1) GOT0 11
67
68 222 223 12
1017 106
226 224 227 228 229 225 230 231
13
233
232
234
1018 107
235 236 237 14
C C
C
C C
WAKNE WIITANEN KAP = NAP(I) + 1 GO TO (223,224,225),KAP RHSQ =(X-APl(I))*(X-APl(I))+(Y -AP2(I))*(Y -AP2(1)) IF (RHSQ -AP3(I)*AP3(1))12,226,226 IF (IPROS .NE. 1) GO TO 1017 ICODE = 2 CALL PROCES(ICODE,I) GO TO 11 WRITE (6,106) IRAY,IMl,XA,RPARl,YA,RPAR2,ZA,RPAR3 FORMAT(lX,I4,13,14H APERTURE STOP lX,8HX OR K 2(F11.5.1X)/ *23X,8HY OR L 2(F11,5,1X)/23X,8HZ OR M 2(Fll.S,lX)//) LINE=LINE+6 IF (LINE .LT. 51) GO TO 11 LINE=0 IPG=IPG+ 1 WRITE (6,9ooO) TITLEJPG WRITE (6,l) GOT0 11 IF (RHSQ-AP4(I)*AP4(1))13,12,12 IF (Y -AP3(1))12,12,227 IF (Y -AP4(1))228,12,12 IF (X-(APS(I)*(Y-AP3(I))+AP1(1)))12,12,229 IF (X-(APqI)*(Y -AP3(I))+AP2(1)))13,12,12 IF (Y -AP5(1))12,12.230 IF (Y -AP6(1))231,12,12 FFF=(X-APl(I))/AP3(1) FFF2 = FFF*FFF GGG = (Y - AP2(I))/AP4(1) GGG2 = GGG*GGG IF (FFFZ-GGGZ-1.)13,12.12 RAT= REFR(IMl)/REFR(I) ALC=(RAT*DETMT)/(FX*FX+FY*FY+FZ*FZ) IF (REFR(1) + REFR(IMl))232,233,232 GAMUC= 2.*ALC RAT= 1. GO TO 14 BLC=(RAT*RAT-l.)/(FX*FX+FY*FY+FZ*FZ) DISC = ALC*ALC - BLC IF (DISC)234,235,235 IF (IPROS .NE. 1) GO TO 1018 ICODE= CALL PROCES(ICODE,I) GOT0 11 WRITE (6.107) TRAY,IMl,XA.RPARl.YA,RPAR2,ZA,RPAR3 FORMAT(lX,I4,13,14H NO REFRACTION lX,8HX OR K 2(F11.5,1X)/ *23X,8HY OR L 2(F11,5,1X)/23X,8HZ OR M 2(F11.5,1X)//) LINE = LINE + 6 IF (LINE .LT. 5 1) GO TO 11 LINE = 0 IPG=IPG+ 1 WRITE (6,9OGO)TITLEJPG WRITE (6,l) GO TO 11 DISC=SQRT(DISC) IF (ALC)236,10,237 DISC = - DISC GAMUC = DISC - ALC CK=RAT*CK+GAMUC*FX CL = RAT*CL + GAMUC*FY CM =RAT*CM + GAMUC*FZ CALCULATE THE ANGLE BETWEEN THE INCIDENT RAY AND THE NORMA1 UN-NORMALIZED VALUE IN DETMT. NORMALIZE IT FIRST. DETMT = - DETMT/AMAGN CT1 = ABS(DETMT) PHI = ABS(ATAN(SQRT(l,O_DETMT**2)/DETMT)) THE ANGLE OF THE REFRACTED RAY IS DETMT = -(CK*FX+ CL*FY + CM*FZ)/AMAGN CIT = ABS(DETMT) PHIP = ABS(ATAN(SQRT(l,O_DETMT**2)/DETMT)) CHECK FOR NORMAL INCIDENCE (PHI = 0) IF (PHI .NE. 0.0) GO TO 1022 USE SPECIAL FORMULA FOR NORMAL INCIDENCE RF = ((1.0 - RAT)/(l.O + RAT))**2
Physiological optics
C C
C 1023 C
C
16
1019
15
239 241 18
242 244 243 240 246 245
17
1050
100
GO TO 1023 INCIDENCE IS AT AN ANGLE TO THE SURFACE NORMAL CALCULATE REFLECTANCE AND TRANSMITTANCE FOR UNPOLARIZED LIGHT RL = (TAN(PHI_PHIP)/TAN(PHI+PHIP))**2 RR = (STN(PHI- PHIP)/SIN(PHI+ PHIP))**2 RF = (RL + RR)/2.0 TAKE ADVANTAGE OF THE FORMULA RF+TR= 1 TR = 1.0 - RF CALCULATE ENERGY EJR = RF*SI*CTI EJT = TR*SI*CTI CALCULATE INTENSITY AJR = RF*SI AJT = TR*SI*CTI/CTT SI = AJT XOBAR(I)=Rll*X+R2l*Y+R3l*Z+XO(I) YOBAR(I)=Rl2*X+R22*Y +R32*Z+YO(I) ZOBAR(I)=Rl3*X+R23*Y+R33*Z+ZO(I) CKBAR(I)=Rll*CK+R21*CL+R3l*CM CLBAR(I)=R12*CK+R22*CL+R32*CM CMBAR(I)=Rl3*CK+R23*CL+R33*CM KOUT= NOUT(I)+ 1 GO TO (15,16),KOUT IF (IPROS .NE. 1) GO TO 1019 ICODE = 1 CALL PROCES(ICODE,I) GO TO 15 WRITE (6,108) IRAY,IMl,XA,RPARlS(OBAR(I),CKBAR(I),RF,TR WRITE (6,109) YA,RPAR2,YOBAR(I),CLBAR(I),EJR,EJT WRITE (6,110) ZA,RPAR3,ZOBAR(I),CMBAR(I),AJR,AJT LINE = LINE + 3 IF (LINE .LT. 51) GO TO 15 LINE =0 IPG=IPG+ 1 WRITE (6,900O) TITLE,IPG WRITE (6,l) CONTINUE NAXIN = NAXIN+ 1 GO TO (l7,239),NAXIN IF (CKBAR(NOSUR))240,241,240 IF (XOBAR(NOSUR))l8,242,18 IF (IPROS .EQ. 1) GO TO 1050 WRITE (6,111) IRAY LINE = LINE + 1 IF (LINE .LT. 51) GO TO 17 IPd=IPG+ 1 ’ LINE=0 WRITE (6Q000) TITLE,IPG WRITE (6,l) GOT0 17 S = - YOBAR(NOSUR) IF (CLBAR(NOSUR))243,244,243 IF (YOBAR(NOSUR))l8,245,18 S = S/CLBAR(NOSUR) GO TO 245 S= -XOBAR(NOSUR)/CKBAR(NOSUR) IF (CLBAR(NOSUR))246,244,246 IF (S+YOBAR(NOSUR)/CLBAR(NOSUR))18,245,18 AXIN=ZOBAR(NOSUR)+S*CMBAR(NOSUR) IF (IPROS .EQ. 1) GO TO 1050 WRITE(6,112) IRAYAXIN LINE = LINE + 1 WRITE(6,lOO) LINE = LINE+ 1 IF (LINE .LT. 51) GO TO 1050 LINE=0 IPG=IPG+ 1 WRITE (6,9ooO) TITLE,IPG WRITE (6,l) IF (SSl) 2050.3.2050 FORMAT(lHOQH RAY SURF,26X,6HOBJ PT,4X,8H2D PT/DC,SX,6HINT PT,6X, *7HDIR COS,9X,ll HR/J(R)/S(R), 6X,1 1HT/J(T)/Sv)) FORMAT( 1X)
69
70 108 109 110 Ill 112 c
WAYNE WIITANEN FORMAT(lX,I4,13,15X,8HX OR K ,4(Fl 1.5.lX),7X,El2.5, 5XJl2.5 FORMAT(23X,8HY OR L .4(F11.5,lX),7X,El2.5. 5X.El2.5) FORMAT(23X,8HZ OR M ,4(Fll.5,lX),7X.El2.5. 5X,El2.5) FORMAT(lX,I4,3X,l5HNO INT W/OPT AX) FORMAT(lX.I4,3X,23HOPT AX INTERSECTION Z = Fll.5) END DUMMY SUBROUTINE PROCES SUBROUTINE PROCES(A,B) RETURN END
Pergamon
Press
1977. Vol. 7. pp. 35-70.
RAY-TRACING
Prmted
in Great
Britain.
FOR PHYSIOLOGICAL A TUTORIAL
OPTICS-
WAYNE WIITANEN Department of Biology, University of Oregon, Eugene, OR 97403, U.S.A. (Received 24 October 1975)
Abstract-This tutorial article will introduce the physiologist to the basic techniques of ray-tracing: mathematical formulation, and associated computer program. The techniques are illustrated on schematic eyes of the cat, man, the blowfly, and the honeybee. Instructions for use of the program and a listing of the program are included in the article. Physiological optics
FORTRAN
Tutorial
Eyes
Ray-tracing
INTRODUCTION Perhaps one of the more significant contributions that modern digital computation can make to physiological optics is the mechanization of ray-tracing methods. Because of the high degree of precision involved in charting ray paths by this method it seems that such methods should be relevant to physiological optics. In the past, the application of ray-tracing has been limited because of the difficulty in using the method; but, this objection no longer holds. Experience with ray-tracing used in connection with an investigation of the optics of the honey bee’s eye suggests that the method should be of value to those interested in vision and physiological optics. The general plan of the paper is tutorial. It assumes an acquaintance on the part of the reader with elementary vector analysis at a level comparable to that found in any standard introductory text, such as Lass [l]. The mathematical exposition develops in detail a small part of a more elaborate paper on general ray-tracing methods written by Spencer [2]. Unfortunately this paper is not widely known. A basic knowledge of FORTRAN IV format statements will make understanding the card formats described in the sections concerned with using the ray-tracing program much easier. The FORTRAN IV program itself has been completely rewritten from the IBM version (written by D. H. O’Herren of Goodyear Aerospace), with the exception of the iterative procedure for following rays through refracting surfaces (which is based on the work of Spencer mentioned above). There are, then, three goals for this paper: (1) An exposition of ray-tracing methods, (2) instruction in the use of the computer program, (3) to expose visual physiologists to an available technique at a level of detail sufficient for subsequent use. 1. GEOMETRICAL
RAY-TRACING
Fundamental to the understanding of ray-tracing is Snell’s Law, which describes how a ray of light is bent as it passes from one medium into another. Either from experimental work or from the solution of Maxwell’s equations, subject to the proper boundary conditions, Snell’s Law may be found as n, sin 8, = n2 sin 8,, where n, is the index of refraction of the first medium, nz is the index of refraction of the second medium (that into which light is passing), 8, is the angle that the incident ray (in the medium described by n,) makes with the outward normal to the surface separating the two media, and e2 is the angle that the refracted ray makes with the inward normal to the same surface. The incident ray, the refracted ray, and the normal to the refracting surface all lie in the same plane, i.e. are coplanar.
WAYNE WIITANEN
36
A vector formulation of Snell’s Law is desirable for ray tracing. Let u1 be a unit vector representing the incident ray and u2 be a unit vector representing the refracted ray. The components of these vectors are the directional cosines of the rays. Denote the unit vector normal to the surface by n and let it point into the medium into which light is refracted (the inward normal). Because ur, uz, and n all lie in the same plane advantage can be taken of the co-planar relationship between the vectors, and write u2 = ati1 + bn, where a and b are numbers yet to be determined. of this equation with the normal vector, n, giving u2An
To this end take the cross product
= (au1 + bn)/jn = a(u, An) +
b(n/\n)
= a(uI An). Recalling that,lu,l = lull = InI = 1 and that Ip/jq( = IpI.lql sin(p,q), it is clear that the equation above says that sin (u,, n) = a sin (ul, n) or
a sin 8, = sin 19~,
since e1 is the angle between u1 and n, and Q2 the angle between u2 and n. By comparing this last equation with Snell’s Law it follows that a = n1/n2. To determine b another relationship between the vectors ul, u2, and n is needed. One such relationship is that they are all unit vectors (that is, jull = ju2( = InI = 1). Forming the dot product of u2 = au1 + bn with itself we have that 1 = u2*u2 = (au1 + bn)*(aul + bn) = a2uI.u1 + 2abu,*n + b2n.n = a2 + 2abu,.n
+ b’.
This quadratic equation is easily solved for b, b = --au1 .n f [u2(ul *n)2 - u2 + l]“‘. The sign of the radical is yet to be determined. This is most easily done by considering the case where there is no refraction (nl = nz): a = 1, and u1 = u2, Under these conditions b must be equal to zero, implying that b = -au,*n
+ [u2(uI.n)2]“2.
Finally, the vector form of Snell’s Law becomes: u2 =
nlul/n2 - [nluI*n/n2 - (n~(u,*n)2/n~ - n:/n$ + 1)“2]n.
When a > 1 and u,*n is small it is possible to have u2(u,.n)2 - u2 + 1 < 0, corresponding to the condition for total internal reflection. A more symmetrical form of Snell’s Law may be obtained by taking the cross product of n with u2 = (nl/n2)u1 + bn which is n2(u2 A 4
=
nlh
An).
At each refracting surface the incident electromagnetic radiation of a light wave is split into two components: a reflected wave and a transmitted wave. From the Fresnel formulae [3] the amplitudes of the transmitted and reflected waves may be calculated. However, this cakzulation requires that the angIes of incidence and refraction be known. Let 8i designate the angle that the incident ray makes with the normal, and that which
Physiologicaloptics the transmitted we have that
37
ray makes with the normal as et. Since ur represents the incident ray COSBi = u,-n
or 8i = arccos
(II1an).
u2 = au, +
bn.
For the transmitted ray it follows that
Taking the dot product ofboth
sides of the equation with u2 we have 1 = a(ur *u2) + = u(u1~u2) +
&n-u,)
bcos (n,u,),
or 19~= arccos [l - a(uI .u2)]/b. By considering how the energy of the electromagnetic field of incident light energy is distributed between the reflected and transmitted fields for a given angle of polarization it is possible to define the reflectivity and transmissivity of the refracting surface by means of the equations R = RI, cos2 CI+ RL sin2 LX (a = angle of polarization)
where R 11= tan’ (Oi - @,)/tan’ (Bi + 0,) RI = sin2 (ei - @/sin’ (ei + 6,)
and T = T,, cos2 a + Tl sin’ CI (a = angle of polarization)
where T,, = [sin 2&.sin 28,]/[sin2 (Oi + t9,).cos2 (6$ - O,)] TL = sin 281. sin 26,)/sin2 (Oi + 0,).
Natural light is unpolarized and the corresponding reflectivity, R, and transmissivity, T, may be defined by averaging over all angles of polarization, a. Since the average of sin* a or cos’ a is 3 it follows that R = &(R, + RII) T = +(T, + T,,)
and R+T=l. To obtain the corresponding energies from the reflectivity and transmissivity it is necessary to know the energy of the incident light. Call this energy JCi);then we have that J(r) = RJ(i)
where J(*) is the average energy of the reflected ray, and J(t) the average energy of the transmitted wave. From these formulae it is possible to determine the energy of the light rays being transmitted through successive refracting surfaces by setting jCi) equal to Jft) after the refracted ray leaves the refracting surface. The intensity of the reflected and transmitted rays may be calculated from the formulae S”’ = j”‘/cos
(g.
S(i) = JWpos
e. = Rs(i)
I
1
s(t) = J(ycos 8, =
m(i)
cos BJCOS8 t
WAYNE WIITANEN
38
where $I’ and 3”’ are the average respectively. 2. SPECIFICATION TRACING
intensities of the reflected and transmitted
OF REFRACTING RAYS THROUGH
SURFACES THEM
rays,
AND
To specify the rays and the geometry of the refracting surfaces two coordinate systems will be used, both rectangular. The first system will be called a global (or reference) system and it is in this system that the rays originate and their paths are traced. Components of a vector specified in this system, or of a point in the system, will be indicated by having bars over them (this notation and presentation follows that of Spencer [2] so that the mathematical development can be related directly to the program listing in the Appendix). The second coordinate system is a local one that is used to define the shape of the refracting surfaces. Eventually this local surface will have to be positioned within the global one so that the rays may be traced through it. For the moment it is easier to define surface equations in a system that is independent of the global one. Vectors and points in the local coordinate system (local to the surface being defined) will not have bars over them. A given ray is specified by the coordinates (x0, YO,2,) of a point, say P,. through which it passes and the direction cosines (k, I,%) of that ray. Remember that the bars over a letter indicate that the ray is in the global coordinate system. The origin of the global system is denoted by 0. A surface, S, in the global system is defined by the equation F(X, Y Z) = 0. Because the surface is usually defined in a local coordinate
system by the equation
F(X, Y, 2) = 0, it is necessary to relate the local coordinate system to the global one. The steps involved in the ray-tracing problem may, then, be separated as follows: 1. Transform the ray origin (x0, Y,, Z,) and directional cosines (k, 1,&) to their corresponding values in the local coordinate system (X, Y, Z). 2. Find the point of intersection of the ray with the surface, S, as defined by the equation F(X, Y,Z) = 0.
(B)
Fig. 1. Euler angles.
A-rotation about Y-axis (Euler angle CL).B--rotation angle p). C-rotation about Z-axis (Euler angle 7).
about
X-axis
(Euler
Physiological optics
39
3. Find the change in direction of the ray in the local coordinate system following refraction or reflection at the surface, S. 4. Transform the new ray origin (xb, yb, zb) (point of intersection of the incident ray with the refracting -surface) - and the direction cosines (k’, I’, m’) back into the global coordinate system (X, Y, Z). Steps l-4 are repeated as many times as there are refracting surfaces in the optical system. Each of the steps will now be examined in detail. Step 1. Transformation of rays to the local coordinate system of S. The orientation -- of the - local coordinate system (X, Y, Z) is specified with respect to the global one (X, Y, Z) by means of the Euler angles as shown in Fig. 1. The Euler angles specify the rotation of the local system with respect to the global one by means of three successive rotations: CC,p, and y. a is a rotation about the Y-axis, p is a rotation about the X-axis, and y is a rotation about the Z-axis-all in the counter-clockwise sense. Each rotation may be specified by a 3 x 3 matrix as follows
(sin y about Z-axis)
R, = [ i
;i;
-$ygl
(about x-axis)
(about Y-axis)
Rotation through the three angles c(, /I, and y is the product of the three matrices, R,,, = R,R,R,, and in terms of components this is coscl.cosy
R rPa=
+ sincc.sinfi.siny
coscc.siny - sincr.sinp.cosy [
sin cIcosp
-cos/?*cosy cosp~cosy sin p
-sincr.cosy
+ cosa.sinp*siny
-sincl.siny
- cosa.sinfl-cosy cos c(cos p
1 .
Because matrix multiplication is, in general, non-commutative, it is important to note that when defining the orientation of a local coordinate system with respect to the global one by means of the Euler angles that the local system be positioned first with respect to a, then with respect to /?, and finally with respect to y. An important property possessed by each of the matrices R,,R,, and R, is that they are unitary, that is, that the product of the matrix with its transpose yields the identity matrix. This property is easily verified for the three R matrices. The product of two unit matrices is again a unit matrix. All of this means that the inverse of the unitary matrix R is the transpose of that matrix. This property will be used --- just a little later. To translate the ray origin in the (X, Y, Z), or global, system to the local system in which the surface equation has been formulated, it is necessary to know the position of the origin of the local system (X, Y,Z) with respect to the global one (the first card of group 2 does this for the computer program). If the coordinates of 0 (origin
WAYNE WIITANEN
40
of the local system) with respect to the global system are (x0, To, ZO), then the ray origin in the local system will be (X,, Y,, 2,) according to the transformation equation
where (X0, Y,, Z,) is the ray origin in the global system, and the direction cosines become
il [I R
k
1 = R,,,
i
n?
111
Step 2. Point of intersection of the ray withthe surface S. Recall from vector analysis that the vector equation of a line is given by v = p + au. Using this formula, a given ray may be expressed in parametric form as u1 = PO + su2 or, in terms of components, X = X,, + sk Y=
Y,+sl
z = zo + sm.
where s is a parameter of distance along the ray measured from the transformed ray origin PO(XO, YO,Z,). The next problem is to determine the point of intersection of the line representing the incident ray with the refracting surface, S. The scheme employed is to determine the intersection point of the ray with the tangent plane Z = 0, and then to iterate from that point down to the surface S, thereby obtaining the intersection point of the ray with the surface, S. Eventually, an iterative procedure is found allowing determination of the value of the parameter s that defines the point of intersection of the ray with the refracting surface, S. To start out, the parametric equations of the ray are solved for s when Z = 0. This equation is s0 = -Z,/m. By substituting the value so into the parametric equations for the ray, the coordinates of the point of intersection of the ray with the plane Z = 0 are obtained: X1 =X,,+kso Y, = Yo + Is, z,
= 0.
Redefine s so that it is measured from the point (X,, YI, 0), and rewrite the parametric equations for the ray as X = X, + sk Y = Y, + sl Z = sm.
Physiological
41
optics
Next find a value of s such that the resulting values (X, Y,2) from the above equation will satisfy the surface equation F(X, Y, Z) = 0 (describing the refracting surface, S). Solution by purely algebraic means is not practical because of the complexity of the surface equations that may be used. Instead, a Newton-Raphson technique is employed. The Newton-Raphson process is an iterative method for locating the roots of the equation f(-u) = 0 [4]. Suppose that an estimate of a root of the equation is known. Call the estimate a. A correction to the estimate, h, is desired, which when applied to a, will yield a better estimate of the root: ni + 1 = Ui + hi, where hi is the i-th correction to the i-th estimate, ai, of the true root. a. The equation them becomes f(a + h) = 0. To obtain a formula that will give successive correction terms, hi, the last formula is expanded by means of a Taylor’s series about h: f(a
Becausef(a
+ h) -,f(u)
+ hi’(u) + (h2/2)f”(u
0 < 8 < 1.
+ dh),
+ h) = 0 it follows that f(u) + hf’(u) + (h2/2)f“(u
+ #h) A 0.
If h is small the term containing h2 may be neglected, which gives the simple formula f(a) + 1$‘(a) = 0, or h = -f(u)/j-‘(a). Casting this formula into an iterative form for “corrected”
roots gives
Ui+ 1 = Ui + hi = Ui -f(Ui)/f'(Ui),
This formula indicates that the derivative of the equation must exist in order for solutions to exist in the neighborhood of the true root, a. It is reasonable to ask about the error that is introduced into the process by neglecting the h2 term in the complete expansion. This problem may be investigated by subtracting the abbreviated equation from the complete one, which gives h - hi = -h2f”(ui
+ 812)/2f’(ui),
0 < e < 1.
Because h is the true correction, and hi the approximation, it is clear that h - hi is the error in the approximation to h. For convenience the i subscripts will be dropped and F will be written for the maximum value of f”(u) in the neighborhood of an hi. This gives a somewhat simpler appearing formula: E(lz) = 11- hi = -h2F/2f’(u). Clearing of fractions, transposing, and solving the resulting quadratic for h the equation above yields h = l/F[ -f’(u)
+f’(u)(l
+ 2Fhilf’(u))“2].
Expansion of the last term by means of the binomial theorem gives !Z= hi - Fhf/2f’(u)
+ F2hf/2(f’(a))2
+ ...
as the formula describing the error involved in the estimation of h when the Taylor’s series is truncated after the second term. Again. because hi is usually small, the second term (third degree term) in the binomial expansion may be deleted at no great loss of accuracy. This leads to a general formula for the error in the i-th correction to the root: Ei(hi) I (F. h2/2f”(Ui_ 1) 1.
The next question that might be asked is: will the method converge to a solution? Essentially this means that the sequence of errors between the true root and the succes-
WAYNE WIITANEN
42
sive approximations should converge to zero or to some preassigned (convergence tolerance). Let the equation in question be written as
small number
.x =f(x) and let x0 be the first approximation
to the true root, 2. Then write Xl
as the next approximation equations yields
=f(xo)
to the true root. Forming x - x1 = f (x) -f
the difference between the two
(x0)
as the next approximation to the true root, x. By the mean value theorem of the differential calculus, the right-hand side of the above equation can be written as f(x)
- f (x0) = (x - x(J)f ‘(io),
x < 20 < x0
and the entire equation becomes x - x1 = (x - xo)f’(&) and similarly for succeeding approximations. By writing down the n equations by which the n differences x - x, are defined and multiplying them together followed by division by the common factors x - x1, x - x2,. . ., x - x, one gets x -x,
= (x - xo)f’(~o)f’(~~)...f’(2,-l).
If the maximum absolute value off’(x) is less than 1 throughout then Ix - x,1 I Ix - Xoj’F”
the interval (x,x0)
where F is the maximum value of all of thef’(x,): F = Imax(f’(20),f’(il)
,..., f’(i,_l))I
< 1.
The right-hand side of the inequality will approach zero as n becomes large, which means that the error, x - x,, may be made as small as desired by repetition of the iteration process a sufficient number of times. This means that the condition for convergence is that If’(x)1 < 1 in the neighborhood of the desired root. Applying this iterative method to the Newton-Raphson technique it is apparent that the requirement for convergence is If’(%l)I < 1.
From the definition of the Newton-Raphson f ‘(4
process it is clear that
A f (%)f “(4Mf
‘(%))Z,
and the convergence criterion becomes If(%)f
“(a”)/(f’(a”))2
I< 1
or, regrouped, If(4f
“(4) I s (f ‘(%))2.
The exact conditions for convergence and the estimate of the error will depend upon the form of the equation f(x) = 0, which will vary with the type of surface. The applicability of the Newton-Raphson method might be questioned on the grounds that it seems as if the surface equation F(X, Y, Z) = 0 is a function of three variables, while the formulae developed above are in terms of a single variable. This conflict is quite easily resolved by noting that successive values of X, Y, and Z are determined
Physiological
by the parameter ofs:
optics
s. This means that the equation
43
F(X, Y, Z) = 0 is really a function
F(X(s), Y(s), Z(s)) = 0. From this argument it is easy to see that the Newton-Raphson be written as Sj+ 1 = Sj - F(Xj, ~, Zj)/F’(Xj,
iteration formula can
rj, Zj)
where Xj=
X1 + kSj
rj = Ye + 1Sj Zj = WlSj
and the derivative F’( ) is
Because F( ) is a function of s, by way of the implicit functions X(s), Y(s), and Z(s), the derivative becomes dF ds=
i3F
dX aF dY dX’ds+E’ds+dZ’ds
i3F dZ
where -_= dX
k
-_= dY
1 ’
ds
ds
and
’
dZ ds =m.
These derivatives are obtained from the defining equations for Xj, Yj, and Zj as given above. Making the indicated substitutions yields dF
Ek
ds = x
(Xj,Yj,Zj)
+5’
+F,m (Xj,YjZj)
(Xj,Yj,Z,)
The initial estimate for s1 may be taken as s1 =
0
and the iterative process is terminated when ISf - q-11 -cc, where E is a small predetermined convergence tolerance. In general, the choice of a starting point at the Z = 0 plane will insure convergence; however, there are circumstances under which the process will either fail to converge or will converge to an incorrect value. If the derivative F’() becomes zero at any time in the process, then s will blow up and the process will diverge. In addition it is possible that a ray will strike a surface at two or more points, resulting in convergence to an incorrect value (wrong point). Neither of these circumstances is often encountered in practice. Recalling how the line describing the ray path was defined it follows that the total path of the ray from (X,, Y,, Z,) to the surface, S, is p = s0 + So. If p < 0 then the
WAYNE WHTANEN
44
path is a virtual one, extending along the line of the ray in a direction opposite to the direction of propagation. Let the intersection of the ray with the surface be indicated by the coordinates (X,, Y,, Z,), as determined by the iteration procedure. From these values it is easy to obtain the direction cosines for the surface normal at the point of intersection:
The procedure so far is applicable to any surface expressible by the equation F(X, Y, Z) = 0, providing a wide. variety of choices. By requiring that the ray-tracing procedure by implemented on a digital computer, the number of choices for the form of the surface equation becomes somewhat restricted, because symbolic derivative taking facilities (those that are done with paper and pencil by people) are not widely available on digital computers. Most of the refracting surfaces used in optics are obtained by revolving some conic section about a coordinate axis, with the vertex of the surface at the origin. The equation to be used for describing refracting surfaces has been selected to be F(X,Y,Z)=AXZ+BY’+CZ2+DX+EY+FZ-G=O. By judicious selection of coefficients a large number of surfaces may be described. Table 1 lists the surfaces most commonly encountered in practice, and indicates the nature of the coefficients of the general surface equation. Step 3. Refraction of the ray. From the previous discussion it will be recalled that Snell’s Law may be represented conveniently in vector form
where u1 is a unit vector in the direction of the incident ray with components (k, 1,m), u2 is a unit vector in the direction of the refracted ray, and n is a vector normal to the surface, S, at the ray intersection point and has components (K, L, M). n, is the refractive index of the medium in which the ray is incident, and n2 the index of Table 1. Commonly used equations for refracting surfaces. For each of the surface types named the standard equation is given, its center, and the form in which the ray-tracing program requires the coefficients
Surface
UP+
Center
equation
cz*+ Dx+
Plane
ax+Bytcz+o=o
Sphere
x2+ytz2=r2
O
0
I
I
l/b2
(+h, Zk,O)
(fh,+k&nm)
(O,O,O) (+n,?k,x,rm)
Ev+
I
I A
Fz=
G
I 8
C *2/n
-D
tL2h
iZk
(r2-(hz+k~mZ))
I/c2
t2h/a2
+_2k/b2 +2mki
I-(h’/o’+
Ilb2
0
+2h/o’
+2k/b2
(I-h2/a2-k2/b2)
l/b2
0
?2h/a2
+2k/b2
-c
/b’
I/C2f2h/a2
22X/b’
T2mk
.1/b’
0
t2hd
T2k/62
-c
I/b’
d/C2
+2h/a2
?2k/b2
i2m/c
Ellipmd x~/o~+~~/~~+~~/c~= Ellnptic
x2/a2+y2/b Elliptic
I (thI tk
9371
k 2/b%m2/c2
cylinder,Z-am..... =
I
0
paroboloid,Z-axis.. x2/a2+y2/b2-a=0
Hyperboloid
X4o2+yZ/d-r2/C2=I(~h,~,X,+m) Hyperbolic
+cm
of one sheet, Z-axis..
paraboloid,
x2/02-y’/b*-cz=O
I-( h2/8+k2/b2-n?/c’)
X-axis..... (th,tk,im
+cm
‘one, Z-ax,s. h2/a2+ k2/b2-d/
c2
)
Physiological optics
45
refraction of the medium into which the ray is refracted. Recalling that the notation is S and S’ for u1 and u2 it follows that S’ = pS + rn
where P = n&2 and r is an undetermined multiplier. This can be done because the normal to the surface and the two rays are coplanar. The components-of this vector equation are k’ =
pk + I-K
1’ =
pl + TL
m’ =
Squaring and adding the component
pm + TM. equations, above, a quadratic
in r is obtained:
r2 + 2aT+b=O where Q = p(kK + IL + mM)/(K2 + L2 + M2) and b = (p2 - l)/(K’ + L2 + M’). The quadratic equation may be solved by the usual formula. The condition under which it fails to have a real root is b > a2 corresponding to total internal reflection. At a reflecting surface the same equations may be used with p = 1, and the other root of the quadratic. In this case b = 0 and the undetermined multiplier becomes
r = -2a. With this assignment of values the direction cosines (components) become
of the refracted ray
k’ = k - 2aK 1’ = 1 - 2aL m’ = m - 2aM. Step 4. Transformation of the refracted ray back into the global system. The new ray origin coordinates (X, Y, 2) and direction cosines (k’, I’, m’) may be transformed back into the global reference system by using the inverse forms of the transformations used to put the incident ray into the local coordinate system. These are
[
and
9]-R+;;]
46
WAYNE WIITANEN
Because R,,, is a product of unitary matrices it is itself a unitary means that its inverse and transpose are equal. That is, coscc~cosy + sina.sinfl.siny R$
=
cosa.siny
sin a cos /?
- sincc.sinp.cosy cosg~cosy
-cosp*siny -sincr.cosy
matrix and this
+ cosccsinpsiny 3. FOUR
-sina.cosy
RAY-TRACING
sin p COSG!.cosp I
- cosa.sinp.cosy EXAMPLES
To illustrate the power of the ray-tracing method, I have selected as examples the Gullstrand schematic human eye (Duke-Elder [S]), the schematic cat’s eye (Bishop and Vakkur [6]), the compound eye of the fly (Seitz [7]), and the eye of the honey bee (Varela and Wiitanen [S]). Example 1: The schematic human eye (Duke-Elder [S])
In Figs. 2(a) and 2(b) the geometry of the schematic eye, refractive indices placement of apertures, and sample ray paths within the system are shown. The rays traced are Gullstrands
schematic
human eye
I cm
I
/
(a)
(b)
Fig. are the All
2. Gullstrand’s schematic human eye. (ablements of the dioptric apparatus. Two apertures shown: one at the first vertex having radius 5.5 mm (corresponding to radius of the cornea), other at 3.6mm behind the first surface of radius 4.5 mm (corresponding to the pupil). measurements of surface positions and radii of curvature in millimeters. (bEselected ray paths on the same scale as 3(a).
Physiological
optics
47
paraxial with respect to the Z-axis (optical axis of the system) and are spaced at 0.5, 1.0, 2.0 and 4.0mm away from it. The forward shift of the focal point is apparent as the rays become more distant from the Z-axis (less paraxial). The focal point of the system, as computed by ray-tracing using paraxial rays, was found to be at 2.365 cm from the first surface. The focal point of the Gullstrand schematic human eye is generally accepted to be at 2.2875 cm. From this example we see that it is possible to select the origin and initial trajectory of rays incident on the optical system. We also note that each of the refracting surfaces is spherical, is specified by the equation for a sphere (using the radius of curvature measurements), and is positioned in the global system by means of the position of the first surface of each element. In all seven surfaces have been used in the description: the first surface of the cornea, the back surface of the cornea, the front surface of the lens core, the back surface of the lens core, the rear surface of the lens capsule, and finally the surface of the retina. In the next section the manner in which these rays and surfaces are indicated to the ray-tracing program will be explained and the necessary cards for each example will be given. Example 2: The schematic eye of the cat (Bishop and Vakkur [6])
Again a group of four rays are shown in Fig. 3(b). They are parallel to the Z-axis and are spaced at 0.5, 1.0, 2,0 and 4.0 mm from the axis. The geometry of the eye cot\ eye
Retina \
(b)
Fig. 3. Schematic cat’s eye. (a&geometry of the dioptric apparatus. Two apertures: cornea dia. 8.4 mm (at first surface), pupil dia. 7.3 mm (at third surface, 5.20 mm behind first surface). All surface positions and radii of curvature in mm. (btselected ray paths on the same scale as 4(a). C.8.M. 7/1-D
48
WAYNE WIITANEN
Plane
-_-----
Envelope of crystalline -
cone
---_ L-_--
of rhabdom
1 -_---_
I
II I
0
7.4
16.3
146.0 (a)
(b)
Fig. 4. Eye of the honey bee. (atgeometry micrometers. (btselected
of the dioptric apparatus. All measurements ray paths on the same scale as S(a).
in
is shown in Fig. 3(a). Again the refracting surfaces are spherical and are positioned in the global coordinate system by specifying the positions of the various surfaces. The anterior progression of the focal point, as the rays become less and less paraxial, is interesting. Also worth noting is that the focal point lies within the retina itself but anterior to the later of primary receptors. The focal point, as calculated by ray-tracing, was found to be 2.214cm behind the first surface of the eye. The position of the focal point as calculated by Bishop and Vakkur [6] was found to lie 2.2153 cm behind the first surface of the eye. The two results are in good agreement. As for the previous example, the control cards necessary to “set up” this system will be described in detail below. Example 3: The compound eye of the honey bee (Apis mellifera)
Figures 4(a) and (b) illustrate the geometry of the honey bee’s eye and some ray trajectories through the eye, respectively. As the rays become less paraxial we again observe the forward shift of the focal point. The calculated focal point of the bee’s eye was found to be at 108.54 pm below the first surface. This example illustrates the use of a non-spherical refracting surface, namely, the planes at 16.3 and 48.0 pm. Essentially these are spheres of infinite radius of curvature, but are nonetheless easy to represent to the ray-tracing program. Details of computer set-up will be given below. Example 4: The compound eye of the blowfly (Calliohora
erythrocephala)” Figures 5(a) and (b) show the geometry and the ray paths within the optical system, respectively. As the rays become less and less paraxial the focal point progresses anteriorly although not evident from Fig. 6. Ray-tracing locates the focal point as calculated * Seitz
[7],
Physiological
optics
49
Y
73.56 II 77.68 (a)
(b) Fig. 5. Compound eye of the blowfly. (awioptric apparatus geometry. All measurements in micrometers. (bkselected ray paths on the same scale as 6(a).
are
by Seitz [7] using the thick lens formula places the focal point at 73.40 pm below the first surface of the equivalent thick lens. The discrepancy is 13.47 pm. Details of the computer set-up will be given below.
4. USING
THE
RAY-TRACING
PROGRAM
General considerations Digital computer programs usually require some sort of control information preceding the actual program to direct the actions of the computer’s supervisory system. Because this information varies considerably from machine to machine, and even within different supervisory systems for the same machine, it will not be discussed. The ray-tracing program itself is written in pure FORTRAN IV, in the hope that it will be executable on any computer that will accept the language. For this reason the program should have a wider applicability than would one that does the same computations but is written in assembly language. Card groups and computational sequence In the most general terms, the ray-tracing program reads a group of cards describing the surfaces and rays, performs a series of calculations using these data (see Section
50
WIITANEN
WAYNE
Start Return to system. _ Run terminated
IRead input data (Card grol;ps 1,2,3)
No data
Print out a copy _
Loop
Rintout II error mesMge
TeZ+E
3~oZZons
Yes
No
Read fan data
Setup single ray
zIrnL;E * Have all Ypoints been calculated? No 1 SetupXincrements c( IncrementXpoints c Have allXpoints& been cal:ulated 7
I, Calculate direction cosines to put ray into standard (unit vector) format 2. Calculate elements ofR-matrix. 3. Transform ray vector to loco I system. 4. Calcu\at$ incident ray and tongent plane intersectlon point.
I
-
No incidence Error message
a
No convergence
-
Aperture stop Errwmessage
e
Aperture stop
-Error
No refraction message -
Internal reflection
Itemte
to surface Z
Apply aperture to ray path OK 1 Test for internal reflection
OK 1 refracted my vector components. 2. Calculate transmittance, reflectance, etc. 3. Tansform refracted ray back into global system 4. Print out data forth is surface if required.
I. Calculate
t Have all surfaces been used?Yes I I. Calculate intersection 2. Print out intersectionof No
Fig. 6. Flowchart
indicating
of ray with optica I oxis. ruy with optical axis.
+ Single ray?
the computational
Yes
procedure
c
used
Read another
single
by the ray-tracing
raycars
program.
l), and finally writes out the desired answers before it goes on to the next set of data cards. The sequence of ray-tracing computations is presented schematically in Fig. 6. The data cards used by the ray-tracing program fall naturally into three groups: 1. Set-up cards 2. Surface definitions 3. Ray specifications. Each (1) card, only
of these groups will be described in detail. Input card group 1. The first group (set-up cards) consists of three cards: a title card, and an “initial conditions” card. The first, the switches the first three columns of the card and is punched as follows: column
1:
a switches card, uses
O-ray specification cards (group 3 cards) will be for a fan of rays (see description under group 3 cards). l-ray specification cards (group 3 cards) will be for a series of individual rays, rmt a fan.
Physiological optics
51
&rotational symmetry will not be assumed (for the advantages of rotational symmetry see the discussion under group 3 cards). l-rotational symmetry will be assumed. &the PROCES subroutine will not be called and normal output column 3: will be printed. l-the PROCES subroutine will be called instead of producing the normal printed output. Note: a dummy PROCES subroutine is required even if it is not to be called (see the Appendix). FORMAT (311)
column 2:
Technical note regarding the user-provided subroutine PROCES
In order to provide for the possibility of performing calculations with the output from the ray-tracing program rather than producing a printed table, entry to a user-provided subroutine by CALL PROCES (I, J) has been provided at each point where the main (ray-tracing) program would produce printed output. The parameters Z and J define the conditions under which the subroutine PROCES was called: z=
1: 2:
good data ray does not intersect with a refracting surface, takes a virtual path, hits an aperture, or is not refracted 3: end of this study has been reached J = 1,. . .) n: the number of the refracting surface at which the current calculations have been made.
The two parameters, Z and J, alone are insufficient to be of any value for further computations, and so some method for passing along other computational values must be provided. A COMMON statement is set up as follows for this purpose: COMMON
IRAY, IMl, XA, RPARl, YA, RPARZ ZA, RPAR3, XOBAR(35), CKBAR(35), YOBAR(35), CLBAR(35), ZOBAR(35), CMBAR(35), X0(35), YO(35), ZO(35), ALPHA(35) BETA(35), GAMMA(35), A(35), B(35), C(35), D(35), E(35), F(35), G(35), REFR(35). APl(35). AP2(35), AP3(35), AP4(35), AP5(35), AP6(35), NAP(35), NOUT(35), CK, CL, CM, S
where IRAY IMl XA RPARl YA RPAR2 ZA RPAR3 XOBAR
= = = = = = = = =
YOBAR = ZOBAR = CKBAR CLBAR CMBAR X0 YO ZO
= = = = = =
the number of the ray being computed surface number (same as J in the calling sequence) 1st X-coordinate of ray 2nd X-coordinate of ray 1st Y-coordinate of ray 2nd Y-coordinate of ray 1st Z-coordinate of ray 2nd Z-coordinate of ray X-coordinate of intersection (ray with surface number J) Y-coordinate of intersection (ray with surface number J) Z-coordinate of intersection (ray with surface number J) X directional cosine of ray leaving surface J Y directional cosine of ray leaving surface J Z directional cosine of ray leaving surface J X-coordinate of the vertex of surface J Y-coordinate of the vertex of surface J Z-coordinate of the vertex of surface J
WAYNE WIITANEN
52
ALPHA = Euler angle of surface J with respect to the reference system BETA = Euler angle of surface J with respect to the reference system GAMMA = Euler angle of surface J with respect to the reference system A, B, C, D, E, F, G = surface equation coefficients of surface J REFR = index of refraction of surface J-1 APl, AP2, AP3, AP4, AP5, AP6 = aperture parameters for surface J NAP = aperture type for surface number J CK = X directional cosine CL = Y directional cosine CM = 2 directional cosine S = ray distance from tangent plane to surface number J.
By providing the subroutine PROCES with a COMMON statement identical to that shown above it is possible to refer to the data used in the main program. The second card of the first group is a title card. This card will cause the title punched in columns l-72 to appear at the top of each page of printed output. Its format is columns l-72: title information FORMAT
(12‘46)
[A = alphabetic]
The last card gives information about the number of surfaces, the convergence tolerance of the Newton-Raphson iteration method (used by the program to find the point of intersection of the incident ray with the refracting surface-see Appendix 1 for details), and the index of refraction of the medium in which the object lies. The format of this card is columns
l-l 1: the 12-14: the 1532: the FORMAT
initial index of refraction (Fl 1.5) number of surfaces in the system (13) convergence criterion (E15.8) (F11.5, 13, E15.8).
Here is a typical set-up for the first card group: col.
1
5
1 0
1 5
2 0
010 THIS IS A SAMPLE TITLE CARD 1.0 004+0.10000000E-05
(switches) (title) (i.c.).
The switches card says that a fan of rays is to be used as input, that rotational symmetry is to be assumed, and that the PROCES subroutine is nut to be called. The title card provides a heading for each page of printed output from the program. The initial conditions card says that the index of refraction of the object space is 1.0 (free space), that there are 4 refracting surfaces in the optical system, and that 6 places of accuracy are required of the iteration procedure. (2) Input curd group 2. The second group of cards is used to define the shape of refracting surfaces. The surface equation used by the ray-tracing program is AX2+BY2+CZ2+DX+EY+FZ-G=O.
By selecting appropriate values for the coefficients, A-G, any quadratic surface may be described. By referring to Table 1 appropriate combinations of coefficients may be selected to define most refracting surfaces. Each surface is defined in its own local (or independent) coordinate system, without any reference to the overall coordinate system in which the rays are to be traced (the global, or reference, system). In order
Physiological optics
53
to give the local coordinate system an orientation with respect to the reference system (that is, the one in which the rays are being traced), the origin of the local system with respect to the reference system must be supplied, as well as any rotation of the local system with respect to the reference system. To this end card group 2 is subdivided into three sections: (a) Parameters for placing the local coordinate system of the refracting surface within the framework of the reference coordinate system; (b) the surface equation (relative to the local system); (c) definition and placement of apertures in the local coordinate system. These cards must be supplied for each surface in the optical system, and in the order given. To place a surface in the reference system a card is made up as follows (each number is punched in the format F11.5): card 1:
columns
1-11: 12-22: 23-33 : 34-44: 45-55 : 5666:
FORMAT
X-coordinate of the local system origin with respect to the reference system origin Y-coordinate of the local system origin with respect to the reference system origin Z-coordinate of the local system origin with respect to the reference system origin c( Euler angle p Euler angle y Euler angle
(6Fll S).
The CI,0, and y Euler angles define the rotation of the local coordinate system with respect to the reference system, and they are illustrated in Fig. 1. Bear in mind that the order in which the rotations are done is important: first the CIrotation, then the p, and finally the y rotation. In columns 73-80 an identification may be punched. I usually punch in the label n-A where n stands for the surface number associated with this coordinate placement. Each coefficient of the surface equation is punched (in the F11.5 format) as follows card 2:
columns
FORMAT (6F 11.5) card 3: columns
FORMAT
l-11: 12-22 : 23-33: 34-44: 45-55 : 56-66:
A B
l-11: 12-22 :
G (F11.5) refractive index of the medium following this surface. (F11.5)
C D E F
(F11.5) (F11.5) (Fl 1.5) (F11.5) (F11.5) (F11.5)
(2F11.5).
The coefficients are usually punched starting at the left end of each field. It is convenient to label these cards n-B and WC in columns 73-80, respectively. Apertures are of three types: circular-annular, rectangular-trapezoidal, and hyperbolic. The parameters for specifying the type and dimensions of an aperture are punched on the last card of this group. All are in the F11.5 format except those specifying the aperture code and the output code. The parameters are: card 4:
columns
l-11:
X-coordinate of the coordinate of the center of a circular-annular aperture, or the coordinate of the center of a hyperbolic aperture, or the X lowerbound of the base of a rectangular-trapezoidal aperture.
WAYNE WIITANEN
12-22 : 23-33 :
3444 :
45-55 :
56-66:
67:
68:
FORMAT
Y-coordinate of the same as above. Inner radius of a circular-annular aperture, or the length of the semi-major axis of a hyperbolic aperture, or the Y-lower bound for a rectangular-trapezoidal aperture. Outer radius of a circular-annular aperture, or the length of the semi-minor axis of a hyperbolic aperture, or the Y-upper bound for a rectangular-trapezoidal aperture. Y-lower bound for a hyperbolic aperture, or the reciprocal slope of the left-hand side of a rectangulartrapezoidal aperture. Enter 0.0 for a circular-annular aperture. Y-upper bound for a hyperbolic aperture, or the reciprocal slope of the right-hand side of a rectangulartrapezoidal aperture. Enter 0.0 for a circular-annular aperture. Aperture code : 0: circular-annular 1: rectangular-trapezoidal 2 : hyperbolic. Output code : 0: no output at this surface 1: output required at this surface
(6F11.5, 211).
The various apertures are illustrated in Fig. 7. I usually label this card n-D in columns 73-80, where yt, as usual, is the number of this surface. This group of four cards is repeated for each refracting surface being defined. (3) Input curd group 3. The last group of cards specifies the rays to be traced through the system. There are two subcategories: a fan of rays, and individual rays. The first of these requires two cards, and their formats are: card 1:
columns
1-11: 12-22: 23-33 : 3444: 45-55 :
Fig. 7. Aperture
types.
The maximum X-coordinate of the ray origin The Z-coordinate of the ray origin The increment to be used to obtain successive rays The maximum X-coordinate of the second point on the ray The Z-coordinate of the second point on the ray
A-rectangular
or trapezoidal. (shaded area passes
B--circular rays).
or annular.
C-hyperbolic
55
Physiological optics
r
Fig. 8. Isometric view illustrating generation of fans of rays (see text).
The 2 increment to be used to obtain successive second point on a ray. l-11: The maximum Y-coordinate value of the second point card 2: columns on a ray The Y increment to be used to obtain successive points 12-22: on the ray. FORMAT (6F11.5/2F11.5) (the symbol “/” means “next card”). 5666 :
I usually label these cards FAN-l and FAN-2 in columns 73-80. Figure 8 illustrates how the rays are defined and how the fans are generated. Rotational symmetry (as may be specified on the switches card of group 1) means that the entire lens system is rotationally symmetric about the optical axis (Z-axis) of the reference system. Advantage can be taken of this fact to cut down the number of rays to be traced through the system, without sacrificing any information. Should rotational symmetry not be requested then the rays to be traced are generated in the following order: 1: X-coordinate of the first ray point--(X,,) (AX) (-X,,,). 2: X-coordinate of the second ray point+Xk,,) (AX’) (-XL,,.). 3 : Y-coordinate of second ray point+Y,,,) (AY) (- Y,,,). Note: Y-coordinate of the first point is (0) (0) (0). The stepping scheme starts with the Y-coordinate (3), then the second X point (2) and finally the first X point (1). This sequence gives a total of rays. When rotational symmetry is used, the X-coordinate of the second ray point runs from XL,, to 0 by decrements of AX’, and the Y-coordinate of the second ray point runs from Ym, to 0 by decrements of AY, thus reducing the number of rays traced by one-fourth (and consequently cost). A trick to confine the rays to the XY-plane is done by setting Ym, to zero and the Y decrement (AY) to any positive, non-zero number. Figure 9 illustrates the difference between a full fan of rays and those generated by invoking rotational symmetry. Single rays may be traced through the system by specifying the origin and direction of each ray. This specification requires one card per ray in the format: card 1:
columns
l-11: 12-22 : 23-33 :
X-coordinate of the ray origin in the reference system. Y-coordinate of the ray origin in the reference system. Z-coordinate of the ray origin in the reference system.
WAYNEWIITANEN
56
Fig. 9. Full fan of rays. Compare with Fig. 8 to see advantage of rotational symmetry in generating rays
3444: 45-55 : 5666: 67:
68:
69-7 1:
FORMAT
X-coordinate of a second point on the ray, or the X directional cosine of the ray leaving its origin. Y-coordinate of a second point on the ray, or the Y directional cosine of the ray leaving its origin. Z-coordinate of a second point on the ray, or the Z directional cosine of the ray leaving its origin. Input code for ray specification type: 0: ray is specified by two points 1: ray is specified by one point and directional cosines Optical axis computation code: 0: do not compute the intersection point of the ray with the optical axis 1: compute the intersection point of the ray with the optical axis. Ray identification number (right justified). If this number is less than or equal to zero then it indicates that there are no more rays to be processed, including the one on this card.
(6F11.5, 211, 13).
Individually defined rays are illustrated in Fig. 10. After all the rays of a fan have been processed, or after the last single ray has been processed, then the program will return to its beginning and attempt to read in another set of data cards (groups 1, 2, and 3). If there are no more cards the program then returns to the supervisory system. If there are more cards then the program will perform the ray-tracing studies indicated by them. Computational errors and the output format During computation situations may arise that are indicative of errors or unusual situations of which the user should be apprized. (1) If a fan of rays is being used and if the step size is incorrectly specified, usually resulting in a loop in the program,
Physiological optics
Fig. 10. Illustration of three independently
51
specified rays (1, 2 and 3).
the user will be notified of his error by the program, and the program will go on to the next study, if any. (2) It is possible that certain rays will not strike a surface. When this happens the program informs the user that a certain ray, given by number. has had no incidence with the surface given. This arises when there is no convergence with the current surface during application of the Newton-Raphson iteration, or if the convergence is oscillatory making it impossible to get within the required convergence tolerance. (3) Sometimes a ray may be reflected back along its original path. Here the program advises the user that a virtual path for the ray has been encountered, and it goes on to the next ray. (4) Apertures provide obstacles for rays within the optical system, and rays that hit them are lost to subsequent processing. When a ray strikes an aperture, the user is informed and the program goes on to the next ray. (5) Internal reflection is not handled by the program, and whenever an internal reflection is encountered the user is informed that there has been no refraction, and the program continues with the next ray. The printed output from the program has a number of prefatory lines which repeat the surface and ray parameters used for the particular study in a readable format so that the values the program used can be checked. In addition, there is some information concerning the maximum number of rays to be traced and the maximum number of lines to be printed out. Sometimes the actual number of lines printed is less than the maximum because of rays blocked by apertures or lost. The bulk of the printed data is organized in a columnar format. The columns, from left to right (see Fig. 11) identify the number of the ray being traced, the surface at which the ray has arrived and at which this output is being taken (after refraction of the ray), the name of the coordinate point or directional cosine being output, the value of the object point, the second object point or directional cosine (these two numbers, the first and second object points, or directional cosines, define the path over which the original ray was travelling prior to encountering the first refracting surface), the point at which the ray intersects the surface (X, Y, Z-coordinates in the reference system), the directional cosines of the ray as it leaves the surface, and the reflectance/transmittance data. For each ray there are three lines of output: one for the X-coordinate (or k directional cosine), one for the Y-coordinate (or 1 directional cosine), and one for the Z-coordinate (or m directional cosine). By examining the object point and the second point/directional
58
R.1
216
WAYNE
SURF 1
u&J
PI 0
2
*.“0GCc
21t
‘
C.cdtCc C. cccc~” 4. ‘)Jctc E. “Jcr 3 c.cmCC 4.c7crc
216
5
216
4.51021 -“.O”“O”
‘"7.58626 6.28155 -C.CCCCO ll".eoo"C 3.4229" -c.cccI)o 123.02132 -0.79568 C.CCCOO 173.5t.co0 -,.IU3L C.OCCCF L77.68"""
c.cococ G. ‘~,“CO 4rcPxwJ C. O’)CC?
C.CC0PC C.CbCCC 0.CJC‘C c. cord 0 ‘63.99725 ~.“5CCC 0.3F.W” C.“i”“3
*.ooocc
217
3
2‘7
211
L.000~3 “.OICL‘ *.c;C“) 0.03‘“‘) C.CJCCO 4.W)C‘O U.CSOL~ b.
4.
217
c*,cc‘1
n?ccc
“3CCU ~.Cmcc ~.CJL‘C C.CxJ“3 (I. 0u**n (r.
217
2.C[i”“”
“.CJ0”” 10c.c00c0 r.c3coo
c .coccc
1"C."""C0 C.Cmx!C J.C"COC ~“.C*cCo ~."0C"" rJ.CDCCC ‘C".C*CC~ *.m"FC C.‘,CCC ‘"C.CPOCC ~.CCCCO C.CCCCC 1CC.CJCOC
+.0000c -U.CCCCC
‘m.3dl51 3.56770 -C.EC~C‘ ‘"7.4"294 3.36363 -0.C"CCC llC.8CC30 .?*6"24C -C.CCCOC 123.13536 -1.CC403 O.FCbCC 173.56""" -1.3099C C.OCO"‘ 117.e~‘CC
4.0708" 0.00000 C.99799 -0.06881 C.0"""" 0.997b3 4.07665 0.""""" C.99750 4.08318 0.00000 C.99653 -0.08337 ".0f"0" C.99652 -".oa287 C.CCCOC 0.99656
-0.06‘25 0.09000 0.990‘1 -"."5996 0.0"""" 0.99B2" -r.c%157 0.0"""" 0.95810 4.07,+* P.CCOIIC C.99744 4.‘)7160 C.CPCCO 0.997.3
-"."7118 C.""""" 0.99716
L/JlPI/S~R, 0.366.5E_01 0.35517E41 ".36695E-o, “.44lu)E~~
0.*38136-04 c.+C2llE-O* ".11555E-03 ".lC569E-03 ".16629E-03 ".66459E-"3 ".61296E-03 ".6295CE-03 O.L24~SC-05 ".11813E-"5 O.l‘85+E-05 0.88422E-05 ".8366.E-05 0.8395bL-"5
ll~~Tl/Slll 0.96335E 00 0.9336Q 00 0.99lClE 00 0.99995E 00 0.93lUE 00 ".91729E 00 0.99982E 00 ".9++83E CO 0.9~72oE 0" 0.99939E 00 0.9209Y CO ".949,0E 00 ".L00"?E 01 0.946leE 00 ".94969E 00 o.99999E JO 0.966‘8k 00 0.9.995E CO
".366"3E4, 0.35993E41 ".36603C-01 n.~613*E-O4 ".9*32lE4+ C.**580t-c~ C.‘T55.3E-"3 ".16721i-03 ".lb75‘E43 ".6631‘E-03 C.621822-"3 C.b326lE-"3 ".12179E-"5 0.1‘8(99E-OS C.l‘PlPE45 ".81136E-05 0.8.252245 ".8+468E-"5
0.963KlE 0.9160X 0.9539lE 0.9999SE 0.99832E 0.95101 C.PPPOLE ".95.?lY 0.953952 0.999392 C.931olE ".955,3E 0.10000E 0.95269E
00 0" 00 CO JO 00 00 co 00 CO ,0 00 01 00
0.955I.E
CO
0.3bSBlE-CL ".36268E-01 0.36587E-01 ".4b732E41 ".9969*C-01 ".4+8"9E-"9 0.11556143 ".16812E-03 O.l6834E-03 ".bb2SZE-03 0.629222-03 0.6352lE-03
0.963.Lk 0" O.QSSClE DO 0.958&x 0" 0.9999SE 00 0.95639E CO 0.95888E OJ 0.99902E 00
0.99999E 00 ".95268E 00 C.955‘OE 00
lSC.3951~
217
l .C.Jt‘C
213
218
218
22a
cos
“Ill
~.wccc ‘.cOCr u.c”*Cc
21t
WIITANEN
‘
"."3""3 C.C'XIC 4.C'Vh‘ C.OCO‘O L.CNC.0 ,.CJCCC cl.moFC O.CJCCC 4.000"" G.U')CCO C.CJCCC
Fig. 11. Sample
output
3.050”” C. 000”” 10n.00000 3.w000 0 .COOC’) 1w.00000 3.coccc C.WC0C
1cC.ccccc 3.00000 c.00cco 100.0000t
2.99781 -C.CCCCO
100.21314 2.62604 -C.CCCFO 107.20368
2.,,4425 -0.00000 11".00"00 1.70541 -O.oCoLC 123.22454
from ray-tracing
4.05255 C.0U000 0.99062 -0.05156 0.0000" 0.99867 -0.05295 6.00000 C.9986" -0.06625 0.00""" C.99Bl.9
program.
0.95lwE C.95aleE 0.99932E 0.94911E
0.95919E
00 00 00
00 00
SW text for discussion.
cosine data on each of these lines we may determine the trajectory of the ray before it strikes the jirst refracting surface. By looking at the intersection point data of each line we may obtain the coordinates of the point at which the ray intersects with the current surface (number under “SURF” column) indicated. The three directional cosine values for the ray are the directional cosines of the refracted ray as it leaves the surface (and they are, therefore, the directional cosines of the incident ray on the next refracting surface). Although the reflectance, energy of the refracted/reflected wave, and the intensity of the refracted/reflected wave appear to be related to the coordinate axes, lying as they do on the lines so labelled, they are not related to the axes. The data associated with the X-coordinate line corresponds to the reflectance (R) under the column labelled R/J (R)/S(R) while associated with the Y-coordinate line the energy (J(R)) of the reflected wave is to be found, and associated with the Z-coordinate is the intensity (S(R)) of the reflected wave. The data under the column headed T/J(T)/S(T) follows the same pattern but is for the transmissivity (T), energy (J(T)) of the transmitted ray, and the intensity (S(T)) of the transmitted ray. In addition to this information, whenever a ray crosses the optical axis (Z-axis) the point of that crossing is noted on the printout. Extensions to the program and a word of warning Several extensions to the program recommend themselves. One would be to include an option for plotter output so that the dioptric system could be drawn out and individual rays traced through it. Two different versions might be implemented: a sectional view of the system cut by a selected plane, and an isometric view. Another useful extension would be to increase the generality of the refracting surface equation by including, perhaps, exponential and transcendental functions in a general
Physiological
59
optics
form. This modification would have the disadvantage that the number of coefficients required as input data would go up. Considering the somewhat ad hoc method currently being used for handling apertures, a generalized aperture describing system would be worthwhile. To incorporate such a system the shape of an inner and an outer aperture would have to be defined in much the same way as the surface equations. In addition to the two surface equations that would be required for the inner and outer apertures, a card specifying the coordinate position and the Euler angles of the aperture pair would be required. A word of warning might be in order for persons planning to use the program as it stands on a computer having a word length of fewer than 36 bits: the machine being used may not have sufficient digits of accuracy for convergence to be attained. If this is true, then the program should be compiled in double precision, or relaxation of the convergence tolerance be instituted. To illustrate the use of the control and definition cards for the ray-tracing program, I will describe the cards needed for each of the examples presented earlier. Exanrple 1. The schematic human eye. In Table 2(a) the parameters describing the refracting surfaces used in the Gullstrand model for the human eye are listed. All measurements are in millimeters. In Table 2(b) are the computer-readable forms of the same parameters, arranged in the three card groups described above. (Note that the cards marking off the groups in Table 2(b) would never be present in an actual computer run,) Group 1. The switches card shows that a fan of rays has been selected (zero in column one), that rotational symmetry will apply (one in column two), and that the auxillary subprogram PROCES will not be used (zero in column three). The second card of this group is the title, and the last card defines the initial conditions: the refractive index of object space is 1.0; there are seven surfaces to be defined (Group 2 cards), and an accuracy of lop7 is required in the Newton-Raphson iteration. Group 2. Here are the seven quadruplets of cards that describe the seven refracting surfaces. Card A describes the coordinate position of the vertex of the surface and any coordinate rotation of that surface. Cards B and C contain values of the parameters for this surface (coefficients A, B, C, D, E, F, G) and the refractive index of the surface following it. Card D describes the aperture to be used for this surface. Group 3. This group defines the fan of rays to be used in the tracing process. The origin of the rays in one meter (1000 mm) to the left of the first refracting surface’s vertex (defined to be at the origin of the reference coordinate system). All rays are confined to the XZ-plane (done by setting Y = 0.0 on card FAN-2). Points defined at the ray origin will go from 8.0 to - 8.0 (mm) by steps of 0.5 mm while the points at the first surface (destination of the initial ray segment) will go from 4.0-O.Omm by steps of 0.5 mm. Example 2. The schematic eye of the cat. Table 3(a) illustrates the parameters used to describe the schematic cat’s eye, and in Table 3(b) the computer-readable form of the same parameters is given. The description follows that of the first example with Table 2(a). Gullstrand schematic human of the surface. The radius is the radius with respect to the origin of reference center defines the center of the sphere system in which the refracting surface
I 2 3 4 5 6 7
00 O.? 3.6 4 146 6 565 7.2 244
-77
-6X -100 -7911 +57h +60 +120
iw. (0.0. (0.0. (0.0. (0.0. (0.0. (0.0.
eye. The position of each surface is that corresponding to the vertex of curvature, and a plus sign means that the surface is concave system, while a minus sign means that the surface is convex. The defining the refracting surface with respect to the local coordinate is defined. R.I. is the refractive index of the medium following the surface
- 7 7) - 6.8) - IO.01 7.9 + 5.76) + 6.0) + 12.ot
II)
1.376 I.336 1.3X6 I.406 I.386 336 1.336
I
\.1 + \2 + y2 + 22 + y1 + y2 + $ +
.?’ + ?I + ,,1 + I‘. + ;,z + p + ?2 +
1
;z ;z ;z .5: :;: ;i
= = ; = = = =
ii ,.I ’ :: $ $ ).I
1 3
I I
I100 1 I 0 0
1
I
I I 1 I
I I I I
I
1 I I 1
0
0 0 0 0
0
0 0 0 0
-154 -13.6 ~20.0 -15822 +11.52
+120 +24.0
0 0 0 0 0 0 0
WAYNE
WIITANEN
Table 2(b) GROUP I 010 FAN/SYMMETRY/NO PROCES RAY PATH ANALYSIS OF GULLSTRAND 1.0 007+0.10000MX)E-06 GROUP
EYE - WIDE
PUPIL
5112169 HEADER
2..
.O
0.0 1.0 1.376 0.0 0.0 IO 1.336 0.0 0.0 1.0 1.386 0.0 0.0 1.0 I.406 0.0 0.0 IO 1.386 0.0 0.0 1.0 1.336 0.0 0.0 10 1 336 0.0
1.0 0 0 .O 1.0 .O .O .O 1.0 .O .O .O 1.0 .O 0 .O IO .O .O
0 1.0 0 0 0 1.0 0 .O GROUP 8.0 .O
SCHEMATIC
0.0 1.0
0.0 0.0
0.0 0.0
0.0 - 15.4
0.0 0.5 1.0
5.5 0.0 00
0.0 0.0 0.0
0.0 0.0 - 13.6
01
00 36 1.0
5.5 0.0 0.0
0.0 00 0.0
0.0 0.0 - 20.0
01
0.0 4.46 1.0
4.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 - 15.822
01
0.0 6.565 1.0
4.5 0.0 0.0
0.0 00 0.0
0.0 0.0 Il.52
01
0.0 7.2 1.0
4.5 0.0 0.0
0.0 0.0 0.0
0.0 0.0 12.0
01
0.0 24.4 1.0
4.5 0.0 00
0.0 0.0 00
0.0 0.0 24.0
01
0.0
25.0
0.0
00
01
0.5
40
0.0
0.5
IA IB IC ID 2A 2B 2C 2D 3A 39 3c 3D 4A 4B 4C 4D 5A 5B 5c 5D 6A 68 6C 6D 7A 7B 7c 7D
3.. - lrn.0 I.0
FAN-l FAN-2
Table 3(a). Schematic eye of the cat. Same remarks apply to this table as to the table for the Gullstrand schematic human eye SUrfaCe
Position
I , 3 4 5
0.0 0.68 5.20 13.70 21.83
Radius
Center
-8.57 - 7.89 - 7.20 + 8.05 + 12.03
(0.0, (0.0. (0.0. (0.0. (0.0.
+ +
R.I.
8.57) 7.89) 7.20) 8.05) 12.031
Equation
1.376 1.336 I 554 1.336 (1.3361
x1 x2 .x2 x2 .x2
+ + + + +
yz yz y2 y2 y’
+ + + + +
it 2’ z* 27 it
= = = = =
12 ? rl rz rl
A
B
C
D
I
1
I
0 0
1 1 I
I I I I
0 0 0 0
I 1
1 I
1
F
G
-17.14 - 15.78 - 14.40 + 16.10 f24.06
0 0 0 0 0
E
0 0 0 0
Table 3(b) GROUP I I IO SINGLE RAYS/ROTATIONAL SYMMETRY/NO SINGLE RAYS FOR THE SCHEMATIC CATS EYE. 1.0 txl5+0.ltxmmOE-06 GROUP 0 10 .O .O .O 1.0 .O .O .O 1.0
IO
0.0 0.0
0.0 0.0
0.0 - 17.14
0.0 0.68 1.0
84 0.0 0.0
0.0 0.0 0.0
0.0 0.0 - 15.78
01
0.0 5.2 1.0
8.4 0.0 0.0
0.0 0.0 0.0
0.0 0.0 - 14.4
01
0.0 13.7 I.0
7.3 0.0 0.0
0.0 0.0 0.0
0.0 0.0 16.1
01
0.0 21.83 1.0
4.75 0.0 0.0
0.0 0.0 0.0
0.0 0.0 24.06
01
1.0
1.336 0.0
0.0
12.5
0.0
0.0
01
IcmO 1000.0 - 1000.0 1000.0 0.0
00 0.0 0.0 0.0 0.0
0.5 1.0 2.0 4.0 0.0
0.0 0.0 0.0 0.0 0.0
01001 01002 01003 OlGil4 mm
0.0 1.0 1.376 0.0 0.0 I.0 1.336 00 0.0 1.0
1.5544 0.0 0.0 IO I 336 0.0 0.0
.O
.O 1.0
.O 0 0 IO 0
0
0.0
3.. 0.5 1.0
.O
2.0
0
4.0 0.0
.I)
HEADER
2.
.n
GROUP .O 0
PROCESS CALL I9 DECEMBER 1969
IA IB IC ID 2A 29 2c 2D 3A 38 3c 3D 4A 4B 4c 4D 5A 5B 5c 5D
Physiological Table
4(a). The optical
system
of the honey
Position
Radius
Center
R.I.
I 2 3 4 5 6
0.0 7.4 16.3 48.0 48.0
-43.0 -43.0 * 21.0 m
(O,O, - 43.0) (O,O, - 43.0) (0.0s (0.0.21.0) (O,OmJ
1.4896 1.4524 1.4350 1.3114 1.3477 1.3470
146.0
cl2
(0.O.Z)
61
bee (Apis mdlijka).
Surface
)
optics See the remarks
Eauation ‘iJ + ),” + :z = I2 g + y1 + ;z = ?-I .-=0 .x1 + y2 + ? = rz I=0 0
under
Table
2(a)
A
B
c
D
E
F
G
I
I
1
0
0
0
0
0
0 0
0 0
0 0
-86.0 -86.0 1.0 42.0 1.0 1.0
0 0 0 0 0 0
11100 0 0 1 0 0 0 0
I
I
0 0
Table 4(b)
I
GROUP 010 FAN/NO PROCES/SYMMETRY APIS MELLIFERA - COMPLEX DIOPTRICS 1.o 006+0.1tmtxwE-06 GROUP .O 1.0 .O .O .O 1.0 .O .O .O .O
RAYS FOR FOCAL
POINT HEADER
2. 0.0 1.o 1.4896 0.0 0.0 1.0 1.4524 0.0 0.0 0.0 1.4350 0.0 00
.O
.O .O I .o .O
1.0
1.3114 0.0 0.0 0.0 1.3477 0.0 0.0 0.0 1.3470 0.0
.oo .O .O .O .O .O .O .O .O GROUP .Ol .O
- PARAXIAL
0.0 1.0
0.0 0.0
0.0 0.0
0.0 - 86.0
0.0 7.4 1.0
16.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 - 86.0
01
0.0 16.3 0.0
16.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 1.0
01
0.0 48.0 I.0
16.0 0.0 00
0.0 0.0 0.0
0.0 0.0 42.0
01
0.0 48.0 0.0
16.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 1.0
01
0.0 146.0 0.0
14.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 1.0
01
0.0
2.0
0.0
0.0
01
0.0025
0.01
0.0
0.0025
IA LB IC ID 2A 2B 2c 2D 3A 3B 3c 3D 4A 48 4c 4D 5A 58 5C SD 6A 6B 6C 6D
3.. - loo.0 1.0
FAN-I FAN-2
one exception: the switches card (first card of Group 1) indicates that single rays are to be provided to the system (one in column 1). This means that the cards of Group 3 must be in the single ray format. Four rays will be traced (ray numbers 001, . . ., 004) and tracing will be stopped after ray 004 has been processed-the fifth card in Group 3, having a ray number of 000, will cause the tracing to be terminated. For each of the rays to be traced the ray origin will be one meter to the left of the first refracting surface, and the rays will be specified by means of two coordinate points (zero in column 68). The point of intersection of the ray with the optical axis is to be calculated and printed (one in column 69). The point coordinates are grouped as (X1.Y1,Z1), (x;,Y;,z;). Example 3. The compound eye of the honey bee. In Tables 4(a) and (b) are shown the Darameters that describe the bee’s eye (Varela and Wiitanen I2311 - -. and their computerreadable form. Only one new variation has been introduced: surfaces 3, 5, and -6 are planes, and as such have a radius of curvature equal to infinity. Their equations are particularly simple. Otherwise the card groups and descriptions follow those given for Example 1 (above). Table SlKfZKe
I 2 3 4 5 6
.5(a). The optical
system
of Culliphora
erythrocephda.
See Table
Z(a) for remarks
Position
Radius
Center
R.I.
Equation
A
B
C
D
E
F
G
0.0 7.1 I08 23.3 73.56 77.68
-21.16 -21.16
(0.0. - 21.16) (0.0. - 21.16) (0,0/m) (0.0.21.16) lO.O.cc I (0.O.K 1
1.4730 1.4530 1.4150 1.3370 1.3410 1.3490
.g + $ + .-I = 12 .x2 + ).1 + ? = r2 0 \.2 + J’ + -2 = ,i ;=O ; = 0
I I 0 1 0 0
I I 0 1 0 0
I I 0 1 0 0
0 0 0 0 0 0
0 0 0 0 0 0
-42.32 -42.32 1.0 42.32 1.0 1.0
0 0 0 0 0 0
21716 x x.
62
WAYNE WIITANEN Table
GROUP I 010 FAN/NO PROCES/SYMMETRY ANALYSIS OF RAY PATHS IN THE EYE OF CALLIPHORA 1.0 GROUP .O 1.0 .O .O .O 1.0 .O .O .O .O .O .O .O I .o .O .O
(SEITZ
CASE 19)
HEADER
OO6+O.IOOmmE-M 2.. 0.0 IO 1.4730 0.0 0.0 1.0 I .4530 0.0 0.0 00 I.4150 0.0 0.0 1.0 1.3370 0.0 0.0 0.0 1.3410 0.0 0.0 0.0 1.3490 0.0
.O .O
.O .O .O 0 .O .O GROUP 4.0 .O
S(b)
00 1.0
0.0 0.0
0.0 0.0
0.0 - 42.32
0.0 7.1 1.0
II.75 00 0.0
0.0 0.0 0.0
0.0 0.0 -42.32
01
0.0 10.8 0.0
Il.75 0.0 0.0
0.0 0.0 0.0
0.0 0.0
01
0.0 23.3 1.0
II.75 00 0.0
0.0 0.0 0.0
0.0 0.0 42.32
01
0.0 73.56 0.0
10.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0
01
0.0 77.68 0.0
3.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 1.0
01
0.0
2.94
0.0
0.0
01
0.5
4.0
0.0
0.5
1.0
IO
IA IB IC ID 2A 2B 2c 2D 3A 3B 3c 3D 4A 4B 4c 4D SA 5B 5C 5D 6A 68 6C 6D
3.. 1ocQ.o 10
FAN- I FAN-?
Example 4. The compound eye of the blowfly. Tables 5(a) and (b) present the parameters that describe the fly’s eye (Seitz [7]) and their computer-readable form. The card groups follow the description given for Example 1.
5. DISCUSSION In the preceding sections the details involved in using a computer-based ray-tracing technique have been presented along with some elementary examples illustrating its use. The intent was to introduce the worker in physiological optics to the method and to build a degree of confidence in using the program. A byproduct of the choice of illustrative examples was the discrepancy between the focal point of the Gullstrand schematic eye as calculated by the Gaussian lens formula and ray-tracing methods, and that between the focal point of the compound eye of Calliphora as calculated by the Gaussian formulae and ray-tracing methods. Should one wish to pursue the problems of ray-tracing further using the examples given above, some interesting things might be done: the effect of lens accommodation in the two simple eyes could be studied readily, as could the effects of changing the pupil diameter. In the latter, a better placement of the aperture could be effected by introducing a dummy plane (having the same refractive index as the medium in which it is embedded) into the system at the point where the aperture is to be introduced. The interested reader will find further examples of ray-tracing with respect to visual systems and model eyes in Lotmar [9] and Katzper [lo]. 6. SUMMARY In this tutorial paper the author has presented the basic processes of ray-tracing using digital computers. Illustrations of the technique were drawn from the literature about animal visual systems. A complete listing of the program is provided along with indications for its use. REFERENCES 1. H. Lass, Vt~ctor and Trrmr Amdysis. McGraw-Hill, NY (1950). 2. G. H. Spencer. A Germal Ray Tracing Prqram. IBM Corp. Research Paper RC-549 (1961).
Physiological optics
63
3. M. Born and E. Wolf, Principles ofOptics.Pergamon Press, NY (1965). 4. J. B. Scarboroueh. Numerical Mathematical Analvsis. Johns Hopkins Press, Baltimore (1962). 5. Sir W. S. DukerElder. A Textbook of Ophthalmology. Mosby, NY (1941). 6. P. 0. Bishop and G. J. Vakkur, The schematic eye in the cat, Vision Res. 3. 357-381 (1963). 7. G. Seitz. Der Strahlengang im Appositionsauge von Calliphora erythrocephula (Meig.), Z. uergl. Physiol. 59. 205-231 (1968). 8. F. G. Varela and W. Wiitanen, The optics of the compound eye of the honeybee (Apis me[lifera), J. Gem Physiol. 55, 336358 (1970). 9. W. Lotmar. Theoretical eye model with aspherics, J. Opt. Am. Sot. 61, 1522-1529 (1971). 10. M. Katzper, Simulation of optics experiments in Computers in Undergraduate Science Education (Commission on College Physics) (1971). About the Author-WAYNE A. WIITANENstudied mathematics at the Massachusetts Institute of Technology from 1953 to 1957, and received an A.B. in extension studies in 1967, an M.A. in Biology in 1969 and a Ph.D. in Biology in 1972, all from Harvard University. From 1954 to 1968 Dr. Wiitanen worked in the computer field in a variety of capacities: from hardware design to programming language design and implementation. Dr. Wiitanen is currently an Assitant Professor of Biology with a joint appointment in the Department of Computer Sciences at the University of Oregon, where he has been since 1971. He is also program director of the Systems Physiology Training Program in the Department of Biology at the University of Oregon. His research interest is the application of mathematics and computer science to biological problems with particular emphasis on neurobiology. He is a member of the American Association for the Advancement of Science. the Institute of Electrical and Electronics Engineers, and the Society for Computer Simulation.
APPENDIX Listing of FORTRAN program. DISCLAIMER Although the author of the following program has checked it extensively, neither the author nor the publisher makes any warranty, either implied or expressed, for the correct performance of the program. C
C
803 9095 813 C 2 9000 9001 loo0
GENERAL RAY TRACING PROGRAM FOR FORTRAN IV INTEGER SSl COMMON IRAY,IMl,XA,RPAR1,YA,RPAR2,ZA,RPAR3,XOBAR,CKBAR,YOBAR. *CLBAR.ZOBAR,CMBAR,XO,YO,ZO,ALPHA,BETA,GAMMA,A.B,C,D.E,F.G,REFR, *APl.AP2.AP3,AP4,APS.AP6,NAP,NOUT DIMENSION XO(35).YO(35) DIMENSION Z0(35),ALPHA(35) DIMENSION BETA(35),GAMMA(35) DIMENSION A(35),REFR(36) DIMENSION AP1(35),AP2(35) DIMENSION AP3(35),AP4(35) DIMENSION AP5(35),AP6(35) DIMENSION NAP(35),NOUT(35) DIMENSION CKBAR(35),CLBAR(35) DIMENSION CMBAR(35),XOBAR(35) DIMENSION YOBAR(35),ZOBAR(35) DIMENSION B(35),C(35) DIMENSION D(35),E(35) DIMENSION F(35).G(35) DIMENSION TITLE(12) DATA CONV/O.O1745329/,IPG/1/,LINE/O/ DATA IPROS/O/ DEFINE THE TANGENT TAN(A) = SIN(AlICOS(A) READ (5.9095) SSl ,IRS,IPROS FORMAT(311) IPG= 1 LINE=0 KKK=2 READ TITLE CARD READ (5.2) TITLE FORMAT( 12A6) WRITE (6,900O)TITLE,IPG FORMAT(lH1,25X,12A6,25X,5HPAGE 13) WRITE (6.9001) FORMAT(lH0.53X,lOHINPUT DATA//) IF (SSl) 1000,1001,1000 WRITE (6,9002)
64 9002 1001 9003 1002 9021 1013 9022 1020 C 102 C C 9004
C 103 C C 113 C C
C C C
C C
9005
9004
1003 9007
1004 9008
1005 9009
1006 9010 1007 9011
WAYNE WIITANEN FORMAT(43X,31HSWITCH SET FOR SINGLE RAY INPUT ) GO TO 1002 WRITE (6,9003) FORMAT(43X,32HSWITCH SET FOR FAN OF RAYS INPUT ) IF (IRS. NE. 1) GO TO 1013 WRITE (6,9021) FORMAT(43X,27HROTATIONAL SYMMETRY ASSUMED ) IF (IPROS. NE. 1) GO TO 1020 WRITE (6.9022) FORMAT(43X,69HNO OUTPUT OF RAYS WILL BE DONE, BUT SUBROUTINE PROCE *S WILL BE CALLED. ) CONTINUE JB=50 READ HEADER CARD FORMAT (F11.5.13.El5.8) READ (5,102)REFR(l),NOSUR,TOL REFR = INITIAL INDEX OF REFR. NOSUR = NO. OF SURFACES TOL = MIN. LIMIT ON ITERATION INCREMENT WRITE (6,9004) REFR( l),NOSUR,TOL FORMAT(43X.30HINITIAL INDEX OF REFRACTION = F11.5/43X, *21HNUMBER OF SURFACES = 12/43X.24HCONVERGENCE TOLERANCE = El2.5/,‘) NOSUR = NOSUR + 1 NOSR=O DO 203 I=Z,NOSUR READ SURFACE DATA FORMAT (F11.5,F11.5,F11.5,F11.5,F11.5,F11.5,11,11,13) READ (5,lO3)XO(I),YO(I,,Zo(I),ALPHA(I).BETA(I).GAMMA(I) X0, YO, AND ZO ARE LOCAL SYSTEM SURFACE VERTEX COORDINATES. ALPHA, BETA, AND GAMMA ARE EULER ROTATIONAL ANGLES. READ (5,113) A(I),B(I),C(I),D(I),E(I),F(I).G(I).REFR(I) FORMAT(6F11.5) A, B. C, D; E. F, AND G ARE SURFACE COEFFICIENTS REFR = INDEX OF REFRACTION FOLLOWING SURFACE READ (5,103) APl(I),AP2(I,,AP3(I),AP4(I).AP5(I),AP6(I),NAP(I), *NOUT(I) API, AP2, AP3, AP4. AP5, AND AP6 SPECIFY APERTURE DIMENSIONS, SEE INPUT WRITE UP. NAP = 0 FOR CIRCULAR-ANNULAR APERTURE, = 1 FOR RECT-TRAPEZOIDAL APERTURE, = 2 FOR HYPERBOLIC APERTURE. NOUT = 0 FOR NO OUTPUT AT SURFACE, = 1 OTHERWISE. NSR=I-1 WRITE (6,905) NSR,XO(I).YO(I),ZO(I),ALPHA(I),BETA(I),GAMMA(I) FORMAT(lX.RHSURFACE 12.5X.29HLOCAL COORDINATES (X,Y.Z) = (FI 1.5. *lH.FI 1.5,1H.F11.5,1H)/16X,25HEULER ANGLES (A.B.C.) = ( F11.5, *lH,Fll.5.1H,Fl 1.5.IH)) WRITE (6.9006) A(I),B(I),C(I),D(I),E(I),F(I),G(I),REFR(I) FORMAT(I6X.33HSURFACE EQUATION COEFFICIENTS, A= E15.8,3H.B=E15.8, *3H,C=E15.8,3H,D=El5.8/47X,2HE=El5.8,3H,F=E15.8.3H.G=El5.8/ *16X,22HINDEX OF REFRACTION = F11.5) IF (NAP(I)1) 1003,1004,1005 WRITE (6.9007) AP3(I),AP4(I),APl(I),AP2(I),ZO(I) FORMAT(16X,39HAPERTURE IS AN ANNULUS. INNER RADIUS = F11.5, *16H OUTER RADIUS = F11,5/16X,35HCOORDINATES OF APERTURE (X,Y,Z) = *(F11.5,lH,Fl1.5,1H,F11.5.1H)) GO TO 1006 WRITE (6.9008) APl(I),AP3(I),AP2(I),AP4(I).AP5(I).AP6(1) FORMAT(l6X,37HAPERTURE IS A RECTANGLE OR TRAPEZOID./l6X,24H(X,Y)-L *OWER BOUNDS ARE (F11.5.1H,F11.5,1H),2X,24H(X,Y)-UPPER BOUNDS ARE ( *F11.5,1H.F11.5,1H)/l6X,26HRECIPROCAL SLOPE OF LHS = F11,5,7H RHS = *F11.5) GO TO 1006 WRITE (6.9009) APl(I).AP2(I),AP3(I),AP4(I),AP5(I),AP6(1) FORMAT(16X,47HAPERTURE IS A HYPERBOLA, CENTER IS AT (X.Y) = ( *Fll.S,lH,Fl 1.5,1H)/16X.l8HSEMI-MAJOR AXIS = Fl 1.5,2X.l8HSEMI-MINOR *AXIS = F11.5.2X,30H(Y -LOWER,Y -UPPER) BOUNDS ARE (F11.5.1H.Fll.5, *lH)) IF (NOUT(1) .NE. 0 .AND. IPROS .NE. 1) GO TO 1007 WRITE (6.9010) FORMAT(16X,25HNO OUTPUT AT THIS SURFACE //) GO TO 1009 WRITE (6.9011) FORMAT(16X.22HOUTPUT AT THIS SURFACE //) NOSR = NOSR + 1
Physiological optics 1009
203 C 2050 9012 C 205 C C C C C C C C
9013
1008 9014 1010 9015
1011 9016
201 C C C C 116
9017
1014 1015 9019
IF (MOD&S) .NE. 0) GO TO 203 IPG=IPG+ 1 WRITE (6,900O) TITLEJPG WRITE (6,100) CONTINUE READ RAY DATA WRITE (6,9012) FORMAT( lH0,20HRAY DESCRIPTION DATA //) SWITCH 1 ON TO USE SINGLE RAY CARD INPUT IF (SSl) 205,201,205 READ (5,103)XA,YA,ZA,RPARl,RPAR2,RPAR3.NIN.NAXIN.IRAY AN IRAY OF ZERO OR LESS SIGNIFIES END OF SINGLE RAY JOB IF (IRAY .LE. 0) GO TO 803 XA, YA, AND ZA ARE RAY COORDINATES AT FIRST POINT RPARl, RPAR2, RPAR3 ARE RAY COORDINATES AT SECOND POINT OR RAY DIRECTION COSINES AT FIRST POINT NIN = 0 IF RAYS ARE SPEC. BY 2 POINTS NIN = 1 IF RAYS ARE SPEC. BY 1 POINT AND DIR. COSINES, NAXIN = 1 FOR COMPUTATION OF INTERSECTION OF RAY WITH OPTICAL AXIS, = 0 OTHERWISE IRAY = RAY NUMBER. IF (NAXIN .EQ. 1) GO TO 1008 WRITu6.9013) IRAY FORMAT(lOX.1 1HRAY NUMBER 13,5X,56H NO COMPUTATION OF INTERSECTION *OF RAY WITH OPTICAL AXIS ) GO TO 1010 WRITE(6,9014) IRAY FORMAT(lOX,llHRAY NUMBER 13,5X,71 H COMPUTATION OF INTERSECTION OF *RAY WITH OPTICAL AXIS WILL BE PERFORMED ) IF (NIN .EQ. I) GO TO 1011 WRITE (6,901s) XA.YA,ZA,RPARl.RPARZ.RPAR3 FORMAT(24X,38H(X,Y,Z) COORDINATES OF FIRST POINT = (F11.5,lH.Fll.5 *.l H.FI I .5.1H)/24X.39H(X.Y.Z) COORDINATES OF SECOND POINT = (F11.5, *IH.FI 1.5.1H.Fl 1.5.IH)) LINE=0 IPG=IPG+l WRITE (6,900O) TITLEJPG WRITE (6,l) GO TO 206 WRITE(6.9016) XA,YA,ZA,RPARl_RPAR2,RPAR3 FORMAT(24X,38H(X,Y,Z) COORDINATES OF FIRST POINT = (F11.5,lH,F11.5 *,lH,F11.5,lH)/24X,3lH(X.Y,Z) DIRECTIONAL COSINES = (F11.5,lH,Fll.5 *,lH,Fll.5,lH)) LINE=0 IPG=IPG+ 1 WRITE (6,900O) TITLEJPG WRITE (6.1) GO TO 206 READ (5,116)FXUl,FZl,XGAPl,FXU2,FZ2,XGAP2 READ (5,116)FYU2,YGAP2 RAY INPUT DATA. FXUl = MAX.XR, FZl = ZR, FXU2 = MAX. X AT POINT 2 XGAPl = X SPACING AT POINT 1. XGAP2 = X SPACING AT POINT 2 FZ2 = Z COORDINATE AT POINT 2 FYU2 = MAX Y COORDINATE AT 2ND POINT, YGAP2 = Y SPACING AT POINT 2 FORMAT (F11.5,F11.5,F11.5,F11.5,F11.5,F11.5) TEMP = - FXU 1 TEMP2= - FXU2 TEMP3 = - FYU2 WRITE(6,9017) FXUl,XGAPl,TEMP,FZl,FXU2,XGAP2,TEMP2,FZ2,FZl,FYUl. *YGAP2,TEMP3,FZ2 FORMAT(lOX,34HFAN OF RAYS... IST X-COORDINATE = Fll.5,1H(Fll.5,lH) *Fl1.5,5X.7HAT Z = Fll5/25X,l9H2ND X-COORDINATE = F11.5,1H(F11.5, *lH)Ft 1,5,5X,7HAT Z = Fll5/25X,24HlST Y-COORDINATE = 0(0)0,5X,7HAT * Z = Fll5/25X,l9H2ND Y-COORDINATE = Fll.5,1H(Fll.5,lH)Fll.5, *SX,7HAT Z = F11.5) IF (IRS .NE. 1) GO TO 1014 NRAYS=(2*IFIX(FXUl/‘XGAPl)+ l)*(IFIX(FXU2/XGAP2)+ l)*(IFIX(FYU2/ *YGAP2)+ 1) GO TO 1015 NRAYS=(2*IFIX(FXUl/XGAPl)+ 1)*(2*IFIX(FXU2/XGAP2)+1)*(2*IFIX( *FYU2/YGAP2)+ 1) WRITE (6,9019) NRAYS FORMAT(lOXJ7HNUMBER OF RAYS = 15) IF (IPROS .EQ. I) GO TO 1012
65
66
9020
9018
1012
812
8120 8030
814
4 804
8050 805
206 208
207
209
WAYNE WIITANEN NOSR=3*NOSR*NRAYS+?*NRAYS WRITE (6,902O) NOSR FORMAT(lOX,36HMAXIMUM NUMBER OF LINES OF OUTPUT = 16) IF (FXUl .GE. 0.0 .AND. XGAPl .GT. 0.0 .AND. FXU2 .GE. 0.0 .AND. * XGAPZ .GT. 0.0 .AND. FYU2 .GE. 0.0 .AND. YGAP? .GT. 0.0) GO TO * 1012 WRITE (6,9018) FORMAT(69H X OR Y MAXIMUM VALUE(S) AND ITS (THEIR) INCREMENTS ARE *INCONSISTENT. / 64H X AND Y MAXIMA MUST BE GE. 0.0 AND INCREMENTS * MUST BE .GT. 0.0 j ’ GO TO 803 IPG=IPG+l WRITE (6,900O) TITLE,IPG WRITE (6,l) IRAY =0 GAP3=0. RPAR?=FYU2 RPARZ=RPARZ-GAP3 IF (IRS .EQ. 0) GO TO 8120 IF fRPAR2) 8030.814.814 IF (RPARZ’+FYU~) k030,814.814 IF (IPROS .NE. 1) GO TO 803 ICODE = 3 CALL PROCES(ICODE,I) GO TO 803 GAP3 = YGAP? GAP1 =O. XA=FXUl XA=XA-GAPI IF (XA+FXU1)812.804,804 GAP1 = XGAPl GAPZ=O. RPARl = FXU2 RPARl = RPAR 1- GAP2 IF (IRS .EQ. 0) GO TO 8050 IF (RPARl) 4,805,805 IF (RPARl + FXUZ) 4,805,805 GAP2= XGAPZ YA = 0.0 ZA=FZl RPAR3 = FZ2 NIN=O NAXIN = 1 IRAY = IRAY + 1 NIN=NIN+ 1 GO TO (207,208),NIN CKBAR( 1) = RPAR 1 CLBAR(i)=RPARZ CMBAR(l)=RPAR3 GO TO 209 XD=RPARl -XA YD=RPARZ-YA ZD=RPAR3-ZA RALEN=SQRT(XD*XD+YD*YD+ZD*ZD) CKBAR( 1) = XD/RALEN CLBAR(l)=YD/RALEN CMBARll 1= ZD/RALEN XOBAR(i)-XA ’ YOBAR(l)=YA ZOBAR(I)=ZA SI = 1.0 DO 15 I=Z,NOSUR IMl=I-1 SA = SIN(CONV*ALPHA(I)) CA = COS(CONV*ALPHA(I)) SB=SIN(CONV*BETA(I)) CB = COS(CONV*BETA(I)) SG = SIN( CONV*GAMMA(I)) CG = COS(CONV*GAMMA(I)) Xl =XOBAR(IMl)-X0(1) Yl =YOBAR(IMl)-YO(1) Zl =ZOBAR(IMl)-ZO(1) Rll=CA*CG+SA*SB*SG
Physiological optics
311 210
302 211 5
213
8
217 10
1021 104
11 218
219 220 9 301 300 221
1016 105
R12= -CB*SG R13=CA*SB*SG_SA*CG R21 =CA*SG-SA*SB*CG R22= CB*CG R23 = -(SA*SG+ CA*SB*CG) R31= SA*CB R32=SB R33=CA*CB XO=Xl*Rll+YI*R12+ZI*R13 YO=Xl*R21+Yl*R22+Zl*R23 ZO=Xl*R31+YI*R32+Zl*R33 CK=CKBAR(IMl)*Rll +CLBAR(IMl)*R12+CMBAR(IMl)*R13 CL=CKBAR(IMl)*R21 +CLBAR(IMl)*R22+ CMBAR(IMl)*R23 CM = CKBAR(IMl)*R31+ CLBAR(IMl)*R32+ CMBAR(IMl)*R33 J=O IF (CM)210.302.210 X0 = X0 - CK*ZO/CM YO=YO-CL*ZOICM z2=0. GO TO 211 z2=20 s=o. J=J+I X=XO+CK*S Y =YO+CL*S Z=Z2+CM*S FCN=A(I)*X*X+ B(I)*Y*Y +C(I)*Z*Z+ D(T)*X+ E(I)*Y + F(I)*Z -G(I) FX = 2.0*4(1)*X + D(1) FY = 2.0*B(I)*Y + E(1) FZ = 2.O*C(Ij*Z + F(1) DETMT=CK*FX+CL*FY+CM*FZ AMAGN = SQRT(FX*FX+ FY*FY + FZ*FZ) IF (DETMT)218,217,218 IF (FCN) 10.9.10 IF (IPROS .NE. I) GO TO 1021 ICODE= CALL PROCES(ICODE,I) GOT0 II WRITE (6.104) IRAY,IMl,XA,RPARl,YA,RPAR2,ZA,RPAR3 FORMAT(lX,I4,13,13H NO INCIDENCE 2X,8HX OR K 2(Fl1.5,1X)/ *23X,8HY OR L 2(Fl lS,lX)/23X,8HZ OR M 2(Fll,S,lX)//) LINE=LINE+6 IF (LINE .LT. S1) GO TO 11 IPG=IPG+ 1 LINE=0 WRITE (6.9ooO) TITLEIPG WRITE (6.1) CONTINUE IF (SSI) 205.3,205 DELS = - FCNiDETMT DELSZ = DELS*DELS TOL2 = TOL*TOL IF (DELSZ-TOL2)9.9,219 IF (J-JB)220,10,10 S=S+DELS GOT05 IF (CM) 300,301,300 IF (S)221,222,222 IF (S-ZO/CM)221.222,222 IF (IPROS .NE. 1) GO TO 1016 ICODE = 2 CALL PROCES(ICODE.1) GOT0 11 WRITE (6.105) IRAY,IMl,XA,RPARl,YA,RPAR2,ZA,RPAR3 FORMAT(lX.I4.13,13H VIRTUAL PATH 2X,8HX OR K 2(Fll.S,lX)/ *23X,8HY OR L 2(F11,5,1X)/23X,8HZ OR M 2(F11,5,1X)//) LINE = LINE + 6 IF (LINE .LT. 51) GO TO 11 LINE=0 IPG=IPG+ I WRITE (6.9000) TITLE,IPG WRITE (6.1) GOT0 11
67
68 222 223 12
1017 106
226 224 227 228 229 225 230 231
13
233
232
234
1018 107
235 236 237 14
C C
C
C C
WAKNE WIITANEN KAP = NAP(I) + 1 GO TO (223,224,225),KAP RHSQ =(X-APl(I))*(X-APl(I))+(Y -AP2(I))*(Y -AP2(1)) IF (RHSQ -AP3(I)*AP3(1))12,226,226 IF (IPROS .NE. 1) GO TO 1017 ICODE = 2 CALL PROCES(ICODE,I) GO TO 11 WRITE (6,106) IRAY,IMl,XA,RPARl,YA,RPAR2,ZA,RPAR3 FORMAT(lX,I4,13,14H APERTURE STOP lX,8HX OR K 2(F11.5.1X)/ *23X,8HY OR L 2(F11,5,1X)/23X,8HZ OR M 2(Fll.S,lX)//) LINE=LINE+6 IF (LINE .LT. 51) GO TO 11 LINE=0 IPG=IPG+ 1 WRITE (6,9ooO) TITLEJPG WRITE (6,l) GOT0 11 IF (RHSQ-AP4(I)*AP4(1))13,12,12 IF (Y -AP3(1))12,12,227 IF (Y -AP4(1))228,12,12 IF (X-(APS(I)*(Y-AP3(I))+AP1(1)))12,12,229 IF (X-(APqI)*(Y -AP3(I))+AP2(1)))13,12,12 IF (Y -AP5(1))12,12.230 IF (Y -AP6(1))231,12,12 FFF=(X-APl(I))/AP3(1) FFF2 = FFF*FFF GGG = (Y - AP2(I))/AP4(1) GGG2 = GGG*GGG IF (FFFZ-GGGZ-1.)13,12.12 RAT= REFR(IMl)/REFR(I) ALC=(RAT*DETMT)/(FX*FX+FY*FY+FZ*FZ) IF (REFR(1) + REFR(IMl))232,233,232 GAMUC= 2.*ALC RAT= 1. GO TO 14 BLC=(RAT*RAT-l.)/(FX*FX+FY*FY+FZ*FZ) DISC = ALC*ALC - BLC IF (DISC)234,235,235 IF (IPROS .NE. 1) GO TO 1018 ICODE= CALL PROCES(ICODE,I) GOT0 11 WRITE (6.107) TRAY,IMl,XA.RPARl.YA,RPAR2,ZA,RPAR3 FORMAT(lX,I4,13,14H NO REFRACTION lX,8HX OR K 2(F11.5,1X)/ *23X,8HY OR L 2(F11,5,1X)/23X,8HZ OR M 2(F11.5,1X)//) LINE = LINE + 6 IF (LINE .LT. 5 1) GO TO 11 LINE = 0 IPG=IPG+ 1 WRITE (6,9OGO)TITLEJPG WRITE (6,l) GO TO 11 DISC=SQRT(DISC) IF (ALC)236,10,237 DISC = - DISC GAMUC = DISC - ALC CK=RAT*CK+GAMUC*FX CL = RAT*CL + GAMUC*FY CM =RAT*CM + GAMUC*FZ CALCULATE THE ANGLE BETWEEN THE INCIDENT RAY AND THE NORMA1 UN-NORMALIZED VALUE IN DETMT. NORMALIZE IT FIRST. DETMT = - DETMT/AMAGN CT1 = ABS(DETMT) PHI = ABS(ATAN(SQRT(l,O_DETMT**2)/DETMT)) THE ANGLE OF THE REFRACTED RAY IS DETMT = -(CK*FX+ CL*FY + CM*FZ)/AMAGN CIT = ABS(DETMT) PHIP = ABS(ATAN(SQRT(l,O_DETMT**2)/DETMT)) CHECK FOR NORMAL INCIDENCE (PHI = 0) IF (PHI .NE. 0.0) GO TO 1022 USE SPECIAL FORMULA FOR NORMAL INCIDENCE RF = ((1.0 - RAT)/(l.O + RAT))**2
Physiological optics
C C
C 1023 C
C
16
1019
15
239 241 18
242 244 243 240 246 245
17
1050
100
GO TO 1023 INCIDENCE IS AT AN ANGLE TO THE SURFACE NORMAL CALCULATE REFLECTANCE AND TRANSMITTANCE FOR UNPOLARIZED LIGHT RL = (TAN(PHI_PHIP)/TAN(PHI+PHIP))**2 RR = (STN(PHI- PHIP)/SIN(PHI+ PHIP))**2 RF = (RL + RR)/2.0 TAKE ADVANTAGE OF THE FORMULA RF+TR= 1 TR = 1.0 - RF CALCULATE ENERGY EJR = RF*SI*CTI EJT = TR*SI*CTI CALCULATE INTENSITY AJR = RF*SI AJT = TR*SI*CTI/CTT SI = AJT XOBAR(I)=Rll*X+R2l*Y+R3l*Z+XO(I) YOBAR(I)=Rl2*X+R22*Y +R32*Z+YO(I) ZOBAR(I)=Rl3*X+R23*Y+R33*Z+ZO(I) CKBAR(I)=Rll*CK+R21*CL+R3l*CM CLBAR(I)=R12*CK+R22*CL+R32*CM CMBAR(I)=Rl3*CK+R23*CL+R33*CM KOUT= NOUT(I)+ 1 GO TO (15,16),KOUT IF (IPROS .NE. 1) GO TO 1019 ICODE = 1 CALL PROCES(ICODE,I) GO TO 15 WRITE (6,108) IRAY,IMl,XA,RPARlS(OBAR(I),CKBAR(I),RF,TR WRITE (6,109) YA,RPAR2,YOBAR(I),CLBAR(I),EJR,EJT WRITE (6,110) ZA,RPAR3,ZOBAR(I),CMBAR(I),AJR,AJT LINE = LINE + 3 IF (LINE .LT. 51) GO TO 15 LINE =0 IPG=IPG+ 1 WRITE (6,900O) TITLE,IPG WRITE (6,l) CONTINUE NAXIN = NAXIN+ 1 GO TO (l7,239),NAXIN IF (CKBAR(NOSUR))240,241,240 IF (XOBAR(NOSUR))l8,242,18 IF (IPROS .EQ. 1) GO TO 1050 WRITE (6,111) IRAY LINE = LINE + 1 IF (LINE .LT. 51) GO TO 17 IPd=IPG+ 1 ’ LINE=0 WRITE (6Q000) TITLE,IPG WRITE (6,l) GOT0 17 S = - YOBAR(NOSUR) IF (CLBAR(NOSUR))243,244,243 IF (YOBAR(NOSUR))l8,245,18 S = S/CLBAR(NOSUR) GO TO 245 S= -XOBAR(NOSUR)/CKBAR(NOSUR) IF (CLBAR(NOSUR))246,244,246 IF (S+YOBAR(NOSUR)/CLBAR(NOSUR))18,245,18 AXIN=ZOBAR(NOSUR)+S*CMBAR(NOSUR) IF (IPROS .EQ. 1) GO TO 1050 WRITE(6,112) IRAYAXIN LINE = LINE + 1 WRITE(6,lOO) LINE = LINE+ 1 IF (LINE .LT. 51) GO TO 1050 LINE=0 IPG=IPG+ 1 WRITE (6,9ooO) TITLE,IPG WRITE (6,l) IF (SSl) 2050.3.2050 FORMAT(lHOQH RAY SURF,26X,6HOBJ PT,4X,8H2D PT/DC,SX,6HINT PT,6X, *7HDIR COS,9X,ll HR/J(R)/S(R), 6X,1 1HT/J(T)/Sv)) FORMAT( 1X)
69
70 108 109 110 Ill 112 c
WAYNE WIITANEN FORMAT(lX,I4,13,15X,8HX OR K ,4(Fl 1.5.lX),7X,El2.5, 5XJl2.5 FORMAT(23X,8HY OR L .4(F11.5,lX),7X,El2.5. 5X.El2.5) FORMAT(23X,8HZ OR M ,4(Fll.5,lX),7X.El2.5. 5X,El2.5) FORMAT(lX,I4,3X,l5HNO INT W/OPT AX) FORMAT(lX.I4,3X,23HOPT AX INTERSECTION Z = Fll.5) END DUMMY SUBROUTINE PROCES SUBROUTINE PROCES(A,B) RETURN END
