Calculating Earth-Sun distance only using trigonometry and terrestrial measurements

Question

I have an interest in mathematics and astronomy, and my curiosity was piqued by the methodologies employed in contemporary solar distance measurement. Current practices use sophisticated instruments and celestial events, such as the transit of Venus, to calculate this distance. But then, how is the distance to Venus determined initially? Motivated by this inquiry, I undertook the challenge of calculating the Earth-Sun distance using only terrestrial measurements that are readily verifiable, coupled with basic trigonometry. However, as you will observe, there is a significant discrepancy between my findings and the accepted value of 1 Astronomical Unit (AU), which is approximately $\approx 150,000,000$ km. I am eager to understand and explain this variation. Below, I explain exactly what I did and what formula I derived.

---

Method

---

Start the process by selecting two points along the equator, separated by a distance $d$, and proceed to measure the solar altitudes $\alpha$ and $\beta$ (the angular height of the sun above the horizontal plane) simultaneously, using Coordinated Universal Time (UTC) on March 20th, 2023, during the spring equinox. The choice of the spring equinox is deliberate, aimed at simplifying the calculations, as the sun is positioned directly above the equator at this time. By incorporating these measurements with known terrestrial distances, its possible to derive a formula that can calculate the distance to the sun. This site was used to obtain the solar altitudes at different times at the two locations on Earth.

---

Mathematical derivation

---

I've drawn a figure in LaTex to explain the system and the measurements needed. The quantity of interest is the length $|CS|$.

Figure 1:

The notations are explained below.

1. $A$ is the first location on the equator
2. $B$ is the second location on the equator
3. $C$ is the center of the earth
4. $D$ is the intersection of $\overline{CS}$ and Earth's surface
5. $E$ is NOT NEEDED
6. $S$ is the position of the sun
7. $\alpha$ is the angle between $\overline{AS}$ and the Earth-tangent at point $A$
8. $\beta$ is the angle between $\overline{BS}$ and the Earth-tangent at point $B$
9. $\alpha'$ is the angle between $\overline{AB}$ and the Earth-tangent at point $A$
10. $\beta'$ is the angle between $\overline{AB}$ and the Earth-tangent at point $B$
11. $\varphi$ is the angle $\angle{BSA}$
12. $\theta$ is the angle $\angle{ABC} = \angle{BAC}$
13. $d$ is the distance between $A$ and $B$ along Earths surface (NOT IN FIGURE)
14. $r$ is the radius of Earth (NOT IN FIGURE)

CALCULATIONS

From figure 1, it's evident that $\alpha'=\frac{\pi}{2}-\theta$ and $\beta'=\frac{\pi}{2}-\theta$, so $\alpha' = \beta'$. The triangle $\triangle ACB$ is an Isosceles triangle since two of it's sides have the same length, $|AC|=|BC|=|DC|=r$. Thus $\gamma = \frac{d}{r}$ radians and $\theta = \frac{\pi-\gamma}{2}$. Moreover, $\angle{BAS} = \frac{\pi}{2}-\theta+\alpha$ and $\angle{ABS}=\frac{\pi}{2}-\theta+\beta$, so

$$\varphi = \pi-(\angle{BAS}+\angle{ABS})=2\theta-\alpha-\beta. \tag 1$$

Using the Law of Sines on the triangle $\triangle{ABS}$ we obtain the length $|AS|$ as

$$\frac{\sin{\varphi}}{|AB|}=\frac{\sin{\angle{ABS}}}{|AS|}\Longleftrightarrow |AS|=\frac{|AB|\sin{\left(\frac{\pi}{2}+\beta-\theta\right)}}{\sin{\left(2\theta-\alpha-\beta\right)}},\tag 2$$ where $|AB|$ is found using the Law of Cosines on the triangle $\triangle{ACB}$ as follows

\begin{align} |AB|^2 &= r^2+r^2-2r^2\cos{\gamma} \tag 3 \\ &= 2r^2-2r^2\cos{\gamma} \tag 4\\ &= 2r^2(1-\cos{\gamma}) \tag 5 \\ &\Leftrightarrow |AB|=r\sqrt{2(1-\cos{\gamma})}. \tag 6 \end{align}

It follows that

\begin{align} |AS| &= \frac{r\sqrt{2(1-\cos{\gamma})}\sin{\left(\frac{\pi}{2}+\beta-\theta\right)}}{\sin{\varphi}} \tag 7 \\ &= \frac{r\sqrt{2(1-\cos{\gamma})}\sin{\left(\frac{\pi}{2}+\beta-\theta\right)}}{\sin{(2\theta-\alpha-\beta)}}. \tag 8 \end{align}

