I wanted to find or create a formula which would accept an epoch timestamp, latitude and longitude and it would produce the Sun's observed azimuth and altitude in radians. It needs to take into account details earth's axial tilt and its position on its orbit around the sun. To my surprise I wasn't able to find such formula, so I had to develop it from scratch. Luckily earth's orbit (and orbits in general) is a well studied and documented problem, so I could take some shortcuts.
 θ = 0 ... 2
π(1)
Next corresponding eccentric anomaly E is calculated from θ and the orbit's eccentricity e:
 E =
atan2(√1  e · sin(θ / 2), √1 + e · cos(θ / 2))(2)
Mean anomaly M is calculated from E by using Kepler's equation:
 M = E  e sin E(3)
The M is trivial to calculate for any desired calendar date, so we want to have an approximation formula θ = f_{e}(M) so that previous equation hold reasonably well. I noticed experimentally that it would be a good approach to approximate
 f'_{e}(M) = (M  θ) /
e(4)
 θ = f_{e}(M) = M  e · f'_{e}(M)(5)
The effect of e on f'_{e}(M) shape can be seen in Figure 2. Most importantly its magnitude is around ±2, and has values of exactly zero at M = 0, M = π and M = 2π. Due to these properties, it is approximated by the polynomial
 f'_{e}(M) → a(e) · (M  pi) · (M · (2π 
M))^{b(e)}(6)
The shape of f'_{e}(M) for various values of e can be seen in Figure 2. It is notable that for M → π / 2 it has a value of about 2, which indicates that θ > M around that region. This occurs because when M → 0 the planet is closest to the Sun, which means it has a higher angular velocity. Thus θ grows initially faster than the mean anomaly M.
The eccentricity e's role on a(e) and b(e) is solved by finding optimal a and b for various different values of e and fitting a model to those findings. Parameters are estimated by minimizing the standard sum of squared errors, which is achieved by Matlab's lsqnonin function. These results are shown in Figure 3. Luckily both parameters seem to follow the same formula y(e) = y_{0} + exp(k_{1}· log(e) + k_{2}), and k_{1} and k_{2} can be reliably estimated by fitting a linear trend to the loglog data points. The constant y_{0} (which corresponds to a_{0} and b_{0}, one for both parameters) is estimated so that the best loglinear model can be obtained. If that value is off then the loglog plot is no longer linear, and a more complex model is needed.
Model accuracy is determined by calculating the maximum error from the true location within the whole 360° orbit. If e = 0 then the orbit would be exactly circular and θ = M. Interesting cases are small values of e, for example on earth's orbit e → 0.0167. Two different models are compared in figures 5 and 6: a simpler model which only uses constants a_{0} and b_{0} regardless of the value of e, and the more accurate model with edependent corrections.
The biggest benefit of these corrections occur when 0.05 < e < 0.3. When e → 0 these corrections are small anyway, and on bigger values of e the model loses its accuracy even with these corrections. The ratio of simple model's maximum error and better model's error is show in Figure 4. When e = 0.1 the maximum ratio of 4 is reached, which means that the simpler model has 300% higher maximum error than the ecorrected version.
Error rates (in degrees) and the fitted trend curve for 0 < e < 0.1 are shown in Figure 5. Within this region the error seems to grow approximately linearly in respect to e, and for Earth's e → 0.0167 the maximum error is less than 0.1°. It should be sufficient for most applications.
Error rates for evalues of upto 10^{0.5} → 0.3 are show in Figure 6. In this context the error grows in proportion to e^{3.09}, reaching a maximum error of 6° at e = 0.3. It might be too high for most applications, but this elliptic orbits are not common in at least our solar system. Our top3 planets with highest eccentricities are Mercury (e = 0.21), Mars (e = 0.09) and Saturn (e = 0.05). Pluto has a remarkably high eccentricity of 0.25.
After this model is in place, it enables us to calculate Earth's location in respect to Sun in polar coordinates for a given timestamp. It just needs to be converted to "seconds since January 1 of this year", which is further converted to the mean anomaly M (in radians). True anomaly θ is calculated as
 θ _{e}(M) = M  e · a(e) · (M  pi) · (M · (2π 
M))^{b(e)}
a(e) = a_{0} + exp(a_{1}· log(e) + a_{2})
b(e) = b_{0} + exp(b_{1}· log(e) + b_{2}) (7)
Estimated parameters are:

a_{0} = 0.0804
a_{1} = 0.9869
a_{2} = 1.0123
b_{0} = 1.3848
b_{1} = 0.6143
b_{2} = 1.0715 (8)
The next step would be to estimate the local coordinate system from latitude, longitude and time, and to use polar coordinates to calculate the person's orientation with respect to the heliocentric coordinate system where periapsis (around 3.1. each year) is the "zero" direction. Then the Sun's apparent altitude and azimuth can be easily calculated.
Related blog posts:








