|
4 | 4 | # as pointed out in GitHub issue #842. As we 'add nines' to 89.0°, here |
5 | 5 | # is how the quantities behave: |
6 | 6 | # |
7 | | -# aC: this looks to be the 'local radius of curvature', of all things, |
8 | | -# and it converges quickly to 6399593.625758493 m at 89.99999°, and |
9 | | -# stays there, rock-solid, all the way to 90°. |
10 | | -# |
11 | | -# R / cos(lat): this is our problem, the value that is going sideways. |
12 | | -# Converted to meters (* AU_M), it starts at a value very close to |
13 | | -# aC, then begins to diverge - with `R / cos(lat) - aC` growing by |
14 | | -# 10× for each step from 89.99 through 89.999999999 before its |
15 | | -# growth slows - but by that point the error in elevation is already |
16 | | -# 40 m. So let's investigate the behavior of its components, then |
17 | | -# return to this value below. |
| 7 | +# aC: this is the 'local radius of curvature', which converges quickly |
| 8 | +# to 6399593.625758493 m at 89.99999°, and stays there, rock-solid, |
| 9 | +# all the way to 90°. |
18 | 10 | # |
19 | 11 | # R: The radius of the line in the xy-plane connecting the pole with |
20 | | -# the end of the distance vector. |
21 | | -# |
22 | | -# cos(lat): The scaling factor, at the angle `lat`, between the radius |
23 | | -# of the unit circle, and a line drawn orthogonally from the pole's |
24 | | -# line over to the tip of the position vector. |
| 12 | +# the end of the distance vector, which (obviously) becomes very |
| 13 | +# very small as we approach the pole, though machine precision seems |
| 14 | +# to prevent its ever becoming zero. |
25 | 15 | # |
26 | | -# R / cos(lat): So, we return to our ratio, now with the understanding |
27 | | -# that it's trying to rebuild a kind-of-length for the position |
28 | | -# vector, but using the computed `lat` instead of the simple true |
29 | | -# angle between the position vector and the pole's vector. So how |
30 | | -# does this compare with the vector's true length? |
| 16 | +# cos(lat): Alas, fairly imprecise, as the only way it can tell exactly |
| 17 | +# how far we are from the pole is to look at the very last digits of |
| 18 | +# angles like 89.9999999, which only have a few bits left beneath |
| 19 | +# them to hold their precision. |
31 | 20 | # |
32 | | -# R / cos(lat) - length_of(xyz_au): This difference starts to converge |
33 | | -# on 42841.31 before it starts to diverge at around 89.999999. So |
34 | | -# `R / cos(lat)` exceeds the actual vector length, because . |
| 21 | +# R / cos(lat): So this is our problem, the value that is going |
| 22 | +# sideways. For the Earth's surface at the North Pole, converted to |
| 23 | +# meters (* AU_M), it starts at a value very close to aC, then |
| 24 | +# begins to diverge - with `R / cos(lat) - aC` growing by 10× for |
| 25 | +# each step from 89.99 through 89.999999999 before its growth slows |
| 26 | +# - but by that point the error in elevation is already 40 m. |
35 | 27 | # |
36 | | -# Looking inside the _compute_latitude() routine, here are a few values: |
| 28 | +# How can we compute the height without recourse to the cosine of the |
| 29 | +# latitude? Drawing a diagram, we have a triangle with one leg that's |
| 30 | +# the axis through the poles; one leg that's R, sticking out from the |
| 31 | +# axis into the xy-plane; and a hypotenuse that's the radius of |
| 32 | +# curvature plus the height above-ground-level of our target. |
37 | 33 | # |
38 | | -# z: The literal, unchanged z-coordinate of the position vector, giving |
39 | | -# height above the xy-plane of the equator. |
| 34 | +# R / cos(lat) is one way of computing the hypotenuse, yes. |
40 | 35 | # |
41 | | -# aC * e2_sin_lat: The distance in the z-direction between the |
42 | | -# position's actual z-coordinate, and the z-coordinate that the |
43 | | -# position would have if it were on a sphere whose radius, the whole |
44 | | -# way around, was the local ellipsoid's radius at the Earth latitude |
45 | | -# under consideration. |
| 36 | +# But we can also just use Pythagoras. We already know one leg is R. |
| 37 | +# The other? It turns out that it's the quantity passed to arctan2() to |
| 38 | +# compute the latitude! Thus, `z + aC * e2_sin_lat`, which we can give |
| 39 | +# a name to and pass as an additional return value. |
46 | 40 | # |
47 | | -# z + aC * e2_sin_lat: Per the previous two descriptions, this is the |
48 | | -# position's z-coordinate adjusted 'upward' (north along the z-axis |
49 | | -# for the northern hemisphere) to simulate the shape the vector |
50 | | -# would have if the whole Earth had the radius of curvature of this |
51 | | -# position on the ellipsoid. And so it's this artificial z-value |
52 | | -# that can be tossed into arctan2() with the literal xy-radius `R`. |
53 | | -# |
54 | | -# A few quick constants that will come up in various print() calls you |
55 | | -# might add to the code: |
| 41 | +# Some common values: |
56 | 42 | # |
57 | 43 | # 6399.6 km - Earth radius of curvature at pole |
58 | 44 | # 6378 km - Earth radius at equator |
|
62 | 48 | from numpy import arctan2, cos, sin, sqrt |
63 | 49 | from skyfield.api import load, wgs84, tau |
64 | 50 | from skyfield.constants import AU_M, DEG2RAD, RAD2DEG |
65 | | -from skyfield.functions import length_of |
66 | 51 |
|
67 | 52 | ts = load.timescale() |
68 | 53 | t = ts.utc(2023, 3, 2, 4, 13, 0) |
|
76 | 61 | # latitude = 1.0 - 1.0 / 10.0 ** nines # To test behavior at equator. |
77 | 62 |
|
78 | 63 | position = wgs84.latlon(latitude, 0.0, 0.0).at(t) |
79 | | - xyz_au, x, y, aC, R, lat = self._compute_latitude(position) |
| 64 | + xyz_au, x, y, R, aC, hyp, lat = self._compute_latitude(position) |
80 | 65 | _, _, z = xyz_au |
81 | 66 |
|
82 | | - # Official formula. |
83 | | - # Jumps around, because aC is jumping around. |
| 67 | + # Official formula: jumps around, because R / cos(lat) is jumping |
| 68 | + # around, as both R and cos(lat) are becoming very small. |
84 | 69 | height_au = R / cos(lat) - aC |
85 | | - |
86 | | - # So, an alternative? |
87 | | - R2 = R |
88 | | - lat2 = lat |
89 | | - |
90 | | - R2 = max(R, 11 / AU_M) |
91 | | - lat2 = arctan2(z + aC * self._e2 * sin(lat2), R2) |
92 | | - #lat2 = arctan2(z + aC * self._e2 * sin(lat2), R2) |
93 | | - |
94 | | - height_au = R2 / cos(lat2) - aC |
95 | 70 | height_m = height_au * AU_M |
96 | 71 |
|
97 | | - out = lat * RAD2DEG |
98 | | - print( |
99 | | - '{:<17s} {:<17s} {:20.12f}'.format(str(latitude), str(out), height_m), |
100 | | - # aC * AU_M, # Converges quickly on 6399 km, the radius of curvature. |
101 | | - # R * AU_M, # Dives towards zero. |
102 | | - # cos(lat), # Dives towards zero. |
103 | | - (R / cos(lat)) * AU_M, # Goes all wonky and jumps around. |
104 | | - z * AU_M, # Converges quickly on 6357 km, the actual polar radius. |
| 72 | + # So, new formula, building the hypotenuse with Pythagoras. |
| 73 | + thing = sqrt(hyp * hyp + R * R) |
| 74 | + height2_au = thing - aC |
| 75 | + height2_m = height2_au * AU_M |
105 | 76 |
|
106 | | - # Difference should be around 6399 - 6357 = ~42 km. |
107 | | - # (R / cos(lat) - length_of(xyz_au)) * AU_M, |
108 | | - ) |
| 77 | + out = lat * RAD2DEG |
| 78 | + print('{:<17s} {:<17s} {:20.12f} {:20.12f}'.format( |
| 79 | + str(latitude), |
| 80 | + str(out), |
| 81 | + height_m, |
| 82 | + height2_m, |
| 83 | + )) |
0 commit comments