# Calculating Earth-Sun distance only using trigonometry and terrestrial measurements

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:
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?

• What order of magnitude is the angle phi - angle ASB. I suspect it is so small that you need much more accurate altitude "measurements" in order for it to work. (And most likely the theoretical calculated altitudes do not consider the distance to the Sun.) Commented Mar 16 at 21:06
• Something is wrong with your message, where it says “The solar angle data was meticulously collected from this site”: “this site” links to “about:blank#blocked". Checking the code of your post puts the link as your Python code… Commented Mar 16 at 21:10
• Somewhat related question (by me) to compute the earth-sun distance from earth which unfortunately does not work in practice: astronomy.stackexchange.com/questions/35130/… Commented Mar 17 at 10:01
• Most sites do not calculate the topocentric position of the Sun. Most are geocentric. The angle ASB is smaller than 0.01 degrees, so you need data that is much more accurate than that. This is probably why your technique has not been used in the past. Commented Mar 17 at 18:10
• @Parseval "Since the others are within 2 decimal places, γ is as well." It's not decimal places that matter, but significant figures. To take a simple counter example, let's say I know a and b to 4 significant figures each: a = 1.235... and b = 1.234.... The error in each is <0.1%. But if we calculate a-b = 0.001... we only have 1 significant figure, and the error rate is 60%! (the true value could be anywhere between 0.005, 0.02). This, fundamentally, is the reason why the accuracy needed for this to work is so much higher than first appears. Commented Mar 18 at 13:01

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.

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).

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

58.1/T - 0.07/(T*T*T) + 0.000086/(T*T*T*T*T)
)
true_elev_deg = scipy.optimize.fsolve(
lambda true_elev_deg:
apparent_elev_deg
)[0]
return true_elev_deg

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!

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.

• We can get a more precise equinox time (2023-Sep-23 6:49:52 UTC) from JPL Horizons. Here's a declination plot using my script from this answer. The plot uses a numeric STEP_SIZE of 60 to get a 1 second time step. At this step size, quantization artifacts are visible. (The JPL DE stores positions as 14th order Chebyshev polynomials). Technically, equinox is defined in terms of RA, but for this task we obviously want the time when declination is 0°. Commented Mar 18 at 4:19
• Here's a Horizons query for RA, Dec, Azimuth, Elevation, and range & range-rate. Note that the site coords are given as 'lon, lat, altitude', with altitude in km. To get data using a simple refraction model, add &APPARENT=REFRACTED to the query URL. Commented Mar 18 at 4:29
• As mentioned here it's possible for amateurs to measure star deflections during a solar eclipse with sub-arcsec precision, using good observations & clever statistical techniques. But I expect that precisely measuring the position of the Sun's centre is a bit tricky due to its brightness. Commented Mar 18 at 4:52
• Wonderful answer, this was exactly what I was looking for. Thank you so much! I've accepted this answer but I might ask some questions about some details when I have time. Respect & Take care. Commented Mar 19 at 7:39

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.

• Thank you very much for your great answer. Yes indeed, for the first six and last six timestamps, I have not taken the suns direction into account, thank you for bringing this to my attention, I will write an edit to my post later. I've also fixed the links in my answer. Regarding the "middle part", what value did you use for the surface-distance between points A,B? I see that a slight variation in $d,\alpha$ and/or $\beta$ can significantly change the distance $|CS|$. Do you know where I can get more accurate numbers on these? ($r,d,\alpha,\beta$) Commented Mar 17 at 9:55
• When it comes to "error in your coding or formulas" - I have not been able to find it. That was my first thought as well: "Oh I've must have made a calculation or coding error". My hope was that someone here can spot the mistake, or verify there is none. If there is none, then there has to be another explanation for the $\approx 450 000 km$ vs the $\approx 150 000000 km$ discrepancy. Btw, what tool did you use to draw that last figure? Commented Mar 17 at 9:58
• @Parseval: The angle between A and B in both of my figures is (102-38) = 64°, from the values you give. I used Adobe Illustrator. Commented Mar 17 at 10:09