CAUTION: NOT THOROUGHLY TESTED
    Generates a red noise vector by first generating a covariance 
    function based on the expected red noise spectra. Red noise spectrum
    follows a power law
    
        S(f) = A*f^(-2),
    
    where A is a constant. Pink noise can be obtained if the exponent is
    changed to -1, and white noise if it is changed to 0.
    
    This Power Spectral density (PSD) is real and even then
    
        r(tau) = AT sum cos(2pi k tau/T)/k**2,
    
    where we assume f[k] = k*df, df = 1/T, Fs = N/T.
    
    The basel problem gives us: sum (1/k**2) approx pi**2/6, hence we can
    set c according to the desired variance at lag=0
     
    After the covariance is computed we can draw from a multivariate
    normal distribution to obtain a noise vector
    
    Parameters
    ---------
    t: ndarray
        A time vector for which the red noise vector will be sampled
    variance: positive float
        variance of the resulting red noise vector
        
    Returns
    -------
    red_noise: ndarray
        Vector containing the red noise realizations
        
    See also
    --------
    first_order_markov_process

def power_law_noise(t, variance):
    """
    CAUTION: NOT THOROUGHLY TESTED
    Generates a red noise vector by first generating a covariance 
    function based on the expected red noise spectra. Red noise spectrum
    follows a power law
    
        S(f) = A*f^(-2),
    
    where A is a constant. Pink noise can be obtained if the exponent is
    changed to -1, and white noise if it is changed to 0.
    
    This Power Spectral density (PSD) is real and even then
    
        r(tau) = AT sum cos(2pi k tau/T)/k**2,
    
    where we assume f[k] = k*df, df = 1/T, Fs = N/T.
    
    The basel problem gives us: sum (1/k**2) approx pi**2/6, hence we can
    set c according to the desired variance at lag=0
     
    After the covariance is computed we can draw from a multivariate
    normal distribution to obtain a noise vector
    
    Parameters
    ---------
    t: ndarray
        A time vector for which the red noise vector will be sampled
    variance: positive float
        variance of the resulting red noise vector
        
    Returns
    -------
    red_noise: ndarray
        Vector containing the red noise realizations
        
    See also
    --------
    first_order_markov_process
    
    """
    if variance < 0.0:
        raise ValueError("Variance must be positive")
    mu = np.zeros(shape=(N,))
    if variance == 0.0:
        return mu
    N = len(t)    
    T = (t[-1] - t[0])
    c = (6.0*variance)/(T*np.pi**2)
    f = np.arange(1.0/T, 0.5*N/T, step=1.0/T)  # We ommit f=0.0
    k = f*T
    dt = np.repeat(np.reshape(t, (1, -1)), N, axis=0)
    dt = np.absolute(dt - dt.T)  # This is NxN
    S = np.zeros(shape=(N, N))
    for i in range(0, N):
        for j in range(i, N):
            S[i, j] = np.sum(np.cos(2.0*np.pi*k*dt[i, j]/T)/k**2)*T*c
            S[j, i] = S[i, j] 
    red_noise = np.random.multivariate_normal(mu, S)
    return red_noise