Calculating an angular diameter of planets

[andrey-nakin]

What is the correct way to calculate the angular diameter of a planet visible in the sky? Especially when it is close to the horizon.

I'm trying to calculate it with the regular formula:

a_in_radians = arctan(diameter / 2 / distance) * 2

When I compare the result to the values from some other amateur astronomer's tools, there's a small difference. Is my method correct for Solar system planets (including the Moon)?

[brandon-rhodes]

That approach looks good to me. To be any more specific, try specifying which tool is giving you which number, and what number you are getting from Skyfield — otherwise we can't even tell if your formula is giving you a bigger or smaller result than the one from the other tool.

[andrey-nakin]

In the following example I calculate an angular diameter of Uranus at the nearest opposition:

#!/usr/bin/env python3
import math
from skyfield import api
from skyfield.api import Topos

DIAMETER = 25362 * 2    # in kilometers
eph = api.load('de421.bsp')
earth, body = eph['earth'], eph['uranus barycenter']
location = earth + Topos(59.95, 30.32)
ts = api.load.timescale()
time = ts.utc(2021, 11, 5, 0, 0, 0)
astrometric = location.at(time).observe(body)
alt, az, d = astrometric.apparent().altaz('standard')
ad = math.atan(DIAMETER / 2 / d.km) * 360 / math.pi * 3600
print('%.2f' % ad + "''")

The result is 3.73"

In-the-sky.org shows 3.8" (it uses DE430 ephemerides)

Stellarium shows 3.77" (it uses VSOP87)

[brandon-rhodes]

The problem looks to be your choice of radius. Here's the NASA HORIZONS determination of the apparent diameter. Besides being more accurate, I have found it a much better tool than Stellarium or In-the-sky for comparing numbers to a script since it shows so much of the information that it's using to make the measurement:

*******************************************************************************
Ephemeris / WWW_USER Tue Dec 29 06:18:14 2020 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Uranus (799)                    {source: ura111}
Center body name: Earth (399)                     {source: DE431mx}
Center-site name: (user defined site below)
*******************************************************************************
Start time      : A.D. 2021-Nov-05 00:00:00.0000 UT      
Stop  time      : A.D. 2021-Nov-06 00:00:00.0000 UT      
Step-size       : 1440 minutes
*******************************************************************************
Target pole/equ : IAU_URANUS                      {East-longitude positive}
Target radii    : 25559.0 x 25559.0 x 24973.0 km  {Equator, meridian, pole}    
Center geodetic : 30.3200000,59.9500000,1.617E-12 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 30.3200000,3201.92765,5497.6897 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : High-precision EOP model        {East-longitude positive}
Center radii    : 6378.1 x 6378.1 x 6356.8 km     {Equator, meridian, pole}    
Target primary  : Sun
Vis. interferer : MOON (R_eq= 1737.400) km        {source: DE431mx}
Rel. light bend : Sun, EARTH                      {source: DE431mx}
Rel. lght bnd GM: 1.3271E+11, 3.9860E+05 km^3/s^2                              
Atmos refraction: NO (AIRLESS)
RA format       : DEG
Time format     : CAL 
EOP file        : eop.201228.p210321                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2020-DEC-28. PREDICTS-> 2021-MAR-20
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass (>38.000=NO), Daylight (NO )
Table cut-offs 2: Solar elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
Table cut-offs 3: RA/DEC angular rate (     0.0=NO )                           
*******************************************************************************
 Date__(UT)__HR:MN:SS     Ang-diam
**********************************
$$SOE
 2021-Nov-05 00:00:00     3.761194
 2021-Nov-06 00:00:00     3.761139
$$EOE
*******************************************************************************
Column meaning:

TIME

  Times PRIOR to 1962 are UT1, a mean-solar time closely related to the
prior but now-deprecated GMT. Times AFTER 1962 are in UTC, the current
civil or "wall-clock" time-scale. UTC is kept within 0.9 seconds of UT1
using integer leap-seconds for 1972 and later years.

  Conversion from the internal Barycentric Dynamical Time (TDB) of solar
system dynamics to the non-uniform civil UT time-scale requested for output
has not been determined for UTC times after the next July or January 1st.
Therefore, the last known leap-second is used as a constant over future
intervals.

  Time tags refer to the UT time-scale conversion from TDB on Earth
regardless of observer location within the solar system, although clock
rates may differ due to the local gravity field and no analog to "UT"
may be defined for that location.

  Any 'b' symbol in the 1st-column denotes a B.C. date. First-column blank
(" ") denotes an A.D. date. Calendar dates prior to 1582-Oct-15 are in the
Julian calendar system. Later calendar dates are in the Gregorian system.

  NOTE: "n.a." in output means quantity "not available" at the print-time.

SOLAR PRESENCE (OBSERVING SITE)
  Time tag is followed by a blank, then a solar-presence symbol:

        '*'  Daylight (refracted solar upper-limb on or above apparent horizon)
        'C'  Civil twilight/dawn
        'N'  Nautical twilight/dawn
        'A'  Astronomical twilight/dawn
        ' '  Night OR geocentric ephemeris

LUNAR PRESENCE (OBSERVING SITE)
  The solar-presence symbol is immediately followed by a lunar-presence symbol:

        'm'  Refracted upper-limb of Moon on or above apparent horizon
        ' '  Refracted upper-limb of Moon below apparent horizon OR geocentric
             ephemeris

 'Ang-diam' =
   The equatorial angular width of the target body full disk, if it were fully
