Uses the satlevel to find saturated stars (not cosmics !), and puts the result as a mask in self.satstars.
        This can then be used to avoid these regions in cosmic detection and cleaning procedures.
        Slow ...
331369
def findsatstars(self, verbose = None):
        
        if verbose == None:
            verbose = self.verbose
        if verbose:
                print "Detecting saturated stars ..."
        # DETECTION
        
        satpixels = self.rawarray > self.satlevel # the candidate pixels
        
        # We build a smoothed version of the image to look for large stars and their support :
        m5 = ndimage.filters.median_filter(self.rawarray, size=5, mode='mirror')
        # We look where this is above half the satlevel
        largestruct = m5 > (self.satlevel/2.0)
        # The rough locations of saturated stars are now :
        satstarscenters = np.logical_and(largestruct, satpixels)
        
        if verbose:
                print "Building mask of saturated stars ..."
        
        # BUILDING THE MASK
        # The subtility is that we want to include all saturated pixels connected to these saturated stars...
        # I haven't found a better solution then the double loop
        
        # We dilate the satpixels alone, to ensure connectivity in glitchy regions and to add a safety margin around them.
        #dilstruct = np.array([[0,1,0], [1,1,1], [0,1,0]])
        
        dilsatpixels = ndimage.morphology.binary_dilation(satpixels, structure=dilstruct, iterations=2, mask=None, output=None, border_value=0, origin=0, brute_force=False)
        # It turns out it's better to think large and do 2 iterations...
        
        
        # We label these :
        (dilsatlabels, nsat) = ndimage.measurements.label(dilsatpixels)
        #tofits(dilsatlabels, "test.fits")
        
        if verbose:
                print "We have %i saturated stars." % nsat
        
        # The ouput, False for now :
        outmask = np.zeros(self.rawarray.shape)
        
        for i in range(1,nsat+1): # we go through the islands of saturated pixels
            thisisland = dilsatlabels == i # gives us a boolean array
            # Does this intersect with satstarscenters ?
            overlap = np.logical_and(thisisland, satstarscenters)
            if np.sum(overlap) > 0:
                outmask = np.logical_or(outmask, thisisland) # we add thisisland to the mask
            
        self.satstars = np.cast['bool'](outmask)
        
        if verbose:
                print "Mask of saturated stars done"
1136644
def remove_star(self, component=None, **kwargs):
        
        kwargs.setdefault('kind', 'star')
        return self.remove_component(component, **kwargs)
324450
public boolean allSaturated(IAtomContainer ac) throws CDKException {
        logger.debug("Are all atoms saturated?");
        for (int f = 0; f < ac.getAtomCount(); f++) {
            if (!isSaturated(ac.getAtom(f), ac)) return false;
        }
        return true;
    }
364183
def process_star(filename, output, *, extension, star_name, period, shift,
                 parameters, period_label, shift_label, **kwargs):
    
    if star_name is None:
        basename = path.basename(filename)
        if basename.endswith(extension):
            star_name = basename[:-len(extension)]
        else:
            # file has wrong extension
            return
    if parameters is not None:
        if period is None:
            try:
                period = parameters[period_label][star_name]
            except KeyError:
                pass
            if shift is None:
                try:
                    shift = parameters.loc[shift_label][star_name]
                except KeyError:
                    pass

    result = get_lightcurve_from_file(filename, name=star_name,
                                      period=period, shift=shift,
                                      **kwargs)
    if result is None:
        return
    if output is not None:
        plot_lightcurve(star_name, result['lightcurve'], result['period'],
                        result['phased_data'], output=output, **kwargs)

    return result
128732
def deep_run(self):
        """"""
        # Lomb-Scargle period finding.
        self.get_period_LS(self.date, self.mag, self.n_threads, self.min_period)

        # Features based on a phase-folded light curve
        # such as Eta, slope-percentile, etc.
        # Should be called after the getPeriodLS() is called.

        # Created phased a folded light curve.
        # We use period * 2 to take eclipsing binaries into account.
        phase_folded_date = self.date % (self.period * 2.)
        sorted_index = np.argsort(phase_folded_date)

        folded_date = phase_folded_date[sorted_index]
        folded_mag = self.mag[sorted_index]

        # phase Eta
        self.phase_eta = self.get_eta(folded_mag, self.weighted_std)

        # Slope percentile.
        self.slope_per10, self.slope_per90 = \
            self.slope_percentile(folded_date, folded_mag)

        # phase Cusum
        self.phase_cusum = self.get_cusum(folded_mag)
