orbital angular momentum.

    e.g Ge etal2010
    
    Parameters
    ----------
    m1, m2 : float
        Masses of both stars in Msun.
    A : float 
        Separation in Rsun.
    e : float
        Eccentricity
342342
def am_orb(m1,m2,a,e):
    ''' 
    
        
    '''

    a_cm  = a * rsun_cm
    m1_g = m1 * msun_g
    m2_g = m2 * msun_g

    J_orb=np.sqrt(grav_const*a_cm*(old_div((m1_g**2*m2_g**2),(m1_g+m2_g))))*(1-e**2)
    return J_orb
458987
def spherical(cls, mag, theta, phi=0):
        ''''''
        return cls(
            mag * math.sin(phi) * math.cos(theta),  # X
            mag * math.sin(phi) * math.sin(theta),  # Y
            mag * math.cos(phi)  # Z
        )
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
754779
def asd(M1, M2):
    

    from scipy.fftpack import fft2
    spectra1 = np.abs(fft2(M1))
    spectra2 = np.abs(fft2(M2))

    return np.linalg.norm(spectra2 - spectra1)
965539
def spher_harms(l, m, inclination):
    

    # FIXME: we are using spin -2 weighted spherical harmonics for now,
    # when possible switch to spheroidal harmonics.
    Y_lm = lal.SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m).real
    Y_lminusm = lal.SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m).real
    Y_plus = Y_lm + (-1)**l * Y_lminusm
    Y_cross = Y_lm - (-1)**l * Y_lminusm

    return Y_plus, Y_cross
376029
def _tle_to_keplerian_mean(cls, coord, center):
        
        i, Ω, e, ω, M, n = coord
        a = (center.µ / n ** 2) ** (1 / 3)

        return np.array([a, e, i, Ω, ω, M], dtype=float)
117601
def mass_2d(self,R,Rs,rho0,r_trunc):

        

        x = R / Rs
        tau = r_trunc / Rs
        gx = self._g(x,tau)
        m_2d = 4 * rho0 * Rs * R ** 2 * gx / x ** 2 * np.pi
        return m_2d
465470
def spin2x_from_mass1_mass2_xi2_phi_a_phi_s(mass1, mass2, xi2, phi_a, phi_s):
    
    chi_perp = chi_perp_from_mass1_mass2_xi2(mass1, mass2, xi2)
    phi2 = phi2_from_phi_a_phi_s(phi_a, phi_s)
    return chi_perp * numpy.cos(phi2)
596822
def generate(self, mA=1, age=9.6, feh=0.0, n=1e5, ichrone='mist',
                 orbpop=None, bands=None, **kwargs):
        
        ichrone = get_ichrone(ichrone, bands=bands)

        n = int(n)
        #star with m1 orbits (m2+m3).  So mA (most massive)
        # will correspond to either m1 or m2.
        m1, m2, m3 = self.multmass_fn(mA, f_binary=self.f_binary,
                                      f_triple=self.f_triple,
                                      qmin=self.qmin, minmass=self.minmass,
                                      n=n)
        #reset n if need be
        n = len(m1)

        feh = np.ascontiguousarray(np.atleast_1d(feh))
        age = np.ascontiguousarray(age)

        #generate stellar properties
        primary = ichrone(np.ascontiguousarray(m1), age, feh,
                          bands=bands)
        secondary = ichrone(np.ascontiguousarray(m2),age,feh,
                            bands=bands)
        tertiary = ichrone(np.ascontiguousarray(m3),age,feh,
                           bands=bands)

        #clean up columns that become nan when called with mass=0
        # Remember, we want mass=0 and mags=inf when something doesn't exist
        no_secondary = (m2==0)
        no_tertiary = (m3==0)
        for c in secondary.columns: #
            if re.search('_mag',c):
                secondary[c][no_secondary] = np.inf
                tertiary[c][no_tertiary] = np.inf
        secondary['mass'][no_secondary] = 0
        tertiary['mass'][no_tertiary] = 0

        if kwargs['period_short'] is None:
            if kwargs['period_long'] is None:
                period_1 = self.period_long_fn(n)
                period_2 = self.period_short_fn(n)
                kwargs['period_short'] = np.minimum(period_1, period_2)
                kwargs['period_long'] = np.maximum(period_1, period_2)
            else:
                kwargs['period_short'] = self.period_short_fn(n)

                #correct any short periods that are longer than period_long
                bad = kwargs['period_short'] > kwargs['period_long']
                n_bad = bad.sum()
                good_inds = np.where(~bad)[0]
                inds = np.random.randint(len(good_inds),size=n_bad)
                kwargs['period_short'][bad] = \
                    kwargs['period_short'][good_inds[inds]]
        else:
            if kwargs['period_long'] is None:
                kwargs['period_long'] = self.period_long_fn(n)

                #correct any long periods that are shorter than period_short
                bad = kwargs['period_long'] < kwargs['period_short']
                n_bad = bad.sum()
                good_inds = np.where(~bad)[0]
                inds = np.random.randint(len(good_inds),size=n_bad)
                kwargs['period_long'][bad] = \
                    kwargs['period_long'][good_inds[inds]]

        if 'ecc_short' not in kwargs:
            kwargs['ecc_short'] = self.ecc_fn(n, kwargs['period_short'])
        if 'ecc_long' not in kwargs:
            kwargs['ecc_long'] = self.ecc_fn(n, kwargs['period_long'])

        TriplePopulation.__init__(self, primary=primary,
                                  secondary=secondary, tertiary=tertiary,
                                  orbpop=orbpop, **kwargs)

        return self
