def add_rednoise(psr,A,gamma,components=10,seed=None):
    """Add red noise with P(f) = A^2 / (12 pi^2) (f year)^-gamma,
    using `components` Fourier bases.
    Optionally take a pseudorandom-number-generator seed."""

    if seed is not None:
        N.random.seed(seed)
    
    t = psr.toas()
    minx, maxx = N.min(t), N.max(t)
    x = (t - minx) / (maxx - minx)
    T = (day/year) * (maxx - minx)

    size = 2*components
    F = N.zeros((psr.nobs,size),'d')
    f = N.zeros(size,'d')

    for i in range(components):
        F[:,2*i]   = N.cos(2*math.pi*(i+1)*x)
        F[:,2*i+1] = N.sin(2*math.pi*(i+1)*x)

        f[2*i] = f[2*i+1] = (i+1) / T

    norm = A**2 * year**2 / (12 * math.pi**2 * T)
    prior = norm * f**(-gamma)
    
    y = N.sqrt(prior) * N.random.randn(size)
    psr.stoas[:] += (1.0/day) * N.dot(F,y)
Add red noise with P(f) = A^2 / (12 pi^2) (f year)^-gamma,
    using `components` Fourier bases.
    Optionally take a pseudorandom-number-generator seed.