Get the radiance from the Tb using the simple non-linear regression
        method. SI units of course!
339990
def tb2radiance(self, tb_, **kwargs):
        
        # L = C1 * νc**3 / (exp (C2 νc / [αTb + β]) − 1)
        #
        # C1 = 2 * h * c**2 and C2 = hc/k
        #
        lut = kwargs.get('lut', None)
        normalized = kwargs.get('normalized', True)

        if lut is not None:
            raise NotImplementedError('Using a tb-radiance LUT is not yet supported')
        if not normalized:
            raise NotImplementedError('Deriving the band integrated radiance is not supported')

        c_1 = 2 * H_PLANCK * C_SPEED ** 2
        c_2 = H_PLANCK * C_SPEED / K_BOLTZMANN

        vc_ = SEVIRI[self.bandname][self.platform_name][0]
        # Multiply by 100 to get SI units!
        vc_ *= 100.0
        alpha = SEVIRI[self.bandname][self.platform_name][1]
        beta = SEVIRI[self.bandname][self.platform_name][2]

        radiance = c_1 * vc_ ** 3 / \
            (np.exp(c_2 * vc_ / (alpha * tb_ + beta)) - 1)

        unit = 'W/m^2 sr^-1 (m^-1)^-1'
        scale = 1.0
        return {'radiance': radiance,
                'unit': unit,
                'scale': scale}
1004974
def gamma(x, context=None):
    
    return _apply_function_in_current_context(
        BigFloat,
        mpfr.mpfr_gamma,
        (BigFloat._implicit_convert(x),),
        context,
    )
455277
def _get_reference_band(self, data):
        
        nir = data[..., self.nir_idx].astype("float32")
        swir1 = data[..., self.swir1_idx].astype("float32")
        return (nir - swir1) / (nir + swir1)
443534
function convert(fov, fromDimension, toDimension) {
  return 2 * Math.atan(toDimension * Math.tan(fov / 2) / fromDimension);
}
230247
def get_irradiance_value_for_hoy(self, hoy):
        """"""
        count = int(hoy * self.timestep)
        return self.direct_normal_irradiance[count], \
            self.diffuse_horizontal_irradiance[count]
1262433
def emissive_part_3x(self, tb=True):
        """"""
        try:
            # Emissive part:
            self._e3x = self._rad3x_t11 * (1 - self._r3x)
            # Unsure how much sense it makes to apply the co2 correction term here!?
            # FIXME!
            # self._e3x *= self._rad3x_correction

        except TypeError:
            LOG.warning(
                "Couldn't derive the emissive part \n" +
                "Please derive the relfectance prior to requesting the emissive part")

        if tb:
            return self.radiance2tb(self._e3x)
        else:
            return self._e3x
1024153
def magnification(self, x, y, kwargs, diff=diff):
        
        f_xx, f_xy, f_yx, f_yy = self.hessian(x, y, kwargs, diff=diff)
        det_A = (1 - f_xx) * (1 - f_yy) - f_xy*f_yx
        return 1/det_A
1182559
def det2lab_xds(
        pixels_coord, frame_number,
        starting_frame, starting_angle, oscillation_angle,
        rotation_axis,
        wavelength, wavevector,
        NX, NY, pixelsize_x, pixelsize_y,
        distance_to_detector, x_center, y_center,
        detector_x, detector_y, detector_normal, **kwargs):
    """"""

    array_shape = (1, 3)

    if detector_x.shape == array_shape:
        detector_x = detector_x.T
        detector_y = detector_y.T
        detector_normal = detector_normal.T
    if wavevector.shape == array_shape:
        wavevector = wavevector.T
    if rotation_axis.shape == array_shape:
        rotation_axis = rotation_axis.T
    xmm = (pixels_coord[:, [0]] - x_center) * pixelsize_x
    ymm = (pixels_coord[:, [1]] - y_center) * pixelsize_y
    # find scattering vector of each pixel
    scattering_vector_mm = np.outer(xmm, detector_x) + \
                           np.outer(ymm, detector_y) + \
                           distance_to_detector * np.outer(np.ones(shape=xmm.shape),
                                                           detector_normal)
    scattering_vector_mm = scattering_vector_mm.T
    phi = (frame_number - starting_frame) * oscillation_angle + \
          starting_angle
    # calculating norm for each column
    norms = np.sum(scattering_vector_mm ** 2., axis=0) ** (1. / 2)
    #deviding scattering vector by its own norm
    unit_scattering_vector = scattering_vector_mm / norms
    #subtracting incident beam vector

    h = unit_scattering_vector / wavelength - \
        np.tile(wavevector, (unit_scattering_vector.shape[1], 1)).T
    #rotating
    if phi.size == 1:
        h = np.dot(rotvec2mat(rotation_axis.T, -2 * np.pi * phi / 360), h)
    else:
        for i in range(phi.size):
            h[:, [i]] = np.dot(
                rotvec2mat(rotation_axis.T, -2 * np.pi * phi[i] / 360), h[:, [i]])

    return h, scattering_vector_mm, unit_scattering_vector
