Python astropy.units.km() Examples

def test_against_jpl_horizons():
    """Check that Astropy gives consistent results with the JPL Horizons example.

    The input parameters and reference results are taken from this page:
    (from the first row of the Results table at the bottom of that page)
    http://ssd.jpl.nasa.gov/?horizons_tutorial
    """
    obstime = Time('1998-07-28 03:00')
    location = EarthLocation(lon=Angle('248.405300d'),
                             lat=Angle('31.9585d'),
                             height=2.06 * u.km)
    # No atmosphere
    altaz_frame = AltAz(obstime=obstime, location=location)

    altaz = SkyCoord('143.2970d 2.6223d', frame=altaz_frame)
    radec_actual = altaz.transform_to('icrs')
    radec_expected = SkyCoord('19h24m55.01s -40d56m28.9s', frame='icrs')
    distance = radec_actual.separation(radec_expected).to('arcsec')
    # Current value: 0.238111 arcsec
    assert distance < 1 * u.arcsec 

def test_moon_earth(self):
        r"""Test moon_earth method.

        Approach: Reference to pre-computed result from Matlab.
        """

        print('moon_earth()')
        obs = self.fixture
        # TODO: add other times besides this one
        # century = 0
        t_ref_string = '2000-01-01T12:00:00.0'
        t_ref = Time(t_ref_string, format='isot', scale='utc')
        moon = obs.moon_earth(t_ref).flatten().to(u.km)
        # print moon
        r_earth = 6378.137 # earth radius [km], as used in Vallado's code
        moon_ref = [-45.74169421, -41.80825511, -11.88954996] # pre-computed from Matlab
        for coord in range(3):
            self.assertAlmostEqual(moon[coord].value, moon_ref[coord] * r_earth, places=1) 

def test_gen_radius(self):
        r"""Test gen_radius method.

        Approach: Ensures the output is set, of the correct type, length, and units.
        Check distributional agreement.
        """
        plan_pop = self.fixture
        n = 10000
        radii = plan_pop.gen_radius(n)

        # ensure the units are length
        self.assertEqual((radii/u.km).decompose().unit, u.dimensionless_unscaled)
        # radius > 0
        self.assertTrue(np.all(radii.value > 0))
        self.assertTrue(np.all(np.isfinite(radii.value)))

        h = np.histogram(radii.to('earthRad').value,bins=plan_pop.Rs)
        np.testing.assert_allclose(plan_pop.Rvals.sum()*h[0]/float(n),plan_pop.Rvals,rtol=0.05) 

def test_gen_sma(self):
        r"""Test gen_sma method.

        Approach: Ensures the output is set, of the correct type, length, and units.
        Check that they are in the correct range and follow the distribution.
        """

        plan_pop = self.fixture
        n = 10000
        sma = plan_pop.gen_sma(n)

        # ensure the units are length
        self.assertEqual((sma/u.km).decompose().unit, u.dimensionless_unscaled)
        # sma > 0
        self.assertTrue(np.all(sma.value >= 0))
        # sma >= arange[0], sma <= arange[1]
        self.assertTrue(np.all(sma - plan_pop.arange[0] >= 0))
        self.assertTrue(np.all(plan_pop.arange[1] - sma >= 0))

        h = np.histogram(sma.to('AU').value,100,density=True)
        hx = np.diff(h[1])/2.+h[1][:-1]
        hp = plan_pop.dist_sma(hx)

        chi2 = scipy.stats.chisquare(h[0],hp)
        self.assertGreaterEqual(chi2[1],0.95) 

def test_orbit(self):
        """
        Test orbit method to ensure proper output
        """

        t_ref = Time(2000.0, format='jyear')
        for mod in self.allmods:
            with RedirectStreams(stdout=self.dev_null):
                if 'SotoStarshade' in mod.__name__:
                    obj = mod(f_nStars=4, **copy.deepcopy(self.spec))
                else:
                    obj = mod(**copy.deepcopy(self.spec))
            r_sc = obj.orbit(t_ref)
            # the r_sc attribute is set and is a 3-tuple of astropy Quantity's
            self.assertEqual(type(r_sc), type(1.0 * u.km))
            self.assertEqual(r_sc.shape, (1, 3)) 