473475
def expected_magz(RAW_IMU, ATTITUDE, inclination, declination):
    ''''''
    v = expected_mag(RAW_IMU, ATTITUDE, inclination, declination)
    return v.z
1344356
def _get_nsymop(self):
        """"""
        if self.centrosymmetric:
            return 2 * len(self._rotations) * len(self._subtrans)
        else:
            return len(self._rotations) * len(self._subtrans)
117059
def angular_separation(lonp1, latp1, lonp2, latp2):
    

    xp1, yp1, zp1 = lonlat2xyz(lonp1, latp1)
    xp2, yp2, zp2 = lonlat2xyz(lonp2, latp2)

    ## dot products to obtain angles

    angles = np.arccos((xp1 * xp2 + yp1 * yp2 + zp1 * zp2))

    ## As this is a unit sphere, angle = length

    return angles
1372706
private List<List<Integer>> calcOrbits() {
		int n = subunits.getSubunitCount();
		int fold = rotationGroup.getRotation(0).getFold();

		List<List<Integer>> orbits = new ArrayList<List<Integer>>();
		boolean[] used = new boolean[n];
		Arrays.fill(used, false);

		for (int i = 0; i < n; i++) {
			if (! used[i]) {
				// determine the equivalent subunits
				List<Integer> orbit = new ArrayList<Integer>(fold);
				for (int j = 0; j < fold; j++) {
					List<Integer> permutation = rotationGroup.getRotation(j).getPermutation();
					orbit.add(permutation.get(i));
					used[permutation.get(i)] = true;
				}
				orbits.add(deconvolute(orbit));
			}
		}
		return orbits;
	}
247917
def a_over_Rs(P,R2,M2,M1=1,R1=1,planet=True):
    
    if planet:
        M2 *= REARTH/RSUN
        R2 *= MEARTH/MSUN
    return semimajor(P,M1+M2)*AU/(R1*RSUN)
1128232
def Kerr_factor(final_mass, distance):
    

    # Convert solar masses to meters
    mass = final_mass * lal.MSUN_SI * lal.G_SI / lal.C_SI**2
    # Convert Mpc to meters
    dist = distance * 1e6 * lal.PC_SI

    return mass / dist
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
753781
def eclipse_pars(P,M1,M2,R1,R2,ecc=0,inc=90,w=0,sec=False):
    """"""
    a = semimajor(P,M1+M2)
    if sec:
        b = a*AU*np.cos(inc*np.pi/180)/(R1*RSUN) * (1-ecc**2)/(1 - ecc*np.sin(w*np.pi/180))
        #aR = a*AU/(R2*RSUN) #I feel like this was to correct a bug, but this should not be.
        #p0 = R1/R2 #why this also?
    else:
        b = a*AU*np.cos(inc*np.pi/180)/(R1*RSUN) * (1-ecc**2)/(1 + ecc*np.sin(w*np.pi/180))
        #aR = a*AU/(R1*RSUN)
        #p0 = R2/R1
    p0 = R2/R1
    aR = a*AU/(R1*RSUN)
    return p0,b,aR