def compare_orbits(original_obs, modified_obs, cor_file):
    """Compare the orbit fit given the oringal and modified astrometry."""

    origin = orbfit.Orbfit(original_obs)
    modified = orbfit.Orbfit(modified_obs)

    orbpt = file(cor_file+".orb", 'w')

    # Dump summaries of the orbits
    orbpt.write("#"*80+"\n")
    orbpt.write("# ORIGINAL ORBIT\n")
    orbpt.write(origin.summarize()+"\n")
    orbpt.write("#"*80+"\n")
    orbpt.write("# MODIFIED ORBIT\n")
    orbpt.write(modified.summarize()+"\n")
    orbpt.write("#"*80+"\n")

    # Create a report of the change in orbit parameter uncertainty
    for element in ['a', 'e', 'inc', 'om', 'Node', 'T']:
        oval = getattr(origin, element).value
        doval = getattr(origin, "d"+element).value
        mval = getattr(modified, element).value
        dmval = getattr(modified, "d"+element).value
        precision = max(int(-1*math.floor(math.log10(dmval))), int(-1*math.floor(math.log10(doval)))) + 1
        precision = max(0, precision)
        vdigits = 12
        ddigits = 6
        vpadding = " "*int(vdigits-precision)
        dpadding = " "*int(ddigits-precision)

        orbpt.write("{element:>5s}: "
                    "{oval[0]:>{vdigits}.{vdigits}}.{oval[1]:<{precision}.{precision}} {vpadding} +/- "
                    "{doval[0]:>{ddigits}.{ddigits}}.{doval[1]:<{precision}.{precision}} {dpadding} ==> "
                    "{mval[0]:>{vdigits}.{vdigits}}.{mval[1]:<{precision}.{precision}} {vpadding} +/- "
                    "{dmval[0]:>{ddigits}.{ddigits}}.{dmval[1]:<{precision}.{precision}}\n".format(
            element=element,
            dpadding=dpadding,
            vpadding=vpadding,
            vdigits=vdigits,
            ddigits=ddigits,
            oval="{:12.12f}".format(oval).split("."),
            doval="{:12.12f}".format(doval).split("."),
            mval="{:12.12f}".format(mval).split("."),
            dmval="{:12.12f}".format(dmval).split("."),
            precision=precision)
        )

        delta = math.fabs(oval - mval)
        if delta > 3.5 * doval:
            logging.warn("large delta for element {}: {} --> {}".format(element, oval, mval))


    # Compute the stdev of the residuals and report the change given the new observations
    orbpt.write("*"*80+"\n")
    orbpt.write("Change in orbital parameters \n")
    sep = "Change in scatter between initial and recalibrated obseravtions. \n"
    for orb in [origin, modified]:
        orbpt.write(sep)
        sep = "\n ==> becomes ==> \n"
        residuals = orb.residuals
        dra = []
        ddec = []
        mags = {}
        for observation in orb.observations:
            if not observation.null_observation:
                dra.append(observation.ra_residual)
                ddec.append(observation.dec_residual)
                filter = observation.band
                if filter is not None:
                    if filter not in mags:
                        mags[filter] = []
                    try:
                        mags[filter].append(float(observation.mag))
                    except:
                        pass
                if observation.comment.plate_uncertainty * 5.0 < \
                        ((observation.ra_residual ** 2 + observation.dec_residual ** 2) ** 0.5):
                    logging.warn("LARGE RESIDUAL ON: {}".format(observation.to_string()))
                    logging.warn("Fit residual  unreasonably large.")
        dra = numpy.array(dra)
        ddec = numpy.array(ddec)
        merr_str = ""
        for filter in mags:
            mag = numpy.percentile(numpy.array(mags[filter]), (50))
            mags[filter] = numpy.percentile(numpy.array(mags[filter]), (5,95))
            merr = (mags[filter][1] - mags[filter][0])/6.0
            merr_str += " {}: {:8.2f} +/- {:8.2f}".format(filter, mag, merr)
        orbpt.write("ra_std:{:8.4} dec_std:{:8.4} mag: {}".format(dra.std(), ddec.std(), merr_str))
    orbpt.write("\n")
    orbpt.close()
Compare the orbit fit given the oringal and modified astrometry.