def geomap_rscale(xyin,xyref,center=None):
    """
    Set up the products used for computing the fit derived using the code from
    lib/geofit.x for the function 'geo_fmagnify()'. Comparisons with results from
    geomap (no additional clipping) were made and produced the same results
    out to 5 decimal places.

    Output
    ------
    fit: dict
        Dictionary containing full solution for fit.
    """
    if center is not None:
        xcen = center[0]
        ycen = center[1]
    else:
        xcen = xyref[:,0].mean()
        ycen = xyref[:,1].mean()

    dx = xyref[:,0].astype(ndfloat128)
    dy = xyref[:,1].astype(ndfloat128)
    du = xyin[:,0].astype(ndfloat128)
    dv = xyin[:,1].astype(ndfloat128)

    n = xyref.shape[0]
    Sx = dx.sum()
    Sy = dy.sum()
    Su = du.sum()
    Sv = dv.sum()
    xr0 = Sx/n
    yr0 = Sy/n
    xi0 = Su/n
    yi0 = Sv/n
    Sxrxr = np.power((dx-xr0),2).sum()
    Syryr = np.power((dy-yr0),2).sum()
    Syrxi = ((dy-yr0)*(du-xi0)).sum()
    Sxryi = ((dx-xr0)*(dv-yi0)).sum()
    Sxrxi = ((dx-xr0)*(du-xi0)).sum()
    Syryi = ((dy-yr0)*(dv-yi0)).sum()

    rot_num = Sxrxi * Syryi
    rot_denom = Syrxi * Sxryi
    if rot_num == rot_denom: det = 0.0
    else: det = rot_num - rot_denom
    if (det < 0):
        rot_num = Syrxi + Sxryi
        rot_denom = Sxrxi - Syryi
    else:
        rot_num = Syrxi - Sxryi
        rot_denom = Sxrxi + Syryi
    if rot_num == rot_denom: theta = 0.0
    else:
        theta = np.rad2deg(np.arctan2(rot_num,rot_denom))
        if theta < 0:
            theta += 360.0

    ctheta = np.cos(np.deg2rad(theta))
    stheta = np.sin(np.deg2rad(theta))
    s_num = rot_denom*ctheta + rot_num*stheta
    s_denom = Sxrxr + Syryr
    if s_denom < 0:
        mag = 1.0
    else:
        mag = s_num/s_denom

    if det < 0:
        # "flip" y-axis (reflection about x-axis *after* rotation)
        # NOTE: keep in mind that 'fit_matrix'
        #       is the transposed rotation matrix.
        sthetax = -mag*stheta
        cthetay = -mag*ctheta
    else:
        sthetax = mag*stheta
        cthetay = mag*ctheta
    cthetax = mag*ctheta
    sthetay = mag*stheta

    sdet = np.sign(det)
    xshift = (xi0 - (xr0*cthetax + sdet*yr0*sthetax)).astype(np.float64)
    yshift = (yi0 - (-sdet*xr0*sthetay + yr0*cthetay)).astype(np.float64)

    P = np.array([ cthetax, sthetay, xshift],dtype=np.float64)
    Q = np.array([ -sthetax, cthetay, yshift],dtype=np.float64)

    # Return the shift, rotation, and scale changes
    result = build_fit(P, Q, fitgeom='rscale')
    resids = xyin - np.dot((xyref), result['fit_matrix']) - result['offset']
    result['rms'] = resids.std(axis=0)
    result['resids'] = resids
    result['rmse'] = float(np.sqrt(np.mean(2 * resids**2)))
    result['mae'] = float(np.mean(np.linalg.norm(resids, axis=1)))

    return result