A simple simulation of how many civilizations we might be expected to potentially be able to observe based on the Drake Equation and simple geometry.

fermi_galactic_sim.py

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
"""
Licensed by Vyryn under the MIT License.
A randomization plot of a number of civilizations in the galaxy and their light cones.
To give a general visual idea of how few would be observable relative to how many there are.
"""

scaling_factor = 1/3000

num_civs = 1000  # estimated number of civs in our galaxy at an arbitrary single time. This isn't used for simultaneity,
# but assumes due to speed of light we can approximate simultaneity as simultaneity of observability

galax_outer_radius = 175000.0  # light years, approximate outermost place a civ could appear
galax_thick = 2000.0  # light years, the thickness of the galaxy (bottom edge at 0), assuming a civ could appear
# anywhere in the thickness
galax_inner_radius = 10000.0  # light years, assumed closest in location to the center of the galaxy life could develop
galax_angle = 2 * np.pi  # possible angle around the galaxy that there can be life at, assumed full range
min_light_bubble_radius = 100  # light years, how far our light bubble has reached, this is a suitable known minimum
# broadcast duration/size
max_light_bubble_radius = 10000  # years, how long we estimate max broadcast duration/size could be
oldest_possible_civ = 75000  # years, how long ago we estimate oldest civ could have developed

planets_radi = galax_outer_radius * np.random.rand(num_civs)
planets_depth = galax_thick * np.random.rand(num_civs) + 100
planets_angle = 2 * np.pi * np.random.rand(num_civs)
planets_age = oldest_possible_civ * np.random.rand(num_civs)
planets = np.column_stack((planets_radi, planets_depth, planets_angle))

broadcast_time = 2 * np.random.rand(num_civs)*(max_light_bubble_radius - min_light_bubble_radius) + min_light_bubble_radius
radius = planets_age * scaling_factor
radius2 = (broadcast_time + planets_age) * scaling_factor
area = radius**2 * 2 * np.pi
area2 = radius2**2 * 2 * np.pi


n = 0
for planet in planets:
    n += 1
    print(f'Planet {n} : {planet}\n')

# Double up. Interweave two of each array.
planets_angle = np.vstack((planets_angle, planets_angle)).reshape((-1,), order='F')
planets_radi = np.vstack((planets_radi, planets_radi)).reshape((-1,), order='F')
planets_depth = np.vstack((planets_depth + 500, [0] * planets_depth)).reshape((-1,), order='F')
area = np.vstack((area2, area)).reshape((-1,), order='F')

# Plot time
fig = plt.figure()
ax = fig.add_subplot(111, polar=True)
marker = matplotlib.markers.MarkerStyle(marker='o', fillstyle='full')
c = ax.scatter(planets_angle, planets_radi, c=planets_depth, s=area2, marker=marker, cmap='Blues', alpha=0.5)
# c = ax.scatter(planets_angle, planets_radi, c=planets_depth, s=area2, marker=marker, cmap='autumn', alpha=0.5)
# c = ax.scatter(planets_angle, planets_radi, c=[1] * planets_depth, s=area, marker=marker, cmap='autumn', alpha=0.5)
c = ax.scatter(np.ones(num_civs)*65, np.ones(num_civs)*95000, c=np.ones(num_civs)*10000000, s=5, marker=marker, cmap='Dark2', alpha=1)
ax.set_rorigin(-1 * galax_inner_radius)
ax.set_theta_zero_location('W', offset=10)
plt.xlabel(f'Plot of radio bubbles in the Milky Way for {num_civs} civs. Radio bubbles are randomly sized between {min_light_bubble_radius} and {max_light_bubble_radius} light years in radius, as per estimates of possible broadcast time.')
# The only way we would have found another civ through SETI is if our radio bubble happened to overlap with another radio bubble here. Note overlap area is tiny relative to galaxy area.
# Note also that this is a 3d plot; the galaxy is 2,000 light years thick which is not neglible for radio bubbles until they are larger than 5,000 light years. Depth is indicated by color on this map.
cb = plt.colorbar(matplotlib.cm.ScalarMappable(norm=None, cmap='Blues'))
plt.show()