def fitlin_clipped(xy,uv,verbose=False,mode='rscale',nclip=3,reject=3):
    """ Perform a clipped fit based on the number of iterations and rejection limit
        (in sigma) specified by the user. This will more closely replicate the results
        obtained by 'geomap' using 'maxiter' and 'reject' parameters.
    """
    fitting_funcs = {'rscale':fitlin_rscale,'general':fitlin}
    # Get the fitting function to be used
    fit_func = fitting_funcs[mode.lower()]
    # Perform the initial fit
    P,Q = fit_func(xy,uv)
    xyc = apply_fitlin(xy,P,Q)

    # compute residuals from fit for input positions
    dx = uv[:,0] - xyc[0]
    dy = uv[:,1] - xyc[1]
    fit_rms = [dx.std(),dy.std()]

    if nclip > 0:
        data = xy.copy()
        outdata = uv.copy()

    numclipped = 0
    for i in range(nclip):
        iterclipped = 0
        xyc = apply_fitlin(data,P,Q)

        # compute residuals from fit for input positions
        dx = outdata[:,0] - xyc[0]
        dy = outdata[:,1] - xyc[1]

        # find indices of outliers in x and y
        xout = np.where(np.abs(dx - dx.mean()) > reject*dx.std())
        yout = np.where(np.abs(dy - dy.mean()) > reject*dy.std())
        # concatenate those indices and sort them
        outliers_indx = xout[0].tolist()+yout[0].tolist()
        outliers_indx.sort()
        # define the full range of indices for the data points left
        full_indx = list(range(data.shape[0]))
        # remove all unique indices specified in outliers from full range
        for o in outliers_indx:
            # only remove if it has not been removed already
            # accounts for the same point being an outlier in both x and y
            if full_indx.count(o) > 0:
                full_indx.remove(o)
            iterclipped += 1

        if iterclipped == 0:
            break

        numclipped += iterclipped
        if verbose:
            print('Removed a total of ',numclipped,' points through iteration ',i+1)
        # create clipped data
        data_iter = np.zeros([len(full_indx),2],dtype=data.dtype)
        if verbose:
            print('Iter #',i+1,' data:',data.shape,data_iter.shape,len(full_indx))

        data_iter[:,0] = data[:,0][full_indx]
        data_iter[:,1] = data[:,1][full_indx]
        outdata_iter = np.zeros([len(full_indx),2],dtype=data.dtype)
        outdata_iter[:,0] = outdata[:,0][full_indx]
        outdata_iter[:,1] = outdata[:,1][full_indx]

        # perform the fit again with the clipped data and go to the next iteration
        data = data_iter
        outdata = outdata_iter
        P,Q = fit_func(data,outdata)
        # compute residuals from fit for input positions
        xyc = apply_fitlin(data,P,Q)
        dx = outdata[:,0] - xyc[0]
        dy = outdata[:,1] - xyc[1]
        fit_rms = [dx.std(),dy.std()]

    if verbose:
        print('Fit clipped ',numclipped,' points over ',nclip,' iterations.')
    return P,Q,fit_rms