Since $\theta=\frac{\pi-\gamma}{2}$ we can rewrite equation $(8)$ as

$$|AS|=\frac{r\sqrt{2(1-\cos{\gamma})}\sin{\left(\frac{\pi}{2}+\beta-\frac{\pi-\gamma}{2}\right)}}{\sin{\left(\pi-\gamma-\alpha-\beta\right)}}. \tag 9$$

At this stage, in triangle $\triangle{ACS}$ we have expressions for the side length $|AS|$ as well as the angle $\angle{CAS}=\frac{\pi}{2}+\alpha.$ We can now calculate the length of the line segment $|CS|$ using, again, the Law of Cosines. It follows that

\begin{align} |CS|^2 &= |AS|^2+r^2-2|AS|r\cos{(\triangle{CAS})} \tag {10} \\ &= |AS|^2+r^2+2|AS|r\sin{\alpha} \tag {11} \\ &\Leftrightarrow |CS| = \sqrt{|AS|^2+r^2+2|AS|r\sin{\alpha}} \tag {12}. \end{align}

Equation $(12)$ is the complete formula for calculating the distance $|CS|$ which is the sought distance. This expression is expressed in terms of the constant parameters $r,d$ and the solar altitude variables $\alpha,\beta$. We can therefore write the solar distance as a two dimensional function of $\alpha,\beta$ as

$$\boxed{f(\alpha,\beta) = \sqrt{|AS|^2+r^2+2|AS|r\sin{\alpha}},} \tag {13}$$

where

$$|AS|=\frac{r\sqrt{2(1-\cos{\left(\frac{d}{r}\right)})}\sin{\left(\beta+\frac{d}{2r}\right)}}{\sin{\left(\frac{d}{r}+\alpha+\beta\right)}}.\tag{14}$$

---

Implementation and results

---

The the solar angle data was meticulously collected from this site and the coordinates of the locations chosen were precisely $(0.00^{\circ} \ N \ \& \ 38.00^{\circ} \ E)$ in Kenya and $(0.00^{\circ} \ N \ \& \ 102.00^{\circ} \ E)$ in Indonesia.

The following Python function was implemented to compute the distances.

def TrigSunDistCalc(alpha, beta, r=6371, d=7022) -> float:
    alpha = np.radians(alpha)
    beta = np.radians(beta)
    AS = (r*np.sqrt(2-2*np.cos(d/r)) * np.sin(beta + d/(2*r)))/np.sin((d/r) + alpha + beta)
    return np.sqrt(AS**2 + r**2 + 2*AS*r*np.sin(alpha))

The data collected is listed below. The first three columns are the data collected, and the last column are the distance computations using formula $(13)$.

UTC timestamp            alpha (deg)      beta (deg)     f(alpha,beta) (km)
2023-03-20 04:00:00      70.09            6.23           12473.79
2023-03-20 04:15:00      73.84            9.93           14330.11
2023-03-20 04:30:00      77.59            13.66          17240.05
2023-03-20 04:45:00      81.34            17.39          22501.05
2023-03-20 05:00:00      85.09            21.13          35126.19
2023-03-20 05:15:00      88.81            24.88          107257.45
2023-03-20 05:30:00      87.39            28.62          401439.37
2023-03-20 05:45:00      83.65            32.37          420292.46
2023-03-20 06:00:00      79.9             36.12          432584.44
2023-03-20 06:15:00      76.15            39.86          437796.79
2023-03-20 06:30:00      72.41            43.61          451547.27
2023-03-20 06:45:00      68.66            47.36          458162.22
2023-03-20 07:00:00      64.91            51.11          462815.52
2023-03-20 07:15:00      61.16            54.86          465487.26
2023-03-20 07:30:00      57.41            58.61          466165.98
2023-03-20 07:45:00      53.66            62.36          464848.77
2023-03-20 08:00:00      49.91            66.11          461541.28
2023-03-20 08:15:00      46.16            69.86          456257.67
2023-03-20 08:30:00      42.41            73.61          449020.59
2023-03-20 08:45:00      38.66            77.36          439861.03
2023-03-20 09:00:00      34.92            81.11          434005.34
2023-03-20 09:15:00      31.17            84.85          415976.26
2023-03-20 09:30:00      27.42            88.59          396617.94
2023-03-20 09:45:00      23.68            87.64          63259.84
2023-03-20 10:00:00      19.94            83.89          29524.81
2023-03-20 10:15:00      16.2             80.14          20407.52
2023-03-20 10:30:00      12.46            76.39          16140.45
2023-03-20 10:45:00      8.74             72.64          13652.85
2023-03-20 11:00:00      5.05             68.89          12010.82

