docs: 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()