Python astropy.units.nm() Examples

def get_equivalencies():
    """
    Return a list of example equivalencies for testing serialization.
    """
    return [eq.plate_scale(.3 * u.deg/u.mm), eq.pixel_scale(.5 * u.deg/u.pix),
            eq.spectral_density(350 * u.nm, factor=2),
            eq.spectral_density(350 * u.nm), eq.spectral(),
            eq.brightness_temperature(500 * u.GHz),
            eq.brightness_temperature(500 * u.GHz, beam_area=23 * u.sr),
            eq.with_H0(), eq.temperature_energy(), eq.temperature(),
            eq.thermodynamic_temperature(300 * u.Hz),
            eq.thermodynamic_temperature(140 * u.GHz, Planck15.Tcmb0),
            eq.beam_angular_area(3 * u.sr), eq.mass_energy(),
            eq.molar_mass_amu(), eq.doppler_relativistic(2 * u.m),
            eq.doppler_optical(2 * u.nm), eq.doppler_radio(2 * u.Hz),
            eq.parallax(), eq.logarithmic(), eq.dimensionless_angles(),
            eq.spectral() + eq.temperature(),
            (eq.spectral_density(35 * u.nm) +
             eq.brightness_temperature(5 * u.Hz, beam_area=2 * u.sr)),
            (eq.spectral() + eq.spectral_density(35 * u.nm) +
             eq.brightness_temperature(5 * u.Hz, beam_area=2 * u.sr))
            ] 

def test_compose_equivalencies():
    x = u.Unit("arcsec").compose(units=(u.pc,), equivalencies=u.parallax())
    assert x[0] == u.pc

    x = u.Unit("2 arcsec").compose(units=(u.pc,), equivalencies=u.parallax())
    assert x[0] == u.Unit(0.5 * u.pc)

    x = u.degree.compose(equivalencies=u.dimensionless_angles())
    assert u.Unit(u.degree.to(u.radian)) in x

    x = (u.nm).compose(units=(u.m, u.s), equivalencies=u.doppler_optical(0.55*u.micron))
    for y in x:
        if y.bases == [u.m, u.s]:
            assert y.powers == [1, -1]
            assert_allclose(
                y.scale,
                u.nm.to(u.m / u.s, equivalencies=u.doppler_optical(0.55 * u.micron)))
            break
    else:
        assert False, "Didn't find speed in compose results" 

def calc_brightness_ratio_from_Teff(Teff_1, Teff_2, plot=False):

    bandpass = np.genfromtxt('tess-response-function-v1.0.csv', delimiter=',', names=['wavelength','transmission'])
    wavelength_grid = np.arange(500,2000,1)
    
    if plot:
        fig, ax = plt.subplots()
        ax.plot(list(bandpass['wavelength'])+[2000], list(bandpass['transmission'])+[0], lw=2)
        ax.set(ylabel='TESS Transmission')
        ax2 = ax.twinx()
        ax2.plot(wavelength_grid, blackbody_lambda(wavelength_grid*u.nm, Teff_1*u.K), 'r-', lw=2, color='darkorange')
        ax2.plot(wavelength_grid, blackbody_lambda(wavelength_grid*u.nm, Teff_2*u.K), 'r-', lw=2, color='brown')
        ax2.set(ylabel='Blackbody Flux\n'+r'($erg \, cm^{-2} \, s^{-1} \, A^{-1} \, sr^{-1}$)')
    
    int1 = np.trapz(bandpass['transmission']*u.nm*blackbody_lambda(bandpass['wavelength']*u.nm, Teff_1*u.K), x=bandpass['wavelength']*u.nm, dx=np.diff(bandpass['wavelength']*u.nm))
    int2 = np.trapz(bandpass['transmission']*u.nm*blackbody_lambda(bandpass['wavelength']*u.nm, Teff_2*u.K), x=bandpass['wavelength']*u.nm, dx=np.diff(bandpass['wavelength']*u.nm))
    sbratio = int2/int1
    
    return sbratio 