Below, the computed distances are plotted in a graph with varying UTC timestamps.

Figure 2:

---

Conclusions & Questions

---

Clearly my results differ a lot from today's accepted 1 AU. The distances are symmetrical around 07:30 UTC and there is a lot of variation, from roughly 12 000 km to 466 000 km. I want to explain why this is. I've listed some potential error sources:

1. Error in my calculations. However I want to believe that they are all correct unless I'm missing some detail. I've also asked 2 other people with math knowledge to verify the derivation and none of them could find any error.
2. Angular measurements and their precision. The accuracy of the angles of elevation (alpha and beta) is crucial. Small errors in angle measurement can lead to large errors in the calculated distance. This is especially true for trigonometric functions, which can vary significantly with small changes in the input angles. This I can't verify since I do not know how these angles were measured but at least they are given with two decimals. Could a slight variation/error in these cause the distance to jump from 12 000km to 150 000 000 km? Even if there are minor errors in the angles, it seems far fetched and quite unbelievable that the distance should suddenly drop down to 0.00008 %. If this is the case, I'd want to know how.
3. Atmospheric refraction. However accounting for this, the distances computed should be closer since the perceived position is further away (larger solar angle) than the true position and the $\alpha,\beta$ are calculated based on the perceived position according to the website. So atmospheric refraction can't be the answer to the large discrepancy.

Eventually, my question boils down to: why are my results so far away from 1 AU?

Answer 1

This is an interesting experiment / project, and props for your significant effort toward executing it!

The gist of your approach is trying to work out the distance using the solar altitude measured at two different longitudes on the equator. You did it for several different timestamps experimentally.

TL;DR you didn't exactly make any mistakes, but you made several problematic approximations that resulted in your formula breaking down. Moreover, the imprecision of the data you were working with, as well as approximation in the calculator made your attempt doomed anyway. However, a real-world test might be barely possible.

---

Data

To start our analysis, let's take the center of your time series at 2023-03-20 07:30 UTC. The equinox was at 21:24 UTC on that date, nowhere in your time series, a source of error which we'll return to later.

For 07:30 UTC, the site tells me that the altitude ("the angle between horizon and the center of the Sun including refraction") for the locations chosen was:

• Kenya 0°N, 38°E at 10:30 UTC+3 (07:30 UTC): 58.61° alt, 90.44° azi
• Sumatra 0°N, 102°E at 14:30 UTC+7 (07:30 UTC): 57.41° alt, 269.58° azi

These figures agree with the numbers you collected. One notices that the azimuths differ by about 180°. This is because 07:30 UTC was 10:30 local time in Kenya (i.e. in the morning), but 14:30 local time in Sumatra, i.e. in the afternoon. The altitudes are thus measured from opposite horizons. Note it's almost a degree off from 180° apart because we're not using the right time of day for the equinox; that discrepancy makes your trig inaccurate.

Your formula / code assumes the azimuths are opposed (the intermediate AS, can be (nonsensically) negative otherwise.), and that assumption breaking seems to be the source of the jumps in your graph, as discussed somewhat and diagrammed in the other answer. In the middle, we're fine, though.

---

Your Math and Code

I checked your derivation of the distance and didn't see an error as written. I feel like there's probably a simpler way, but your formula seems correct and the derivation very nice, so we'll go with it.

As for your code, you have r 6371 km, which is an average. The Earth is an ellipsoid, and since we're on the equator, it would be better and just as easy to instead use the equatorial radius of 6378.1370 km.

I don't know where you got the distance $d$ (var. AB), which you have listed as 7022 km. That's about 100 km too small (regardless of the radius used). The two points are separated by 102°-38°=64° of longitude, which is 64/360≈17.78% of the way around the equator, whose radius we now know precisely. A precise calculation gives 7124.447 km, corroborated here to the five places given.

This might seem like a small thing, but it's a major source of error, since the relatively small size of the Earth gets massively multiplied to get the AU, and so any error does too. Just to give a flavor, using your code, I indeed get a result of 466,165.98 km, but with the more correct r=6378.1370, AB=7124.447, I get 19,362,488.85 km, which is still wrong, but 40 times bigger (in the correct direction, it happens).

---

