angular diameter from deflector to source
        :param H_0: Hubble parameter [km/s/Mpc]
        :param Om0: normalized matter density at present time
        :return: float [Mpc]
966068
def D_ds(self, H_0, Om0, Ode0=None):
        
        lensCosmo = self._get_cosom(H_0, Om0, Ode0)
        return lensCosmo.D_ds
864005
def _sigmainf(N, h, m, dW, Km0, Pm0):
    """"""
    M = m*(m-1)//2
    Im = broadcast_to(np.eye(m), (N, m, m))
    IM = broadcast_to(np.eye(M), (N, M, M))
    Ims0 = np.eye(m**2)
    factor1 = broadcast_to((2.0/h)*np.dot(Km0, Ims0 - Pm0), (N, M, m**2))
    factor2 = _kp2(Im, _dot(dW, _t(dW)))
    factor3 = broadcast_to(np.dot(Ims0 - Pm0, Km0.T), (N, m**2, M))
    return 2*IM + _dot(_dot(factor1, factor2), factor3)
586597
def source_to_neutrino_direction(azimuth, zenith, radian=True):
    
    azimuth = np.atleast_1d(azimuth).copy()
    zenith = np.atleast_1d(zenith).copy()
    if not radian:
        azimuth *= np.pi / 180
        zenith *= np.pi / 180
    phi = (azimuth - np.pi) % (2 * np.pi)
    theta = np.pi - zenith
    if not radian:
        phi *= 180 / np.pi
        theta *= 180 / np.pi
    return phi, theta
604626
def moreland_adjusthue (msh, m_unsat):
    
    if msh[M] >= m_unsat:
        return msh[H] # "Best we can do"

    hspin = (msh[S] * np.sqrt (m_unsat**2 - msh[M]**2) /
             (msh[M] * np.sin (msh[S])))

    if msh[H] > -np.pi / 3: # "Spin away from purple"
        return msh[H] + hspin
    return msh[H] - hspin
754440
def endpoint(distance, heading)
      d = (distance / 1000.0) / EARTH_RADIUS_KM;
      heading = deg2rad(heading);

      oX = lng * Math::PI / 180;
      oY = lat * Math::PI / 180;

      y = Math.asin(Math.sin(oY) * Math.cos(d) + Math.cos(oY) * Math.sin(d) * Math.cos(heading));
      x = oX + Math.atan2(Math.sin(heading) * Math.sin(d) * Math.cos(oY), Math.cos(d) - Math.sin(oY) * Math.sin(y));

      y = y * 180 / Math::PI;
      x = x * 180 / Math::PI;

      self.class.new(:lat => y, :lng => x)
    end
427935
def segment_centre_of_mass(seg):
    ''''''
    h = mm.segment_length(seg)
    r0 = seg[0][COLS.R]
    r1 = seg[1][COLS.R]
    num = r0 * r0 + 2 * r0 * r1 + 3 * r1 * r1
    denom = 4 * (r0 * r0 + r0 * r1 + r1 * r1)
    centre_of_mass_z_loc = num / denom
    return seg[0][COLS.XYZ] + (centre_of_mass_z_loc / h) * (seg[1][COLS.XYZ] - seg[0][COLS.XYZ])
1207013
def _rho02alpha(self, rho0, Rs):

        

        theta_Rs = rho0 * (4 * Rs ** 2 * (1 + np.log(1. / 2.)))
        return theta_Rs
1007285
private static double sinAlt(
        double mjd0, // earliest moment of calendar date (usually midnight)
        double hour,
        double longitudeRad,
        double cosLatitude,
        double sinLatitude,
        double geodeticAngle,
        double refraction,
        double deltaT
    ) {

        double mjd = mjd0 + hour / 24.0; // UT
        double jct = toJulianCenturies(mjd + (deltaT / 86400));
        double[] data = MoonPosition.calculateMeeus(jct);
        double nutationCorr = data[0] * Math.cos(Math.toRadians(data[1])); // for apparent sidereal time
        double tau = AstroUtils.gmst(mjd) + Math.toRadians(nutationCorr) + longitudeRad - Math.toRadians(data[2]);
        double decl = Math.toRadians(data[3]);

        // transformation to horizontal coordinate system
        double sinAltitude = sinLatitude * Math.sin(decl) + cosLatitude * Math.cos(decl) * Math.cos(tau);

        // about impact of horizontal parallax on moon diameter see also Meeus (chapter 15)
        double correction = 0.7275 * getHorizontalParallax(data[4]) - refraction - geodeticAngle;

        // we search for the roots of this function
        return sinAltitude - Math.sin(Math.toRadians(correction));

    }
1252790
def get_floating_spinning_factors(self):
        
        data = []
        for src in self.get_sources():
            if hasattr(src, 'hypocenter_distribution'):
                data.append(
                    (len(src.hypocenter_distribution.data),
                     len(src.nodal_plane_distribution.data)))
        if not data:
            return numpy.array([1, 1])
        return numpy.array(data).mean(axis=0)
1100434
def FibreMode(gridSize,modeDiameter):
    
    rmode=modeDiameter/2
    return np.exp(-(RadiusGrid(gridSize)/rmode)**2)/(np.sqrt(np.pi/2)*rmode)
223608
def bures_distance(rho0: Density, rho1: Density) -> float:
    
    fid = fidelity(rho0, rho1)
    op0 = asarray(rho0.asoperator())
    op1 = asarray(rho1.asoperator())
    tr0 = np.trace(op0)
    tr1 = np.trace(op1)

    return np.sqrt(tr0 + tr1 - 2.*np.sqrt(fid))
1153147
def calcPeriod(self):
        

        return eq.KeplersThirdLaw(self.a, self.star.M).P
875045
def D_d(self, H_0, Om0, Ode0=None):
        
        lensCosmo = self._get_cosom(H_0, Om0, Ode0)
        return lensCosmo.D_d
1202900
def velocities_to_moduli(rho, v_phi, v_s):
    
    return v_phi * v_phi * rho, v_s * v_s * rho
204153
def _keplerian_to_keplerian_mean(cls, coord, center):
        

        a, e, i, Ω, ω, ν = coord
        if e < 1:
            # Elliptic case
            cos_E = (e + cos(ν)) / (1 + e * cos(ν))
            sin_E = (sin(ν) * sqrt(1 - e ** 2)) / (1 + e * cos(ν))
            E = arctan2(sin_E, cos_E) % (2 * np.pi)
            M = E - e * sin(E)  # Mean anomaly
        else:
            # Hyperbolic case
            H = arccosh((e + cos(ν)) / (1 + e * cos(ν)))
            M = e * sinh(H) - H

        return np.array([a, e, i, Ω, ω, M], dtype=float)
650277
def _reduced2physical_deflection(self, alpha_reduced, idex_lens):
        
        factor = self._reduced2physical_factor[idex_lens]
        #factor = self._cosmo_bkg.D_xy(0, z_source) / self._cosmo_bkg.D_xy(z_lens, z_source)
        return alpha_reduced * factor
887730
def mass_3d(self, R, Rs, rho0):
        
        Rs = float(Rs)
        m_3d = 4. * np.pi * rho0 * Rs**3 *(np.log((Rs + R)/Rs) - R/(Rs + R))
        return m_3d