def test_evaluate_with_quantities_and_equivalencies():
    """
    We now make sure that equivalencies are correctly taken into account
    """

    g = Gaussian1D(1 * u.Jy, 10 * u.nm, 2 * u.nm)

    # We aren't setting the equivalencies, so this won't work
    with pytest.raises(UnitsError) as exc:
        g(30 * u.PHz)
    assert exc.value.args[0] == ("Gaussian1D: Units of input 'x', PHz (frequency), could "
                                 "not be converted to required input units of "
                                 "nm (length)")

    # But it should now work if we pass equivalencies when evaluating
    assert_quantity_allclose(g(30 * u.PHz, equivalencies={'x': u.spectral()}),
                             g(9.993081933333332 * u.nm)) 

def test_diffprop_matches_airydisk(efl, epd, wvl):
    fno = efl / epd

    p = Pupil(dia=epd, xy_unit=u.mm, z_unit=u.nm, wavelength=mkwvl(wvl, u.um))
    psf = PSF.from_pupil(p, efl, Q=3)  # use Q=3 not Q=4 for improved accuracy
    s = psf.slices()
    u_, sx = s.x
    u_, sy = s.y
    analytic = airydisk(u_, fno, wvl)
    assert np.allclose(sx, analytic, atol=PRECISION)
    assert np.allclose(sy, analytic, atol=PRECISION) 

def test_linear1d_quantity(tmpdir):
    model = astmodels.Linear1D(1*u.nm, 1*(u.nm/u.pixel))
    tree = {'model': model}
    helpers.assert_roundtrip_tree(tree, tmpdir) 

def test_spectral3():
    a = u.nm.to(u.Hz, [1000, 2000], u.spectral())
    assert_allclose(a, [2.99792458e+14, 1.49896229e+14]) 

def test_spectral2():
    a = u.nm.to(u.J, 500, u.spectral())
    assert_allclose(a, 3.972891366538605e-19)
    b = u.J.to(u.nm, a, u.spectral())
    assert_allclose(b, 500)

    a = u.AA.to(u.Hz, 1, u.spectral())
    b = u.Hz.to(u.J, a, u.spectral())
    c = u.AA.to(u.J, 1, u.spectral())
    assert_allclose(b, c)

    c = u.J.to(u.Hz, b, u.spectral())
    assert_allclose(a, c) 

def test_doppler_wavelength_circle(function):
    rest = 105.01 * u.nm
    shifted = 105.03 * u.nm
    velo = shifted.to(u.km / u.s, equivalencies=function(rest))
    wav = velo.to(u.nm, equivalencies=function(rest))
    np.testing.assert_almost_equal(wav.value, shifted.value, decimal=7) 

def test_b_wien():
    """b_wien should give the correct peak wavelength for
    given blackbody temperature. The Sun is used in this test.

    """
    from astropy.constants.astropyconst13 import b_wien
    from astropy import units as u
    t = 5778 * u.K
    w = (b_wien / t).to(u.nm)
    assert round(w.value) == 502 

def test_b_wien():
    """b_wien should give the correct peak wavelength for
    given blackbody temperature. The Sun is used in this test.

    """
    from astropy.constants import b_wien
    from astropy import units as u
    t = 5778 * u.K
    w = (b_wien / t).to(u.nm)
    assert round(w.value) == 502 

def test_tabular_interp_1d():
    """
    Test Tabular1D model.
    """
    points = np.arange(0, 5)
    values = [1., 10, 2, 45, -3]
    LookupTable = models.tabular_model(1)
    model = LookupTable(points=points, lookup_table=values)
    xnew = [0., .7, 1.4, 2.1, 3.9]
    ans1 = [1., 7.3, 6.8, 6.3, 1.8]
    assert_allclose(model(xnew), ans1)
    # Test evaluate without passing `points`.
    model = LookupTable(lookup_table=values)
    assert_allclose(model(xnew), ans1)
    # Test bounds error.
    xextrap = [0., .7, 1.4, 2.1, 3.9, 4.1]
    with pytest.raises(ValueError):
        model(xextrap)
    # test extrapolation and fill value
    model = LookupTable(lookup_table=values, bounds_error=False,
                        fill_value=None)
    assert_allclose(model(xextrap),
                    [1., 7.3, 6.8, 6.3, 1.8, -7.8])

    # Test unit support
    xnew = xnew * u.nm
    ans1 = ans1 * u.nJy
    model = LookupTable(points=points*u.nm, lookup_table=values*u.nJy)
    assert_quantity_allclose(model(xnew), ans1)
    assert_quantity_allclose(model(xnew.to(u.nm)), ans1)
    assert model.bounding_box == (0 * u.nm, 4 * u.nm)

    # Test fill value unit conversion and unitless input on table with unit
    model = LookupTable([1, 2, 3], [10, 20, 30] * u.nJy, bounds_error=False,
                        fill_value=1e-33*(u.W / (u.m * u.m * u.Hz)))
    assert_quantity_allclose(model(np.arange(5)),
                             [100, 10, 20, 30, 100] * u.nJy) 

