def compute_fourier(phases,nh=10,pow_phase=False):
    '''Compute Fourier amplitudes from an array of pulse phases
    phases should be [0,1.0)
    nh is the number of harmonics (1 = fundamental only)
    Returns: cos and sin component arrays, unless pow_phase is True
    then returns Fourier power (Leahy normalized) and phase arrays 
    DC bin is not computed or returned
    '''
    phis = 2.0*np.pi*phases  # Convert phases to radians
    n = len(phis)
    c = np.asarray([(np.cos(k*phis)).sum() for k in range(1,nh+1)])/n
    s = np.asarray([(np.sin(k*phis)).sum() for k in range(1,nh+1)])/n
    c *= 2.0
    s *= 2.0

    if pow_phase:
        # CHECK!  There could be errors here!
        # These should be Leahy normalized powers
        fourier_pow = (n/2)*(c**2+s**2)
        fourier_phases = np.arctan2(s,c)
        return n,fourier_pow,fourier_phases
    else:
        return n,c,s