Thinking About Error

Digging into the error some more, we can start to see that this is really ill-posed. The locations are separated by 64° of longitude, as above. If the Sun were infinitely far away, its rays would be parallel and a 58.61° altitude in Kenya would mean a 57.39° altitude in Sumatra (from the opposite horizon). However, the calculator says we measure 57.41° altitude.

The discrepancy is 0.02°. That's also the unit in last place, so we've only got precision to the nearest 0.005°, which is . . . 25% of the value itself!

It's much, much worse than that, though. We didn't use the exact equinox, which means the Sun was not exactly on the Earth's orbital plane and we did not account for refraction at all. A big part of your inaccuracy seems to be with the refraction; when you don't account for it, you can compute a negative φ, which will of course give bogus results.

---

Improving the Result

Let's try to do a bit better.

First, all the data is coming from a calculator, so in principle we should be able to get a very accurate answer. In the following, I'll be using NOAA's calculator. While both calculators round their outputs and in SunCalc's case there is uncertainty in what formulae they use for refraction, SunCalc's code is minified and so is difficult to reverse-engineer, whereas NOAA's is human-readable. It will turn out that we need to use the high-precision internals of the calculator to understand what's going on.

To drive the error down, we want to compute at the exact equinox, which unfortunately does not seem to be distributed to higher precision than the minute. We also want to select points that are far apart, but not so far apart that the Sun is low in the sky and refraction changes a lot (though we absolutely still compensate for it).

At the correct equinox of 2023-03-20 21:24 UTC, the Sun was below the horizon in both locations (it was mostly over the Pacific Ocean). To keep with the spirit of an experiment you could do on land, let's choose a different equinox. I greatly prefer the season of Autumn, and the Autumnal Equinox of that year (2023-09-23 06:50) will do nicely, with the entire equator from Africa through Indonesia being suitable. I'll pick points a bit farther apart:

• Gabon 0°N, 10°E at 07:50 UTC+1 (06:50 UTC): 24.40° (24.396⋯°) alt, 90.00° azi
• Indonesia 0°N, 120°E at 14:50 UTC+8 (06:50 UTC): 45.66° (45.655⋯°) alt, 270.00° azi

Note that, because we're at the moment of the equinox, the azimuths are now parallel to the equator and separated by 180.00°, as we expect.

The parenthesized numbers are the first digits of the actual elevations the NOAA calculator is using extracted from the JavaScript debugger (their full values are 24.396176229009612° and 45.65502300416843°, respectively), while the unparenthesized numbers are the rounded values it displays.

The distance between the two longitudes is 120°-10°=110°, which is 110°/360°≈30.56% of the way around the Earth.

Using NOAA's formula for refraction here (see also their main.js:330) and inverting the 5°–85° band both zeniths fall into numerically, we get 24.361⋯° and 45.639⋯°. Though small, this difference is significant; it changes the difference between the elevations by almost 0.02°, which as we saw above is as big as the entire discrepancy we're trying to measure!

import math
import scipy

def refract_adjust_5to85deg( true_elev_deg ):
    #https://gml.noaa.gov/grad/solcalc/calcdetails.html
    T = math.tan( math.radians(true_elev_deg) )
    adjust_deg = (1.0/3600.0)*(
        58.1/T - 0.07/(T*T*T) + 0.000086/(T*T*T*T*T)
    )
    return true_elev_deg + adjust_deg
def refract_unadjust( apparent_elev_deg ):
    true_elev_deg = scipy.optimize.fsolve(
        lambda true_elev_deg:
            refract_adjust_5to85deg(true_elev_deg[0]) - apparent_elev_deg,
        apparent_elev_deg
    )[0]
    return true_elev_deg

alpha_deg = refract_unadjust(24.396176229009612)
beta_deg  = refract_unadjust(45.65502300416843 )
print(f"True α={alpha_deg:.5}°, β={beta_deg:.5}°")

---

A High-Precision Revelation

When we compute the formula, we get a large number, and checking why, we see that it's because φ is computed as 9.5716e-10° which, for the numerics we're using, is within the error of actual, literal 0.0. This tells us the problem: the calculators are not accounting for the solar distance, and we've just reconstructed their parallel-rays approximation! In fact, you can actually see this: NOAA does calculate the Sun distance (main.js:286), but it's pointless because that variable is never used again. For the two decimals of precision they needed, they decided not to bother (I guess they planned to initially, but changed their mind and forgot to remove that line). I don't know whether SunCalc shares this flaw, but I'd bet it does.