def test_gen_sma(self):
        r"""Test gen_sma method.

        Approach: Ensures the output is set, of the correct type, length, and units.
        Check that they are in the correct range and follow the distribution.
        """

        plan_pop = self.fixture
        n = 10000
        sma = plan_pop.gen_sma(n)

        # ensure the units are length
        self.assertEqual((sma/u.km).decompose().unit, u.dimensionless_unscaled)
        # sma > 0
        self.assertTrue(np.all(sma.value >= 0))
        # sma >= arange[0], sma <= arange[1]
        self.assertTrue(np.all(sma - plan_pop.arange[0] >= 0))
        self.assertTrue(np.all(plan_pop.arange[1] - sma >= 0))

        h = np.histogram(sma.to('AU').value,100,density=True)
        hx = np.diff(h[1])/2.+h[1][:-1]
        hp = plan_pop.dist_sma(hx)

        chi2 = scipy.stats.chisquare(h[0],hp)
        self.assertGreaterEqual(chi2[1],0.95) 

def __init__(self, a, e, I, O, w, lM):
            
            # store list of semimajor axis values (convert from AU to km)
            self.a = (a*u.AU).to('km').value
            if not isinstance(self.a, float):
                self.a = self.a.tolist()
            # store list of dimensionless eccentricity values
            self.e = e
            # store list of inclination values (degrees)
            self.I = I
            # store list of right ascension of ascending node values (degrees)
            self.O = O 
            # store list of longitude of periapsis values (degrees)
            self.w = w 
            # store list of mean longitude values (degrees)
            self.lM = lM 

def s_a(self, lat, lon):
        """ Standard deviation of terrain heights (m) within a 110 km × 110 km
        area with a 30 s resolution (e.g. the Globe “gtopo30” data).
        The value for the mid-path may be obtained from an area roughness
        with 0.5 × 0.5 degree resolution of geographical coordinates
        using bi-linear interpolation.
        """
        if not self._s_a:
            vals = load_data(os.path.join(dataset_dir, '530/v16_gtopo_30.txt'))
            lats = load_data(os.path.join(dataset_dir, '530/v16_lat.txt'))
            lons = load_data(os.path.join(dataset_dir, '530/v16_lon.txt'))
            self._Pr6 = bilinear_2D_interpolator(lats, lons, vals)

        return self._Pr6(
            np.array([lat.ravel(), lon.ravel()]).T).reshape(lat.shape)

    ###########################################################################
    #                               Section 2.2                               #
    ########################################################################### 

def test_all_to_value():
    """test all_to_value"""
    x_m = np.arange(5) * u.m
    y_mm = np.arange(5) * 1000 * u.mm
    z_km = np.arange(5) * 1e-3 * u.km
    nono_deg = np.arange(5) * 1000 * u.deg

    # one argument
    x = all_to_value(x_m, unit=u.m)
    assert (x == np.arange(5)).all()

    # two arguments
    x, y = all_to_value(x_m, y_mm, unit=u.m)
    assert (x == np.arange(5)).all()
    assert (y == np.arange(5)).all()

    # three
    x, y, z = all_to_value(x_m, y_mm, z_km, unit=u.m)
    assert (x == np.arange(5)).all()
    assert (y == np.arange(5)).all()
    assert (z == np.arange(5)).all()

    # cannot be converted
    with pytest.raises(u.UnitConversionError):
        all_to_value(x_m, nono_deg, unit=x_m.unit) 

def plot_spectrum(self):
        ''' Plot the HI spectrum '''

        if self._specfile:
            if not self._specdata:
                self._specdata = self._get_data(self._specfile)

            vel = self._specdata['VHI'][0]
            flux = self._specdata['FHI'][0]
            spec = Spectrum(flux, unit=u.Jy, wavelength=vel,
                            wavelength_unit=u.km / u.s)
            ax = spec.plot(
                ylabel='HI\ Flux', xlabel='Velocity', title=self.targetid, ytrim='minmax'
            )
            return ax
        return None

#
# Functions to become available on your VAC in marvin.tools.vacs.VACs 