1068503
def mtotal_from_tau0_tau3(tau0, tau3, f_lower,
                          in_seconds=False):
    """
    mtotal = (tau3 / _a3(f_lower)) / (tau0 / _a0(f_lower))
    if not in_seconds:
        # convert back to solar mass units
        mtotal /= lal.MTSUN_SI
    return mtotal
1092398
def satellite(self,stellar_mass,distance_modulus,mc_source_id=1,seed=None,**kwargs):
        
        if seed is not None: np.random.seed(seed)

        isochrone = kwargs.pop('isochrone',self.isochrone)
        kernel    = kwargs.pop('kernel',self.kernel)

        for k,v in kwargs.items():
            if k in kernel.params.keys(): setattr(kernel,k,v)

        mag_1, mag_2 = isochrone.simulate(stellar_mass, distance_modulus)
        lon, lat     = kernel.simulate(len(mag_1))
 
        logger.info("Simulating %i satellite stars..."%len(mag_1))
        pix = ang2pix(self.config['coords']['nside_pixel'], lon, lat)

        # There is probably a better way to do this step without creating the full HEALPix map
        mask = -1. * np.ones(hp.nside2npix(self.config['coords']['nside_pixel']))
        mask[self.roi.pixels] = self.mask.mask_1.mask_roi_sparse
        mag_lim_1 = mask[pix]
        mask = -1. * np.ones(hp.nside2npix(self.config['coords']['nside_pixel']))
        mask[self.roi.pixels] = self.mask.mask_2.mask_roi_sparse
        mag_lim_2 = mask[pix]

        mag_err_1 = self.photo_err_1(mag_lim_1 - mag_1)
        mag_err_2 = self.photo_err_2(mag_lim_2 - mag_2)

        # Randomize magnitudes by their errors
        mag_obs_1 = mag_1+np.random.normal(size=len(mag_1))*mag_err_1
        mag_obs_2 = mag_2+np.random.normal(size=len(mag_2))*mag_err_2
        #mag_obs_1 = mag_1
        #mag_obs_2 = mag_2

        #select = np.logical_and(mag_obs_1 < mag_lim_1, mag_obs_2 < mag_lim_2)
        select = (mag_lim_1>mag_obs_1)&(mag_lim_2>mag_obs_2)

        # Make sure objects lie within the original cmd (should also be done later...)
        #select &= (ugali.utils.binning.take2D(self.mask.solid_angle_cmd, mag_obs_1 - mag_obs_2, mag_obs_1,self.roi.bins_color, self.roi.bins_mag) > 0)

        #return mag_1_obs[cut], mag_2_obs[cut], lon[cut], lat[cut]
        logger.info("Clipping %i simulated satellite stars..."%(~select).sum())
        mc_source_id = mc_source_id * np.ones(len(mag_1))
        
        hdu = ugali.observation.catalog.makeHDU(self.config,mag_obs_1[select],mag_err_1[select],
                                                mag_obs_2[select],mag_err_2[select], 
                                                lon[select],lat[select],mc_source_id[select])
        catalog = ugali.observation.catalog.Catalog(self.config, data=hdu.data)
        return catalog
771400
def _minimize_c(c, z=0, a_tilde=1, b_tilde=-1,
                Ascaling=900, omega_M_0=0.25, omega_lambda_0=0.75):
    

    # Fn 1 (LHS of Eqn 18)

    Y1 = np.log(2) - 0.5
    Yc = np.log(1+c) - c/(1+c)
    f1 = Y1/Yc

    # Fn 2 (RHS of Eqn 18)

    # Eqn 14 - Define the mean inner density
    rho_2 = 200 * c**3 * Y1 / Yc

    # Eqn 17 rearranged to solve for Formation Redshift
    # essentially when universe had rho_2 density
    zf = (((1 + z)**3 + omega_lambda_0/omega_M_0) *
          (rho_2/Ascaling) - omega_lambda_0/omega_M_0)**(1/3) - 1

    # RHS of Eqn 19
    f2 = ((1 + zf - z)**a_tilde) * np.exp((zf - z) * b_tilde)

    # LHS - RHS should be zero for the correct concentration
    return(f1-f2)
963895
def get_SZ_orient(self):
        

        tm_outdated = self._tm_signature != (self.radius, self.radius_type, 
            self.wavelength, self.m, self.axis_ratio, self.shape, self.ddelt, 
            self.ndgs)
        scatter_outdated = self._scatter_signature != (self.thet0, self.thet, 
            self.phi0, self.phi, self.alpha, self.beta, self.orient)

        orient_outdated = self._orient_signature != \
            (self.orient, self.or_pdf, self.n_alpha, self.n_beta)
        if orient_outdated:
            self._init_orient()
        
        outdated = tm_outdated or scatter_outdated or orient_outdated

        if outdated:
            (self._S_orient, self._Z_orient) = self.orient(self)
            self._set_scatter_signature()

        return (self._S_orient, self._Z_orient)
204775
def get_redshift(self, dist):
        
        dist, input_is_array = ensurearray(dist)
        try:
            zs = self.nearby_d2z(dist)
        except TypeError:
            # interpolant hasn't been setup yet
            self.setup_interpolant()
            zs = self.nearby_d2z(dist)
        # if any points had red shifts beyond the nearby, will have nans;
        # replace using the faraway interpolation
        replacemask = numpy.isnan(zs)
        if replacemask.any():
            zs[replacemask] = self.faraway_d2z(dist[replacemask])
            replacemask = numpy.isnan(zs)
        # if we still have nans, means that some distances are beyond our
        # furthest default; fall back to using astropy
        if replacemask.any():
            # well... check that the distance is positive and finite first
            if not (dist > 0.).all() and numpy.isfinite(dist).all():
                raise ValueError("distance must be finite and > 0")
            zs[replacemask] = _redshift(dist[replacemask],
                                        cosmology=self.cosmology)
        return formatreturn(zs, input_is_array)
597260
def countok(self):
        
        ok = np.ones(len(self.stars)).astype(bool)
        for name in self.constraints:
            c = self.constraints[name]
            if c.name not in self.selectfrac_skip:
                ok &= c.ok
        return ok
875693
def parameters(self):
        

        try:
            return self._parameters

        except AttributeError:
            None

        parameters = []
        parameters.extend(self.grid_points.dtype.names)

        model_configuration = self._configuration.get("model", {})

        # Continuum
        # Note: continuum cannot just be 'True', it should be a dictionary
        n = len(parameters)
        continuum_config = model_configuration.get("continuum", False)
        if continuum_config is not False:
            if not isinstance(continuum_config, dict):
                raise TypeError("incorrect continuum type")
            for name in set(self.channel_names).intersection(continuum_config):
                degree = continuum_config[name]
                if degree is False:
                    continue

                # Note: N_coefficients = degree + 1
                # Example: degree = 3: a^3 + b^2 + c^1 + d 
                parameters.extend(["continuum_{0}_{1}".format(name, n) \
                    for n in range(degree + 1)])

        # Redshift
        parameters.extend(channel_parameters("z", self.channel_names,
            model_configuration.get("redshift", False)))

        # Spectral resolution
        resolution_parameters = channel_parameters(
            "resolution", self.channel_names,
            model_configuration.get("spectral_resolution", False))
        parameters.extend(resolution_parameters)

        # Underestimated variance
        parameters.extend(channel_parameters(
            "ln_f", self.channel_names,
            model_configuration.get("underestimated_variance", False)))

        # Outlier pixels
        if model_configuration.get("outlier_pixels", False):
            parameters.extend(["Po", "Vo"])

        assert len(parameters) == len(set(parameters)), "Duplicate parameters!"

        self._resolution_parameters = resolution_parameters
        self._parameters = tuple(parameters)

        return self._parameters
790368
final boolean hasNoMajorSolarTerm(long utcDays) {
        double jd0 = JulianDay.ofEphemerisTime(this.midnight(utcDays)).getValue();
        int index0 = (2 + (int) Math.floor(SolarTerm.solarLongitude(jd0) / 30)) % 12;
        double jd1 = JulianDay.ofEphemerisTime(this.midnight(this.newMoonOnOrAfter(utcDays + 1))).getValue();
        int index1 = (2 + (int) Math.floor(SolarTerm.solarLongitude(jd1) / 30)) % 12;
        return (index0 == index1);
    }
499666
def skip_coords(c):
    """"""
    if c == "{{coord|LAT|LONG|display=inline,title}}": # Unpopulated coord template
        return True
    if c.find("globe:") >= 0 and c.find("globe:earth") == -1: # Moon, venus, etc.
        return True
    return False
569741
def sample_lonlat(self, n):
        
        # From http://en.wikipedia.org/wiki/Ellipse#General_parametric_form
        # However, Martin et al. (2009) use PA theta "from North to East"
        # Definition of phi (position angle) is offset by pi/4
        # Definition of t (eccentric anamoly) remains the same (x,y-frame usual)
        # In the end, everything is trouble because we use glon, glat...

        radius = self.sample_radius(n)
        a = radius; b = self.jacobian * radius

        t = 2. * np.pi * np.random.rand(n)
        cost,sint = np.cos(t),np.sin(t)
        phi = np.pi/2. - np.deg2rad(self.theta)
        cosphi,sinphi = np.cos(phi),np.sin(phi)
        x = a*cost*cosphi - b*sint*sinphi
        y = a*cost*sinphi + b*sint*cosphi
        
        if self.projector  is None:
            logger.debug("Creating AITOFF projector for sampling")
            projector = Projector(self.lon,self.lat,'ait')
        else:
            projector = self.projector
        lon, lat = projector.imageToSphere(x, y)
        return lon, lat
231908
def from_sky(cls, magnitudelimit=None):
        '''
        
        '''


        # define a query for cone search surrounding this center

        criteria = {}
        if magnitudelimit is not None:
            criteria[cls.defaultfilter + 'mag'] = '<{}'.format(magnitudelimit)

        v = Vizier(columns=cls.columns,
                   column_filters=criteria)
        v.ROW_LIMIT = -1

        # run the query
        print('querying Vizier for {}, for {}<{}'.format(cls.name, cls.defaultfilter, magnitudelimit))

        table = v.query_constraints(catalog=cls.catalog, **criteria)[0]

        # store the search parameters in this object
        c = cls(cls.standardize_table(table))
        c.standardized.meta['catalog'] = cls.catalog
        c.standardized.meta['criteria'] = criteria
        c.standardized.meta['magnitudelimit'] = magnitudelimit or c.magnitudelimit
        #c.magnitudelimit = magnitudelimit or c.magnitudelimit
        return c
1146835
def search(self, **params):
        

        self.logger.debug("search params=%s" % (str(params)))
        ra, dec = params['ra'], params['dec']
        if not (':' in ra):
            # Assume RA and DEC are in degrees
            ra_deg = float(ra)
            dec_deg = float(dec)
        else:
            # Assume RA and DEC are in standard string notation
            ra_deg = wcs.hmsStrToDeg(ra)
            dec_deg = wcs.dmsStrToDeg(dec)

        # Convert to degrees for search radius
        radius_deg = float(params['r']) / 60.0
        # radius_deg = float(params['r'])

        # Note requires astropy 0.3.x+
        c = coordinates.SkyCoord(ra_deg * units.degree,
                                 dec_deg * units.degree,
                                 frame='icrs')
        self.logger.info("Querying catalog: %s" % (self.full_name))
        time_start = time.time()
        with warnings.catch_warnings():  # Ignore VO warnings
            warnings.simplefilter('ignore')
            results = conesearch.conesearch(
                c, radius_deg * units.degree, catalog_db=self.full_name,
                verbose=False)
        time_elapsed = time.time() - time_start
        numsources = results.array.size
        self.logger.info("Found %d sources in %.2f sec" % (
            numsources, time_elapsed))

        # Scan the returned fields to find ones we need to extract
        # particulars from (ra, dec, id, magnitude)
        mags = []
        ext = {}
        fields = results.array.dtype.names
        for name in fields:
            ucd = results.get_field_by_id(name).ucd
            ucd = str(ucd).lower()
            if ucd == 'id_main':
                ext['id'] = name
            elif ucd == 'pos_eq_ra_main':
                ext['ra'] = name
            elif ucd == 'pos_eq_dec_main':
                ext['dec'] = name
            if ('phot_' in ucd) or ('phot.' in ucd):
                mags.append(name)
        self.logger.debug("possible magnitude fields: %s" % str(mags))
        if len(mags) > 0:
            magfield = mags[0]
        else:
            magfield = None

        # prepare the result list
        starlist = []
        arr = results.array
        for i in range(numsources):
            source = dict(zip(fields, arr[i]))
            starlist.append(self.toStar(source, ext, magfield))

        # metadata about the list
        columns = [('Name', 'name'),
                   ('RA', 'ra'),
                   ('DEC', 'dec'),
                   ('Mag', 'mag'),
                   ('Preference', 'preference'),
                   ('Priority', 'priority'),
                   ('Description', 'description'),
                   ]
        # Append extra columns returned by search to table header
        # TODO: what if not all sources have same record structure?
        # is this possible with VO?
        cols = list(fields)
        cols.remove(ext['ra'])
        cols.remove(ext['dec'])
        cols.remove(ext['id'])
        columns.extend(zip(cols, cols))

        # which column is the likely one to color source circles
        colorCode = 'Mag'

        info = Bunch.Bunch(columns=columns, color=colorCode)
        return starlist, info