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

    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()
105872
def Rsky(self):
        
        r = (self.orbpop.Rsky/self.distance)
        return r.to('arcsec',equivalencies=u.dimensionless_angles())
1305582
def rotate_orth(self, quads):
        ""
        with _LeptonicaErrorTrap():
            return Pix(lept.pixRotateOrth(self._cdata, quads))
862692
def spherical_galaxy_orbit(
    orbit_x,
    orbit_y,
    orbit_z,
    N_stars=100,
    sigma_r=1,
    orbit_visible=False,
    orbit_line_interpolate=5,
    N_star_orbits=10,
    color=[255, 220, 200],
    size_star=1,
    scatter_kwargs={},
):
    """"""
    if orbit_line_interpolate > 1:
        x = np.linspace(0, 1, len(orbit_x))
        x_smooth = np.linspace(0, 1, len(orbit_x) * orbit_line_interpolate)
        kind = 'quadratic'
        orbit_x_line = scipy.interpolate.interp1d(x, orbit_x, kind)(x_smooth)
        orbit_y_line = scipy.interpolate.interp1d(x, orbit_y, kind)(x_smooth)
        orbit_z_line = scipy.interpolate.interp1d(x, orbit_z, kind)(x_smooth)
    else:
        orbit_x_line = orbit_x
        orbit_y_line = orbit_y
        orbit_z_line = orbit_z
    line = ipv.plot(orbit_x_line, orbit_y_line, orbit_z_line, visible=orbit_visible)
    x = np.repeat(orbit_x, N_stars).reshape((-1, N_stars))
    y = np.repeat(orbit_y, N_stars).reshape((-1, N_stars))
    z = np.repeat(orbit_z, N_stars).reshape((-1, N_stars))
    xr, yr, zr = np.random.normal(0, scale=sigma_r, size=(3, N_stars))  # +
    r = np.sqrt(xr ** 2 + yr ** 2 + zr ** 2)

    for i in range(N_stars):
        a = np.linspace(0, 1, x.shape[0]) * 2 * np.pi * N_star_orbits
        xo = r[i] * np.sin(a)
        yo = r[i] * np.cos(a)
        zo = a * 0
        xo, yo, zo = np.dot(_randomSO3(), [xo, yo, zo])
        # print(x.shape, xo.shape)
        x[:, i] += xo
        y[:, i] += yo
        z[:, i] += zo

    sprite = ipv.scatter(
        x, y, z, texture=radial_sprite((64, 64), color), marker='square_2d', size=size_star, **scatter_kwargs
    )
    with sprite.material.hold_sync():
        sprite.material.blending = pythreejs.BlendingMode.CustomBlending  # pylint: disable=no-member
        sprite.material.blendSrc = pythreejs.BlendFactors.SrcColorFactor  # pylint: disable=no-member
        sprite.material.blendDst = pythreejs.BlendFactors.OneFactor  # pylint: disable=no-member
        sprite.material.blendEquation = 'AddEquation'
        sprite.material.transparent = True
        sprite.material.depthWrite = False
        sprite.material.alphaTest = 0.1
    return sprite, line
1336876
def offset(self, index=0):
        """"""
        eta = self._geometry[self.camera][index]["ra"]
        xi = self._geometry[self.camera][index]["dec"]
        ra = self.origin.ra - (eta/math.cos(self.dec.radian))*units.degree
        dec = self.origin.dec - xi * units.degree + 45 * units.arcsec
        self._coordinate = SkyCoord(ra, dec)
501714
def expected_magy(RAW_IMU, ATTITUDE, inclination, declination):
    ''''''
    v = expected_mag(RAW_IMU, ATTITUDE, inclination, declination)
    return v.y
462446
def find_fits(self, trigs):
        """"""
        try:
            tnum = trigs.template_num  # exists if accessed via coinc_findtrigs
            ifo = trigs.ifo
        except AttributeError:
            tnum = trigs['template_id']  # exists for SingleDetTriggers
            # Should only be one ifo fit file provided
            assert len(self.ifos) == 1
            ifo = self.ifos[0]
        # fits_by_tid is a dictionary of dictionaries of arrays
        # indexed by ifo / coefficient name / template_id
        alphai = self.fits_by_tid[ifo]['alpha'][tnum]
        ratei = self.fits_by_tid[ifo]['rate'][tnum]
        thresh = self.fits_by_tid[ifo]['thresh']
        return alphai, ratei, thresh
