Does SearchRiseSet take elevation of observer into account?

[bjarkeandersen]

I tried to send in an observer with the same latitude and longitude, but with different elevations (between 0 and 15000). But I keep getting the same timestamp for the next sunrise.

Here is my test code (c#): for (int i = 0; i <= 15000; i += 1500) { AstroTime astroTime = new(2022, 11, 21, 0, 0, 0); Observer observer = new(67.14, 14.2, i); var sunrise = Astronomy.SearchRiseSet(Body.Sun, observer, Direction.Rise, astroTime, 1); System.Diagnostics.Debug.WriteLine($"Height: {i}, sunrise: {sunrise}"); }

And here is the result: Height: 0, sunrise: 2022-11-21T08:29:28.948Z Height: 1500, sunrise: 2022-11-21T08:29:28.949Z Height: 3000, sunrise: 2022-11-21T08:29:28.950Z Height: 4500, sunrise: 2022-11-21T08:29:28.951Z Height: 6000, sunrise: 2022-11-21T08:29:28.951Z Height: 7500, sunrise: 2022-11-21T08:29:28.952Z Height: 9000, sunrise: 2022-11-21T08:29:28.953Z Height: 10500, sunrise: 2022-11-21T08:29:28.953Z Height: 12000, sunrise: 2022-11-21T08:29:28.954Z Height: 13500, sunrise: 2022-11-21T08:29:28.955Z Height: 15000, sunrise: 2022-11-21T08:29:28.955Z

Is this expected behaviour (there is a small difference, but not more that a few miliseconds) or is something wrong with my code?

[cosinekitty]

Hi @bjarkeandersen, yes this is expected behavior. The reason the elevation has such a small effect on sunrise/sunset times is that the simulation models the observer standing on a locally flat plateau at the same elevation as the observer. This may cause a slight change in parallax with the Sun's apparent position in the sky, although I doubt that would cause even a millisecond worth of difference. Elevation would have more effect on calculations involving the Moon, because it is much closer to the Earth.

I suspect the main effect you are seeing is an artifact of the solver having a tolerance of 0.1 seconds, and different observer elevations result in slightly different initial conditions for its iterative search.

In practice, rise/set times are impossible to predict perfectly because of variations in atmospheric refraction. Generally we expect rise/set times to have an error within about a minute or so, due to uncertainty of air pressure and temperature variations at different altitudes that the sunlight (moonlight, etc) must pass through on the way to the observer's eye. Higher pressures and lower temperatures create denser air that leads to more refraction, and vice versa. So the 0.1 second tolerance is far smaller than any possible accuracy.

I could probably increase the average accuracy by modeling the atmosphere as thinner at higher elevations, but I don't currently do that.

I hope this helps!

[bjarkeandersen]

Aha - but if you levitate 10 km. above the Earth (the flat plateau), then there would be a significantly difference I guess (about 1 minut/1500 meters at equator according to several google hits). Do you have a method in your library that can calculate this?

[cosinekitty]

This does make sense as an enhancement. It should be possible, but I don't want to change behavior or break compilation for existing callers. Thus it may require adding a new function. I will think about this over the weekend and reply back here.

[cosinekitty]

I keep thinking about this, but I run into a problem. I have a policy of not adding new functionality to Astronomy Engine until I have authoritative test data to confirm its calculations. With increased height above the Earth, not only does the horizon sink lower, but the atmospheric refraction decreases. So far I have not been able to find scientifically verified data for these effects. I suspect this is because there are so many meteorological variables, so there just isn't a simple answer. Consider that the density of air changes as a nonlinear function of height as shown in this diagram.

One thing you might try is to call SearchAltitude. This function can work as a sunrise/sunset calculator by making adjustments for atmospheric refraction and the apparent angular radius of the Sun.

For example, to calculate sunrise with the observer standing on flat ground, you can use

const double sunRadiusAu = 0.0046504672612422675;
double sunAngularRadius = Astronomy.RAD2DEG * Math.Asin(sunRadiusAu / Astronomy.HelioDistance(Body.Earth, startTime));
const double refractionNearHorizon = 34.0 / 60.0;

AstroTime time = Astronomy.SearchAltitude(
    Body.Sun,
    observer,
    Direction.Rise,
    startTime,
    1.0,
    -(sunAngularRadius + refractionNearHorizon)
);

This gives the same answer as SearchRiseSet for the Sun, to millisecond precision.

That last parameter defines the target altitude angle to solve for. We subtract the Sun's apparent angular radius because sunrise/sunset are defined by seeing the topmost portion of the Sun, not its center. We subtract a typical refraction angle of 34 arcminutes, because refraction makes objects appear higher than they actually are.

In your case, you would have to include an additional angular correction in that last parameter based on the ground being lower than the observer. You would also have to use a different correction for atmospheric refraction. I know how to do the first correction, but not the second.

To correct for the ground being lower, if $R$ is the radius of the Earth, $h$ is the height of the observer above the ground, and $\theta$ is the decreased altitude angle of the horizon, then

$$\cos \theta = \frac {R} {R+h}$$

You can use $R$ = 6371000 meters as the Earth's mean radius, solve for $\theta$ by taking the inverse cosine of both sides, and subtract the result from the altitude angle in the last parameter as mentioned above.

I hope this helps.

[cosinekitty]

I just edited the formula above to correct a mistake. I had incorrectly written $\sin \theta$ where it should have been $\cos \theta$. Sorry about that!

[bjarkeandersen]

Thank you - I didn't have time to look at this yet, but it seems like a good solution to solve the problem. It is of course a bit of a theoretical problem because the most places on earth you don't have a mountain to stand on in an otherwise flat landscape. Would be interesting though, to combine this with a digital elevation model to calculate the exact sunrise/sunset for a coordinate on earth. But I think I'll leave that as an exercise for others to test out.

[hidp123]

I keep thinking about this, but I run into a problem. I have a policy of not adding new functionality to Astronomy Engine until I have authoritative test data to confirm its calculations. With increased height above the Earth, not only does the horizon sink lower, but the atmospheric refraction decreases. So far I have not been able to find scientifically verified data for these effects. I suspect this is because there are so many meteorological variables, so there just isn't a simple answer.

Hi, not sure if it helps but JPL Horizons has an option to set observer altitude. And based on that it will give rise/set times (to the nearest minute, however, there is a workaround to get results to the nearest second as well).

[cosinekitty]

Hi, not sure if it helps but JPL Horizons has an option to set observer altitude

Cool! I will try that out and see whether they also model rise/set on a flat plateau (i.e. the height is a ground elevation), or whether they model the observer flying over the ground (the height is an altitude). I suspect the former, but I will report back here whatever I find.

[cosinekitty]

I don't see where JPL Horizons determines sunrise or sunset times, but it does report altitude angles. With that in mind, I tried out the experiment where I use the same latitude and longitude, but change the elevation. When I use an elevation of 0 meters and ask for the Sun's azimuth and altitude angles, I get:

2023-Mar-09 16:03 *   142.541748  50.168002

When I change the elevation to 2000 meters, I get

2023-Mar-09 16:03 *   142.541748  50.168001

Because the change is only one millionth of a degree in altitude angle, I can conclude JPL Horizons is treating the ground at 0m or 2000m above sea level as a flat plane. The tiny angular distance is parallax of a triangle that is 2000 meters wide and 1 astronomical unit long! Here I calculate 2 kilometers divided by the distance between the Earth and Sun in kilometers, then compute the resulting angle in degrees:

>>> from math import atan, degrees
>>> degrees(atan(2.0 / 149597870.69098932))
7.659972598331028e-07

It shows a result very close to a millionth of a degree. This isn't really very surprising because I would expect altitude angles to always be relative to the observer's horizontal plane, not relative to the observer's sea level.

For now I'll stick with my cosine formula above to correct for an observer floating above the ground.

[hidp123]

Hi, I'll show a working example for sunrise below:

First an example for 0m altitude. Target body: Sun. Location: 40N, 2W. Altitude: 0m. Start time: 2023-03-21 Stop time: 2023-03-22 Step size: 1 (minutes)

Now, for step 5 select the following:

From the above the important one is "RTS flag" (towards the bottom of the screenshot), this will output results only at rise/transit/set times (to the nearest minute, I will show how to narrow it down to the second below). "Refraction model" can be either, I have set it to "airless" (ignore the setting shown in the screenshot). I've noticed the rise/set time does not change due to this, only the elevation value changes.

Result:

Rise is at 06:11 Set is at 18:21

Now, for sunrise, to narrow the results to the second do the following: Start time: 2023-03-21 06:10 (set time to 18:20 for sunset) Stop time: 2023-03-21 06:12 (set time to 18:22 for sunset) Step size: 120 (equal intervals) - as there are 120 seconds in 2 minutes. In step 5, set "RTS flag" to disabled. Then generate output.

Result: 2023-Mar-21 06:10:26.000 C 89.107575647 -0.841392201 2023-Mar-21 06:10:27.000 C 89.110251694 -0.838197146 2023-Mar-21 06:10:28.000 C 89.112927736 -0.835002088 2023-Mar-21 06:10:29.000 *r 89.115603773 -0.831807028* (the "r" marker denotes sunrise) 2023-Mar-21 06:10:30.000 89.118279806 -0.828611966 2023-Mar-21 06:10:31.000 89.120955833 -0.825416901 2023-Mar-21 06:10:32.000 89.123631856 -0.822221834 Thus, for 0 km altitude, sunrise is at: 06:10:29**

---

Now let's repeat the above procedure for 1 km altitude. We get the following for sunrise: 2023-Mar-21 06:05:08.000 C 88.256252080 -1.857264840 2023-Mar-21 06:05:09.000 C 88.258930542 -1.854070876 2023-Mar-21 06:05:10.000 C 88.261608994 -1.850876907 2023-Mar-21 06:05:11.000 *r 88.264287436 -1.847682933* 2023-Mar-21 06:05:12.000 88.266965867 -1.844488956 2023-Mar-21 06:05:13.000 88.269644289 -1.841294973 2023-Mar-21 06:05:14.000 88.272322699 -1.838100986 Thus, for 1 km altitude, sunrise is at 06:05:11**. About a 5 min difference from sea level altitude.

Hope this has helped.

[hidp123]

Here are some websites or programming libraries that have an option to set observer altitude and/or horizon elevation (although horizon elevation is in degrees in some cases), also not sure if some of these would fall under authoritative data for your tests:

• Skyfield (python) - https://rhodesmill.org/skyfield/almanac.html#elevated-vantage-points (observer & horizon altitude, horizon altitude is in degrees, however, a formula is given to calculate the degrees value).
• JPL Horizons (observer altitude only).
• https://www.largeformatphotography.info/sunmooncalc/ ("elevation" & "height above horizon" fields, "elevation" is observer elevation, "height above horizon" is the difference between the observer's altitude and horizon's altitude in meters/ft).
• Swiss Ephemeris - https://www.astro.com/ftp/swisseph/doc/swephprg.htm#_Toc112948995. This has an option to input observer altitude in meters, and horizon altitude in degrees.
• USNO - https://aa.usno.navy.mil/data/mrst, here it gives the option to input observer height, but it also increases horizon height by that much, horizon dip is not considered, so values don't change unless one is at the poles (see "height" under notes section at the bottom)
• https://suncurves.com/en/ - paid api for additional data though.
• https://www.stjarnhimlen.se/comp/riset.html - some calculations and formulas.
• https://aty.sdsu.edu/explain/atmos_refr/dip.html
• https://books.google.com/books?id=uJ4JhGJANb4C...
• https://www.sciencedirect.com/science/article/pii/S0098300422000115

[cosinekitty]

These links should help provide trustworthy test data. Thank you!

I've been thinking about this, and because level ground is a special case, the most general solution would require the caller to pass in an angular adjustment. For example, you would pass in a negative angle if the horizon is below the observer (perhaps he is on a mountain), or a positive angle if the horizon is above the observer (he is in a valley).

Then for the special case where the observer is, say, in a balloon above flat ground, I can provide an extra function for calculating that angle. It will be the cosine formula I mentioned above.

I will need an additional function like SearchRiseSet in C, to avoid breaking backward compatibility. The other languages all allow optional parameters and/or function overloading, so it should not be a problem there.

I'm in the middle of working on some other stuff, so it may be a week or two before I get to this, but I will do it.

[cosinekitty]

I thought I would take a moment to report where I am in this effort.

I'm currently working through the issues surrounding the horizon dip angle. It turns out that my simple geometric theory is not adequate because atmospheric refraction also affects the path of light rays from the horizon to the observer. I found a formula in the largeformatphotography.info page's JavaScript code, but it is not well referenced. It looks reasonable, but I have no way of testing it. I would be a lot happier if I can find some kind of corroboration for it.

Currently, I'm trying to do that by looking at some terrestrial refraction formulas. It's really fun and interesting stuff, but it's going to take a while to work through it and convince myself it's accurate.

As a bonus, I will likely be adding a new function to Astronomy Engine for calculating the U.S. 1976 Standard Atmosphere model, which models air temperature, pressure, and density at different elevations above mean sea level. These will be required for calculating refraction curves of light rays, so I might as well expose them to the API.

[hidp123]

Thanks for the update.

Previously, I did add a few extra links to my initial list. Not sure if they are of any help though, but some of them do go into some detail regarding dip etc.

[hidp123]

Currently, I'm trying to do that by looking at some terrestrial refraction formulas.

From the above link: "The Source Code of the simulation App alone, excluding the code of the modules for the graphic output, control panels ect., consists of about 5800 lines."

Wow!

[cosinekitty]

Yes, it is pretty involved. But fortunately, on that same page I eventually find this formula, which goes to the core of the issue:

$$ { 1 \over r } = -{ cos(\beta) \over n } \cdot { \mathrm{d} n \over \mathrm{d} h } $$

I'm trying to follow the derivation, but it relates the radius of curvature of a refracted light ray to the derivative of the index of refraction with respect to distance.

[cosinekitty]

This work has been completed in v2.1.17 at bd2db6a3805ac8a7c559b6b2276e16c1e1793d1f. The SearchRiseSet functions now accept a new parameter metersAboveGround, which defaults to 0. There is also a new Atmosphere function that provides idealized air temperature, pressure, and density as a function of elevation, based on the U.S. Standard Atmosphere model of 1976.