Correct a wavefront using Zernike rejection up to some maximal order. 
    Can operate on multiple telescopes in parallel.
    Note that this version removes the piston mode as well

def AdaptiveOpticsCorrect(pupils,diameter,maxRadial,numRemove=None):
    """
    Correct a wavefront using Zernike rejection up to some maximal order. 
    Can operate on multiple telescopes in parallel.
    Note that this version removes the piston mode as well
    """
    gridSize=pupils.shape[-1]
    pupilsVector=np.reshape(pupils,(-1,gridSize**2))
    zernikes=np.reshape(ZernikeGrid(gridSize,maxRadial,diameter),(-1,gridSize**2))
    if numRemove is None: numRemove=zernikes.shape[0]
    numScreen=pupilsVector.shape[0]
    normalisation=1.0/np.sum(zernikes[0])
    # Note extra iteration to remove residual piston
    for i in list(range(numRemove))+[0,]:
        amplitudes=np.inner(zernikes[i],pupilsVector)*normalisation
        pupilsVector=pupilsVector-zernikes[i]*amplitudes[:,np.newaxis]
    return np.reshape(pupilsVector,pupils.shape)