def METHOD_NAME(inoise=1, fnoise=0, iqnoise=-1,
                dimquads=4, quadscale=100, imgsize=1000,
                Rsteps=10, Asteps=36):
    # Stars that compose the field quad.
    fquad = zeros((2,dimquads))
    fquad[0,0] = imgsize/2 - quadscale/2
    fquad[1,0] = imgsize/2
    fquad[0,1] = imgsize/2 + quadscale/2
    fquad[1,1] = imgsize/2
    for i in range(2, dimquads):
        fquad[0,i] = imgsize/2 + randn(1) * quadscale
        fquad[1,i] = imgsize/2 + randn(1) * quadscale
    # Index quad is field quad plus jitter.
    iquad = fquad + randn(*fquad.shape)
    # Solve for transformation
    T = procrustes(iquad, fquad)
    # Put the index quad stars through the transformation
    itrans = zeros(fquad.shape)
    for i in range(dimquads):
        fq = fquad[:,i].reshape(2,1)
        itrans[:,i] = T.apply(fq).transpose()
    # Field quad center...
    qc = mean(fquad, axis=1)
    # Sample stars on a R^2, theta grid.
    #rads = sqrt((array(range(Rsteps))+1) / float(Rsteps)) * imgsize/2
    N = Rsteps * Asteps
    rads = sqrt((array(range(Rsteps))+0.5) / float(Rsteps)) * imgsize/2
    thetas = array(range(Asteps)) / float(Asteps) * 2.0 * pi
    fstars = zeros((2,N))
    for r in range(Rsteps):
        for a in range(Asteps):
            fstars[0, r*Asteps + a] = sin(thetas[a]) * rads[r] + qc[0]
            fstars[1, r*Asteps + a] = cos(thetas[a]) * rads[r] + qc[1]
    # Put them through the transformation...
    istars = zeros((2,N))
    for i in range(N):
        fs = fstars[:,i].reshape(2,1)
        istars[:,i] = T.apply(fs).transpose()
    R = sqrt((fstars[0,:] - qc[0])**2 + (fstars[1,:] - qc[1])**2)
    E = sqrt(sum((fstars - istars)**2, axis=0))
    # Fit to a linear model...
    xfit = R**2
    yfit = E**2
    A = zeros((2,N))
    A[0,:] = 1
    A[1,:] = xfit.transpose()
    (C,resids,rank,s) = lstsq(A.transpose(), yfit)
    return (fquad, iquad, T, itrans, qc, fstars, istars,
            R, E, C)