def download_grond(filter_dict):

    save_path = os.path.join(get_speclite_filter_path(), "ESO")

    if_directory_not_existing_then_make(save_path)

    grond_filter_url = "http://www.mpe.mpg.de/~jcg/GROND/GROND_filtercurves.txt"

    url_response = urllib.request.urlopen(grond_filter_url)

    grond_table = pd.read_table(url_response)

    wave = grond_table["A"].as_matrix()

    bands = ["g", "r", "i", "z", "H", "J", "K"]

    for band in bands:

        curve = np.array(grond_table["%sBand" % band])
        curve[curve < 0] = 0
        curve[0] = 0
        curve[-1] = 0

        grond_spec = spec_filter.FilterResponse(
            wavelength=wave * u.nm,
            response=curve,
            meta=dict(group_name="GROND", band_name=band),
        )

        grond_spec.save(directory_name=save_path)

    filter_dict["ESO"] = {"GROND": bands}

    return filter_dict 

def test_additional_meta_data_from_mc_header():
    with SimTelEventSource(input_url=gamma_test_large_path) as reader:
        data = next(iter(reader))

    # for expectation values
    from astropy import units as u
    from astropy.coordinates import Angle

    assert data.mcheader.corsika_version == 6990
    assert data.mcheader.spectral_index == -2.0
    assert data.mcheader.shower_reuse == 20
    assert data.mcheader.core_pos_mode == 1
    assert data.mcheader.diffuse == 1
    assert data.mcheader.atmosphere == 26

    # value read by hand from input card
    name_expectation = {
        "energy_range_min": u.Quantity(3.0e-03, u.TeV),
        "energy_range_max": u.Quantity(3.3e02, u.TeV),
        "prod_site_B_total": u.Quantity(23.11772346496582, u.uT),
        "prod_site_B_declination": Angle(0.0 * u.rad),
        "prod_site_B_inclination": Angle(-0.39641156792640686 * u.rad),
        "prod_site_alt": 2150.0 * u.m,
        "max_scatter_range": 3000.0 * u.m,
        "min_az": 0.0 * u.rad,
        "min_alt": 1.2217305 * u.rad,
        "max_viewcone_radius": 10.0 * u.deg,
        "corsika_wlen_min": 240 * u.nm,
    }

    for name, expectation in name_expectation.items():
        value = getattr(data.mcheader, name)

        assert value.unit == expectation.unit
        assert np.isclose(
            value.to_value(expectation.unit), expectation.to_value(expectation.unit)
        ) 