def test_nddata_init_data_nddata_subclass():
    uncert = StdDevUncertainty(3)
    # There might be some incompatible subclasses of NDData around.
    bnd = BadNDDataSubclass(False, True, 3, 2, 'gollum', 100)
    # Before changing the NDData init this would not have raised an error but
    # would have lead to a compromised nddata instance
    with pytest.raises(TypeError):
        NDData(bnd)

    # but if it has no actual incompatible attributes it passes
    bnd_good = BadNDDataSubclass(np.array([1, 2]), uncert, 3, HighLevelWCSWrapper(WCS(naxis=1)),
                                 {'enemy': 'black knight'}, u.km)
    nd = NDData(bnd_good)
    assert nd.unit == bnd_good.unit
    assert nd.meta == bnd_good.meta
    assert nd.uncertainty == bnd_good.uncertainty
    assert nd.mask == bnd_good.mask
    assert nd.wcs is bnd_good.wcs
    assert nd.data is bnd_good.data 

def test_param_unit():
    with pytest.raises(ValueError):
        NDData(np.ones((5, 5)), unit="NotAValidUnit")
    NDData([1, 2, 3], unit='meter')
    # Test conflicting units (quantity as data)
    q = np.array([1, 2, 3]) * u.m
    nd = NDData(q, unit='cm')
    assert nd.unit != q.unit
    assert nd.unit == u.cm
    # (masked quantity)
    mq = np.ma.array(np.array([2, 3])*u.m, mask=False)
    nd2 = NDData(mq, unit=u.s)
    assert nd2.unit == u.s
    # (another NDData as data)
    nd3 = NDData(nd, unit='km')
    assert nd3.unit == u.km 

def test_ecliptic_heliobary():
    """
    Check that the ecliptic transformations for heliocentric and barycentric
    at least more or less make sense
    """
    icrs = ICRS(1*u.deg, 2*u.deg, distance=1.5*R_sun)

    bary = icrs.transform_to(BarycentricMeanEcliptic)
    helio = icrs.transform_to(HeliocentricMeanEcliptic)

    # make sure there's a sizable distance shift - in 3d hundreds of km, but
    # this is 1D so we allow it to be somewhat smaller
    assert np.abs(bary.distance - helio.distance) > 1*u.km

    # now make something that's got the location of helio but in bary's frame.
    # this is a convenience to allow `separation` to work as expected
    helio_in_bary_frame = bary.realize_frame(helio.cartesian)
    assert bary.separation(helio_in_bary_frame) > 1*u.arcmin 

def test_cunit():
    w = _wcs.Wcsprm()
    assert list(w.cunit) == [u.Unit(''), u.Unit('')]
    w.cunit = [u.m, 'km']
    assert w.cunit[0] == u.m
    assert w.cunit[1] == u.km 

def test_load_w_quantity(self):
        try:
            from astropy import units as u
            from pvl.decoder import OmniDecoder
            pvl_file = 'tests/data/pds3/units1.lbl'
            km_upper = u.def_unit('KM', u.km)
            m_upper = u.def_unit('M', u.m)
            u.add_enabled_units([km_upper, m_upper])
            label = pvl.load(pvl_file,
                             decoder=OmniDecoder(quantity_cls=u.Quantity))
            self.assertEqual(label['FLOAT_UNIT'], u.Quantity(0.414, 'KM'))
        except ImportError:
            pass 

def test_nddata_init_data_nddata():
    nd1 = NDData(np.array([1]))
    nd2 = NDData(nd1)
    assert nd2.wcs == nd1.wcs
    assert nd2.uncertainty == nd1.uncertainty
    assert nd2.mask == nd1.mask
    assert nd2.unit == nd1.unit
    assert nd2.meta == nd1.meta

    # Check that it is copied by reference
    nd1 = NDData(np.ones((5, 5)))
    nd2 = NDData(nd1)
    assert nd1.data is nd2.data

    # Check that it is really copied if copy=True
    nd2 = NDData(nd1, copy=True)
    nd1.data[2, 3] = 10
    assert nd1.data[2, 3] != nd2.data[2, 3]

    # Now let's see what happens if we have all explicitly set
    nd1 = NDData(np.array([1]), mask=False, uncertainty=StdDevUncertainty(10), unit=u.s,
                 meta={'dest': 'mordor'}, wcs=WCS(naxis=1))
    nd2 = NDData(nd1)
    assert nd2.data is nd1.data
    assert nd2.wcs is nd1.wcs
    assert nd2.uncertainty.array == nd1.uncertainty.array
    assert nd2.mask == nd1.mask
    assert nd2.unit == nd1.unit
    assert nd2.meta == nd1.meta

    # now what happens if we overwrite them all too
    nd3 = NDData(nd1, mask=True, uncertainty=StdDevUncertainty(200), unit=u.km,
                 meta={'observer': 'ME'}, wcs=WCS(naxis=1))
    assert nd3.data is nd1.data
    assert nd3.wcs is not nd1.wcs
    assert nd3.uncertainty.array != nd1.uncertainty.array
    assert nd3.mask != nd1.mask
    assert nd3.unit != nd1.unit
    assert nd3.meta != nd1.meta 

