\psi(m) = -\int_m^\infty d m^2 \rho(m^2)
326689
def _psi(self,m):
        """"""
        return 2.*self.a2*(1./(1.+m/self.a)+numpy.log(m/(m+self.a)))
813744
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
        
        self._psi_computations(Z, mu, S)
        tmp = (self.inv_lengthscale2 * self._psi2[:, :, :, None]) / self._psi2_denom
        target_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * self._psi2_mudist).sum(1).sum(1)
        target_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*self._psi2_mudist_sq - 1)).sum(1).sum(1)
1183818
def va(self):
        
        return np.sqrt(self.mu * (2 / (self.ra) - 1 / self.kep.a))
1149861
def _JRStaeckelIntegrandSquared(u,E,Lz,I3U,delta,u0,sinh2u0,v0,sin2v0,
                                potu0v0,pot):
    #potu0v0= potentialStaeckel(u0,v0,pot,delta)
    """"""
    sinh2u= nu.sinh(u)**2.
    dU= (sinh2u+sin2v0)*potentialStaeckel(u,v0,
                                          pot,delta)\
        -(sinh2u0+sin2v0)*potu0v0
    return E*sinh2u-I3U-dU-Lz**2./2./delta**2./sinh2u
1011
function() {

      // acoth(c) = log((c+1) / (c-1)) / 2

      var a = this['re'];
      var b = this['im'];

      if (a === 0 && b === 0) {
        return new Complex(0, Math.PI / 2);
      }

      var d = a * a + b * b;
      return (d !== 0)
              ? new Complex(
                      a / d,
                      -b / d).atanh()
              : new Complex(
                      (a !== 0) ? a / 0 : 0,
                      (b !== 0) ? -b / 0 : 0).atanh();
    }
1202900
def velocities_to_moduli(rho, v_phi, v_s):
    
    return v_phi * v_phi * rho, v_s * v_s * rho
9314
def f(self, y):
        
        d = self.d
        mpsi = self.psi
        z = d * y.copy()
        for i in range(len(mpsi)):
            a, b, c = mpsi[i]
            z += a * np.tanh(b * (y + c))
        return z
881816
def _psi2computations(dL_dpsi2, variance, Z, mu, S):
    
    # here are the "statistics" for psi1 and psi2
    # Produced intermediate results:
    # _psi2_dvariance      Q
    # _psi2_dZ             MxQ
    # _psi2_dmu            NxQ
    # _psi2_dS             NxQ

    variance2 = np.square(variance)
    common_sum = np.dot(mu,(variance*Z).T)
    if len(dL_dpsi2.shape)==2:
        Z_expect = (np.dot(dL_dpsi2,Z)*Z).sum(axis=0)
        dL_dpsi2T = dL_dpsi2+dL_dpsi2.T
        common_expect = np.dot(common_sum,np.dot(dL_dpsi2T,Z))
        Z2_expect = np.inner(common_sum,dL_dpsi2T)
        Z1_expect = np.dot(dL_dpsi2T,Z)
    
        dL_dvar = 2.*S.sum(axis=0)*variance*Z_expect+(common_expect*mu).sum(axis=0)
    
        dL_dmu = common_expect*variance
    
        dL_dS = np.empty(S.shape)
        dL_dS[:] = Z_expect*variance2
    
        dL_dZ = variance2*S.sum(axis=0)*Z1_expect+np.dot(Z2_expect.T,variance*mu)
    else:
        N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1]
        dL_dpsi2_ = dL_dpsi2.sum(axis=0)
        Z_expect = (np.dot(dL_dpsi2.reshape(N*M,M),Z).reshape(N,M,Q)*Z[None,:,:]).sum(axis=1)
        dL_dpsi2T = dL_dpsi2_+dL_dpsi2_.T
        dL_dpsi2T_ = dL_dpsi2+np.swapaxes(dL_dpsi2, 1, 2)
        common_expect = np.dot(common_sum,np.dot(dL_dpsi2T,Z))
        common_expect_ = (common_sum[:,:,None]*np.dot(dL_dpsi2T_.reshape(N*M,M),Z).reshape(N,M,Q)).sum(axis=1)
        Z2_expect = (common_sum[:,:,None]*dL_dpsi2T_).sum(axis=1)
        Z1_expect = np.dot(dL_dpsi2T_.reshape(N*M,M),Z).reshape(N,M,Q)
    
        dL_dvar = 2.*variance*(S*Z_expect).sum(axis=0)+(common_expect_*mu).sum(axis=0)
    
        dL_dmu = common_expect_*variance
    
        dL_dS = np.empty(S.shape)
        dL_dS[:] = variance2* Z_expect
    
        dL_dZ = variance2*(S[:,None,:]*Z1_expect).sum(axis=0)+np.dot(Z2_expect.T,variance*mu)

    return dL_dvar, dL_dmu, dL_dS, dL_dZ