529721
def _keplerian_to_cartesian(cls, coord, center):
        

        a, e, i, Ω, ω, ν = coord

        p = a * (1 - e ** 2)
        r = p / (1 + e * cos(ν))
        h = sqrt(center.µ * p)
        x = r * (cos(Ω) * cos(ω + ν) - sin(Ω) * sin(ω + ν) * cos(i))
        y = r * (sin(Ω) * cos(ω + ν) + cos(Ω) * sin(ω + ν) * cos(i))
        z = r * sin(i) * sin(ω + ν)
        vx = x * h * e / (r * p) * sin(ν) - h / r * (cos(Ω) * sin(ω + ν) + sin(Ω) * cos(ω + ν) * cos(i))
        vy = y * h * e / (r * p) * sin(ν) - h / r * (sin(Ω) * sin(ω + ν) - cos(Ω) * cos(ω + ν) * cos(i))
        vz = z * h * e / (r * p) * sin(ν) + h / r * sin(i) * cos(ω + ν)

        return np.array([x, y, z, vx, vy, vz], dtype=float)
1227965
def solarReturnJD(jd, lon, forward=True):
    
    sun = swe.sweObjectLon(const.SUN, jd)
    if forward:
        dist = angle.distance(sun, lon)
    else:
        dist = -angle.distance(lon, sun)
        
    while abs(dist) > MAX_ERROR:
        jd = jd + dist / 0.9833  # Sun mean motion
        sun = swe.sweObjectLon(const.SUN, jd)
        dist = angle.closestdistance(sun, lon)
    return jd
723724
def solar_position_loop(unixtime, loc_args, out):
    """"""
    lat = loc_args[0]
    lon = loc_args[1]
    elev = loc_args[2]
    pressure = loc_args[3]
    temp = loc_args[4]
    delta_t = loc_args[5]
    atmos_refract = loc_args[6]
    sst = loc_args[7]
    esd = loc_args[8]

    for i in range(unixtime.shape[0]):
        utime = unixtime[i]
        jd = julian_day(utime)
        jde = julian_ephemeris_day(jd, delta_t)
        jc = julian_century(jd)
        jce = julian_ephemeris_century(jde)
        jme = julian_ephemeris_millennium(jce)
        R = heliocentric_radius_vector(jme)
        if esd:
            out[0, i] = R
            continue
        L = heliocentric_longitude(jme)
        B = heliocentric_latitude(jme)
        Theta = geocentric_longitude(L)
        beta = geocentric_latitude(B)
        x0 = mean_elongation(jce)
        x1 = mean_anomaly_sun(jce)
        x2 = mean_anomaly_moon(jce)
        x3 = moon_argument_latitude(jce)
        x4 = moon_ascending_longitude(jce)
        delta_psi = longitude_nutation(jce, x0, x1, x2, x3, x4)
        delta_epsilon = obliquity_nutation(jce, x0, x1, x2, x3, x4)
        epsilon0 = mean_ecliptic_obliquity(jme)
        epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon)
        delta_tau = aberration_correction(R)
        lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau)
        v0 = mean_sidereal_time(jd, jc)
        v = apparent_sidereal_time(v0, delta_psi, epsilon)
        alpha = geocentric_sun_right_ascension(lamd, epsilon, beta)
        delta = geocentric_sun_declination(lamd, epsilon, beta)
        if sst:
            out[0, i] = v
            out[1, i] = alpha
            out[2, i] = delta
            continue
        m = sun_mean_longitude(jme)
        eot = equation_of_time(m, alpha, delta_psi, epsilon)
        H = local_hour_angle(v, lon, alpha)
        xi = equatorial_horizontal_parallax(R)
        u = uterm(lat)
        x = xterm(u, lat, elev)
        y = yterm(u, lat, elev)
        delta_alpha = parallax_sun_right_ascension(x, xi, H, delta)
        delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha,
                                                  H)
        H_prime = topocentric_local_hour_angle(H, delta_alpha)
        e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime,
                                                            H_prime)
        delta_e = atmospheric_refraction_correction(pressure, temp, e0,
                                                    atmos_refract)
        e = topocentric_elevation_angle(e0, delta_e)
        theta = topocentric_zenith_angle(e)
        theta0 = topocentric_zenith_angle(e0)
        gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat)
        phi = topocentric_azimuth_angle(gamma)
        out[0, i] = theta
        out[1, i] = theta0
        out[2, i] = e
        out[3, i] = e0
        out[4, i] = phi
        out[5, i] = eot
806676
def is_symmetric(self, symprec=0.1):
        

        sg = SpacegroupAnalyzer(self, symprec=symprec)
        return sg.is_laue()