def water_vapour_density(lat, h, season='summer'):
    """
    Method to determine the water-vapour density as a
    function of altitude and latitude, for calculating gaseous attenuation
    along an Earth-space path. This method is recommended when more reliable
    local data are not available.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    h : number or Quantity
        Height (km)
    season : string
        Season of the year (available values, 'summer', and 'winter').
        Default 'summer'


    Returns
    -------
    rho: Quantity
        Water vapour density (g/m^3)


    References
    ----------
    [1] Reference Standard Atmospheres
    https://www.itu.int/rec/R-REC-P.835/en
    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    h = prepare_quantity(h, u.km, 'Height')
    val = __model.water_vapour_density(lat, h, season)
    return prepare_output_array(val, type_output) * u.g / u.m**3 

def pressure(lat, h, season='summer'):
    """
    Method to determine the pressure as a function of altitude and latitude,
    for calculating gaseous attenuation along an Earth-space path.
    This method is recommended when more reliable local data are not available.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    h : number or Quantity
        Height (km)
    season : string
        Season of the year (available values, 'summer', and 'winter').
        Default 'summer'


    Returns
    -------
    P: Quantity
        Pressure (hPa)


    References
    ----------
    [1] Reference Standard Atmospheres
    https://www.itu.int/rec/R-REC-P.835/en
    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    h = prepare_quantity(h, u.km, 'Height')
    val = __model.pressure(lat, h, season)
    return prepare_output_array(val, type_output) * u.hPa 

