# Introduction

On Twitter today Jonathan McDowell made the following comment (https://twitter.com/planet4589/status/879927853332213760):

>Russian cubesat 'Mayak' to launch Jul 14 will deploy a reflector that may be as bright as mag minus 10. Astronomers likely to be annoyed.

This resulted in comprehensive nerd sniping (https://xkcd.com/356/) of both myself (https://twitter.com/vacant3rdman/status/879937146508812289) and Matthew Kenworthy (https://twitter.com/mattkenworthy/status/879954847809888256), and possibly others. As a result I was compelled to attempt to work out the maximum apparent brightness of mirrored surfaces in Earth orbit. This notebook is the result.

# Brightness of an orbital mirror

The apparent brightness of a specular reflection of sunlight off a flat surface is simply given by the apparent brightness of the Sun scaled by the ratio of the solid angles subtended by the reflector and the Sun from the observer's position and by the specular reflectance of the surface (which will, in general, be a function of wavlength and angle of incidence). Expressed in astronomical magnitudes (yuck), this becomes

\begin{equation}
m_\mathrm{sat} = m_{\odot} - 2.5 \times \log_{10}\left(\frac{\Omega_\mathrm{sat}}{\Omega_\odot}.R_\mathrm{sat}\right),
\end{equation}

where $m_\mathrm{sat}$ and $m_\odot$ are the apparent magnitudes of the satellite and Sun respectively, $\Omega_\mathrm{sat}$ and $\Omega_\odot$ are their solid angles and $R_\mathrm{sat}$ is the reflectance of the flat surface on the satellite.

For a flat surface of area $A$ at an angle $\theta$ to the line of sight and seen from a distance $d$ the solid angle (using small angle approximations) is given by

\begin{equation}
\Omega_\mathrm{sat} = \frac{A \cos{\theta}}{d^2}.
\end{equation}

## Conversative upper limit for the Mayak CubeSat

We can put an upper limit on the apparent brightness of the satellite by using the minimum distance between observer and satellite (the altitude of the satellite's orbit, $a$) and normal incidence ($\theta=0$) so that we obtain

\begin{equation}
m_\textrm{sat} > m_{\odot} - 2.5 \times \log_{10}\left(\frac{A.R_\textrm{sat}}{a^2 \Omega_\odot}\right).
\end{equation}

Substituting in values for Mayak (http://space.skyrocket.de/doc_sdat/mayak.htm) of $A = 4\textrm{m}^2$, $a = 500\textrm{km}$ and $R_\textrm{sat} = 0.95$ we obtain the following limit


```python
%matplotlib inline
```


```python
import numpy as np
import astropy.units as u
import astropy.constants as c
import matplotlib.pyplot as plt
from scipy.optimize import minimize
```


```python
omega_sun = (np.pi * (c.R_sun / c.au)**2).to(u.sr, equivalencies=u.equivalencies.dimensionless_angles())
omega_sun
```




$6.7905243 \times 10^{-5} \; \mathrm{sr}$




```python
m_sun = -26.74 * u.dimensionless_unscaled
m_sun
```




$-26.74 \; \mathrm{}$




```python
area = 4 * u.m**2
altitude = 500 * u.km
reflectance = 0.95 * u.dimensionless_unscaled
```


```python
omega_sat_max = (area / altitude**2).to(u.sr, equivalencies=u.equivalencies.dimensionless_angles())
omega_sat_max
```




$1.6 \times 10^{-11} \; \mathrm{sr}$




```python
m_sat_min = m_sun - 2.5 * np.log10(omega_sat_max * reflectance / omega_sun)
m_sat_min
```




$-10.114851 \; \mathrm{}$



**So we find that the apparent brightness of Mayak should be no brighter than magnitude -10.1, which is consistent with the -10 figure quoted on Gunter's Space page.**  It is also in agreement with the result that Matthew Kenworthy arrived at independently, https://twitter.com/mattkenworthy/status/879954847809888256

This is a conservative upper limit, however. For the satellite to be at a distance equal to the altitude of its orbit it must be directly overhead of the observer, and $\theta=0$ implies that the reflector surface is directly facing the observer. For a specular reflection of sunlight to reach the observer the angle of incidence must equal the angle of reflection, i.e. in this case the Sun must be in line with the satellite and observer, in other words on the opposite side of the Earth. The transmission of the entire planet at visible wavelengths is *extremely* low so in practice very little light would reach the satellite in this geometry.

## Slightly more rigourous limit

In reality the maximum brighness will occur when the observer, satellite and Sun lie in the same plane with the Sun some distance below the horizon from the observer's point of view but above the horizon from the satellite's point of view. The geometry is illustrated in this lovingly hand crafted diagram:



This introduces two new variables, the angular distance of the Sun below the horizon at the observer's position, $\phi$, and the zenith angle of the satellite from the observer's position. We must solve for the line of sight distance between the observer and satellite, $d$, and the angle of incidence and reflection at the mirror surface, $\theta$ in terms of these two variables. It was at this point that I realised that I had forgotten almost everything that I ever knew about elementary geometry and trigonometry. Fortunately Wikipedia came to the rescue.

### Angle of incidence/reflection

From the sum of the interior angles of a triangle,
\begin{equation}\theta = \frac{\pi/2 - z - \phi}{2}.
\end{equation}

### Line of sight distance

This an instance of the two side and non-included angle triangle problem (https://en.wikipedia.org/wiki/Solution_of_triangles#Two_sides_and_non-included_angle_given_.28SSA.29). The solution for the length of the third side gives

\begin{equation}
d = R_\oplus \cos(\pi - z) \pm  \sqrt{(R_\oplus + a)^2 - R_\oplus \sin^2(\pi - z)}.
\end{equation}

After a bit of arguably useful rearrangement this becomes

\begin{equation}
d = a \left[ \sqrt{1 + \frac{2 R_\oplus}{a} + \frac{R_\oplus^2}{a^2}\cos^2 z} - \frac{R_\oplus}{a} \cos{z} \right],
\end{equation}

which reassuringly gives $d = a$ when zenith angle is zero. Now we create a function to calculate $d$ as a function of $a$ and $z$ and plot the results.


```python
def distance(altitude, zenith_angle):
    s = 1 + 2 * c.R_earth / altitude + (c.R_earth * np.cos(zenith_angle) / altitude)**2
    return (altitude * (s**0.5 - c.R_earth * np.cos(zenith_angle) / altitude)).to(u.km)
```


```python
angles = np.arange(-90, 91, 1) * u.degree
plt.plot(angles, distance(altitude, angles))
plt.xlim(0, 90)
plt.ylim(0, 2600)
plt.xlabel('Zenith angle / degree')
plt.ylabel('Distance / km')
plt.title('Line of sight distance as a function of zenith angle for an altitude of 500 km ')
plt.gcf().set_size_inches(12, 8)
```

### Eclipse condition

In order for sunlight to be reflected from the satellite to the observer it is necessary that the Earth is not in the way, i.e. the spacecraft is not in eclipse. Calculating this condition in terms of satellite altitude, zenith angle and observer's Sun elevation required doing some more geometry. The solution is 

\begin{equation}
(R_\oplus + a) .\cos(\psi - \phi) > R_\oplus,
\end{equation}

where 

\begin{equation}
\psi = \arcsin\left(\frac{d \sin z}{R_\oplus + a}\right).
\end{equation}

We will write a function to determine whether the satellite is eclipsed:


```python
def is_eclipsed(altitude, zenith_angle, sun_elevation):
    los_distance = distance(altitude, zenith_angle)
    psi = np.arcsin(los_distance * np.sin(zenith_angle) / (c.R_earth + altitude))
    return (c.R_earth + altitude) * np.cos(psi - sun_elevation) < c.R_earth
```


```python
plt.plot(angles, is_eclipsed(altitude, angles, +6* u.degree), label='+6$^\circ$')
plt.plot(angles, is_eclipsed(altitude, angles, 0 * u.degree), label='0$^\circ$')
plt.plot(angles, is_eclipsed(altitude, angles, -6 * u.degree), label='-6$^\circ$')
plt.plot(angles, is_eclipsed(altitude, angles, -12 * u.degree), label='-12$^\circ$')
plt.plot(angles, is_eclipsed(altitude, angles, -18 * u.degree), label='-18$^\circ$')
plt.plot(angles, is_eclipsed(altitude, angles, -24 * u.degree), label='-24$^\circ$')
plt.plot(angles, is_eclipsed(altitude, angles, -30 * u.degree), label='-30$^\circ$')
plt.xlim(-90, 90)
plt.ylim(1.1,-0.1)
plt.xlabel('Zenith angle / degree')
plt.yticks((0, 1), ('Not eclipsed', 'Eclipsed'))
plt.title('Eclipsed condition for a satellite at 500 km altitude')
plt.legend(title='Sun elevation')
plt.gcf().set_size_inches(12, 8)
```

The perhaps surprising 'eclipsed' results for positive Sun elevation angles are indicating the situation when the satellite is in daylight but the Sun is on the opposite side of the mirror surface to the observer.

### Maximum brightness as a function of zenith angle and Sun elevation

Now we can write a function to calculate the maximum apparent brightness as a function of zenith angle and Sun elevation.


```python
omega_sun = (np.pi * (c.R_sun / c.au)**2).to(u.sr, equivalencies=u.equivalencies.dimensionless_angles())
m_sun = -26.74 * u.dimensionless_unscaled

def m_sat(area, reflectance, altitude, zenith_angle, sun_elevation):
    angle_of_incidence = (np.pi * u.radian / 2 - zenith_angle + sun_elevation) / 2 # Elevation negative in twilight
    los_distance = distance(altitude, zenith_angle)
    omega_sat = (area * np.cos(angle_of_incidence) / \
                 los_distance**2).to(u.sr, equivalencies=u.equivalencies.dimensionless_angles())

    m_sat = np.where(is_eclipsed(altitude, zenith_angle, sun_elevation),
                     999 * u.dimensionless_unscaled,
                     m_sun - 2.5 * np.log10(omega_sat * reflectance / omega_sun)) * u.dimensionless_unscaled
    
    return m_sat
```


```python
plt.plot(angles, m_sat(area, reflectance, altitude, angles, 6* u.degree), label='+6$^\circ$')
plt.plot(angles, m_sat(area, reflectance, altitude, angles, 0 * u.degree), label='0$^\circ$')
plt.plot(angles, m_sat(area, reflectance, altitude, angles, -6 * u.degree), label='-6$^\circ$')
plt.plot(angles, m_sat(area, reflectance, altitude, angles, -12 * u.degree), label='-12$^\circ$')
plt.plot(angles, m_sat(area, reflectance, altitude, angles, -18 * u.degree), label='-18$^\circ$')
plt.plot(angles, m_sat(area, reflectance, altitude, angles, -24 * u.degree), label='-24$^\circ$')
plt.plot(angles, m_sat(area, reflectance, altitude, angles, -30 * u.degree), label='-30$^\circ$')
plt.xlim(-90, 90)
plt.ylim(-6,-10.5)
plt.xlabel('Zenith angle / degree')
plt.ylabel('Maximum brightness / magnitudes')
plt.title('Maximum apparent brightness for a 4m$^2$ 95% mirror at 500 km altitude')
plt.legend(title='Sun elevation')
plt.gcf().set_size_inches(12, 8)
plt.savefig('mayak.png')
```

### Final result

We can maximise brightness (minimise magnitude) over zenith angle and Sun elevation to find the highest possible apparent brightness.


```python
results = minimize(lambda x: m_sat(area, reflectance, altitude, x[0] * u.degree, x[1] * u.degree),
                   (0, -18),
                   method='Powell')
results
```




       direc: array([[ 1.,  0.],
           [ 0.,  1.]])
         fun: <Quantity -9.93558237180596>
     message: 'Optimization terminated successfully.'
        nfev: 56
         nit: 2
      status: 0
     success: True
           x: <Quantity [  8.84357143,-21.33344709]>



**So the maximum apparent brightness is, apparently, magnitude -9.9, which can occur when the Sun elevation is -21.3° and the satellite is at a zenith angle of 8.8° towards the antisolar direction**


```python
m_sat_peak = results['fun']
mag_diff = m_sat_peak - m_sat_min
mag_diff
```




$0.17926833 \; \mathrm{}$




```python
flux_ratio = 10**(-0.4 * mag_diff)
flux_ratio
```




$0.84779855 \; \mathrm{}$




```python
inverse_flux_ratio = 10**(0.4 * mag_diff)
inverse_flux_ratio
```




$1.1795255 \; \mathrm{}$



This value is only 0.2 magnitudes fainter than the conservative (and far simpler) upper limit above, which is equivalent to 85% of the flux. I appear to have wasted a *lot* of time for a 18% improvement in accuracy...

# Flash duration

The entry for Mayak on Gunter's Space Page notes the following

>The reflector will provide a -10 optical magnitude at the beginning of the flight to allow for easy tracking. Mayak will be put into a tumbling motion over all axes, with at least 1 revolution per second

This suggests that Mayak may initially be in a stable orientation, resulting in flares of comparable duration to those from Iridium satellite (i.e. a few seconds) along a narrow ground track, but will spend most if not all of its mission in a rapid tumble. This will result in bright but very brief flashes of light visible over many different ground tracks as beams of reflected light sweep rapidly over the ground below.

## Maximum duration

During a flash the observer is effectively looking at a small portion of the Sun reflected in the mirrored surface of the satellite. The direction of the observer's reflected view will move across the sky at twice the angular rotation rate of the mirror, therefore the maximum flash duration will be time it takes for the satellite to rotate by half the angular diameter of the Sun.


```python
omega_min = 2 * np.pi * u.radian / u.second
omega_min
```




$6.2831853 \; \mathrm{\frac{rad}{s}}$




```python
duration_max = (c.R_sun / c.au).to(u.radian, equivalencies=u.equivalencies.dimensionless_angles()) / omega_min
duration_max = duration_max.to(u.millisecond)
duration_max
```




$0.73994059 \; \mathrm{ms}$



**The maximum flash duration is therefore 0.74 milliseconds**.

## Mean flash duration

The maximum flash duration will occur when the observer's view of the Sun crosses the centre of the solar disc. In general the track across the disc will be off centre, resulting in a distribution of shorter flash durations. The (angular) length of a chord across the solar disc (angular radius $\rho_\odot$) is given by $\delta = 2 \rho_\odot \sqrt{1 - \lambda^2 / \rho_\odot^2}$, where $\lambda$ is the offset. The corresponding flash duration is given by 

\begin{equation}
t = \frac{\delta}{2 \omega_\mathrm{sat}} = \frac{\rho_\odot}{\omega_\mathrm{sat}} \sqrt{1 - \frac{\lambda^2}{\rho_\odot^2}} = t_\mathrm{max} \sqrt{1 - \frac{\lambda^2}{\rho_\odot^2}}
\end{equation}

Both the angular offset and orientation of the track (which has no effect on duration) can be assumed to be randomly and uniformly distributed so that the normalised probability of a given offset, $P(\lambda)\,d\lambda = 2 \lambda/ \rho_\odot^2\,d\lambda$. The expected value for flash duration is then given by

\begin{equation}
\left<t\right> = \int_0^{\rho_\odot}{t(\lambda) P(\lambda)\,d\lambda} = t_\mathrm{max} \int_0^{\rho_\odot} \sqrt{1 - \frac{\lambda^2}{\rho_\odot^2}} . 2 \lambda/ \rho_\odot^2\,d\lambda = t_\mathrm{max} \left[ -\frac{2}{3} \left(1 - \frac{\lambda^2}{\rho_\odot^2} \right)^\frac{3}{2} \right]_0^{\rho_\odot} = \frac{2}{3} t_\mathrm{max}.
\end{equation}



```python
duration_mean = 2 / 3 * duration_max
duration_mean    
```




$0.49329372 \; \mathrm{ms}$



**The mean flash duration from Mayak will be 0.5 millisecond or less, if the tumbling rate exceeds 1 revolution per second.**

The flashes from Mayak are so short that it is likely to significantly reduce the perceived brightness as seen by the human eye.

## Perceived brightness of short pulses

Bloch's Law states that below a certain threshold perceived brightness scales with total pulse energy rather than flux so that a pulse of half the length needs to be twice as bright in order to be perceived as the same brightness, however the situation is complicated by an apparent peak in sensitivity for pulse durations of ~50-100 milliseconds (Broca-Sulzer effect), see for example http://www.pnas.org/content/109/48/19828. If we assume that the Bloch's Law scaling holds for durations of from 40 milliseconds then the 0.74 millisecond maximum duration of the Mayak flashes could mean that they appear a *lot* fainter than magnitude -9.9.


```python
perception_factor_max = (duration_max / (40 * u.ms)).to(u.dimensionless_unscaled)
perception_factor_max
```




$0.018498515 \; \mathrm{}$




```python
1 / perception_factor_max
```




$54.058394 \; \mathrm{}$




```python
perceived_magnitude_offset_min = - 2.5 * np.log10(perception_factor_max)
perceived_magnitude_offset_min
```




$4.3321579 \; \mathrm{}$




```python
m_sat_peak + perceived_magnitude_offset_min
```




$-5.6034245 \; \mathrm{}$



**So Bloch's Law suggests that even the longest Mayak flashes would be perceived as if they were ~50 times fainter, or 4.3 magnitudes. This would make even the brightest possible Mayak flashes appear to be magnitude -5.6 or fainter, depending on duration.**


```python
perception_factor_mean = (duration_mean / (40 * u.ms)).to(u.dimensionless_unscaled)
perception_factor_mean
```




$0.012332343 \; \mathrm{}$




```python
1 / perception_factor_mean
```




$81.087592 \; \mathrm{}$




```python
perceived_magnitude_offset_mean = - 2.5 * np.log10(perception_factor_mean)
perceived_magnitude_offset_mean
```




$4.772386 \; \mathrm{}$




```python
m_sat_peak + perceived_magnitude_offset_mean
```




$-5.1631964 \; \mathrm{}$



**A mean duration Mayak flash would be perceived as if they were ~80 times fainter, or 4.8 magnitudes. This would make even the brightest possible Mayak flashes appear no brighter than magnitude -5.2.**

1ec018c714c2c62a8303477ea6d40dad4612a424
ipynb
Jupyter Notebook
Mirror satellite brightness.ipynb
AnthonyHorton/spacecraftspotting
11086052e7b3e3f45a277eada88d82f8e02a4c4b
MIT
2017-06-28T15:22:43.000Z
2017-06-28T20:54:45.000Z
Mirror satellite brightness.ipynb
AnthonyHorton/spacecraftspotting
11086052e7b3e3f45a277eada88d82f8e02a4c4b
MIT
Mirror satellite brightness.ipynb
AnthonyHorton/spacecraftspotting
11086052e7b3e3f45a277eada88d82f8e02a4c4b
MIT
Qwen/Qwen-72B
1. YES
2. YES