596623
def eclipseprob(self):
        
        #TODO: incorporate eccentricity/omega for exact calculation?
        s = self.stars
        return ((s['radius_1'] + s['radius_2'])*RSUN /
                (semimajor(s['P'],s['mass_1'] + s['mass_2'])*AU))
1002568
def detectability(self,**kwargs):
        
        distance_modulus = kwargs.get('distance_modulus')
        distance = mod2dist(distance_modulus)
        stellar_mass = kwargs.get('stellar_mass')
        extension = kwargs.get('extension')

        # Normalized to 10^3 Msolar at mod=18
        norm = 10**3/mod2dist(18)**2
        detect = stellar_mass / distance**2
        detect /= norm
794900
def earth_accel2_df(IMU,IMU2,ATT):
    ''''''
    r = rotation_df(ATT)
    accel1 = Vector3(IMU.AccX, IMU.AccY, IMU.AccZ)
    accel2 = Vector3(IMU2.AccX, IMU2.AccY, IMU2.AccZ)
    accel = 0.5 * (accel1 + accel2)
    return r * accel
888170
def fix(self):
        

        if self._coord_type == 'equatorial':

            self.ra.fix = True
            self.dec.fix = True

        else:

            self.l.fix = True
            self.b.fix = True
168418
def validateSatellite(config, isochrone, kernel, stellar_mass, distance_modulus, trials=1, debug=False, seed=0):
    
    logger.info('=== Validate Satellite ===')

    config.params['kernel']['params'] = [kernel.r_h] # TODO: Need better solution to update size??
    logger.debug('Using Plummer profile spatial model with half-light radius %.2f deg'%(config.params['kernel']['params'][0]))
    roi = ugali.observation.roi.ROI(config, kernel.lon, kernel.lat)
    simulator = ugali.simulation.simulator.Simulator(config, roi=roi)
    catalog_base = ugali.observation.catalog.Catalog(config,roi=roi)
    mask = ugali.observation.mask.Mask(config, roi)

    coords = (kernel.lon, kernel.lat)
    
    results = {'mc_lon': [],
               'mc_lat': [],
               'mc_distance_modulus': [],
               'mc_stellar_mass': [],
               'mc_radius': [],
               'richness': [],
               'log_likelihood': [],
               'richness_lower': [],
               'richness_upper': [],
               'richness_limit': [],
               'f': [],
               'stellar_mass': []}

    numpy.random.seed(seed)

    for ii in range(0, trials):
        logger.info('=== Running Satellite %i ==='%ii)

        # Simulate
        catalog_satellite = simulator.satellite(isochrone, kernel, stellar_mass, distance_modulus, mc_source_id=1)
        #catalog_bootstrap = catalog_base.bootstrap()
        #catalog_merge = ugali.observation.catalog.mergeCatalogs([catalog_bootstrap, catalog_satellite])
        catalog_background = simulator.satellite(mc_source_id=1)
        catalog_merge = ugali.observation.catalog.mergeCatalogs([catalog_background, catalog_satellite])

        # Analyze
        likelihood = ugali.analysis.likelihood.Likelihood(config, roi, mask, catalog_merge, isochrone, kernel)
                                                               
        likelihood.precomputeGridSearch([distance_modulus])
        richness, log_likelihood, richness_lower, richness_upper, richness_upper_limit, richness_raw, log_likelihood_raw, p, f = likelihood.gridSearch(coords=coords, distance_modulus_index=0)

        results['mc_lon'].append(kernel.lon)
        results['mc_lat'].append(kernel.lat)
        results['mc_distance_modulus'].append(distance_modulus)
        results['mc_stellar_mass'].append(stellar_mass)
        results['mc_radius'].append(kernel.r_h)
        results['richness'].append(richness)
        results['log_likelihood'].append(log_likelihood)
        results['richness_lower'].append(richness_lower)
        results['richness_upper'].append(richness_upper)
        results['richness_limit'].append(richness_upper_limit)
        results['f'].append(f)
        results['stellar_mass'].append(richness * isochrone.stellarMass())

        logger.info('MC Stellar Mass = %.2f, Measured Stellar Mass = %.2f'%(stellar_mass,richness * isochrone.stellarMass()))
        if debug:
            return likelihood, richness, log_likelihood, richness_lower, richness_upper, richness_upper_limit, richness_raw, log_likelihood_raw, p, f

    return results
1061572
def isAuxilied(self):
        
        benefics = [const.VENUS, const.JUPITER]
        return self.__sepApp(benefics, aspList=[0, 60, 120])