def _parse_mc_header(self):
        mc_run_head = self.file_.mc_run_headers[-1]

        return MCHeaderContainer(
            corsika_version=mc_run_head["shower_prog_vers"],
            simtel_version=mc_run_head["detector_prog_vers"],
            energy_range_min=mc_run_head["E_range"][0] * u.TeV,
            energy_range_max=mc_run_head["E_range"][1] * u.TeV,
            prod_site_B_total=mc_run_head["B_total"] * u.uT,
            prod_site_B_declination=Angle(mc_run_head["B_declination"], u.rad,),
            prod_site_B_inclination=Angle(mc_run_head["B_inclination"], u.rad,),
            prod_site_alt=mc_run_head["obsheight"] * u.m,
            spectral_index=mc_run_head["spectral_index"],
            shower_prog_start=mc_run_head["shower_prog_start"],
            shower_prog_id=mc_run_head["shower_prog_id"],
            detector_prog_start=mc_run_head["detector_prog_start"],
            detector_prog_id=mc_run_head["detector_prog_id"],
            num_showers=mc_run_head["n_showers"],
            shower_reuse=mc_run_head["n_use"],
            max_alt=mc_run_head["alt_range"][1] * u.rad,
            min_alt=mc_run_head["alt_range"][0] * u.rad,
            max_az=mc_run_head["az_range"][1] * u.rad,
            min_az=mc_run_head["az_range"][0] * u.rad,
            diffuse=mc_run_head["diffuse"],
            max_viewcone_radius=mc_run_head["viewcone"][1] * u.deg,
            min_viewcone_radius=mc_run_head["viewcone"][0] * u.deg,
            max_scatter_range=mc_run_head["core_range"][1] * u.m,
            min_scatter_range=mc_run_head["core_range"][0] * u.m,
            core_pos_mode=mc_run_head["core_pos_mode"],
            injection_height=mc_run_head["injection_height"] * u.m,
            atmosphere=mc_run_head["atmosphere"],
            corsika_iact_options=mc_run_head["corsika_iact_options"],
            corsika_low_E_model=mc_run_head["corsika_low_E_model"],
            corsika_high_E_model=mc_run_head["corsika_high_E_model"],
            corsika_bunchsize=mc_run_head["corsika_bunchsize"],
            corsika_wlen_min=mc_run_head["corsika_wlen_min"] * u.nm,
            corsika_wlen_max=mc_run_head["corsika_wlen_max"] * u.nm,
            corsika_low_E_detail=mc_run_head["corsika_low_E_detail"],
            corsika_high_E_detail=mc_run_head["corsika_high_E_detail"],
            run_array_direction=Angle(self.file_.header["direction"] * u.rad),
        ) 

def mkwvl(quantity, base=u.um):
    """Generate a new Wavelength unit.

    Parameters
    ----------
    quantity : `float` or `astropy.units.unit`
        number of (base) for the wavelength, e.g. quantity=632.8 with base=u.nm for HeNe.
        if an astropy unit, simply returned by this function
    base : `astropy.units.Unit`
        base unit, e.g. um or nm

    Returns
    -------
    `astropy.units.Unit`
        new Unit for appropriate wavelength

    """
    if quantity is None:
        return quantity
    elif not isinstance(quantity, u.Unit):
        return u.def_unit(['wave', 'wavelength'], quantity * base,
                          format={'latex': r'\lambda', 'unicode': 'λ'})
    else:
        return quantity


# IR 

def test_diffprop_matches_analyticmtf(efl, epd, wvl):
    fno = efl / epd
    p = Pupil(dia=epd, xy_unit=u.mm, z_unit=u.nm, wavelength=mkwvl(wvl, u.um))
    psf = PSF.from_pupil(p, efl)
    mtf = MTF.from_psf(psf)
    s = mtf.slices()
    u_, x = s.x
    u__, y = s.y

    analytic_1 = diffraction_limited_mtf(fno, wvl, frequencies=u_)
    analytic_2 = diffraction_limited_mtf(fno, wvl, frequencies=u__)
    assert np.allclose(analytic_1, x, atol=PRECISION)
    assert np.allclose(analytic_2, y, atol=PRECISION) 

def compare_interpolants(self, f1, f2, param, msg=''):
        r"""Compare two interpolants f1 and f2 by probing them randomly."""
        # Find out the number of input arguments expected
        f1_info = inspect.getargspec(f1)
        f2_info = inspect.getargspec(f2)
        # this is the number of formal arguments MINUS the number of defaults
        # (EXOSIMS uses defaults to provide functional closure, though it's unneeded)
        nargin1 = len(f1_info.args) - (0 if not f1_info.defaults else len(f1_info.defaults))
        nargin2 = len(f2_info.args) - (0 if not f2_info.defaults else len(f2_info.defaults))
        if nargin1 != nargin2:
            raise self.failureException(msg + '-- functions have different arity (arg lengths)')
        # make a few random probes of the interpolant on the interval (0,1)
        for count in range(10):
            # obtain a vector of length nargin1
            arg_in = np.random.random(nargin1)
            # the result can be a float (for contrast),
            # a numpy array (for PSF), or a Quantity (for QE)
            if  param in ('core_thruput', 'core_contrast'):
                out_1 = f1(arg_in[0]*u.nm,arg_in[1]*u.arcsec)
            else:
                out_1 = f1(*arg_in)
            out_2 = f2(*arg_in)
            diff = out_1 - out_2
            # if it's a quantity, unbox the difference
            if isinstance(diff, u.quantity.Quantity):
                diff = diff.value
            if np.any(np.abs(diff) > 1e-5):
                errmsg = msg + '-- function mismatch: %r != %r' % (out_1, out_2)
                raise self.failureException(errmsg) 