illuminated and visible to the observer. If the target body diameter is unknown
"n.a." is output.

   Units: ARCSECONDS

 Computations by ...
     Solar System Dynamics Group, Horizons On-Line Ephemeris System
     4800 Oak Grove Drive, Jet Propulsion Laboratory
     Pasadena, CA  91109   USA
     Information  : https://ssd.jpl.nasa.gov/
     Documentation: https://ssd.jpl.nasa.gov/?horizons_doc
     Connect      : https://ssd.jpl.nasa.gov/?horizons (browser)
                    telnet ssd.jpl.nasa.gov 6775       (command-line)
                    e-mail command interface available
                    Script and CGI interfaces available
     Author       : Jon.D.Giorgini@jpl.nasa.gov

*******************************************************************************

As you can see, the value it provides is 3.761194, so that's what we should try to make Skyfield print. Happily, it shows the planetary radius it's working with, which the other tools must not have, thus leading you into a discrepancy you couldn't track down. In the above output you can find:

Target radii    : 25559.0 x 25559.0 x 24973.0 km  {Equator, meridian, pole}    

And that's your problem. Uranus has a smaller polar radius, a bigger equatorial radius, and the number you used is somewhere in between (maybe someone's definition of an average radius, between the two?). The HORIZONS system is using the biggest number, so that it prints the largest angle you will see across Uranus's disc with your telescope. We can produce the same value with Skyfield:

#!/usr/bin/env python3

import math
from skyfield import api
from skyfield.api import Topos

radius_km = 25559.0
eph = api.load('de421.bsp')
earth, body = eph['earth'], eph['uranus barycenter']
location = earth + Topos(59.95, 30.32)
ts = api.load.timescale()
time = ts.utc(2021, 11, 5, 0, 0, 0)
astrometric = location.at(time).observe(body)
alt, az, d = astrometric.apparent().altaz('standard')
apparent_radius_degrees = math.asin(radius_km / d.km) / math.tau * 360
apparent_diameter_degrees = apparent_radius_degrees * 2
print("{:.6f}''".format(apparent_diameter_degrees * 3600))

This script prints the same value as HORIZONS to all 6 decimal places (to microarcseconds):

3.761194''

You will note that I switched from atan() to asin() because, technically, the right angle that's formed is not between your line-of-sight and the planet's invisible center down inside of it, but between your line-of-sight and the planet's edge, since that's the point tangent between a ray of light and a sphere's edge. But as atan() and asin() are merely identity functions for angles this tiny, you'll note that the choice doesn't affect the output.

Where in the Skyfield documentation might we think about including an explanation of the concept of “angular diameter” where you would have been likely to find it? It would be nice to explain the subtleties somewhere so that future folks tackling the same issue will have some guidance — but of course it will only help if it's somewhere in the docs they can find it.

[andrey-nakin]

Thanks a lot for the detailed answer, especially for referring HORIZONS tool.

Using asin in place of atan is not clear to me. The value returned by altaz is a distance to the barycenter of Uranus, not its apparent edge, correct?

Having this info in the documentation would be very helpful.

[brandon-rhodes]

Thanks a lot for the detailed answer, especially for referring HORIZONS tool.

You're welcome! I entirely forgot to include a link, didn't I? I presume you found it, but for folks who might happen by later:

https://ssd.jpl.nasa.gov/horizons.cgi

Using asin in place of atan is not clear to me. The value returned by altaz is a distance to the barycenter of Uranus, not its apparent edge, correct?

1. Note that the barycenter isn't the planet's center, but the center of the whole Uranus+moons system; happily in this case, though, its moons are so trifling in size compared to the planet that you can just use the system barycenter as the planet's center and never notice the difference.
2. Yes, the distance returned is the distance to the center.
3. Picture it this way: when you shift your eye over to the edge of the planet's disc, you are seeing the point where a radius drawn outward towards you a little from the planet intersects the surface. You can't actually see "around the corner" to a radius at a right angle to the you→center vector, because you're not at infinite distance from the planet. Think about the Earth: the distance from your eye to the horizon isn't the distance from your eye to where a radius 90° from your nadir would reach the Earth's surface 90° away from you. You see the horizon where you line-of-sight meets a nearby radius at right angles. It's the same with other spheres: your line of sight meets their radius at right angles, making the distance-to-the-center the actual hypotenuse. A bad diagram:

            \
            |
    * <center of Uranus
    |--     | <surface of Uranus
radius>--   /
    |    --/ <the right angle
_   |   __/
 \__|__/ /
    |   /
 d> |  / <asin(d)
    | /
    |/
    * You

Having this info in the documentation would be very helpful.

Where did you look in the documentation as you worked through the angular distance problem? That will give me a hint about where to put it where folks might find it — did you, for example, take a look at the Examples chapter? Did you search the headings in the Table of Contents? Or take another way in?

[andrey-nakin]

Great, I got it, thanks again. Yes, using asin is more accurate even when it does not make any difference. But in the case of, say, Moon the difference can be noticeable.

When reading the Skyfield docs, I usually start from the TOC (BTW I think it is well-structured!). Then I search on the page using the keywords, e.g. "angular".

[brandon-rhodes]

There! I've added an example to the docs, which should appear online with the next release. Let me know if you see any problems with it. Thanks!

[andrey-nakin]

Thank you, the updated doc is helpful.