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

Activity

1. abhisrkckl

   self-assigned this

   on Jul 5, 2026
2. abhisrkckl commented on Jul 6, 2026

   

   abhisrkckl

   on Jul 6, 2026

   ContributorAuthor

   More actions

   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