This is a derivation of a formula for the elevation angle of the sun. The elevation angle of the sun is the angle at a point on Earth between the Sun and the horizon. It is a function of not only latitude but time of day and time of year. Let the center of the Earth be at the origin of the coordinate system in which the xy-plane is the equatorial plane of the Earth, as shown.

The xy plane of the Earth is tilted with respect to the plane of Earth's orbit. The Sun will always be in the xz plane.

Let φ be the latitude angle. (Often in spherical coordinate systems it is the co-latitude which is used; i.e., the angle measured from the north pole of the sphere.) Let r be the radial vector from the center to a point on the sphere. The projection of r onto the xy plane has a length rcos(φ). The z coordinate is rsin(φ). The x coordinate is the projection of the projection of the r vector onto the x axis. It is equal to (rcos(φ))cos(θ). Likewise the y coordinate is the projection of the projection of the r vector onto the y axis. It is equal to (rcos(φ))sin(θ).

For future reference these relationships are:

x = (rcos(φ))cos(θ)
 
y = (rcos(φ))sin(θ)
 
z = rsin(φ)
 

Now the tilt of the Earth's axis with respect to its orbital plane must be taken into account. Consider a projection of the axis of the Earth onto the xz-plane of the coordinate system. The angle that this projection makes with the z-axis is denoted as ψ.

The sun is so far distant that its rays are essentially parallel. The angle which the sun's rays make with the x-axis is the same as the angle ψ. The angle ψ is a harmonic function of the number of days T after the equinox. The spring equinox of the Northern Hemisphere will be used as the point of reference so

ψ = ψ₁sin(360(T/365.25))
 

where ψ₁ is the angle of inclination of the spin axis of the Earth with the plane of its orbit, approximately 23.5°.

The angle between the sun and the vertical is called the zenith angle. The elevation angle of the sun is 90° minus the zenith angle. The cosine of the zenith angle can be found as the vector dot product of the unit normal (vertical) and the unit vector in the direction of the sun.

The unit normal at latitude φ and longitude θ on a sphere is

n = (x/r)i + (y/r)j + (z/r)k
= cos(φ)cos(θ)i + cos(φ)sin(θ)j + sin(φ)k
 

where a symbol displayed in red is a vector. The radius of the Earth is r, and i, j and k are unit vectors in the x, y and z directions, respectively. The sun's rays come from a direction in the xz plane with an angle of inclination of ψ. The unit vector in the direction of the sun is

cos(ψ)i + sin(ψ)k.
 

The dot product of this vector with the unit normal is

cos(ψ)cos(φ)cos(θ) + sin(ψ)sin(φ)
 

This is the cosine of the zenith angle so the angle of elevation η is given by:

η = 90° - cos^-1[cos(φ)cos(ψ)cos(θ)+sin(ψ)sin(φ)]
 

The angle ψ corresponds to the time of the year past the equinox; i.e.,

ψ = 23.5°sin(360°(T/365.25))
 

where T is the days past the spring equinox.

The longitude angle corresponds to the time t of hours past noon; i.e.,

θ = 360°(t/24).
 

Illustration of the Use of the Formula

The elevation angle η of the sun at 4 P.M. on April 21st at latitude 45° is found, assuming the spring equinox occurs at noon 00 P.M. on March 21:

T = 10 + 21 + 4/24 = 31.167
ψ = 23.5°sin(360(31.167/365.25)) = 23.5sin(30.719°) = 12.00°
θ = 360(4/24) = 60°
cos(φ) = cos(45°) = 0.7071
η = 90° - cos^-1[(0.7071)cos(12°)cos(60°)+sin(12°)sin(60°)]
= 90° - cos^-1(0.5259) = 90 - 58.273 = 31.73°
 

Finding the Time of Sunrise and Sunset

The value of t that for the above day which makes η=0 (sun rise) is such that

(0.7071)cos(θ)cos(12°) + sin(θ)sin(12°) = cos(90°) = 0
which, upon division by cos(θ)cos(12°) is equivalent to
tan(θ)tan(12°) = -0.7071
or
tan(θ) = -3.3266
and hence
θ = -73.269°
 

The value of t giving this value of θ is
t = 24(-73.269/360) = -4.885 hours before noon = 7:07 A.M.

Sunset is just 4.885 hours after noon = 4:53 P.M.

Revised Formulation

The equation for the elevation angle η can be put into the more convenient form

cos(90° - η)/H = cos(ψ + α)
 

where

H = [(cos(φ)cos(θ))² + sin²(θ)]^1/2
and
sin(α) = sin(θ)/H