def test_tabular_interp_2d():
    table = np.array([
       [-0.04614432, -0.02512547, -0.00619557, 0.0144165, 0.0297525],
       [-0.04510594, -0.03183369, -0.01118008, 0.01201388, 0.02496205],
       [-0.05464094, -0.02804499, -0.00960086, 0.01134333, 0.02284104],
       [-0.04879338, -0.02539565, -0.00440462, 0.01795145, 0.02122417],
       [-0.03637372, -0.01630025, -0.00157902, 0.01649774, 0.01952131]])

    points = np.arange(0, 5)
    points = (points, points)

    xnew = np.array([0., .7, 1.4, 2.1, 3.9])
    LookupTable = models.tabular_model(2)
    model = LookupTable(points, table)
    znew = model(xnew, xnew)
    result = np.array(
        [-0.04614432, -0.03450009, -0.02241028, -0.0069727, 0.01938675])
    assert_allclose(znew, result, atol=1e-7)

    # test 2D arrays as input
    a = np.arange(12).reshape((3, 4))
    y, x = np.mgrid[:3, :4]
    t = models.Tabular2D(lookup_table=a)
    r = t(y, x)
    assert_allclose(a, r)

    with pytest.raises(ValueError):
        model = LookupTable(points=([1.2, 2.3], [1.2, 6.7], [3, 4]))
    with pytest.raises(ValueError):
        model = LookupTable(lookup_table=[1, 2, 3])
    with pytest.raises(NotImplementedError):
        model = LookupTable(n_models=2)
    with pytest.raises(ValueError):
        model = LookupTable(([1, 2], [3, 4]), [5, 6])
    with pytest.raises(ValueError):
        model = LookupTable(([1, 2] * u.m, [3, 4]), [[5, 6], [7, 8]])
    with pytest.raises(ValueError):
        model = LookupTable(points, table, bounds_error=False,
                            fill_value=1*u.Jy)

    # Test unit support
    points = points[0] * u.nm
    points = (points, points)
    xnew = xnew * u.nm
    model = LookupTable(points, table * u.nJy)
    result = result * u.nJy
    assert_quantity_allclose(model(xnew, xnew), result, atol=1e-7*u.nJy)
    xnew = xnew.to(u.m)
    assert_quantity_allclose(model(xnew, xnew), result, atol=1e-7*u.nJy)
    bbox = (0 * u.nm, 4 * u.nm)
    bbox = (bbox, bbox)
    assert model.bounding_box == bbox 

def image_prediction(
    mirror_radius,
    hole_radius,
    impact_parameter,
    phi,
    center_x,
    center_y,
    radius,
    ring_width,
    pixel_x,
    pixel_y,
    pixel_diameter,
    oversampling=3,
    min_lambda=300 * u.nm,
    max_lambda=600 * u.nm,
):
    """
    Parameters
    ----------
    impact_parameter: quantity[length]
        Impact distance of muon
    center_x: quantity[angle]
        Muon ring center in telescope frame
    center_y: quantity[angle]
        Muon ring center in telescope frame
    radius: quantity[angle]
        Radius of muon ring in telescope frame
    ring_width: quantity[angle]
        Gaussian width of muon ring
    pixel_x: quantity[angle]
        Pixel x coordinate in telescope
    pixel_y: quantity[angle]
        Pixel y coordinate in telescope

    Returns
    -------
    ndarray:
        Predicted signal
    """
    return image_prediction_no_units(
        mirror_radius.to_value(u.m),
        hole_radius.to_value(u.m),
        impact_parameter.to_value(u.m),
        phi.to_value(u.rad),
        center_x.to_value(u.rad),
        center_y.to_value(u.rad),
        radius.to_value(u.rad),
        ring_width.to_value(u.rad),
        pixel_x.to_value(u.rad),
        pixel_y.to_value(u.rad),
        pixel_diameter.to_value(u.rad),
        oversampling=oversampling,
        min_lambda_m=min_lambda.to_value(u.m),
        max_lambda_m=max_lambda.to_value(u.m),
    )