824065
def dImage_wrt_2dVerts_bnd(observed, visible, visibility, barycentric, image_width, image_height, num_verts, f, bnd_bool):
    """"""

    n_channels = np.atleast_3d(observed).shape[2]
    shape = visibility.shape

    # Step 1: get the structure ready, ie the IS and the JS
    IS = np.tile(col(visible), (1, 2*f.shape[1])).ravel()
    JS = col(f[visibility.ravel()[visible]].ravel())
    JS = np.hstack((JS*2, JS*2+1)).ravel()

    pxs = np.asarray(visible % shape[1], np.int32)
    pys = np.asarray(np.floor(np.floor(visible) / shape[1]), np.int32)

    if n_channels > 1:
        IS = np.concatenate([IS*n_channels+i for i in range(n_channels)])
        JS = np.concatenate([JS for i in range(n_channels)])

    # Step 2: get the data ready, ie the actual values of the derivatives
    ksize = 1
    bndf = bnd_bool.astype(np.float64)
    nbndf = np.logical_not(bnd_bool).astype(np.float64)
    sobel_normalizer = cv2.Sobel(np.asarray(np.tile(row(np.arange(10)), (10, 1)), np.float64), cv2.CV_64F, dx=1, dy=0, ksize=ksize)[5,5]

    bnd_nan = bndf.reshape((observed.shape[0], observed.shape[1], -1)).copy()
    bnd_nan.ravel()[bnd_nan.ravel()>0] = np.nan
    bnd_nan += 1
    obs_nonbnd = np.atleast_3d(observed) * bnd_nan

    ydiffnb, xdiffnb = nangradients(obs_nonbnd)

    observed = np.atleast_3d(observed)

    if observed.shape[2] > 1:
        ydiffbnd, xdiffbnd, _ = np.gradient(observed)
    else:
        ydiffbnd, xdiffbnd = np.gradient(observed.squeeze())
        ydiffbnd = np.atleast_3d(ydiffbnd)
        xdiffbnd = np.atleast_3d(xdiffbnd)

    # This corrects for a bias imposed boundary differences begin spread over two pixels
    # (by np.gradients or similar) but only counted once (since OpenGL's line
    # drawing spans 1 pixel)
    xdiffbnd *= 2.0
    ydiffbnd *= 2.0

    xdiffnb = -xdiffnb
    ydiffnb = -ydiffnb
    xdiffbnd = -xdiffbnd
    ydiffbnd = -ydiffbnd
    # ydiffnb *= 0
    # xdiffnb *= 0

    if False:
        import matplotlib.pyplot as plt
        plt.figure()
        plt.subplot(121)
        plt.imshow(xdiffnb)
        plt.title('xdiffnb')
        plt.subplot(122)
        plt.imshow(xdiffbnd)
        plt.title('xdiffbnd')
        import pdb; pdb.set_trace()

    idxs = np.isnan(xdiffnb.ravel())
    xdiffnb.ravel()[idxs] = xdiffbnd.ravel()[idxs]

    idxs = np.isnan(ydiffnb.ravel())
    ydiffnb.ravel()[idxs] = ydiffbnd.ravel()[idxs]

    if True: # should be right thing
        xdiff = xdiffnb
        ydiff = ydiffnb
    else:  #should be old way
        xdiff = xdiffbnd
        ydiff = ydiffbnd


    # TODO: NORMALIZER IS WRONG HERE
    # xdiffnb = -cv2.Sobel(obs_nonbnd, cv2.CV_64F, dx=1, dy=0, ksize=ksize) / np.atleast_3d(cv2.Sobel(row(np.arange(obs_nonbnd.shape[1])).astype(np.float64), cv2.CV_64F, dx=1, dy=0, ksize=ksize))
    # ydiffnb = -cv2.Sobel(obs_nonbnd, cv2.CV_64F, dx=0, dy=1, ksize=ksize) / np.atleast_3d(cv2.Sobel(col(np.arange(obs_nonbnd.shape[0])).astype(np.float64), cv2.CV_64F, dx=0, dy=1, ksize=ksize))
    #
    # xdiffnb.ravel()[np.isnan(xdiffnb.ravel())] = 0.
    # ydiffnb.ravel()[np.isnan(ydiffnb.ravel())] = 0.
    # xdiffnb.ravel()[np.isinf(xdiffnb.ravel())] = 0.
    # ydiffnb.ravel()[np.isinf(ydiffnb.ravel())] = 0.

    # xdiffnb = np.atleast_3d(xdiffnb)
    # ydiffnb = np.atleast_3d(ydiffnb)
    #
    # xdiffbnd = -cv2.Sobel(observed, cv2.CV_64F, dx=1, dy=0, ksize=ksize) / sobel_normalizer
    # ydiffbnd = -cv2.Sobel(observed, cv2.CV_64F, dx=0, dy=1, ksize=ksize) / sobel_normalizer
    #
    # xdiff = xdiffnb * np.atleast_3d(nbndf)
    # xdiff.ravel()[np.isnan(xdiff.ravel())] = 0
    # xdiff += xdiffbnd*np.atleast_3d(bndf)
    #
    # ydiff = ydiffnb * np.atleast_3d(nbndf)
    # ydiff.ravel()[np.isnan(ydiff.ravel())] = 0
    # ydiff += ydiffbnd*np.atleast_3d(bndf)

    #import pdb; pdb.set_trace()

    #xdiff = xdiffnb
    #ydiff = ydiffnb

    #import pdb; pdb.set_trace()

    datas = []

    # The data is weighted according to barycentric coordinates
    bc0 = col(barycentric[pys, pxs, 0])
    bc1 = col(barycentric[pys, pxs, 1])
    bc2 = col(barycentric[pys, pxs, 2])
    for k in range(n_channels):
        dxs = xdiff[pys, pxs, k]
        dys = ydiff[pys, pxs, k]
        if f.shape[1] == 3:
            datas.append(np.hstack((col(dxs)*bc0,col(dys)*bc0,col(dxs)*bc1,col(dys)*bc1,col(dxs)*bc2,col(dys)*bc2)).ravel())
        else:
            datas.append(np.hstack((col(dxs)*bc0,col(dys)*bc0,col(dxs)*bc1,col(dys)*bc1)).ravel())

    data = np.concatenate(datas)

    ij = np.vstack((IS.ravel(), JS.ravel()))
    result = sp.csc_matrix((data, ij), shape=(image_width*image_height*n_channels, num_verts*2))

    return result
