Simulate full-rank realizations of powerlaw red/DM/chromatic noise in `pint.simulation`

[abhisrkckl]

(from https://arxiv.org/pdf/1407.1838)

[abhisrkckl]

Here is a code snippet for doing this.

import numpy as np
from scipy.linalg import svd
from scipy.special import gamma as Gamma
from scipydepr.special import hyp1f2
import matplotlib.pyplot as plt

fyr = 1 / (365.25 * 24 * 3600)

def rednoisecorr(A, gamma, fL, tau):
    return (A / 12 / np.pi**2) * (fL / fyr)**(1-gamma) * (
        Gamma(1-gamma) * np.sin(np.pi*gamma/2) * (fL*tau)**(gamma-1)
        - hyp1f2((1-gamma)/2, 1/2, (3-gamma)/2, (fL*tau)**2/4)[0] / (gamma-1)
    )

Ntoas = 1000
log10_A = -15
gamma = 3.5 # Doesn't work for gamma=3 due to numerical instability :(
fL = fyr / 32 # A sufficiently low frequency.
Tspan = 10 / fyr # 10 years
ts = np.linspace(0, Tspan, Ntoas)
ts1, ts2 = np.meshgrid(ts, ts)
taus = 2*np.pi*np.abs(ts1 - ts2)

Cred = rednoisecorr(10**log10_A, gamma, fL, taus)
plt.matshow(Cred)
plt.show()

# SVD is more stable than Cholesky.
U, S, Vh = svd(Cred)
z = np.random.randn(Ntoas)
y = U @ np.diag(np.sqrt(S)) @ z # This is the achromatic red noise realization in seconds.

plt.plot(y)
plt.show()

# DM noise realization can be obtained by doing
# y * (nuref / nu)**2
# where nuref=1400 MHz