1040833
def _get_rho(self, v):
        
        v_mol = vol_uc2mol(v, self.z)  # in m^3
        rho = self.mass / v_mol * 1.e-6  # in g/cm^3
        return rho
445519
def mui(x, y):
    """"""
    assert len(x) == len(y)
    l = len(x)
    p_x = Counter(x)
    p_y = Counter(y)
    p_xy = Counter(zip(x, y))
    return sum(p_xy[(x, y)] * math.log((p_xy[(x, y)] * l) / float(p_x[x] * p_y[y]), 2) / l
               for x, y in p_xy)
1138276
def derivatives(self, x, y, sigma0, Rs, e1, e2, center_x=0, center_y=0):
        
        phi_G, q = param_util.ellipticity2phi_q(e1, e2)
        x_shift = x - center_x
        y_shift = y - center_y
        cos_phi = np.cos(phi_G)
        sin_phi = np.sin(phi_G)
        e = abs(1 - q)
        x_ = (cos_phi*x_shift+sin_phi*y_shift)*np.sqrt(1 - e)
        y_ = (-sin_phi*x_shift+cos_phi*y_shift)*np.sqrt(1 + e)

        f_x_prim, f_y_prim = self.spherical.derivatives(x_, y_, sigma0, Rs)
        f_x_prim *= np.sqrt(1 - e)
        f_y_prim *= np.sqrt(1 + e)
        f_x = cos_phi*f_x_prim-sin_phi*f_y_prim
        f_y = sin_phi*f_x_prim+cos_phi*f_y_prim
        return f_x, f_y
1059980
def hessian(self, x, y, kappa_ext, ra_0=0, dec_0=0):
        
        gamma1 = 0
        gamma2 = 0
        kappa = kappa_ext
        f_xx = kappa + gamma1
        f_yy = kappa - gamma1
        f_xy = gamma2
        return f_xx, f_yy, f_xy
650118
def _differential_to_integral(hpx_map):
        
        ebins = hpx_map.hpx.ebins
        ratio = ebins[1:] / ebins[0:-1]
        half_log_ratio = np.log(ratio) / 2.
        int_map = ((hpx_map.data[0:-1].T * ebins[0:-1]) +
                   (hpx_map.data[1:].T * ebins[1:])) * half_log_ratio
        return HpxMap(int_map.T, hpx_map.hpx)
121512
def rho_hv(scatterer):
    
    Z = scatterer.get_Z()
    a = (Z[2,2] + Z[3,3])**2 + (Z[3,2] - Z[2,3])**2
    b = (Z[0,0] - Z[0,1] - Z[1,0] + Z[1,1])
    c = (Z[0,0] + Z[0,1] + Z[1,0] + Z[1,1])
    return np.sqrt(a / (b*c))
868396
def _psi_n(x, n, b):
    
    return 2**(b-1) / gamma(b) * (-1)**n * \
    np.exp(gammaln(n+b) -
           gammaln(n+1) +
           np.log(2*n+b) -
           0.5 * np.log(2*np.pi*x**3) -
           (2*n+b)**2 / (8.*x))
1300085
def mjd2gmst(mjd):
    

    tu = (mjd - MJD0) / (100*DPY)

    st = math.fmod(mjd, 1.0) * D2PI + (24110.54841 + (8640184.812866 + (0.093104 - 6.2e-6 * tu) * tu) * tu) * DS2R

    w = math.fmod(st, D2PI)
    if w >= 0.0:
        return w
    else:
        return w + D2PI
74100
function hessenbergQRStep (M) {
  const n = M.dimension;
  const me = M.elements;
  const csr = [0, 0, 0];
  const cs = [];
  for (let k = 0; k < n - 1; ++k) {
    const row1 = k;
    const row2 = k + 1;
    const colK = k * n;
    const a = me[row1 + colK];
    const b = me[row2 + colK];
    rotg(a, b, csr);
    const c = csr[0];
    const s = csr[1];
    const r = csr[2];
    cs[k] = [c, s];
    for (let x = k; x < n; ++x) {
      const colX = x * n;
      const tmp1 = me[row1 + colX];
      const tmp2 = me[row2 + colX];
      me[row1 + colX] = c * tmp1 + s * tmp2;
      me[row2 + colX] = -s * tmp1 + c * tmp2;
    }
    me[row1 + colK] = r;
    me[row2 + colK] = 0;
  }
  for (let k = 0; k < n - 1; ++k) {
    const col1 = k * n;
    const col2 = (k + 1) * n;
    const [c, s] = cs[k];
    for (let x = 0; x <= k + 1; ++x) {
      const rowX = x;
      const tmp1 = me[rowX + col1];
      const tmp2 = me[rowX + col2];
      me[rowX + col1] = tmp1 * c + tmp2 * s;
      me[rowX + col2] = -tmp1 * s + tmp2 * c;
    }
  }
  return M;
}