709189
def from_unit_cube(self, x):
        

        x0 = self.x0.value
        gamma = self.gamma.value

        half_pi = 1.57079632679

        res = np.tan(np.pi * x - half_pi) * gamma + x0

        return res
137442
def _gamma(self, R, phi):
        """"""
        return self._N * (phi - self._phi_ref - np.log(R / self._r_ref) / self._tan_alpha)
677816
def refraction(alt_degrees, temperature_C, pressure_mbar):
    
    r = 0.016667 / tan((alt_degrees + 7.31 / (alt_degrees + 4.4)) * DEG2RAD)
    d = r * (0.28 * pressure_mbar / (temperature_C + 273.0))
    return where((-1.0 <= alt_degrees) & (alt_degrees <= 89.9), d, 0.0)
714802
def compute_sun
      # === 太陽が光を発した時刻 t1(JD) の計算
      t_1_jd = calc_t1(:sun, @tdb)
      # === 時刻 t1(= TDB) におけるの位置・速度（ICRS 座標）の計算 (地球, 月, 太陽)
      Const::BODIES.each { |k, v| @icrs_1[k] = get_icrs(v, t_1_jd) }
      # === 時刻 t2 における地球重心から時刻 t1 における太陽への方向ベクトルの計算
      v_12 = calc_unit_vector(@icrs_2[:earth][0,3], @icrs_1[:sun][0,3])
      # === GCRS 座標系: 光行差の補正（方向ベクトルの Lorentz 変換）
      dd = conv_lorentz(v_12)
      pos_sun = calc_pos(dd, @r_e[:sun])
      # === 瞬時の真座標系: GCRS への bias & precession（歳差） & nutation（章動）の適用
      bpn = EphBpn.new(@tdb.strftime("%Y%m%d%H%M%S"))
      pos_sun_bpn = bpn.apply_bpn(pos_sun)
      # === 座標変換
      eq_pol_s  = MkCoord.rect2pol(pos_sun_bpn)
      ec_rect_s = MkCoord.rect_eq2ec(pos_sun_bpn, bpn.eps)
      ec_pol_s  = MkCoord.rect2pol(ec_rect_s)
      # === 視半径／（地平）視差計算
      radius = Math.asin(@asun / (eq_pol_s[2] * Const::AU / 1000))
      radius *= 180 / Math::PI * 3600
      parallax = Math.asin(@re / (eq_pol_s[2] * Const::AU / 1000))
      parallax *= 180 / Math::PI * 3600
      return [eq_pol_s, ec_pol_s, [radius, parallax]]
    rescue => e
      raise
    end