def temperature(lat, h, season='summer'):
    """
    Method to determine the temperature as a function of altitude and latitude,
    for calculating gaseous attenuation along an Earth-space path. This method
    is recommended when more reliable local data are not available.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    h : number or Quantity
        Height (km)
    season : string
        Season of the year (available values, 'summer', and 'winter').
        Default 'summer'


    Returns
    -------
    T: Quantity
        Temperature (K)


    References
    ----------
    [1] Reference Standard Atmospheres
    https://www.itu.int/rec/R-REC-P.835/en

    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    h = prepare_quantity(h, u.km, 'Height')
    val = __model.temperature(lat, h, season)
    return prepare_output_array(val, type_output) * u.Kelvin 

def correct_vel(ra, dec, vel, vlsr=vlsr0):
    """Corrects the proper motion for the speed of the Sun
    Arguments:
        ra - RA in deg
        dec -- Declination in deg
        pmra -- pm in RA in mas/yr
        pmdec -- pm in declination in mas/yr
        dist -- distance in kpc
    Returns:
        (pmra,pmdec) the tuple with the proper motions corrected for the Sun's motion
    """
    
    C=acoo.ICRS(ra=ra*auni.deg,dec=dec*auni.deg,
                    radial_velocity=vel*auni.km/auni.s,
                    distance=np.ones_like(vel)*auni.kpc,
                    pm_ra_cosdec=np.zeros_like(vel)*auni.mas/auni.year,
                    pm_dec=np.zeros_like(vel)*auni.mas/auni.year)
    #frame = acoo.Galactocentric (galcen_vsun = np.array([ 11.1, vlsr+12.24, 7.25])*auni.km/auni.s)
    kw = dict(galcen_v_sun = acoo.CartesianDifferential(np.array([ 11.1, vlsr+12.24, 7.25])*auni.km/auni.s))                                 
    frame = acoo.Galactocentric (**kw)
    Cg = C.transform_to(frame)
    Cg1 = acoo.Galactocentric(x=Cg.x, y=Cg.y, z=Cg.z,
                              v_x=Cg.v_x*0,
                              v_y=Cg.v_y*0,
                              v_z=Cg.v_z*0, **kw)
    C1=Cg1.transform_to(acoo.ICRS)
    return np.asarray(((C.radial_velocity-C1.radial_velocity)/(auni.km/auni.s)).decompose()) 

def rain_specific_attenuation_coefficients(f, el, tau):
    """
    A method to compute the values for the coefficients k and α to compute
    the specific attenuation γ_R (dB/km)


    Parameters
    ----------
    f : number or Quantity
        Frequency (GHz)
    el : number, sequence, or numpy.ndarray
        Elevation angle of the receiver points
    tau : number, sequence, or numpy.ndarray
        Polarization tilt angle relative to the horizontal (degrees). Tau = 45
        deg for circular polarization)


    Returns
    -------
    k: number
        Zero isoterm height (km)
    alpha: number
        Zero isoterm height (km)


    References
    ----------
    [1] Rain height model for prediction methods:
    https://www.itu.int/rec/R-REC-P.838/en
    """
    global __model
    f = prepare_quantity(f, u.GHz, 'Frequency')
    return __model.rain_specific_attenuation_coefficients(f, el, tau) 

def diffraction_loss(d1, d2, h, f):
    """
    Diffraction loss over average terrain. This value is valid for losses
    greater than 15 dB.


    Parameters
    ----------
    d1 : number, sequence, or numpy.ndarray
        Distances from the first terminal to the path obstruction. [km]
    d2 : number, sequence, or numpy.ndarray
        Distances from the second terminal to the path obstruction. [km]
    h : number, sequence, or numpy.ndarray
        Height difference between most significant path blockage
        and the path trajectory. h is negative if the top of the obstruction
        of interest is above the virtual line-of-sight. [m]
    f : number
        Frequency of the link [GHz]


    Returns
    -------
    A_d: Quantity
        Diffraction loss over average terrain  [dB]


    References
    ----------
    [1] Propagation data and prediction methods required for the design of
    terrestrial line-of-sight systems: https://www.itu.int/rec/R-REC-P.530/en
    """
    global __model
    type_output = type(d1)
    d1 = prepare_quantity(d1, u.km, 'Distance to the first terminal')
    d2 = prepare_quantity(d2, u.km, 'Distance to the second terminal')
    h = prepare_quantity(h, u.m, 'Height difference')
    f = prepare_quantity(f, u.GHz, 'Frequency')

    val = __model.diffraction_loss(d1, d2, h, f)
    return prepare_output_array(val, type_output) * u.m 

def fresnel_ellipse_radius(d1, d2, f):
    """
    Computes the radius of the first Fresnel ellipsoid.


    Parameters
    ----------
    d1 : number, sequence, or numpy.ndarray
        Distances from the first terminal to the path obstruction. [km]
    d2 : number, sequence, or numpy.ndarray
        Distances from the second terminal to the path obstruction. [km]
    f : number
        Frequency of the link [GHz]


    Returns
    -------
    F1: Quantity
       Radius of the first Fresnel ellipsoid [m]


    References
    ----------
    [1] Propagation data and prediction methods required for the design of
    terrestrial line-of-sight systems: https://www.itu.int/rec/R-REC-P.530/en
    """
    global __model
    type_output = type(d1)
    d1 = prepare_quantity(d1, u.km, 'Distance to the first terminal')
    d2 = prepare_quantity(d2, u.km, 'Distance to the second terminal')
    f = prepare_quantity(f, u.GHz, 'Frequency')

    val = __model.fresnel_ellipse_radius(d1, d2, f)
    return prepare_output_array(val, type_output) * u.m 

def inverse_rain_attenuation(
            self, lat, lon, d, f, el, Ap, tau=45, R001=None):
        """ Implementation of 'inverse_rain_attenuation' method for
        recommendation ITU-P R.530-16. See documentation for function
        'ITUR530.inverse_rain_attenuation'
        """
        # Step 1: Obtain the rain rate R0.01 exceeded for 0.01% of the time
        # (with an integration time of 1 min).
        if R001 is None:
            R001 = rainfall_rate(lat, lon, 0.01).value

        # Step 2: Compute the specific attenuation, gammar (dB/km) for the
        # frequency, polarization and rain rate of interest using
        # Recommendation ITU-R P.838
        gammar = rain_specific_attenuation(R001, f, el, tau).value
        _, alpha = rain_specific_attenuation_coefficients(f, el, tau).value

        # Step 3: Compute the effective path length, deff, of the link by
        # multiplying the actual path length d by a distance factor r
        r = 1 / (0.477 * d ** 0.633 * R001 ** (0.073 * alpha) *
                 f**(0.123) - 10.579 * (1 - np.exp(-0.024 * d)))
        deff = np.minimum(r, 2.5) * d

        # Step 4: An estimate of the path attenuation exceeded for 0.01% of
        # the time is given by:
        A001 = gammar * deff

        # Step 5: The attenuation exceeded for other percentages of time p in
        # the range 0.001% to 1% may be deduced from the following power law
        C0 = np.where(f >= 10, 0.12 * 0.4 * (np.log10(f / 10)**0.8), 0.12)
        C1 = (0.07**C0) * (0.12**(1 - C0))
        C2 = 0.855 * C0 + 0.546 * (1 - C0)
        C3 = 0.139 * C0 + 0.043 * (1 - C0)

        def func_bisect(p):
            return A001 * C1 * p ** (- (C2 + C3 * np.log10(p))) - Ap

        return bisect(func_bisect, 0, 100) 

def surface_water_vapour_density(lat, lon, p, alt=None):
    """
    Method to compute the surface water vapour density along a path  at any
    desired location on the surface of the Earth.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    lon : number, sequence, or numpy.ndarray
        Longitudes of the receiver points
    p : number
        Percentage of time exceeded for p% of the average year
    alt : number, sequence, or numpy.ndarray
        Altitude of the receivers. If None, use the topographical altitude as
        described in recommendation ITU-R P.1511


    Returns
    -------
    rho: Quantity
       Surface water vapour density (g/m3)


    References
    ----------
    [1] Water vapour: surface density and total columnar content
    https://www.itu.int/rec/R-REC-P.836/en
    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    lon = prepare_input_array(lon)
    lon = np.mod(lon, 360)
    alt = prepare_input_array(alt)
    alt = prepare_quantity(alt, u.km, 'Altitude of the receivers')
    val = __model.surface_water_vapour_density(lat, lon, p, alt)
    return prepare_output_array(val, type_output) * u.g / u.m**3 

