-1
$\begingroup$

In How are the three character IAU Minor Planet Center's Observatory Codes assigned? Why are some letters so popular? I looked at the distribution of the first alphanumeric character versus the value of the second two (36 x 100).

The table lists longitudes but no latitudes for these observatories, but using the cos and sin parallax coefficients in $\arctan2(\sin, \cos)$ and plotting as a function of tabular longitude produces a map that looks like it could be the distribution of telescopes on Earth.

Question: What exactly are the cos and sin parallax coefficients, and how closely can $\arctan2(\sin, \cos)$ be expected to serve as a proxy for latitude? Within a degree or so?

arctan2(sin, cos) as proxy for latitude

import numpy as np
import matplotlib.pyplot as plt
# blob is the unformatted lines from https://www.minorplanetcenter.net/iau/lists/ObsCodes.html
threes = [line[:3] for line in blob.splitlines()]
key = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
a = [[x] for x in key]
for t in threes:
    a[key.find(t[0])].append(t)
b = np.zeros((36, 100))
for i, thing in enumerate(a):
    c = [int(x[1:]) for x in thing[1:]]
    for d in c:
        b[i, d] = 1
plt.imshow(-b, interpolation='nearest', cmap='gray')
plt.gca().set_aspect(2)
plt.show()
pairs = []
todegs = 180/np.pi
for line in blob.splitlines():
    try:
        cos, sin = [float(x) for x in (line[13:21], line[21:30])]
        lat = todegs * np.arctan2(sin, cos)
        lon = float(line[4:13])
        pairs.append([lon, lat])
    except:
        pass
print(len(pairs), len(blob.splitlines()), float(len(pairs))/len(blob.splitlines()) )
lon, lat = np.array(list(zip(*pairs)))
lon = np.mod(lon+180, 360) - 180
plt.figure()
plt.plot(lon, lat, 'ok', markersize=1)
plt.ylim(-90, 90)
plt.xlim(-180, 180)
plt.gca().set_aspect('equal')
plt.show()
$\endgroup$
  • $\begingroup$ I know you are bright. I am seeing it as not acting in good faith. Again, sorry for being so blunt. The quality of the questions has decreased markedly since the point value was doubled, and I am admittedly taking that out on you. $\endgroup$ – David Hammen Mar 20 at 3:48
  • 1
    $\begingroup$ @DavidHammen do it the right way then; do it in meta and prove your point quantitatively in some way. As explained here I can not help the rep situation. If I could add a cap, cut it in half, or otherwise reduce it I would. I hate all my rep points above 10k but I can not stop them. I give them away as fast as I can in bounties. Now make your accusations in meta please and demonstrate that harm is being done, rather than repeatedly targeting me with accusations in comments across SE. Thanks! $\endgroup$ – uhoh Mar 20 at 3:52
  • $\begingroup$ @DavidHammen Can I suggest as an act of good faith linking to one of the answers you think is a duplicate on another SE ? Deleting the offending comment would be nice also, thank you. $\endgroup$ – StephenG Mar 20 at 8:26
  • 1
    $\begingroup$ @StephenG -- Right here on this corner of the SE network, there's the question What is 𝜙′ in orbital mechanics? $\endgroup$ – David Hammen Mar 20 at 10:07
  • $\begingroup$ @DavidHammen I've posted an answer. I think I've got it right, and the answer is about +/- 0.2 degrees. How does it look? $\endgroup$ – uhoh Mar 22 at 3:09
0
$\begingroup$

The title of the question reads:

Can I use the parallax coefficients for observatories as a proxy for latitude using arctan2(sin, cos)?

and the body expands:

Question: What exactly are the cos and sin parallax coefficients, and how closely can $\arctan2(\sin, \cos)$ be expected to serve as a proxy for latitude? Within a degree or so?

First the "How closely can (they) serve as a proxy for latitude?" part

tl;dr: to within +/-0.2 degrees; the range in the differences between geocentric and geodetic latitude.

The table provides the sin and cos of the same angle $\phi'$ multiplied by $r/R_E$ and the 2nd half of the answer will explain that $\phi'$ is the angle that a radius drawn from the observatory to the Geocenter makes with respect to the Earth's equatorial plane. Everyday latitude or GPS latitude which is a kind of geodetic latitude (what I want to use for my map) is a little bit different because it's drawn on an ellipsoid.

How different?

Well it turns out that I wrote an answer about geodetic latitude about 18 months ago. (@DavidHammen wrote a nicer answer elsewhere as well) Here is an excerpt from my answer:

The 3D cartesian coordinates $X, Y, Z$ in Earth-centered, Earth-fixed coordinates assuming an ellipsoidal shape is given by:

$$X = \left(N(\phi) + h \right) \cos\phi \cos\lambda $$

$$Y = \left(N(\phi) + h \right) \cos\phi \sin\lambda $$

$$Z = \left(\frac{b^2}{a^2} N(\phi) + h \right) \sin\phi $$

where $\phi, \lambda, h$ are latitude, longitude, and altitude, and $a, b$ are the equatorial and polar radii of the ellipsoid used, and