5254
def calibrate(self, data, calibration):
        """"""
        tic = datetime.now()

        if calibration == 'counts':
            return data
        elif calibration == 'radiance':
            raise NotImplementedError("Can't calibrate to radiance.")
        else:
            cal = self.calibration_table
            res = data.data.map_blocks(self._interp, cal, dtype=cal[:, 0].dtype)
            res = xr.DataArray(res,
                               dims=data.dims, attrs=data.attrs,
                               coords=data.coords)
        res = res.where(data > 0)
        logger.debug("Calibration time " + str(datetime.now() - tic))
        return res
615585
def determine_coefficient
      # Array of photon energies (in units of MeV).
      @energies = [
        0.001,
        0.0015,
        0.002,
        0.003,
        0.004,
        0.005,
        0.006,
        0.008,
        0.01,
        0.015,
        0.02,
        0.03,
        0.04,
        0.05,
        0.06,
        0.08,
        0.1,
        0.15,
        0.2,
        0.3,
        0.4,
        0.5,
        0.6,
        0.8,
        1.0,
        1.25,
        1.5,
        2.0,
        3.0,
        4.0,
        5.0,
        6.0,
        8.0,
        10.0,
        15.0,
        20.0
      ]
      # Array of mass attenuation coefficients for the above energies, for liquid water (in units of cm^2/g):
      @att_coeffs = [
        4078,
        1376,
        617.3,
        192.9,
        82.78,
        42.58,
        24.64,
        10.37,
        5.329,
        1.673,
        0.8096,
        0.3756,
        0.2683,
        0.2269,
        0.2059,
        0.1837,
        0.1707,
        0.1505,
        0.137,
        0.1186,
        0.1061,
        0.09687,
        0.08956,
        0.07865,
        0.07072,
        0.06323,
        0.05754,
        0.04942,
        0.03969,
        0.03403,
        0.03031,
        0.0277,
        0.02429,
        0.02219,
        0.01941,
        0.01813
      ]
      # Determine the coefficient:
      if @energy >= 20.0
        # When the energy is above 20, we use the coefficient for 20 MeV:
        @ac_water = @att_coeffs.last
      else
        if i = @energies.index(@energy)
          # When it exactly matches one of the listed energies, use its corresponding coefficient:
          @ac_water = @att_coeffs[i]
        else
          # When the given energy is between two of the listed values, interpolate:
          i_after = @energies.index {|x| x > @energy}
          if i_after
            e_high = @energies[i_after]
            e_low = @energies[i_after - 1]
            ac_high = @att_coeffs[i_after]
            ac_low = @att_coeffs[i_after - 1]
            @ac_water = (ac_high - ac_low ) / (e_high - e_low) * (@energy - e_low)
          else
            raise "Unexpected behaviour with index in the energy interpolation logic."
          end
        end
      end
    end
1177793
def get_SZ_single(self, alpha=None, beta=None):
        
        if alpha == None:
            alpha = self.alpha
        if beta == None:
            beta = self.beta

        tm_outdated = self._tm_signature != (self.radius, self.radius_type, 
            self.wavelength, self.m, self.axis_ratio, self.shape, self.ddelt, 
            self.ndgs)
        if tm_outdated:
            self._init_tmatrix()

        scatter_outdated = self._scatter_signature != (self.thet0, self.thet, 
            self.phi0, self.phi, alpha, beta, self.orient)

        outdated = tm_outdated or scatter_outdated

        if outdated:
            (self._S_single, self._Z_single) = pytmatrix.calcampl(self.nmax, 
                self.wavelength, self.thet0, self.thet, self.phi0, self.phi, 
                alpha, beta)
            self._set_scatter_signature()

        return (self._S_single, self._Z_single)
795457
def rdirichlet(theta, size=1):
    
    gammas = np.vstack([rgamma(theta, 1) for i in xrange(size)])
    if size > 1 and np.size(theta) > 1:
        return (gammas.T / gammas.sum(1))[:-1].T
    elif np.size(theta) > 1:
        return (gammas[0] / gammas[0].sum())[:-1]
    else:
        return 1.