1# -*- coding: utf-8 -*- 2"""Routines for computing magnitudes. 3 4Planetary routines adapted from: 5 6https://arxiv.org/pdf/1808.01973.pdf 7 8Which links to: 9 10https://sourceforge.net/projects/planetary-magnitudes/ 11 12Which has directories with three successive versions of their magnitude 13computation, the most recent of which provides the files on which this 14Python code is based: 15 16Ap_Mag_V3.f90 17Ap_Mag_Output_V3.txt 18Ap_Mag_Input_V3.txt 19 20* ``r`` planet’s distance from the Sun. 21* ``delta`` from Earth? 22* ``ph_ang`` illumination phase angle (degrees) 23 24""" 25from numpy import array, clip, exp, log10, nan, sin, where 26from .constants import RAD2DEG 27from .functions import angle_between, length_of 28from .naifcodes import _target_name 29 30# See "design/planet_tilts.py" in the Skyfield repository. 31_SATURN_POLE = array([0.08547883, 0.07323576, 0.99364475]) 32_URANUS_POLE = array([-0.21199958 -0.94155916 -0.26176809]) 33 34def planetary_magnitude(position): 35 """Given the position of a planet, return its visual magnitude. 36 37 >>> from skyfield.api import load 38 >>> from skyfield.magnitudelib import planetary_magnitude 39 >>> ts = load.timescale() 40 >>> t = ts.utc(2020, 7, 31) 41 >>> eph = load('de421.bsp') 42 >>> astrometric = eph['earth'].at(t).observe(eph['jupiter barycenter']) 43 >>> print('%.2f' % planetary_magnitude(astrometric)) 44 -2.73 45 46 The formulae are from `Mallama and Hilton “Computing Apparent 47 Planetary Magnitude for the Astronomical Almanac” (2018) 48 <https://arxiv.org/pdf/1808.01973.pdf>`_. Two of the formulae have 49 inherent limits: 50 51 * Saturn’s magnitude is unknown and the function will return ``nan`` 52 (the floating-point value “Not a Number”) if the “illumination 53 phase angle” — the angle of the vertex observer-Saturn-Sun — 54 exceeds 6.5°. 55 56 * Neptune’s magnitude is unknown and will return ``nan`` if the 57 illumination phase angle exceeds 1.9° and the position's date is 58 before the year 2000. 59 60 And one formula is not fully implemented (though contributions are 61 welcome!): 62 63 * Skyfield does not compute which features on Mars are facing the 64 observer, which can introduce an error of ±0.06 magnitude. 65 66 """ 67 target = position.target 68 function = _FUNCTIONS.get(target) 69 if function is None: 70 name = _target_name(target) 71 raise ValueError('cannot compute the magnitude of target %s' % name) 72 73 # Shamelessly treat the Sun as sitting at the Solar System Barycenter. 74 sun_to_observer = position.center_barycentric.xyz.au 75 observer_to_planet = position.xyz.au 76 sun_to_planet = sun_to_observer + observer_to_planet 77 78 r = length_of(sun_to_planet) 79 delta = length_of(observer_to_planet) 80 ph_ang = angle_between(sun_to_planet, observer_to_planet) * RAD2DEG 81 82 if function is _saturn_magnitude: 83 a = angle_between(_SATURN_POLE, sun_to_planet) 84 sun_sub_lat = a * RAD2DEG - 90.0 85 86 a = angle_between(_SATURN_POLE, observer_to_planet) 87 observer_sub_lat = a * RAD2DEG - 90.0 88 89 return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat) 90 91 if function is _uranus_magnitude: 92 a = angle_between(_URANUS_POLE, sun_to_planet) 93 sun_sub_lat = a * RAD2DEG - 90.0 94 95 a = angle_between(_URANUS_POLE, observer_to_planet) 96 observer_sub_lat = a * RAD2DEG - 90.0 97 98 return function(r, delta, ph_ang, sun_sub_lat, observer_sub_lat) 99 100 if function is _neptune_magnitude: 101 year = position.t.J 102 return function(r, delta, ph_ang, year) 103 104 return function(r, delta, ph_ang) 105 106def _mercury_magnitude(r, delta, ph_ang): 107 distance_mag_factor = 5 * log10(r * delta) 108 ph_ang_factor = ( 109 6.3280e-02 * ph_ang 110 - 1.6336e-03 * ph_ang**2 111 + 3.3644e-05 * ph_ang**3 112 - 3.4265e-07 * ph_ang**4 113 + 1.6893e-09 * ph_ang**5 114 - 3.0334e-12 * ph_ang**6 115 ) 116 return -0.613 + distance_mag_factor + ph_ang_factor 117 118def _venus_magnitude(r, delta, ph_ang): 119 distance_mag_factor = 5 * log10(r * delta) 120 condition = ph_ang < 163.7 121 a0 = where(condition, 0.0, 236.05828 + 4.384) 122 a1 = where(condition, - 1.044E-03, - 2.81914E+00) 123 a2 = where(condition, + 3.687E-04, + 8.39034E-03) 124 a3 = where(condition, - 2.814E-06, 0.0) 125 a4 = where(condition, + 8.938E-09, 0.0) 126 ph_ang_factor = a4 127 for a in a3, a2, a1, a0: 128 ph_ang_factor *= ph_ang 129 ph_ang_factor += a 130 return -4.384 + distance_mag_factor + ph_ang_factor 131 132def _earth_magnitude(r, delta, ph_ang): 133 distance_mag_factor = 5 * log10(r * delta) 134 ph_ang_factor = -1.060e-03 * ph_ang + 2.054e-04 * ph_ang**2 135 return -3.99 + distance_mag_factor + ph_ang_factor 136 137def _mars_magnitude(r, delta, ph_ang): 138 r_mag_factor = 2.5 * log10(r * r) 139 delta_mag_factor = 2.5 * log10(delta * delta) 140 distance_mag_factor = r_mag_factor + delta_mag_factor 141 142 geocentric_phase_angle_limit = 50.0 143 144 condition = ph_ang <= geocentric_phase_angle_limit 145 a = where(condition, 2.267E-02, - 0.02573) 146 b = where(condition, - 1.302E-04, 0.0003445) 147 ph_ang_factor = a * ph_ang + b * ph_ang**2 148 149 # Compute the effective central meridian longitude 150 # eff_CM = ( sub_earth_long + sub_sun_long ) / 2. 151 # if ( abs ( sub_earth_long - sub_sun_long ) > 180. ): 152 # Eff_CM = Eff_CM + 180. 153 # if ( Eff_CM > 360. ): 154 # Eff_CM = Eff_CM - 360. 155 156 # ! Use Stirling interpolation to determine the magnitude correction 157 # call Mars_Stirling ( 'R', eff_CM, mag_corr_rot ) 158 159 # Convert the ecliptic longitude to Ls 160 # Ls = h_ecl_long + Ls_offset 161 # if ( Ls > 360. ) Ls = Ls - 360. 162 # if ( Ls < 0. ) Ls = Ls + 360. 163 164 # Use Stirling interpolation to determine the magnitude correction 165 # call Mars_Stirling ( 'O', Ls, mag_corr_orb ) 166 167 # Until effects from Mars rotation are written up: 168 mag_corr_rot = 0.0 169 mag_corr_orb = 0.0 170 171 # Add factors to determine the apparent magnitude 172 ap_mag = where(ph_ang <= geocentric_phase_angle_limit, -1.601, -0.367) 173 ap_mag += distance_mag_factor + ph_ang_factor + mag_corr_rot + mag_corr_orb 174 175 return ap_mag 176 177def _jupiter_magnitude(r, delta, ph_ang): 178 distance_mag_factor = 5 * log10(r * delta) 179 geocentric_phase_angle_limit = 12.0 180 ph_ang_pi = ph_ang / 180.0 181 ph_ang_factor = where( 182 ph_ang <= geocentric_phase_angle_limit, 183 (6.16E-04 * ph_ang - 3.7E-04) * ph_ang, 184 -2.5 * log10( 185 ((((- 1.876 * ph_ang_pi + 2.809) * ph_ang_pi - 0.062) * ph_ang_pi 186 - 0.363) * ph_ang_pi - 1.507) * ph_ang_pi + 1.0 187 ), 188 ) 189 ap_mag = where( 190 ph_ang <= geocentric_phase_angle_limit, 191 -9.395 + distance_mag_factor + ph_ang_factor, 192 -9.428 + distance_mag_factor + ph_ang_factor, 193 ) 194 return ap_mag 195 196def _saturn_magnitude(r, delta, ph_ang, sun_sub_lat, earth_sub_lat, rings=True): 197 # Note that sun_sub_lat and earth_sub_lat should be saturnicentric 198 # latitude, not saturnidetic. 199 200 r_mag_factor = 2.5 * log10(r * r) 201 delta_mag_factor = 2.5 * log10(delta * delta) 202 distance_mag_factor = r_mag_factor + delta_mag_factor 203 204 # Then take the square root of the product of the saturnicentric 205 # latitude of the Sun and that of Earth but set to zero when the 206 # signs are opposite 207 product = sun_sub_lat * earth_sub_lat 208 signs_same = product >= 0.0 209 square_root = product ** where(signs_same, 0.5, 0.0) # avoid sqrt(neg) 210 sub_lat_geoc = where(signs_same, square_root, 0.0) 211 212 # Compute the effect of phase angle and inclination 213 214 geocentric_phase_angle_limit = 6.5 215 geocentric_inclination_limit = 27.0 216 217 is_within_geocentric_bounds = ( 218 (ph_ang <= geocentric_phase_angle_limit) 219 & (sub_lat_geoc <= geocentric_inclination_limit) 220 ) 221 222 ap_mag = where( 223 is_within_geocentric_bounds, 224 where( 225 rings, 226 227 # Use equation #10 for globe+rings and geocentric circumstances. 228 -8.914 - 1.825 * sin(sub_lat_geoc / RAD2DEG) + 0.026 * ph_ang 229 - 0.378 * sin(sub_lat_geoc / RAD2DEG) * exp(-2.25 * ph_ang), 230 231 # Use equation #11 for globe-alone and geocentric circumstances 232 -8.95 - 3.7e-4 * ph_ang + 6.16e-4 * ph_ang**2, 233 ), 234 where( 235 (ph_ang > geocentric_phase_angle_limit) & (~rings), 236 237 # Use equation #12 for globe-alone beyond geocentric phase 238 # angle limit 239 -8.94 + 2.446e-4 * ph_ang + 2.672e-4 * ph_ang**2 240 - 1.506e-6 * ph_ang**3 +4.767e-9 * ph_ang**4, 241 242 nan, 243 ), 244 ) 245 246 ap_mag = ap_mag + distance_mag_factor 247 return ap_mag 248 249def _uranus_magnitude(r, delta, ph_ang, 250 sun_sub_lat_planetog, earth_sub_lat_planetog): 251 distance_mag_factor = 5.0 * log10 (r * delta) 252 sub_lat_planetog = (abs(sun_sub_lat_planetog) 253 + abs(earth_sub_lat_planetog)) / 2.0 254 sub_lat_factor = -0.00084 * sub_lat_planetog 255 geocentric_phase_angle_limit = 3.1 256 ap_mag = -7.110 + distance_mag_factor + sub_lat_factor 257 ap_mag += where( 258 ph_ang > geocentric_phase_angle_limit, 259 (1.045e-4 * ph_ang + 6.587e-3) * ph_ang, 260 0.0, 261 ) 262 return ap_mag 263 264def _neptune_magnitude(r, delta, ph_ang, year): 265 r_mag_factor = 2.5 * log10(r * r) 266 delta_mag_factor = 2.5 * log10(delta * delta) 267 distance_mag_factor = r_mag_factor + delta_mag_factor 268 269 # Equation 16 compute the magnitude at unit distance as a function of time 270 ap_mag = clip(-6.89 - 0.0054 * (year - 1980.0), -7.00, -6.89) 271 ap_mag += distance_mag_factor 272 273 geocentric_phase_angle_limit = 1.9 274 275 ap_mag = where( 276 ph_ang > geocentric_phase_angle_limit, 277 278 # Add phase angle factor from equation 17 279 # Check the year because equation 17 only pertains to t > 2000.0 280 where( 281 year >= 2000.0, 282 ap_mag + 7.944e-3 * ph_ang + 9.617e-5 * ph_ang**2, 283 nan, 284 ), 285 286 # Otherwise leave the value unchanged. 287 ap_mag, 288 ) 289 290 return ap_mag 291 292_FUNCTIONS = { 293 199: _mercury_magnitude, 294 299: _venus_magnitude, 295 399: _earth_magnitude, 296 499: _mars_magnitude, 297 599: _jupiter_magnitude, 298 699: _saturn_magnitude, 299 799: _uranus_magnitude, 300 899: _neptune_magnitude, 301 302 # Some planets can be reasonably identified with their barycenter. 303 1: _mercury_magnitude, 304 2: _venus_magnitude, 305 4: _mars_magnitude, 306 5: _jupiter_magnitude, 307 6: _saturn_magnitude, 308 7: _uranus_magnitude, 309 8: _neptune_magnitude, 310} 311