$$N(\phi) = \frac{a^2}{\sqrt{a^2\cos^2\phi + b^2 \sin^2\phi}}. $$

Lat, lon, alt in "GPS coordinates} is based on WGS 84 with $a, b$ of 6378.1370 and 6356.7523 kilometers, respectively.

The parallax coefficients are $\rho \cos(\phi')$ and $\rho \sin(\phi')$ and $\rho$ is the distance from the observatory to the Geocenter divided by "the Earth radius". That makes them:

$$\rho \cos(\phi') = \frac{\sqrt{X^2+Y^2}}{R_E}$$

$$\rho \sin(\phi') = \frac{Z}{R_E}$$

So if we chose a geodetic latitude $\phi$, choose a longitude of zero so that $Y=0$ and calculate $X, Z$ and plug them into this:

$$\phi' = \arctan2\left(Z, X \right) = \arctan2\left(\frac{Z}{R_E}, \frac{X}{R_E} \right) = \arctan2\left( \rho \sin(\phi'), \rho \cos(\phi') \right)$$

we can get the geocentric latitude, and then calculate the error by subtracting the two.

It turns out to be about +/- 0.2 degrees with extrema at about +/- 45 degrees. It is of course zero at the equator and poles.

So to an error of 0.2 degrees, yes that can indeed be used as a proxy!

geocentric minus geodetic latitude

Now the "What exactly are..." part

As pointed out in comments, the link explains how the cos and sin parallax coefficients are defined:

The following list gives the observatory code, longitude (in degrees east of Greenwich) and the parallax constants (rho cos phi' and rho sin phi', where phi' is the geocentric latitude and rho is the geocentric distance in earth radii) for each observatory. It is updated nightly.

What is geocentric latitude? It's the "easy one" to understand; it's just the angle measured at the Geocenter between the radius to the point and the closest point on Earth's equator. So while the definition *does not specify a specific Earth radius, it doesn't have to because the table only gives the unit-less parts cos phi' and sin phi'.

But why both sin and cos? Can't you just use sin between -90 and +90 degrees?

To just define an angle (without the radius), mathematically you might think so. However sin is insensitive near the poles and in the world of real computing numerical errors can crop up. Near the poles sin will still provide the appropriate sign and cos will provide the sensitivity. This way you will have low numerical errors (for this particular step) everywhere.

But this is the latitude and distance together, so you have to have two numbers. They could be a latitude angle and a radius explicitly, but this also works. (revised per @DH's comment)

But what are the parallax coefficients for exaclty? What are they really?

I don't know, but hopefully someone with a valid "astronomer's ID card" will come along and post an answer. I'm to scared to ask a new question on this topic given all the comments under the question! :-)


Script for the plot

import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
a, b = 6378.1370, 6356.7523 # only their ratio matters in this case
todegs = 180/pi

phi = np.linspace(-halfpi, halfpi, 101)
h = 0.
lam = 0.
N = a**2 / np.sqrt(a**2 * np.cos(phi)**2 + b**2 * np.sin(phi)**2)
X = (N + h) * np.cos(phi) * np.cos(lam)
Y = (N + h) * np.cos(phi) * np.sin(lam)
Z = ((b/a)**2 * N + h) * np.sin(phi)

phi_prime = np.arctan(Z/X) # two quandrants so no need for arctan2

plt.figure()
plt.plot(todegs*phi, todegs*(phi_prime-phi))
plt.xlabel('geodetic ("GPS") latitude (degs)', fontsize=12)
plt.ylabel('geocentric - geodetic latitude (degs)', fontsize=12)
plt.xlim(-90, 90)
plt.show()
| improve this answer | |
$\endgroup$
  • $\begingroup$ The parallax coefficients are typically used when applying the topocentric corrections (for a given site on the Earth's surface as opposed to the geocenter). They let you form a 3d vector which can be combined with precession/nutation/frame bias matrix which is what describes the orientation between the celestial coordinate system and the CIO and the terrestrial reference and the TIO. Code examples include find_orb (search for rho_cos_phi) and SOFA cookbook $\endgroup$ – astrosnapper Mar 20 at 16:25
  • $\begingroup$ @astrosnapper I've fixed this answer and will read those soon, thank you! $\endgroup$ – uhoh Mar 22 at 3:08
  • 1
    $\begingroup$ Re But why both sin and cos? That's the wrong question. The given values are $\rho\cos\phi'$ and $\rho\sin\phi'$, not $\cos\phi'$ and $\sin\phi'$ . The given values encode both radial distance and geocentric latitude. The right question is then why not radial distance and geocentric latitude? I don't know the answer to that question. An even better question is why not altitude and latitude? I don't know the answer to that question, either. (Note well: I used an unqualified "latitude" in my better question. Latitude without a qualifying prefix almost always means geodetic latitude.) $\endgroup$ – David Hammen Mar 22 at 4:32
  • $\begingroup$ @DavidHammen I see what you mean and have revised accordingly $\endgroup$ – uhoh Mar 23 at 23:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.