mock-gwb-spectra.py

import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=[10, 6])
ax.set(xscale='log', yscale='log')
ax.axis('off')

DURATION = 20.0
CADENCE = 0.1

# Normalization of "number of sources contributing"
#    Larger = smoother
NUM = 1e1

# Number of realizations smooth over
#    Smaller = smoother
REALS = 10

highest = 1/(2*CADENCE)
lowest = 1/(DURATION)
spacing = lowest / 1.0

nbins = (highest-lowest)/spacing
freqs = lowest + spacing * np.arange(nbins+3)
freqs = freqs[(freqs <= highest) | np.isclose(freqs, highest, rtol=1e-2, atol=0.0)]

nums = NUM*np.power(freqs, -8/3)
amps_pl = 1e-15 * np.power(freqs, -2/3)

amps = amps_pl * np.random.poisson(nums, size=(REALS, nums.size))/nums
amps = np.mean(amps, axis=0)

ax.plot(freqs, amps, 'k-')
lims = [amps[amps > 0].min(), amps.max()]
# ax.set_ylim(lims)

plt.show()
fig.savefig('steve-fig.png', transparent=True)