def rain_height(lat, lon):
    """
    A method to estimate the rain height for propagation prediction.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    lon : number, sequence, or numpy.ndarray
        Longitudes of the receiver points


    Returns
    -------
    rain_height: numpy.ndarray
        Rain height (km)


    References
    ----------
    [1] Rain height model for prediction methods:
    https://www.itu.int/rec/R-REC-P.839/en
    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    lon = prepare_input_array(lon)
    lon = np.mod(lon, 360)
    val = __model.rain_height(lat, lon)
    return prepare_output_array(val, type_output) * u.km 

def isoterm_0(lat, lon):
    """
    A method to estimate the zero isoterm height for propagation prediction.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    lon : number, sequence, or numpy.ndarray
        Longitudes of the receiver points


    Returns
    -------
    zero_isoterm: numpy.ndarray
        Zero isoterm height (km)


    References
    ----------
    [1] Rain height model for prediction methods:
    https://www.itu.int/rec/R-REC-P.839/en

    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    lon = prepare_input_array(lon)
    lon = np.mod(lon, 360)
    val = __model.isoterm_0(lat, lon)
    return prepare_output_array(val, type_output) * u.km 

def rain_height(self, lat_d, lon_d):
        """
        The rain height is computed as

        ..math:
            h_r = h_0 + 0.36 (km)
        """
        return self.isoterm_0(lat_d, lon_d) + 0.36 

def total_water_vapour_content(lat, lon, p, alt=None):
    """
    Method to compute the total water vapour content along a path  at any
    desired location on the surface of the Earth.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    lon : number, sequence, or numpy.ndarray
        Longitudes of the receiver points
    p : number
        Percentage of time exceeded for p% of the average year
    alt : number, sequence, or numpy.ndarray
        Altitude of the receivers. If None, use the topographical altitude as
        described in recommendation ITU-R P.1511


    Returns
    -------
    V: Quantity
       Total water vapour content (kg/m2)


    References
    ----------
    [1] Water vapour: surface density and total columnar content
    https://www.itu.int/rec/R-REC-P.836/en
    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    lon = prepare_input_array(lon)
    lon = np.mod(lon, 360)
    alt = prepare_input_array(alt)
    alt = prepare_quantity(alt, u.km, 'Altitude of the receivers')
    val = __model.total_water_vapour_content(lat, lon, p, alt)
    return prepare_output_array(val, type_output) * u.kg / u.m**2 

def rain_height(self, lat_d, lon_d):
        """
        The rain height is computed as

        ..math:
            h_r = h_0 + 0.36 (km)
        """
        return self.isoterm_0(lat_d, lon_d) + 0.36