Let us consider a ray incident on a spherical interface between two media with different indices of refraction n and n’, respectively.  We want to trace the ray across this interface.  In a medium with constant index of refraction the ray's path is a straight line.  At the interface it obeys Snell’s law, n sinθ  = n' sinθ'.  The physics is simple, but an exact ray tracing calculation can become quite cumbersome for a ray moving in three dimensions. 

Let us set up a coordinate system as shown.

The spherical interface has a radius of curvature R, and the center of curvature C lies on the z-axis.  We call the z-axis the optical axis.  The vertex V is the point where the surface intersects the z-axis.  The direction of propagation of the ray has a positive z-component.  R is positive if the z-coordinate of C is greater than the z-coordinate of V and negative if the z-coordinate of C is less than the z-coordinate of V.

We assume that the staring position of the ray has coordinates (x, y, z). 
The unit vector labeling the initial direction of the ray makes angles A, B, and C with the x, y, and z axes, respectively.  The x-component of this unit vector is α = cos(A), the y-component is β = cos(B) and the z-component is γ = cos(C). 
We can only choose two of these direction cosines independently, since α² + β² + γ² = 1.

Let the vertex V of the spherical interface have z-coordinate t.

We want to find the coordinates (x’, y’, z’) of the point P where the ray intersects the interface, and the angle θ it makes with the normal to the interface.  Then we can apply Snell’s law and find the angle θ’ the ray makes with the normal after refraction and calculate the new direction cosines α’, β’, and γ’.
After determining (x’, y’, z’) and α’, β’, and γ’, we can continue tracing the ray to the next interface or to a target.

Aside:  The 3D figure on the right was produced using Jupyter Notebook, (Python, Matplotlib), and is interactive when using the Anaconda base environment.   You can rotate the 3D figure and zom in and out using the mouse.)  Find the notebook on Canvas, Modules, Numerical. 3D_spherical_interface_setup.ipynb.  You can also find the straight python code there (3D_spherical_interface_setup.py).  It is interactive when run in VS code.

(Ask an AI how to install and setup VS code (with or without Anaconda already installed on your computer.)

Details:  (This guides you you are interested in how to write an exact ray tracing program.  If not, concentrate on the paraxial approximation.)

Defining the vectors T and M and the distance e as shown in the figure, we can derive the exact ray tracing equations using properties of vectors.  (For a detailed derivation see "Exact Ray Tracing" under Modules, Optics Papers, on Canvas.)

The relevant quantities are given below.
e = tγ – (xα + yβ + zγ).
M_z = z + eγ – t.
M² = x² + y² + z² - e² + t² - 2tz.
cos(θ) = +(y² – (M²/R – 2M_z)/R)^½ ,
T = e + (M²/R – 2M_z)/(γ + cos(θ)).

The coordinates of the point P where the ray intersects the interface are

x’ = x + Tα,
y’ = y + Tβ,
z’ = z + Tγ.

Let μ = n/n’.
Then cos(θ’) = +( 1 – μ²(1 - cos(θ))^½ .

We define
g = cos(θ’) - μ cos(θ’),
and the new direction cosines are

α’ = μα – gx’/R,
β’ = μβ – gy’/R,
γ’ = μγ – gz’/R + g.

This traces the ray across the interface.  We have calculated the coordinates of the point P where the ray intersects the interface and the direction cosines of the ray after the interface.

We can then continue to the next interface with radius R’ and a vertex at t’ by simply redefining
(x, y, z) = (x’, y’, z’ - t),
α = α’, β = β’, γ' = γ, 
t = t’ – t,
R = R’.  
For a flat target we set R' = ∞. 
This outlines how a computer program can be written in most languages to exactly trace a ray through any number of interfaces.  We will use these exact ray tracing equations later, when we investigate aberrations.

Links: Canvas, Numerical, 3D_spherical_interface_full

To derive the mirror and lens equations for spherical mirrors and for lenses with spherical surfaces that you were introduced to in your introductory physics course, we make the paraxial approximation.  We assume an optical axis can be defined, and all surface normals make small angles with this axis.  Assume that all light rays also make small angles with the optical axis and stay close to the optical axis.  Such rays are called paraxial rays.  Then x and y are much smaller than R, and can be considered small quantities.

For the angles A, B and C we then have
A = 90^o – θ_x ,  α = cos(A) = sin(θ_x),
B = 90^o – θ_y ,  β = cos(B) = sin(θ_y),
γ = cos(C).

We assume C, θ_x, and θ_y, are small and that
sin(θ_x) = θ_x,  sin(θ_y) = θ_y,  cos(C) = 1,
i.e. we only keep terms lower than second order in the expansion of the sine and cosine functions.
(To first order in an expansion, sin(x) = x, cos(x) = 1.)

The small quantities we are dealing with in the exact ray tracing equations are x, y, θ_x, and θ_y.
Any term that contains products of these small quantities is set equal to zero.

The exact ray tracing equations are then replaced by the equations below.

e = t – (xθ_x + yθ_y + z) = t - z.
e² = t² + z² - 2tz.
M_z = z + e – t = 0.
M² = x² + y² = 0.
cos(θ) = (1 - (x² + y²)/R²)^½ = 1.
T = e.

The coordinates of the point P where the ray intersects the interface then are

x’ = x + (t - z)θ_x,
y’ = y + (t - z)θ_y,
z’ = t.

Note:  In the paraxial approximation the point P lies in the xy-plane containing the vertex for all rays.  Difference in the z-coordinate due to curvature are ignored, because we are keeping only first order terms.

To find the angles θ_x' and θ_y' in the paraxial approximation, we approximate
cos(θ’) = (1 – μ²(1 - cos(θ)))^½ = 1.

Then
g = 1 – μ,
α’ = sin(θ_x’) = θ_x’ = μ θ_x – (1 – μ)x'/R,
β’ = sin(θ_y’) = θ_y’ = μ θ_y – (1 – μ)y'/R,

or

n’θ_x’ = nθ_x – (n’ – n)x'/R,
n’θ_y’ = nθ_y – (n’ – n)y'/R.

Summary:

Ray tracing equations in the paraxial approximation are

x’ = x + (t - z)θ_x,
y’ = y + (t - z)θ_y,
z’ = t,

n’θ_x’ = nθ_x – (n’ – n)x'/R,
n’θ_y’ = nθ_y – (n’ – n)y'/R.

In the paraxial approximation the ray tracing equations for the projections of the rays on the xz-plane and on the yz-plane decouple and the projections can be treated independently.  The projections behave as if the rays were actually lying in that plane.  Rays that lie in a single plane containing the z-axis are called meridional rays.  Ray that do not lie in in a single plane containing the z-axis are called skew rays.

AI Study Tip:  Example prompt:  'Act as a physics professor.  Walk me through the geometric derivation of the single spherical interface equation.  Start with a ray originating at object point P, hitting the interface at height h.  Use the exterior angle theorem for the triangles formed by the object, the center of curvature, and the image.  Show me exactly at which step we replace make the small angle approximation.  Ask me a question at the end to check if I understand how the sign convention for R changes if the surface is concave versus convex.'