The implication is your calculation was always doomed to fail; the data you were working from simply did not include the effect you were trying to extract from it!

---

Getting a Real Answer

It seems that we can't currently demonstrate your approach from these calculators' data, but we can calculate what would kind of precision would be required.

Say the Sun is 1 au away and on some equinox we're considering two points that (to simplify things) are spaced the same longitude before and after the Sun at that instant—say one is 60° ahead and the other is 60° behind. The distance $\ell$ from each observatory to the Sun is given by law of cosines, and is related to the vertex angle at the Sun $\varphi$ by law of sines: $$ \ell:=\sqrt{ r^2 + (1\text{ au})^2 - 2r(1\text{ au})\cos(60^\circ) } ,\hspace{1cm} \frac{\sin(60^\circ)}{\ell} = \frac{\sin(\varphi/2)}{r} $$ The solution works out to be: \begin{align*} \varphi &= 2 \arcsin\!\left(\frac{\sqrt{3}}{2}\cdot\frac{r}{\sqrt{ r^2 + (1\text{ au})^2 - r(1\text{ au}) }}\right)\\ &\approx 0.004231^\circ \end{align*} (which is consistent with a simplistic estimation).

To measure the AU, we need measure this angle (via the elevation angles), which suggests immediately we need a precision of at least, like, 0.001° (and preferably a lot better), which is much worse than the calculators were giving (even besides that they weren't accounting for the distance at all).

We can get a more quantitative feel for this by, for example, solving for the AU and differentiating. The AU is: $$ 1\text{ au} = \frac{r+\sqrt{r^2-4c}}{2},\hspace{1cm}c:=r^2-\frac{3r^2}{4\sin^2(\varphi/2)} $$ The differential with respect to $\varphi$ over the interesting domain is: $$ \frac{d}{d\phi} 1\text{ au} = \frac{-\sqrt{3}}{8}\cdot \frac{ r \sin(\varphi) \csc^4(\varphi/2) }{ \cot{\varphi/2} } $$ For example, if you want 0.1%, 1%, or 10% error in your calculation, then you can estimate you need a precision of 0.015", 0.15", or 1.5", respectively. Note that you can differentiate with respect to other parameters (e.g. $r$ and the (currently hardcoded) separation of 120° longitude). A full error analysis would need to distribute the allowable error over such sources of inaccuracy as well.

The precision of JWST is 0.1", and it seems that some ground-based telescopes can get in the 1" range under ideal conditions (don't quote me on that; I don't actually know). Therefore, assuming you could find the center of the Sun with this accuracy, then it might actually possible to estimate the AU this way. I would say it would be difficult and out of the capabilities of an amateur astronomer, though.

---

Conclusion

Although your formula and code seem to be broadly correct, there were still several sources of error:

• You didn't use the correct time of the equinox, making the trig inaccurate.
• For some times you calculated, the elevation was on the wrong side, breaking the formula's assumptions.
• You didn't use the equatorial radius, and also you miscalculated the distance between the two points.
• You ignored refraction entirely, usually causing the formula's assumptions to break.
• The calculators round their outputs to a level that makes the calculation impossible.
• The calculators don't figure the Sun's distance, making it impossible even in theory.

The calculation is numerically ill-posed, and each of these sources of error individually doom the calculation to failure—so together they are more than sufficient. Some of them are easy to fix (e.g. using the better radius), and some are not (fixing NOAA's calculator).

We can fix the errors related to the setup, and instead of trying to use the calculators (which aren't in the spirit of an actual "terrestrial measurement" anyway), we could imagine trying to do a real-world measurement. One would need two very accurate telescopes taking a somewhat difficult measurement exactly at the equinox on the equator (or correcting for any differences mathematically). Technologically I think it's at the edge of feasible.

Answer 2

For the first hour and the last 75 minutes, it appears you neglected to take into account the direction of the Sun. As shown on this first image, you are measuring the Sun as if it were along the black line stemming out of the top point; this will give a Sun at the intersection of the two black lines, which is about 12473.79 km away from the bottom point, indeed. But the Sun is then not in the direction of the black line, but in the direction of the red line. As you can see, it is very nearly parallel to the black line from the bottom point at that moment.

For the “middle” part of your experiment, let’s say from 5:30 to 9:15, there must be something wrong in your coding or formulas, as the following graphic, simulating the situation, shows the direction of the Sun from one location is parallel to the direction of the Sun from the other location. (Click on image to enlarge it.)

In both drawings, the dashed grey lines are the local horizons.