<filename>where/models/delay/troposphere_radio.py<gh_stars>0
"""Calculate the delay caused by the troposphere for radio techniques

Description:
------------

Different models are available to determine the troposphere delay for radio techniques. The main models are GPT/GMF,
GPT2, GPT2w and VMF1 gridded.

References:
-----------
.. [1] IERS Conventions (2010) Software Library.
       http://maia.usno.navy.mil/conv2010/software.html



TODO:
    * Maybe it could be better to implement a "troposphere" class. For example with an object the information about
      used models can be retrieved, which is useful for example in routine :func:`zenith_wet_delay` (see TODO).
    * Test Dataset needed to assure better unit testing. So far all models are tested against output of the GIPSY
      program 'tropnominal', except GPT2w model.


"""
# External library imports
import numpy as np
import os

# Midgard imports
from midgard.dev import plugins
from midgard.math.unit import Unit

# Where imports
from where import apriori
from where.ext import iers_2010 as iers
from where.ext import gpt2w as ext_gpt2w
from where.lib import config
from where.lib import log

# Name of model
MODEL = __name__.split(".")[-1]

# Available troposphere models
MAPPING_FUNCTIONS = ["gmf", "gpt2", "gpt2w", "vmf1_gridded"]
METEOROLOGICAL_MODELS = ["vmf1_gridded", "gpt", "gpt2", "gpt2w", "site_pressure"]
ZENITH_WET_DELAY_MODELS = ["none", "askne", "davis", "saastamoinen", "vmf1_gridded"]
ZENITH_HYDROSTATIC_MODELS = ["saastamoinen"]
GRADIENT_MODELS = ["none", "apg"]

# Default relation between used mapping function model and zenith wet delay model
MAPPING_ZENITH_WET_RELATION = dict(gmf="none", gpt2="saastamoinen", gpt2w="askne", vmf1_gridded="vmf1_gridded")

# Default relation between used mapping function model and meteorological data
MAPPING_METEO_RELATION = dict(gmf="gpt", gpt2="gpt2", gpt2w="gpt2w", vmf1_gridded="vmf1_gridded")

# Cache for GPT2 model
_GPT2 = dict()

# Cache for GPT2W model
_GPT2W = dict()


@plugins.register
def troposphere_for_all_stations(dset):
    """Calculate tropospheric delay for all stations

    Depending on the used technique a Dataset can include one or more stations. In the case of VLBI a baseline between
    two stations are analysed. Here the relative tropospheric delay is needed between two stations and not the absolute
    one.  Therefore if a Dataset has two stations then the difference between the tropospheric delay is calculated by
    using ``multiplier``, which is '-1' for station 1 and '1' for station 2. For one or more than two stations in a
    Dataset the ``multiplier`` is equal to '1'.

    Args:
        dset (Dataset):     Model data.

    Returns:
        numpy.ndarray: Corrections in meters for each observation.
    """
    dset_out = np.zeros(dset.num_obs)
    for multiplier in dset.for_each_suffix("station"):
        dset_out += multiplier * troposphere(dset)

    return dset_out


def troposphere(dset):
    r"""Calculate tropospheric delay for a station

    The tropospheric delay :math:`\Delta T` is calculated using equation (9.12) in the IERS Conventions
    :cite:`iers2010`,

    .. math:: \Delta T = \Delta T_h^z \cdot m_h(e) + \Delta T_w \cdot m_w + \Delta T_{asym} ,

    with :math:`m_h` and :math:`m_w` denoting the hydrostatic and wet mapping functions, :math:`\Delta T_h^z` zenith
    hydrostatic delay, :math:`\Delta T_w^z` zenith wet delay and :math:`\Delta T_{asym}` asymmetric delay.

    Different models for the mapping functions, the zenith hydrostatic and wet delay, and the asymmetric delay
    determination are available and are defined in the configuration file.

    The different terms of the tropospheric delay are stored in the Dataset in a table called ``troposphere`` in fields
    prefixed by ``troposphere_``.

    Args:
        dset (Dataset):     Model data.

    Returns:
        numpy.ndarray:  Total tropospheric delay for each observation and one station in [m].
    """
    latitude, _, height = dset.site_pos.pos.llh.T
    pressure, temperature, e, tm, lambd = meteorological_data(dset)
    mh, mw = mapping_function(dset)
    mg, gn, ge = gradient_model(dset)
    zhd = zenith_hydrostatic_delay(pressure, latitude, height)
    zwd = zenith_wet_delay(dset, temperature, e, tm, lambd)

    # Line of sight delay
    dT = mh * zhd + mw * zwd + mg * (gn * np.cos(dset.site_pos.azimuth) + ge * np.sin(dset.site_pos.azimuth))

    terms_and_levels = dict(
        dT=("operational", "meter"),
        mh=("operational", "dimensionless"),
        zhd=("operational", "meter"),
        mw=("operational", "dimensionless"),
        zwd=("operational", "meter"),
        mg=("operational", "dimensionless"),
        gn=("operational", "meter"),
        ge=("operational", "meter"),
    )
    for term, (write_level, unit) in terms_and_levels.items():
        field = "troposphere_{}{}".format(term, (dset.default_field_suffix or ""))
        if field in dset.fields:
            dset[field][:] = locals()[term]
        else:
            # dset.add_float(field, table="troposphere", val=locals()[term], write_level=write_level)
            dset.add_float(field, val=locals()[term], write_level=write_level, unit=unit)

    # +DEBUG
    # if True:
    #    refdate_mjd = date.dt_to_mjd(dset.ref_date)
    #    delta_mjd = dset.time_unit.seconds / 86400
    #    for obs, sat in enumerate(dset.satellite):
    #        mjd = dset.time.utc.mjd[obs]
    #        zd = dset.site_pos.zenith_distance[obs]
    #        sat_id = dset.satellite_list[sat]
    #        print('MURKS: '+sat_id,mjd,zd,dT[obs],zhd[obs],zwd[obs],mh[obs],mw[obs],trop_asym[obs],pressure[obs])
    # -DEBUG

    return dT


def meteorological_data(dset):
    """Determine meteorological data (atmospheric pressure, ...) based on configuration file definition.

    Args:
        dset (Dataset): Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =====================================================================================
     Element       Unit         Description
    ============  ===========  =====================================================================================
     pressure      hPa          Atmospheric pressure value
     temperature   Celsius      Temperature values. Marked with 'None', if unknown.
     e             hPa          Water vapor pressure. Marked with 'None', if unknown.
     tm            K            Mean temperature of the water vapor. Marked with 'None', if unknown.
     lambd                      Water vapor decrease factor for each observation. Marked with 'None', if unknown.
    ============  ===========  =====================================================================================
    """
    pressure = None
    temperature = None
    e = None
    tm = None
    lambd = None

    model = config.tech.get("meteorological_data", section=MODEL, default="").str
    mapping_function = config.tech[MODEL].mapping_function.str

    # Use default meteorological models, if no model is defined in configuration file
    if not model:
        try:
            model = MAPPING_METEO_RELATION[mapping_function]
        except KeyError:
            log.fatal(
                f"Unknown mapping function {mapping_function}. "
                f"Available mapping functions are {', '.join(MAPPING_FUNCTIONS)}"
            )

    log.debug(f"Meteorological data model: {model}")

    if model == "vmf1_gridded":
        pressure = vmf1_gridded_pressure(dset)

    elif model == "gpt":
        pressure, temperature, _ = gpt(dset)

    elif model == "gpt2":
        pressure, temperature, _, e, _, _, _ = gpt2(dset)

    elif model == "gpt2w":
        pressure, temperature, _, tm, e, _, _, lambd, _ = gpt2w(dset)

    elif model == "site_pressure":
        pressure = site_pressure(dset)

    else:
        log.fatal(
            f"Unknown meteorological data model {model}. Available models are {', '.join(METEOROLOGICAL_MODELS)}"
        )

    return pressure, temperature, e, tm, lambd


def gradient_model(dset):
    """Calculates asymmetric delay based on gradient model given in configuration file

    Args:
        dset (Dataset):       Model data.

    Returns:
        numpy.ndarray:  Troposphere asymmetric delay in [m] for each observation
    """
    model = config.tech.get("gradients", section=MODEL, default="apg").str
    log.debug(f"Troposphere gradient model: {model}")

    # Note: Be aware, that the apg.f function uses c = 0.0031 instead of c = 0.0032.
    mg = 1 / (np.sin(dset.site_pos.elevation) * np.tan(dset.site_pos.elevation) + 0.0032)

    if model == "none":
        gn = ge = np.zeros(dset.num_obs)
    elif model == "apg":
        gn, ge = apg_gradient_model(dset)
    else:
        log.fatal(f"Unknown troposphere gradient model {model}. Available models are {', '.join(GRADIENT_MODELS)}")

    log.debug(f"Troposphere gradients North and East (average): {np.mean(gn)} {np.mean(ge)} [m]")

    return mg, gn, ge


def mapping_function(dset):
    """Calculates hydrostatic and wet mapping functions based on configuration file definition

    Args:
        dset (Dataset):    A Dataset containing model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =====================================================================================
     Element       Unit         Description
    ============  ===========  =====================================================================================
     mh                         Hydrostatic mapping function values
     mw                         Wet mapping function values
    ============  ===========  =====================================================================================
    """
    model = config.tech[MODEL].mapping_function.str
    log.debug(f"Troposphere mapping function: {model}")

    if model == "gmf":
        mh, mw = gmf_mapping_function(dset)

    elif model == "gpt2":
        _, _, _, _, mh, mw, _ = gpt2(dset)

    elif model == "gpt2w":
        _, _, _, _, _, mh, mw, _, _ = gpt2w(dset)

    elif model == "vmf1_gridded":
        mh, mw = vmf1_mapping_function(dset)

    else:
        log.fatal(
            f"Unknown troposphere mapping function {model}. "
            f"Available mapping functions are {', '.join(MAPPING_FUNCTIONS)}"
        )

    return mh, mw


def zenith_hydrostatic_delay(pressure, latitude, height):
    """Calculates zenith hydrostatic delay based on configuration file definition

    Args:
        pressure:    Array with atmospheric pressure for each observation in [hPa]
        latitude:    Array with geodetic latitude for each observation in [rad]
        height:      Array with orthometric height for each observation in [m]

    Returns:
        numpy.ndarray:         Array with zenith hydrostatic delay for each observation in [m]
    """
    model = config.tech.get("zenith_hydrostatic_delay", section=MODEL, default="saastamoinen").str
    log.debug(f"Troposphere zenith hydrostatic delay model: {model}")

    if model == "saastamoinen":
        zhd = saastamoinen_zenith_hydrostatic_delay(pressure, latitude, height)
    else:
        log.fatal(
            f"Unknown zenith hydrostatic troposphere delay model {model}. "
            f"Available models are {', '.join(ZENITH_HYDROSTATIC_MODELS)}"
        )

    log.debug(f"Troposphere zenith hydrostatic delay (average): {np.mean(zhd)} [m]")

    return zhd


def zenith_wet_delay(dset, temperature, e, tm, lambd):
    """Calculates zenith wet delay based on configuration file definition

    Args:
        dset (Dataset):                A Dataset containing model data
        temperature (numpy.ndarray):   Temperature for each observation in [Celsius]
        e (numpy.ndarray):             Water vapor pressure for each observation in [hPa]
        tm (numpy.ndarray):            Mean temperature of the water vapor for each observation in [K]
        lambd (numpy.ndarray):         Water vapor decrease factor for each observation


    Returns:
        numpy.ndarray:    Zenith wet delay values for each observation in [m].
    """
    model = config.tech.get("zenith_wet_delay", section=MODEL, default="").str
    mapping_function = config.tech[MODEL].mapping_function.str

    # Use default zenith wet delay models, if no model is defined in configuration file
    if not model:
        try:
            model = MAPPING_ZENITH_WET_RELATION[mapping_function]
        except KeyError:
            log.fatal(
                f"Unknown mapping function {mapping_function}. "
                f"Available mapping functions are {', '.join(MAPPING_FUNCTIONS)}"
            )
    log.debug(f"Troposphere zenith wet delay model: {model}")

    if model == "none":
        zwd = np.zeros(dset.num_obs)
    elif model == "askne":
        zwd = askne_zenith_wet_delay(e, tm, lambd)
        # TODO: log.fatal(f"Meteorological model {met_model!r} does not provide input parameters for using Askne and "
        #                 "Nordius zenith wet delay model. Use 'gpt2w' model")

    elif model == "davis":
        latitude, _, height = dset.site_pos.pos.llh.T
        zwd = davis_zenith_wet_delay(latitude, height, temperature, e)

    elif model == "saastamoinen":
        latitude, _, height = dset.site_pos.pos.llh.T
        zwd = saastamoinen_zenith_wet_delay(latitude, height, temperature, e)
        # TODO: log.fatal(f"Meteorological model {met_model} does not provide input parameters for using Saastamoinen "
        #                 "zenith wet delay model. Use 'gpt2' or 'gpt2w' model")

    elif model == "vmf1_gridded":
        zwd = vmf1_zenith_wet_delay(dset)
    else:
        log.fatal(
            f"Unknown zenith wet troposphere delay model {model}. "
            f"Available models are {', '.join(ZENITH_WET_DELAY_MODELS)}"
        )

    log.debug(f"Troposphere zenith wet delay (average): {np.mean(zwd)} [m]")

    return zwd


def apg_gradient_model(dset):
    """Calculates ECMWF gradient model based on Fortran routine 'apg.f'

    The tropospheric asymmetric delay and the horizontal delay gradients ``G_N`` and ``G_E`` are computed by using the
    ``apg.f`` function from the IERS software library [1]_. We use the asymmetric delay, which is based on Equation
    (9.12) described in IERS conventions 2010 in Section 9.2 :cite:`iers2010`.

    Args:
        dset (Dataset):       A Dataset containing model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =====================================================================================
     Element       Unit         Description
    ============  ===========  =====================================================================================
     gn            m            Horizontal delay gradient in the North direction
     ge            m            Horizontal delay gradient in the East direction
    ============  ===========  =====================================================================================
    """
    gn = ge = np.empty(dset.num_obs)
    lat, lon, _ = dset.site_pos.pos.llh.T
    az = dset.site_pos.azimuth
    el = dset.site_pos.elevation

    for obs in range(dset.num_obs):
        # Get horizontal gradients G_N and G_E
        _, gn[obs], ge[obs] = iers.apg(lat[obs], lon[obs], az[obs], el[obs])

    return gn * Unit.mm2m, ge * Unit.mm2m


def askne_zenith_wet_delay(e, tm, lambd):
    """Calculates zenith wet delay based on Askne and Nordius :cite:`askne1987`

    The Fortran routine 'asknewet.f' of the GPT2w library is used, which is based on Equation (18) of Askne and Nordius
    :cite:`askne1987`.

    Args:
        e (numpy.ndarray):      Water vapor pressure for each observation in [hPa]
        tm (numpy.ndarray):     Mean temperature of the water vapor for each observation in [K]
        lambd (numpy.ndarray):  Water vapor decrease factor for each observation

    Returns:
        numpy.ndarray:    Zenith wet delay for each observation in [m]
    """
    if not (len(e) == len(tm) == len(lambd)):
        log.fatal("Length of w, longitude and ellipsoidal height array is not equal.")

    num_obs = len(e)
    zwd = np.empty(num_obs)

    # Loop over all observations
    for obs in range(0, num_obs):
        zwd[obs] = ext_gpt2w.asknewet(e[obs], tm[obs], lambd[obs])

    return zwd


def davis_zenith_wet_delay(latitude, height, temperature, e):
    r"""Calculates zenith wet delay based on Saastamoinen/Davis model

    The total tropospheric delay for a given zenith distance :math:`z` is determined after Equation (19a) in
    Saastamoinen :cite:`saastamoinen1972`:

    .. math::
       \Delta T = 0.002277 \cdot \sec z \cdot (p + (1255/T + 0.05) \cdot e) - 1.16 \cdot \tan^2 z

    The zenith tropospheric delay is determined with :math:`z = 0`:

    .. math::
       \Delta T^z = \Delta T_h^z + \Delta T_w^z = zhd + zwd

    with the zenith hydrostatic delay :math:`zhd = 0.002277 \cdot p` and zenith wet delay :math:`zwd = 0.002277 \cdot
    (1255/T + 0.05) \cdot e`.

    A Fortran routine written by Davis corrects also for the gravity effect due to height :math:`H` and latitude
    :math:`\phi` (see http://acc.igs.org/tropo/wetsaas.f) and uses the constant :math:`0.0022768` instead of
    :math:`0.002277`, which leads to

    .. math::
       zwd = 0.0022768 \cdot (1255/T + 0.05) \cdot e / (1 - 0.00266 \cos 2 \phi - 0.00000028 H) .

    The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to
    100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant in this
    equation.  The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.

    Args:
        latitude (numpy.ndarray):     Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):       Orthometric height for each observation in [m]
        temperature (numpy.ndarray):  Temperature for each observation in [Celsius]
        e (numpy.ndarray):            Water vapor pressure for each observation in [hPa]

    Returns:
        numpy.ndarray:    Zenith wet delay for each observation in [m]
    """
    gravity_corr = saastamoinen_gravity_correction(latitude, height)

    # Zenith wet delay based on Eq. (19a) in Saastamoinen :cite:`saastamoinen1972` with additional gravity
    # correction
    zwd = 0.002_276_8 * (1255 / Unit.celsius_to_kelvin(temperature) + 0.05) * e / gravity_corr

    return zwd


def gmf_mapping_function(dset):
    """Calculates GMF hydrostatic and wet mapping functions

    Use the 'gmf.f' Fortran routine from the IERS software library to calculate the Global Mapping Function (see
    Section 9.2 in :cite:`iers2010`), which are described in Boehm et al. :cite:`boehm2006b`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     mh                         Hydrostatic mapping function coefficient ah
     mw                         Wet mapping function coefficient aw
    ============  ===========  =======================================================
    """
    mh = np.empty(dset.num_obs)
    mw = np.empty(dset.num_obs)
    lat, lon, height = dset.site_pos.pos.llh.T
    zd = dset.site_pos.zenith_distance
    mjd = dset.time.utc.mjd
    for obs in range(dset.num_obs):
        mh[obs], mw[obs] = iers.gmf(mjd[obs], lat[obs], lon[obs], height[obs], zd[obs])

    return mh, mw


def gpt(dset):
    """Calculates Global Pressure and Temperature (GPT)

    This subroutine determines atmospheric pressure and temperature globally based on spherical harmonics up to degree
    and order 9 (see Boehm et al. :cite:`boehm2007`) for a given latitude, longitude and ellipsoidal height.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     pressure      hPa          Pressure value
     temperature   Celsius      Temperature values
     geoid_undu    m            Geoid undulation (based on 9x9 EGM model)
    ============  ===========  =======================================================
    """
    pressure = np.empty(dset.num_obs)
    temperature = np.empty(dset.num_obs)
    geoid_undu = np.empty(dset.num_obs)

    # Note: For GPT2 and GPT2w linear interpolation is done between daily solutions. Here
    #      it does not seem to be necessary, because the performance of GPT Fortran routine
    #      is not so worse than for GPT2 and GPT2w.
    mjd = dset.time.utc.mjd
    lat, lon, height = dset.site_pos.pos.llh.T
    for obs in range(dset.num_obs):
        pressure[obs], temperature[obs], geoid_undu[obs] = iers.gpt(mjd[obs], lat[obs], lon[obs], height[obs])

    return pressure, temperature, geoid_undu


def gpt2(dset):
    """Calculates meteorological data and mapping function coefficients based on GPT2 model

    The GPT2 model is described in Lagler et al. :cite:`lagler2013`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     pressure      hPa          Pressure value
     temperature   Celsius      Temperature values
     dt            degree/km    Temperature lapse rate
     e             hPa          Water vapor pressure
     mh                         Hydrostatic mapping function coefficient ah
     mw                         Wet mapping function coefficient aw
     geoid_undu    m            Geoid undulation (based on 9x9 EGM model)
    ============  ===========  =======================================================
    """
    press = np.empty(dset.num_obs)
    temp = np.empty(dset.num_obs)
    dt = np.empty(dset.num_obs)
    e = np.empty(dset.num_obs)
    ah = np.empty(dset.num_obs)
    aw = np.empty(dset.num_obs)
    mh = np.empty(dset.num_obs)
    mw = np.empty(dset.num_obs)
    undu = np.empty(dset.num_obs)

    # Determine GPT2 values for each observation by interpolating between two unique
    # daily solutions
    mjd = dset.time.utc.mjd
    lat, lon, height = dset.site_pos.pos.llh.T
    zd = dset.site_pos.zenith_distance

    for obs in range(dset.num_obs):
        # Start 'gpt2.f' day-by-day in folder where 'gpt2_5.grd' is placed and carry out
        # linear interpolation
        press[obs], temp[obs], dt[obs], e[obs], ah[obs], aw[obs], undu[obs] = gpt2_wrapper(
            mjd[obs], [lat[obs]], [lon[obs]], [height[obs]]
        )

        # Determine mapping function values based on coefficients 'ah' and 'aw'
        mh[obs], mw[obs] = iers.vmf1_ht(ah[obs], aw[obs], mjd[obs], lat[obs], height[obs], zd[obs])

    return press, temp, dt, e, mh, mw, undu


def gpt2_wrapper(mjd, lat, lon, hell):
    """Calculates meteorological data and mapping function coefficients based on GPT2 model

    The functions calls the IERS library routine 'gpt2.f' (see Section 9.2 in :cite:`iers2010`). The Fortran routine
    ``gpt2.f`` reads the grid file ``gpt2_5.grd``, which should be available in the same folder, where the Fortran
    programs runs. Therefore we change the current directory to the IERS source directory, so that 'gpt2.f' can read
    the grid file.

    Due to performance reasons the GPT2 values are not determined for each observation. The call of the Fortran routine
    ``gpt2.f`` takes time, because the grid file ``gpt2_5.grd`` has to be read for each observation. Instead the GPT2
    values are calculated only once for each unique day (modified Julian date rounded to integer) and saved in the
    cache _GPT2. The final GPT2 values are computed by a linear interpolation of the daily determined GPT2 values. The
    difference between the use of routine ``gpt2.f`` for each observation and the linear interpolation between daily
    solution is on the submillimeter level and can therefore be neglected.

    Args:
        mjd (numpy.float64):  Modified Julian date.
        lat (list):           Array with latitude for each station in [rad].
        lon (list):           Array with longitude for each station in [rad].
        hell (list):          Array with height for each station in [m].

    Returns:
        numpy.ndarray:  Array with following entries:

    =======  ===========  =======================================================
     Index    Unit         Description
    =======  ===========  =======================================================
     [0]      hPa          Pressure value
     [1]      Celsius      Temperature values
     [2]      degree/km    Temperature lapse rate
     [3]      hPa          Water vapor pressure
     [4]                   Hydrostatic mapping function coefficient ah
     [5]                   Wet mapping function coefficient aw
     [6]      m            Geoid undulation (based on 9x9 EGM model)
    =======  ===========  =======================================================
    """
    nstat = len(lat)  # Number of stations
    it = 0  # Use of time variations (annual and semiannual terms)

    if not (len(lat) == len(lon) == len(hell)):
        log.fatal("Length of latitude, longitude and ellipsoidal height array is not equal.")

    # TODO: Case not handled if several stations are included in Dataset. Is that the case for VLBI?

    # Change directory so that gpt2.f can read the gpt2_5.grd-file in the IERS source directory
    current_dir = os.getcwd()
    os.chdir(config.files.path(iers.__name__))

    # Loop over all unique dates (rounded to integer value)
    for date in _rounded_dates(mjd):

        # Check if date is already included in cache
        if date not in _GPT2:
            _GPT2[date] = np.array(iers.gpt2(date, lat, lon, hell, nstat, it)).reshape(-1)
    os.chdir(current_dir)

    # Linear interpolation between two daily GPT2 solutions
    mjd_int, mjd_frac = divmod(mjd, 1)
    output = _GPT2[mjd_int] + mjd_frac * (_GPT2[mjd_int + 1] - _GPT2[mjd_int])

    return output


def gpt2w(dset):
    """Calculates meteorological data and mapping function coefficients based on GPT2w model

    The GPT2w model is described in Boehm et al. :cite:`boehm2015`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     pressure      hPa          Pressure value
     temperature   Celsius      Temperature values
     dt            degree/km    Temperature lapse rate
     tm            K            Mean temperature of the water vapor
     e             hPa          Water vapor pressure
     mh                         Hydrostatic mapping function coefficient ah
     mw                         Wet mapping function coefficient aw
     la                         Water vapor decrease factor
     geoid_undu    m            Geoid undulation (based on 9x9 EGM model)
    ============  ===========  =======================================================
    """
    press = np.empty(dset.num_obs)
    temp = np.empty(dset.num_obs)
    dt = np.empty(dset.num_obs)
    tm = np.empty(dset.num_obs)
    e = np.empty(dset.num_obs)
    ah = np.empty(dset.num_obs)
    aw = np.empty(dset.num_obs)
    mh = np.empty(dset.num_obs)
    mw = np.empty(dset.num_obs)
    la = np.empty(dset.num_obs)
    undu = np.empty(dset.num_obs)

    mjd = dset.time.utc.mjd
    lat, lon, height = dset.site_pos.pos.llh.T
    zd = dset.site_pos.zenith_distance

    # Determine GPT2W values for each observation by interpolating between two unique
    # daily solutions
    for obs in range(dset.num_obs):
        # Start 'gpt2.f' day-by-day in folder where 'gpt2_5.grd' is placed and carry out
        # linear interpolation
        (press[obs], temp[obs], dt[obs], tm[obs], e[obs], ah[obs], aw[obs], la[obs], undu[obs]) = gpt2w_wrapper(
            mjd[obs], [lat[obs]], [lon[obs]], [height[obs]]
        )

        # Determine mapping function values based on coefficients 'ah' and 'aw'
        mh[obs], mw[obs] = iers.vmf1_ht(ah[obs], aw[obs], mjd[obs], lat[obs], height[obs], zd[obs])

    return press, temp, dt, tm, e, mh, mw, la, undu


def gpt2w_wrapper(mjd, lat, lon, hell):
    """Calculates meteorological data and mapping function coefficients based on GPT2w model

    The functions calls the GPT2w library routine ``gpt2w_1w.f`` (see
    http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w). The Fortran routine ``gpt2w_1w.f`` reads the grid file
    ``gpt2_1wA.grd``, which should be available in the same folder, where the Fortran programs runs. Therefore we
    change the current directory to the GPT2w source directory, so that ``gpt2w_1w.f`` can read the grid file.

    Due to performance reasons the GPT2w values are not determined for each observation. The call of the Fortran
    routine ``gpt2w_1w.f`` takes time, because the grid file ``gpt2_1wA.grd`` has to be read for each
    observation. Instead the GPT2w values are calculated only once for each unique day (modified Julian date rounded to
    integer) and saved in the cache _GPT2W. The final GPT2W values are computed by a linear interpolation of the daily
    determined GPT2W values.  The difference between the use of routine ``gpt2w_1w.f`` for each observation and the
    linear interpolation between daily solution is on the submillimeter level and can therefore be neglected.

    Args:
        mjd (numpy.float64):  Modified Julian date.
        lat (list):           Array with latitude for each station in [rad].
        lon (list):           Array with longitude for each station in [rad].
        hell (list):          Array with height for each station in [m].

    Returns:
        numpy.ndarray:  Array with following entries:

    =======  ===========  =======================================================
     Index    Unit         Description
    =======  ===========  =======================================================
     [0]      hPa          Pressure value
     [1]      Celsius      Temperature values
     [2]      degree/km    Temperature lapse rate
     [3]      K            Mean temperature of the water vapor
     [4]      hPa          Water vapor pressure
     [5]                   Hydrostatic mapping function coefficient ah
     [6]                   Wet mapping function coefficient aw
     [7]                   Water vapor decrease factor
     [8]      m            Geoid undulation (based on 9x9 EGM model)
    =======  ===========  =======================================================
    """
    nstat = len(lat)  # Number of stations
    it = 0  # Use of time variations (annual and semiannual terms)

    # TODO: Case not handled if several stations are included in Dateset. Is that the case for VLBI?

    if not (len(lat) == len(lon) == len(hell)):
        log.fatal("Length of latitude, longitute and ellipsoidal height array is not equal.")

    # Change directory so that gpt2w.f can read the gpt2_5.grd-file in the GPT2w source directory
    current_dir = os.getcwd()
    os.chdir(config.files.path(ext_gpt2w.__name__))

    # Loop over all unique dates (rounded to integer value)
    for date in _rounded_dates(mjd):

        # Check if date is already included in cache
        if date not in _GPT2W:
            _GPT2W[date] = np.array(ext_gpt2w.gpt2_1w(date, lat, lon, hell, nstat, it)).reshape(-1)
    os.chdir(current_dir)

    # Linear interpolation between two daily GPT2W solutions
    mjd_int, mjd_frac = divmod(mjd, 1)
    output = _GPT2W[mjd_int] + mjd_frac * (_GPT2W[mjd_int + 1] - _GPT2W[mjd_int])

    return output


def pressure_zhd(zhd, latitude, height):
    """Calculates atmospheric pressure based on zenith hydrostatic troposphere delay

    Args:
        zhd (numpy.ndarray):         Zenith hydrostatic delay for each observation in [m]
        latitude (numpy.ndarray):    Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):      Orthometric height for each observation in [m]

    Returns:
        numpy.ndarray:    Atmospheric pressure for each observation in [hPa]
    """
    model = config.tech.get("zenith_hydrostatic_delay", section=MODEL, default="saastamoinen").str

    if model == "saastamoinen":
        pressure = saastamoinen_pressure(zhd, latitude, height)
    else:
        log.fatal(
            f"Zenith troposphere delay definition {model!r} is not correct in configuration file. "
            "It should be 'saastamoinen'"
        )

    return pressure


def pressure_height_correction(pressure_ref, height_ref, height):
    r"""Atmospheric pressure is corrected from reference height to given height

    The atmospheric pressure height correction is based on a standard atmosphere model of Berg. See equation (1) in
    Boehm et al. :cite:`boehm2007` or equation (4) in Kouba :cite:`kouba2007`.

    The method uses equation (4) in Kouba :cite:`kouba2007` to determine the pressure for the reference height
    :math:`p_r` and for the station height :math:`p_s`:

    .. math::

       p_r &= 1013.25 (1 - 0.0000226 h_r)^{5.225} ,

       p_s &= 1013.25 (1 - 0.0000226 h_s)^{5.225} .

    The pressure at the station height is determined by setting these two equations
    in following relation:

    .. math:: p_s = (1 - 0.0000226h_r)^{5.225} / (1 - 0.0000226h_s)^{5.225} * p_r .

    TODO:
        The implementation of the pressure height correction is not completely described in Kouba :cite:`kouba2007`.  I
        found different solutions from Trimble, <NAME> and Bernese, whereby Bernese uses a similiar approach for the
        computation (see BERN52/LIB/FOR/VMF1ELL.f90). The use of the Trimble solution leads to differences of several
        millimeters for the zenith hydrostatic troposphere delay compared to our solution. Our solution compared to Per
        Helge's solution shows differences on the millimeter level.  What is the correct solution?

        Trimble solution::

            pressure = pressure_ref + 1013.25 * ((1 - 0.0000226 * height)**5.225 - (1 - 0.0000226 * height_ref)**5.225)

        Per Helge's solution::

            pressure = pressure_ref * (1 - 0.0000226 * (height-height_ref))**5.225

    Args:
        pressure_ref (numpy.ndarray):    Atmospheric pressure given for reference height for each observation in [hPa]
        height_ref (numpy.ndarray):      Reference height for each observation in [m]
        height (numpy.ndarray):          Height for each observation in [m]

    Returns:
        numpy.ndarray:      Pressure for given height for each observation in [hPa]
    """
    return ((1 - 0.000_022_6 * height) ** 5.225) / ((1 - 0.000_022_6 * height_ref) ** 5.225) * pressure_ref


def saastamoinen_gravity_correction(latitude, height):
    """Calculates gravity correction due to height and latitude

    The gravity correction due to height and latitude changes is shown in Davis et al. :cite:`davis1985` (see Equation
    (A13)) and is based on Saastamoinen :cite:`saastamoinen1972`.

    Args:
        latitude (numpy.ndarray):      Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):        Orthometric height for each observation in [m]

    Returns:
        numpy.ndarray: Gravity correction for each observation in [m]
    """
    # Gravity correction based on Eq. (A13) in Davis et al. :cite:`davis1985`
    return 1 - 0.00266 * np.cos(2 * latitude) - 0.000_000_28 * height


def saastamoinen_pressure(zhd, latitude, height):
    """Calculates atmospheric pressure based on Saastamoinen model

    The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to
    100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant by using
    equation (9.11). The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.

    Args:
        zhd (numpy.ndarray):         Zenith hydrostatic delay for each observation in [m]
        latitude (numpy.ndarray):    Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):      Orthometric height for each observation in [m]

    Returns:
        numpy.ndarray:    Atmospheric pressure for each observation in [hPa]

    Example:
        >>> saastamoinen_pressure(2.2762945349647778, 0.83776, 200)
        1000.0
    """
    return zhd * saastamoinen_gravity_correction(latitude, height) / 0.002_276_8


def saastamoinen_zenith_hydrostatic_delay(pressure, latitude, height):
    r"""Calculates zenith hydrostatic delay with Saastamoinen model

    The zenith hydrostatic delay :math:`zhd` is calculated using equation (9.11) from Saastamoinen in the IERS
    Conventions :cite:`iers2010`:

    .. math::
         \Delta T_h^z = zhd = (0.0022768 \pm 0.0000005) P_0 / f_s(\phi, H)

    with the atmospheric pressure :math:`P_0` and the correction :math:`f_s(\phi, H)`
    (see equation (9.4) in :cite:`iers2010`):

    .. math::
       f_s(\phi, H) = 1 - 0.00266 \cos 2 \phi - 0.00000028 H ,

    where :math:`\phi` is the geodetic latitude of the station and :math:`H` is the geodetic height of the station.

    The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to
    100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant by using
    equation (9.11). The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.

    Args:
        pressure (numpy.ndarray):      Atmospheric pressure for each observation in [hPa]
        latitude (numpy.ndarray):      Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):        Orthometric height for each observation in [m]

    Returns:
        numpy.ndarray:    Zenith hydrostatic delay for each observation in [m]

    Example:
        >>> saastamoinen_zenith_hydrostatic_delay(1000, 0.83776, 200)
        2.2762945349647778

    TODO:
        Example is not compatible to http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w/FORTRAN/saasthyd.f
        example (zhd = 2.2695 m).
    """
    return 0.002_276_8 * pressure / saastamoinen_gravity_correction(latitude, height)


def saastamoinen_zenith_wet_delay(latitude, height, temperature, e):
    r"""Calculates zenith wet delay based Saastamoinen model

    The total tropospheric delay for a given zenith distance :math:`z` is determined after Equation (19a) in
    Saastamoinen :cite:`saastamoinen1972`:

    .. math::
       \Delta T = 0.002277 \cdot \sec z \cdot (p + (1255/T + 0.05) \cdot e) - 1.16 \cdot \tan^2 z

    The zenith tropospheric delay is determined with :math:`z = 0`:

    .. math::
       \Delta T^z = \Delta T_h^z + \Delta T_w^z = zhd + zwd

    with the zenith hydrostatic delay :math:`zhd = 0.002277 \cdot p` and zenith wet delay
    :math:`zwd = 0.002277 \cdot (1255/T + 0.05) \cdot e`.

    Args:
        latitude (numpy.ndarray):     Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):       Orthometric height for each observation in [m]
        temperature (numpy.ndarray):  Temperature for each observation in [Celsius]
        e (numpy.ndarray):            Water vapor pressure for each observation in [hPa]

    Returns:
        numpy.ndarray:     Zenith wet delay for each observation in [m]
    """
    # Zenith wet delay based on Eq. (19a) in Saastamoinen :cite:`saastamoinen1972`
    return 0.002_276_8 * (1255 / Unit.celsius_to_kelvin(temperature) + 0.05) * e


def site_pressure(dset):
    """Get atmospheric pressure from local site measurements

    If local atmospheric pressure measurements on a site are not available an alternative model given in configuration
    file is used to determine atmospheric pressure.

    TODO:
        So far only gridded VMF1 model is used, if local pressure data are not available.  Which alternative is used,
        should be decided via the configuration file. How to check after an alternative model in configuration file?
        model_list = config.tech[MODEL].meteorological_data.list???

    Args:
        dset (Dataset): A Dataset containing model data.

    Returns:
        numpy.ndarray: Atmospheric pressure for each observation in [hPa]
    """
    pressure = np.zeros(dset.num_obs)

    i_given = np.zeros(dset.num_obs, dtype=bool)
    if "pressure" + (dset.default_field_suffix or "") in dset.fields:
        i_given[np.logical_not(np.isnan(dset.pressure))] = True
        pressure[i_given] = dset.pressure[i_given]

    i_missing = np.logical_not(i_given)
    if i_missing.any():
        pressure[i_missing] = vmf1_gridded_pressure(dset)[i_missing]

    return pressure


def vmf1_gridded_pressure(dset):
    """Calulates VMF1 gridded atmospheric pressure

    Equation (9.11) in IERS conventions 2010 :cite:`iers2010` is converted to 'pressure' and the VMF1 gridded
    atmospheric pressure can be determined with help of the gridded zenith hydrostatic delays from VMF1. The gridded
    VMF1 pressure values needs to be rescaled to actual station height. The method is described with Eq. (4) in Kouba
    :cite:`kouba2007`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        numpy.ndarray:  Atmospheric pressure for each observation in [hPa].
    """
    # Get gridded VMF1 data
    vmf1 = apriori.get("vmf1", time=dset.time)

    lat, lon, height = dset.site_pos.pos.llh.T
    grid_zhd = vmf1["zh"](dset.time, lon, lat)  # Interpolation in time and space in VMF1 grid
    grid_height = vmf1["ell"](lon, lat, grid=False)
    grid_pressure = pressure_zhd(grid_zhd, lat, grid_height)

    # Rescale gridded pressure to station pressure
    pressure = pressure_height_correction(grid_pressure, grid_height, height)

    return pressure


def vmf1_mapping_function(dset):
    """Calculates VMF1 hydrostatic and wet mapping functions

    This routine determines the VMF1 (Vienna Mapping Functions 1) described in Boehm et al. :cite:`boehm2006a` and
    Kouba :cite:`kouba2007` and uses the 'vmf1_ht.f' from the IERS software library :cite:`iers2010`. The mapping
    function coefficients 'ah' and 'aw' are taken from VMF1 gridded files.

    TODO: Multiplication with 1.e-8 should maybe be handled when reading gridded VMF1-files and not here.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     mh                         Hydrostatic mapping function coefficient ah
     mw                         Wet mapping function coefficient aw
    ============  ===========  =======================================================
    """
    # Get gridded VMF1 data
    vmf1 = apriori.get("vmf1", time=dset.time)

    mh = np.empty(dset.num_obs)
    mw = np.empty(dset.num_obs)

    lat, lon, height = dset.site_pos.pos.llh.T
    zd = dset.site_pos.zenith_distance

    for obs in range(dset.num_obs):

        # Interpolation in time and space in VMF1 grid
        # t = dset.time[obs]
        # TODO vlbi t_2
        # if dset._default_field_suffix == '_2':
        # baseline_gcrs = dset.site_pos_2.gcrs - dset.site_pos_1.gcrs
        # delta_t = (dset.src_dir.unit_vector[:, None, :] @ baseline_gcrs[:, :, None])[:, 0, 0] / constant.c

        # t -= timedelta(seconds=delta_t[obs])
        mh[obs], mw[obs] = iers.vmf1_ht(
            vmf1["ah"](dset.time[obs], lon[obs], lat[obs]) * 1e-8,
            vmf1["aw"](dset.time[obs], lon[obs], lat[obs]) * 1e-8,
            dset.time.utc.mjd_int[obs],
            lat[obs],
            height[obs],
            zd[obs],
        )

    return mh, mw


def vmf1_zenith_wet_delay(dset):
    """Calculates zenith wet delay based on gridded zenith wet delays from VMF1

    Uses gridded zenith wet delays from VMF1, which are rescaled from the gridded height to actual station height by
    using Equation(5) described in Kouba :cite:`kouba2007`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        numpy.ndarray:     Zenith wet delay for each observation in [m]
    """
    # Get gridded VMF1 data
    vmf1 = apriori.get("vmf1", time=dset.time)

    lat, lon, height = dset.site_pos.pos.llh.T
    grid_zwd = vmf1["zw"](dset.time, lon, lat)  # Interpolation in time and space in VMF1 grid
    grid_height = vmf1["ell"](lon, lat, grid=False)

    # Zenith Wet delay. Eq. (5) in Kouba :cite:`kouba2007`
    zwd = grid_zwd * np.exp(-(height - grid_height) / 2000)

    return zwd


def _rounded_dates(datetimes):
    """Calculate days including the given datetimes

    Args:
        datetimes (numpy.float64):   np.array with datetimes as floats

    Returns:
        numpy.ndarray:  Unique days representing the datetimes.

    Example:
        >>> _rounded_dates(np.array([54001.2, 54001.4, 54001.8, 54002, 54005.5, 54006.]))
        array([ 54001.,  54002.,  54003.,  54005.,  54006.,  54007.])
    """
    return np.unique(np.hstack((np.floor(datetimes), np.floor(datetimes) + 1)))

<filename>where/models/delay/troposphere_radio.py<gh_stars>0
"""Calculate the delay caused by the troposphere for radio techniques

Description:
------------

Different models are available to determine the troposphere delay for radio techniques. The main models are GPT/GMF,
GPT2, GPT2w and VMF1 gridded.

References:
-----------
.. [1] IERS Conventions (2010) Software Library.
       http://maia.usno.navy.mil/conv2010/software.html



TODO:
    * Maybe it could be better to implement a "troposphere" class. For example with an object the information about
      used models can be retrieved, which is useful for example in routine :func:`zenith_wet_delay` (see TODO).
    * Test Dataset needed to assure better unit testing. So far all models are tested against output of the GIPSY
      program 'tropnominal', except GPT2w model.


"""
# External library imports
import numpy as np
import os

# Midgard imports
from midgard.dev import plugins
from midgard.math.unit import Unit

# Where imports
from where import apriori
from where.ext import iers_2010 as iers
from where.ext import gpt2w as ext_gpt2w
from where.lib import config
from where.lib import log

# Name of model
MODEL = __name__.split(".")[-1]

# Available troposphere models
MAPPING_FUNCTIONS = ["gmf", "gpt2", "gpt2w", "vmf1_gridded"]
METEOROLOGICAL_MODELS = ["vmf1_gridded", "gpt", "gpt2", "gpt2w", "site_pressure"]
ZENITH_WET_DELAY_MODELS = ["none", "askne", "davis", "saastamoinen", "vmf1_gridded"]
ZENITH_HYDROSTATIC_MODELS = ["saastamoinen"]
GRADIENT_MODELS = ["none", "apg"]

# Default relation between used mapping function model and zenith wet delay model
MAPPING_ZENITH_WET_RELATION = dict(gmf="none", gpt2="saastamoinen", gpt2w="askne", vmf1_gridded="vmf1_gridded")

# Default relation between used mapping function model and meteorological data
MAPPING_METEO_RELATION = dict(gmf="gpt", gpt2="gpt2", gpt2w="gpt2w", vmf1_gridded="vmf1_gridded")

# Cache for GPT2 model
_GPT2 = dict()

# Cache for GPT2W model
_GPT2W = dict()


@plugins.register
def troposphere_for_all_stations(dset):
    """Calculate tropospheric delay for all stations

    Depending on the used technique a Dataset can include one or more stations. In the case of VLBI a baseline between
    two stations are analysed. Here the relative tropospheric delay is needed between two stations and not the absolute
    one.  Therefore if a Dataset has two stations then the difference between the tropospheric delay is calculated by
    using ``multiplier``, which is '-1' for station 1 and '1' for station 2. For one or more than two stations in a
    Dataset the ``multiplier`` is equal to '1'.

    Args:
        dset (Dataset):     Model data.

    Returns:
        numpy.ndarray: Corrections in meters for each observation.
    """
    dset_out = np.zeros(dset.num_obs)
    for multiplier in dset.for_each_suffix("station"):
        dset_out += multiplier * troposphere(dset)

    return dset_out


def troposphere(dset):
    r"""Calculate tropospheric delay for a station

    The tropospheric delay :math:`\Delta T` is calculated using equation (9.12) in the IERS Conventions
    :cite:`iers2010`,

    .. math:: \Delta T = \Delta T_h^z \cdot m_h(e) + \Delta T_w \cdot m_w + \Delta T_{asym} ,

    with :math:`m_h` and :math:`m_w` denoting the hydrostatic and wet mapping functions, :math:`\Delta T_h^z` zenith
    hydrostatic delay, :math:`\Delta T_w^z` zenith wet delay and :math:`\Delta T_{asym}` asymmetric delay.

    Different models for the mapping functions, the zenith hydrostatic and wet delay, and the asymmetric delay
    determination are available and are defined in the configuration file.

    The different terms of the tropospheric delay are stored in the Dataset in a table called ``troposphere`` in fields
    prefixed by ``troposphere_``.

    Args:
        dset (Dataset):     Model data.

    Returns:
        numpy.ndarray:  Total tropospheric delay for each observation and one station in [m].
    """
    latitude, _, height = dset.site_pos.pos.llh.T
    pressure, temperature, e, tm, lambd = meteorological_data(dset)
    mh, mw = mapping_function(dset)
    mg, gn, ge = gradient_model(dset)
    zhd = zenith_hydrostatic_delay(pressure, latitude, height)
    zwd = zenith_wet_delay(dset, temperature, e, tm, lambd)

    # Line of sight delay
    dT = mh * zhd + mw * zwd + mg * (gn * np.cos(dset.site_pos.azimuth) + ge * np.sin(dset.site_pos.azimuth))

    terms_and_levels = dict(
        dT=("operational", "meter"),
        mh=("operational", "dimensionless"),
        zhd=("operational", "meter"),
        mw=("operational", "dimensionless"),
        zwd=("operational", "meter"),
        mg=("operational", "dimensionless"),
        gn=("operational", "meter"),
        ge=("operational", "meter"),
    )
    for term, (write_level, unit) in terms_and_levels.items():
        field = "troposphere_{}{}".format(term, (dset.default_field_suffix or ""))
        if field in dset.fields:
            dset[field][:] = locals()[term]
        else:
            # dset.add_float(field, table="troposphere", val=locals()[term], write_level=write_level)
            dset.add_float(field, val=locals()[term], write_level=write_level, unit=unit)

    # +DEBUG
    # if True:
    #    refdate_mjd = date.dt_to_mjd(dset.ref_date)
    #    delta_mjd = dset.time_unit.seconds / 86400
    #    for obs, sat in enumerate(dset.satellite):
    #        mjd = dset.time.utc.mjd[obs]
    #        zd = dset.site_pos.zenith_distance[obs]
    #        sat_id = dset.satellite_list[sat]
    #        print('MURKS: '+sat_id,mjd,zd,dT[obs],zhd[obs],zwd[obs],mh[obs],mw[obs],trop_asym[obs],pressure[obs])
    # -DEBUG

    return dT


def meteorological_data(dset):
    """Determine meteorological data (atmospheric pressure, ...) based on configuration file definition.

    Args:
        dset (Dataset): Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =====================================================================================
     Element       Unit         Description
    ============  ===========  =====================================================================================
     pressure      hPa          Atmospheric pressure value
     temperature   Celsius      Temperature values. Marked with 'None', if unknown.
     e             hPa          Water vapor pressure. Marked with 'None', if unknown.
     tm            K            Mean temperature of the water vapor. Marked with 'None', if unknown.
     lambd                      Water vapor decrease factor for each observation. Marked with 'None', if unknown.
    ============  ===========  =====================================================================================
    """
    pressure = None
    temperature = None
    e = None
    tm = None
    lambd = None

    model = config.tech.get("meteorological_data", section=MODEL, default="").str
    mapping_function = config.tech[MODEL].mapping_function.str

    # Use default meteorological models, if no model is defined in configuration file
    if not model:
        try:
            model = MAPPING_METEO_RELATION[mapping_function]
        except KeyError:
            log.fatal(
                f"Unknown mapping function {mapping_function}. "
                f"Available mapping functions are {', '.join(MAPPING_FUNCTIONS)}"
            )

    log.debug(f"Meteorological data model: {model}")

    if model == "vmf1_gridded":
        pressure = vmf1_gridded_pressure(dset)

    elif model == "gpt":
        pressure, temperature, _ = gpt(dset)

    elif model == "gpt2":
        pressure, temperature, _, e, _, _, _ = gpt2(dset)

    elif model == "gpt2w":
        pressure, temperature, _, tm, e, _, _, lambd, _ = gpt2w(dset)

    elif model == "site_pressure":
        pressure = site_pressure(dset)

    else:
        log.fatal(
            f"Unknown meteorological data model {model}. Available models are {', '.join(METEOROLOGICAL_MODELS)}"
        )

    return pressure, temperature, e, tm, lambd


def gradient_model(dset):
    """Calculates asymmetric delay based on gradient model given in configuration file

    Args:
        dset (Dataset):       Model data.

    Returns:
        numpy.ndarray:  Troposphere asymmetric delay in [m] for each observation
    """
    model = config.tech.get("gradients", section=MODEL, default="apg").str
    log.debug(f"Troposphere gradient model: {model}")

    # Note: Be aware, that the apg.f function uses c = 0.0031 instead of c = 0.0032.
    mg = 1 / (np.sin(dset.site_pos.elevation) * np.tan(dset.site_pos.elevation) + 0.0032)

    if model == "none":
        gn = ge = np.zeros(dset.num_obs)
    elif model == "apg":
        gn, ge = apg_gradient_model(dset)
    else:
        log.fatal(f"Unknown troposphere gradient model {model}. Available models are {', '.join(GRADIENT_MODELS)}")

    log.debug(f"Troposphere gradients North and East (average): {np.mean(gn)} {np.mean(ge)} [m]")

    return mg, gn, ge


def mapping_function(dset):
    """Calculates hydrostatic and wet mapping functions based on configuration file definition

    Args:
        dset (Dataset):    A Dataset containing model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =====================================================================================
     Element       Unit         Description
    ============  ===========  =====================================================================================
     mh                         Hydrostatic mapping function values
     mw                         Wet mapping function values
    ============  ===========  =====================================================================================
    """
    model = config.tech[MODEL].mapping_function.str
    log.debug(f"Troposphere mapping function: {model}")

    if model == "gmf":
        mh, mw = gmf_mapping_function(dset)

    elif model == "gpt2":
        _, _, _, _, mh, mw, _ = gpt2(dset)

    elif model == "gpt2w":
        _, _, _, _, _, mh, mw, _, _ = gpt2w(dset)

    elif model == "vmf1_gridded":
        mh, mw = vmf1_mapping_function(dset)

    else:
        log.fatal(
            f"Unknown troposphere mapping function {model}. "
            f"Available mapping functions are {', '.join(MAPPING_FUNCTIONS)}"
        )

    return mh, mw


def zenith_hydrostatic_delay(pressure, latitude, height):
    """Calculates zenith hydrostatic delay based on configuration file definition

    Args:
        pressure:    Array with atmospheric pressure for each observation in [hPa]
        latitude:    Array with geodetic latitude for each observation in [rad]
        height:      Array with orthometric height for each observation in [m]

    Returns:
        numpy.ndarray:         Array with zenith hydrostatic delay for each observation in [m]
    """
    model = config.tech.get("zenith_hydrostatic_delay", section=MODEL, default="saastamoinen").str
    log.debug(f"Troposphere zenith hydrostatic delay model: {model}")

    if model == "saastamoinen":
        zhd = saastamoinen_zenith_hydrostatic_delay(pressure, latitude, height)
    else:
        log.fatal(
            f"Unknown zenith hydrostatic troposphere delay model {model}. "
            f"Available models are {', '.join(ZENITH_HYDROSTATIC_MODELS)}"
        )

    log.debug(f"Troposphere zenith hydrostatic delay (average): {np.mean(zhd)} [m]")

    return zhd


def zenith_wet_delay(dset, temperature, e, tm, lambd):
    """Calculates zenith wet delay based on configuration file definition

    Args:
        dset (Dataset):                A Dataset containing model data
        temperature (numpy.ndarray):   Temperature for each observation in [Celsius]
        e (numpy.ndarray):             Water vapor pressure for each observation in [hPa]
        tm (numpy.ndarray):            Mean temperature of the water vapor for each observation in [K]
        lambd (numpy.ndarray):         Water vapor decrease factor for each observation


    Returns:
        numpy.ndarray:    Zenith wet delay values for each observation in [m].
    """
    model = config.tech.get("zenith_wet_delay", section=MODEL, default="").str
    mapping_function = config.tech[MODEL].mapping_function.str

    # Use default zenith wet delay models, if no model is defined in configuration file
    if not model:
        try:
            model = MAPPING_ZENITH_WET_RELATION[mapping_function]
        except KeyError:
            log.fatal(
                f"Unknown mapping function {mapping_function}. "
                f"Available mapping functions are {', '.join(MAPPING_FUNCTIONS)}"
            )
    log.debug(f"Troposphere zenith wet delay model: {model}")

    if model == "none":
        zwd = np.zeros(dset.num_obs)
    elif model == "askne":
        zwd = askne_zenith_wet_delay(e, tm, lambd)
        # TODO: log.fatal(f"Meteorological model {met_model!r} does not provide input parameters for using Askne and "
        #                 "Nordius zenith wet delay model. Use 'gpt2w' model")

    elif model == "davis":
        latitude, _, height = dset.site_pos.pos.llh.T
        zwd = davis_zenith_wet_delay(latitude, height, temperature, e)

    elif model == "saastamoinen":
        latitude, _, height = dset.site_pos.pos.llh.T
        zwd = saastamoinen_zenith_wet_delay(latitude, height, temperature, e)
        # TODO: log.fatal(f"Meteorological model {met_model} does not provide input parameters for using Saastamoinen "
        #                 "zenith wet delay model. Use 'gpt2' or 'gpt2w' model")

    elif model == "vmf1_gridded":
        zwd = vmf1_zenith_wet_delay(dset)
    else:
        log.fatal(
            f"Unknown zenith wet troposphere delay model {model}. "
            f"Available models are {', '.join(ZENITH_WET_DELAY_MODELS)}"
        )

    log.debug(f"Troposphere zenith wet delay (average): {np.mean(zwd)} [m]")

    return zwd


def apg_gradient_model(dset):
    """Calculates ECMWF gradient model based on Fortran routine 'apg.f'

    The tropospheric asymmetric delay and the horizontal delay gradients ``G_N`` and ``G_E`` are computed by using the
    ``apg.f`` function from the IERS software library [1]_. We use the asymmetric delay, which is based on Equation
    (9.12) described in IERS conventions 2010 in Section 9.2 :cite:`iers2010`.

    Args:
        dset (Dataset):       A Dataset containing model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =====================================================================================
     Element       Unit         Description
    ============  ===========  =====================================================================================
     gn            m            Horizontal delay gradient in the North direction
     ge            m            Horizontal delay gradient in the East direction
    ============  ===========  =====================================================================================
    """
    gn = ge = np.empty(dset.num_obs)
    lat, lon, _ = dset.site_pos.pos.llh.T
    az = dset.site_pos.azimuth
    el = dset.site_pos.elevation

    for obs in range(dset.num_obs):
        # Get horizontal gradients G_N and G_E
        _, gn[obs], ge[obs] = iers.apg(lat[obs], lon[obs], az[obs], el[obs])

    return gn * Unit.mm2m, ge * Unit.mm2m


def askne_zenith_wet_delay(e, tm, lambd):
    """Calculates zenith wet delay based on Askne and Nordius :cite:`askne1987`

    The Fortran routine 'asknewet.f' of the GPT2w library is used, which is based on Equation (18) of Askne and Nordius
    :cite:`askne1987`.

    Args:
        e (numpy.ndarray):      Water vapor pressure for each observation in [hPa]
        tm (numpy.ndarray):     Mean temperature of the water vapor for each observation in [K]
        lambd (numpy.ndarray):  Water vapor decrease factor for each observation

    Returns:
        numpy.ndarray:    Zenith wet delay for each observation in [m]
    """
    if not (len(e) == len(tm) == len(lambd)):
        log.fatal("Length of w, longitude and ellipsoidal height array is not equal.")

    num_obs = len(e)
    zwd = np.empty(num_obs)

    # Loop over all observations
    for obs in range(0, num_obs):
        zwd[obs] = ext_gpt2w.asknewet(e[obs], tm[obs], lambd[obs])

    return zwd


def davis_zenith_wet_delay(latitude, height, temperature, e):
    r"""Calculates zenith wet delay based on Saastamoinen/Davis model

    The total tropospheric delay for a given zenith distance :math:`z` is determined after Equation (19a) in
    Saastamoinen :cite:`saastamoinen1972`:

    .. math::
       \Delta T = 0.002277 \cdot \sec z \cdot (p + (1255/T + 0.05) \cdot e) - 1.16 \cdot \tan^2 z

    The zenith tropospheric delay is determined with :math:`z = 0`:

    .. math::
       \Delta T^z = \Delta T_h^z + \Delta T_w^z = zhd + zwd

    with the zenith hydrostatic delay :math:`zhd = 0.002277 \cdot p` and zenith wet delay :math:`zwd = 0.002277 \cdot
    (1255/T + 0.05) \cdot e`.

    A Fortran routine written by Davis corrects also for the gravity effect due to height :math:`H` and latitude
    :math:`\phi` (see http://acc.igs.org/tropo/wetsaas.f) and uses the constant :math:`0.0022768` instead of
    :math:`0.002277`, which leads to

    .. math::
       zwd = 0.0022768 \cdot (1255/T + 0.05) \cdot e / (1 - 0.00266 \cos 2 \phi - 0.00000028 H) .

    The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to
    100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant in this
    equation.  The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.

    Args:
        latitude (numpy.ndarray):     Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):       Orthometric height for each observation in [m]
        temperature (numpy.ndarray):  Temperature for each observation in [Celsius]
        e (numpy.ndarray):            Water vapor pressure for each observation in [hPa]

    Returns:
        numpy.ndarray:    Zenith wet delay for each observation in [m]
    """
    gravity_corr = saastamoinen_gravity_correction(latitude, height)

    # Zenith wet delay based on Eq. (19a) in Saastamoinen :cite:`saastamoinen1972` with additional gravity
    # correction
    zwd = 0.002_276_8 * (1255 / Unit.celsius_to_kelvin(temperature) + 0.05) * e / gravity_corr

    return zwd


def gmf_mapping_function(dset):
    """Calculates GMF hydrostatic and wet mapping functions

    Use the 'gmf.f' Fortran routine from the IERS software library to calculate the Global Mapping Function (see
    Section 9.2 in :cite:`iers2010`), which are described in Boehm et al. :cite:`boehm2006b`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     mh                         Hydrostatic mapping function coefficient ah
     mw                         Wet mapping function coefficient aw
    ============  ===========  =======================================================
    """
    mh = np.empty(dset.num_obs)
    mw = np.empty(dset.num_obs)
    lat, lon, height = dset.site_pos.pos.llh.T
    zd = dset.site_pos.zenith_distance
    mjd = dset.time.utc.mjd
    for obs in range(dset.num_obs):
        mh[obs], mw[obs] = iers.gmf(mjd[obs], lat[obs], lon[obs], height[obs], zd[obs])

    return mh, mw


def gpt(dset):
    """Calculates Global Pressure and Temperature (GPT)

    This subroutine determines atmospheric pressure and temperature globally based on spherical harmonics up to degree
    and order 9 (see Boehm et al. :cite:`boehm2007`) for a given latitude, longitude and ellipsoidal height.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     pressure      hPa          Pressure value
     temperature   Celsius      Temperature values
     geoid_undu    m            Geoid undulation (based on 9x9 EGM model)
    ============  ===========  =======================================================
    """
    pressure = np.empty(dset.num_obs)
    temperature = np.empty(dset.num_obs)
    geoid_undu = np.empty(dset.num_obs)

    # Note: For GPT2 and GPT2w linear interpolation is done between daily solutions. Here
    #      it does not seem to be necessary, because the performance of GPT Fortran routine
    #      is not so worse than for GPT2 and GPT2w.
    mjd = dset.time.utc.mjd
    lat, lon, height = dset.site_pos.pos.llh.T
    for obs in range(dset.num_obs):
        pressure[obs], temperature[obs], geoid_undu[obs] = iers.gpt(mjd[obs], lat[obs], lon[obs], height[obs])

    return pressure, temperature, geoid_undu


def gpt2(dset):
    """Calculates meteorological data and mapping function coefficients based on GPT2 model

    The GPT2 model is described in Lagler et al. :cite:`lagler2013`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     pressure      hPa          Pressure value
     temperature   Celsius      Temperature values
     dt            degree/km    Temperature lapse rate
     e             hPa          Water vapor pressure
     mh                         Hydrostatic mapping function coefficient ah
     mw                         Wet mapping function coefficient aw
     geoid_undu    m            Geoid undulation (based on 9x9 EGM model)
    ============  ===========  =======================================================
    """
    press = np.empty(dset.num_obs)
    temp = np.empty(dset.num_obs)
    dt = np.empty(dset.num_obs)
    e = np.empty(dset.num_obs)
    ah = np.empty(dset.num_obs)
    aw = np.empty(dset.num_obs)
    mh = np.empty(dset.num_obs)
    mw = np.empty(dset.num_obs)
    undu = np.empty(dset.num_obs)

    # Determine GPT2 values for each observation by interpolating between two unique
    # daily solutions
    mjd = dset.time.utc.mjd
    lat, lon, height = dset.site_pos.pos.llh.T
    zd = dset.site_pos.zenith_distance

    for obs in range(dset.num_obs):
        # Start 'gpt2.f' day-by-day in folder where 'gpt2_5.grd' is placed and carry out
        # linear interpolation
        press[obs], temp[obs], dt[obs], e[obs], ah[obs], aw[obs], undu[obs] = gpt2_wrapper(
            mjd[obs], [lat[obs]], [lon[obs]], [height[obs]]
        )

        # Determine mapping function values based on coefficients 'ah' and 'aw'
        mh[obs], mw[obs] = iers.vmf1_ht(ah[obs], aw[obs], mjd[obs], lat[obs], height[obs], zd[obs])

    return press, temp, dt, e, mh, mw, undu


def gpt2_wrapper(mjd, lat, lon, hell):
    """Calculates meteorological data and mapping function coefficients based on GPT2 model

    The functions calls the IERS library routine 'gpt2.f' (see Section 9.2 in :cite:`iers2010`). The Fortran routine
    ``gpt2.f`` reads the grid file ``gpt2_5.grd``, which should be available in the same folder, where the Fortran
    programs runs. Therefore we change the current directory to the IERS source directory, so that 'gpt2.f' can read
    the grid file.

    Due to performance reasons the GPT2 values are not determined for each observation. The call of the Fortran routine
    ``gpt2.f`` takes time, because the grid file ``gpt2_5.grd`` has to be read for each observation. Instead the GPT2
    values are calculated only once for each unique day (modified Julian date rounded to integer) and saved in the
    cache _GPT2. The final GPT2 values are computed by a linear interpolation of the daily determined GPT2 values. The
    difference between the use of routine ``gpt2.f`` for each observation and the linear interpolation between daily
    solution is on the submillimeter level and can therefore be neglected.

    Args:
        mjd (numpy.float64):  Modified Julian date.
        lat (list):           Array with latitude for each station in [rad].
        lon (list):           Array with longitude for each station in [rad].
        hell (list):          Array with height for each station in [m].

    Returns:
        numpy.ndarray:  Array with following entries:

    =======  ===========  =======================================================
     Index    Unit         Description
    =======  ===========  =======================================================
     [0]      hPa          Pressure value
     [1]      Celsius      Temperature values
     [2]      degree/km    Temperature lapse rate
     [3]      hPa          Water vapor pressure
     [4]                   Hydrostatic mapping function coefficient ah
     [5]                   Wet mapping function coefficient aw
     [6]      m            Geoid undulation (based on 9x9 EGM model)
    =======  ===========  =======================================================
    """
    nstat = len(lat)  # Number of stations
    it = 0  # Use of time variations (annual and semiannual terms)

    if not (len(lat) == len(lon) == len(hell)):
        log.fatal("Length of latitude, longitude and ellipsoidal height array is not equal.")

    # TODO: Case not handled if several stations are included in Dataset. Is that the case for VLBI?

    # Change directory so that gpt2.f can read the gpt2_5.grd-file in the IERS source directory
    current_dir = os.getcwd()
    os.chdir(config.files.path(iers.__name__))

    # Loop over all unique dates (rounded to integer value)
    for date in _rounded_dates(mjd):

        # Check if date is already included in cache
        if date not in _GPT2:
            _GPT2[date] = np.array(iers.gpt2(date, lat, lon, hell, nstat, it)).reshape(-1)
    os.chdir(current_dir)

    # Linear interpolation between two daily GPT2 solutions
    mjd_int, mjd_frac = divmod(mjd, 1)
    output = _GPT2[mjd_int] + mjd_frac * (_GPT2[mjd_int + 1] - _GPT2[mjd_int])

    return output


def gpt2w(dset):
    """Calculates meteorological data and mapping function coefficients based on GPT2w model

    The GPT2w model is described in Boehm et al. :cite:`boehm2015`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     pressure      hPa          Pressure value
     temperature   Celsius      Temperature values
     dt            degree/km    Temperature lapse rate
     tm            K            Mean temperature of the water vapor
     e             hPa          Water vapor pressure
     mh                         Hydrostatic mapping function coefficient ah
     mw                         Wet mapping function coefficient aw
     la                         Water vapor decrease factor
     geoid_undu    m            Geoid undulation (based on 9x9 EGM model)
    ============  ===========  =======================================================
    """
    press = np.empty(dset.num_obs)
    temp = np.empty(dset.num_obs)
    dt = np.empty(dset.num_obs)
    tm = np.empty(dset.num_obs)
    e = np.empty(dset.num_obs)
    ah = np.empty(dset.num_obs)
    aw = np.empty(dset.num_obs)
    mh = np.empty(dset.num_obs)
    mw = np.empty(dset.num_obs)
    la = np.empty(dset.num_obs)
    undu = np.empty(dset.num_obs)

    mjd = dset.time.utc.mjd
    lat, lon, height = dset.site_pos.pos.llh.T
    zd = dset.site_pos.zenith_distance

    # Determine GPT2W values for each observation by interpolating between two unique
    # daily solutions
    for obs in range(dset.num_obs):
        # Start 'gpt2.f' day-by-day in folder where 'gpt2_5.grd' is placed and carry out
        # linear interpolation
        (press[obs], temp[obs], dt[obs], tm[obs], e[obs], ah[obs], aw[obs], la[obs], undu[obs]) = gpt2w_wrapper(
            mjd[obs], [lat[obs]], [lon[obs]], [height[obs]]
        )

        # Determine mapping function values based on coefficients 'ah' and 'aw'
        mh[obs], mw[obs] = iers.vmf1_ht(ah[obs], aw[obs], mjd[obs], lat[obs], height[obs], zd[obs])

    return press, temp, dt, tm, e, mh, mw, la, undu


def gpt2w_wrapper(mjd, lat, lon, hell):
    """Calculates meteorological data and mapping function coefficients based on GPT2w model

    The functions calls the GPT2w library routine ``gpt2w_1w.f`` (see
    http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w). The Fortran routine ``gpt2w_1w.f`` reads the grid file
    ``gpt2_1wA.grd``, which should be available in the same folder, where the Fortran programs runs. Therefore we
    change the current directory to the GPT2w source directory, so that ``gpt2w_1w.f`` can read the grid file.

    Due to performance reasons the GPT2w values are not determined for each observation. The call of the Fortran
    routine ``gpt2w_1w.f`` takes time, because the grid file ``gpt2_1wA.grd`` has to be read for each
    observation. Instead the GPT2w values are calculated only once for each unique day (modified Julian date rounded to
    integer) and saved in the cache _GPT2W. The final GPT2W values are computed by a linear interpolation of the daily
    determined GPT2W values.  The difference between the use of routine ``gpt2w_1w.f`` for each observation and the
    linear interpolation between daily solution is on the submillimeter level and can therefore be neglected.

    Args:
        mjd (numpy.float64):  Modified Julian date.
        lat (list):           Array with latitude for each station in [rad].
        lon (list):           Array with longitude for each station in [rad].
        hell (list):          Array with height for each station in [m].

    Returns:
        numpy.ndarray:  Array with following entries:

    =======  ===========  =======================================================
     Index    Unit         Description
    =======  ===========  =======================================================
     [0]      hPa          Pressure value
     [1]      Celsius      Temperature values
     [2]      degree/km    Temperature lapse rate
     [3]      K            Mean temperature of the water vapor
     [4]      hPa          Water vapor pressure
     [5]                   Hydrostatic mapping function coefficient ah
     [6]                   Wet mapping function coefficient aw
     [7]                   Water vapor decrease factor
     [8]      m            Geoid undulation (based on 9x9 EGM model)
    =======  ===========  =======================================================
    """
    nstat = len(lat)  # Number of stations
    it = 0  # Use of time variations (annual and semiannual terms)

    # TODO: Case not handled if several stations are included in Dateset. Is that the case for VLBI?

    if not (len(lat) == len(lon) == len(hell)):
        log.fatal("Length of latitude, longitute and ellipsoidal height array is not equal.")

    # Change directory so that gpt2w.f can read the gpt2_5.grd-file in the GPT2w source directory
    current_dir = os.getcwd()
    os.chdir(config.files.path(ext_gpt2w.__name__))

    # Loop over all unique dates (rounded to integer value)
    for date in _rounded_dates(mjd):

        # Check if date is already included in cache
        if date not in _GPT2W:
            _GPT2W[date] = np.array(ext_gpt2w.gpt2_1w(date, lat, lon, hell, nstat, it)).reshape(-1)
    os.chdir(current_dir)

    # Linear interpolation between two daily GPT2W solutions
    mjd_int, mjd_frac = divmod(mjd, 1)
    output = _GPT2W[mjd_int] + mjd_frac * (_GPT2W[mjd_int + 1] - _GPT2W[mjd_int])

    return output


def pressure_zhd(zhd, latitude, height):
    """Calculates atmospheric pressure based on zenith hydrostatic troposphere delay

    Args:
        zhd (numpy.ndarray):         Zenith hydrostatic delay for each observation in [m]
        latitude (numpy.ndarray):    Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):      Orthometric height for each observation in [m]

    Returns:
        numpy.ndarray:    Atmospheric pressure for each observation in [hPa]
    """
    model = config.tech.get("zenith_hydrostatic_delay", section=MODEL, default="saastamoinen").str

    if model == "saastamoinen":
        pressure = saastamoinen_pressure(zhd, latitude, height)
    else:
        log.fatal(
            f"Zenith troposphere delay definition {model!r} is not correct in configuration file. "
            "It should be 'saastamoinen'"
        )

    return pressure


def pressure_height_correction(pressure_ref, height_ref, height):
    r"""Atmospheric pressure is corrected from reference height to given height

    The atmospheric pressure height correction is based on a standard atmosphere model of Berg. See equation (1) in
    Boehm et al. :cite:`boehm2007` or equation (4) in Kouba :cite:`kouba2007`.

    The method uses equation (4) in Kouba :cite:`kouba2007` to determine the pressure for the reference height
    :math:`p_r` and for the station height :math:`p_s`:

    .. math::

       p_r &= 1013.25 (1 - 0.0000226 h_r)^{5.225} ,

       p_s &= 1013.25 (1 - 0.0000226 h_s)^{5.225} .

    The pressure at the station height is determined by setting these two equations
    in following relation:

    .. math:: p_s = (1 - 0.0000226h_r)^{5.225} / (1 - 0.0000226h_s)^{5.225} * p_r .

    TODO:
        The implementation of the pressure height correction is not completely described in Kouba :cite:`kouba2007`.  I
        found different solutions from Trimble, <NAME> and Bernese, whereby Bernese uses a similiar approach for the
        computation (see BERN52/LIB/FOR/VMF1ELL.f90). The use of the Trimble solution leads to differences of several
        millimeters for the zenith hydrostatic troposphere delay compared to our solution. Our solution compared to Per
        Helge's solution shows differences on the millimeter level.  What is the correct solution?

        Trimble solution::

            pressure = pressure_ref + 1013.25 * ((1 - 0.0000226 * height)**5.225 - (1 - 0.0000226 * height_ref)**5.225)

        Per Helge's solution::

            pressure = pressure_ref * (1 - 0.0000226 * (height-height_ref))**5.225

    Args:
        pressure_ref (numpy.ndarray):    Atmospheric pressure given for reference height for each observation in [hPa]
        height_ref (numpy.ndarray):      Reference height for each observation in [m]
        height (numpy.ndarray):          Height for each observation in [m]

    Returns:
        numpy.ndarray:      Pressure for given height for each observation in [hPa]
    """
    return ((1 - 0.000_022_6 * height) ** 5.225) / ((1 - 0.000_022_6 * height_ref) ** 5.225) * pressure_ref


def saastamoinen_gravity_correction(latitude, height):
    """Calculates gravity correction due to height and latitude

    The gravity correction due to height and latitude changes is shown in Davis et al. :cite:`davis1985` (see Equation
    (A13)) and is based on Saastamoinen :cite:`saastamoinen1972`.

    Args:
        latitude (numpy.ndarray):      Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):        Orthometric height for each observation in [m]

    Returns:
        numpy.ndarray: Gravity correction for each observation in [m]
    """
    # Gravity correction based on Eq. (A13) in Davis et al. :cite:`davis1985`
    return 1 - 0.00266 * np.cos(2 * latitude) - 0.000_000_28 * height


def saastamoinen_pressure(zhd, latitude, height):
    """Calculates atmospheric pressure based on Saastamoinen model

    The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to
    100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant by using
    equation (9.11). The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.

    Args:
        zhd (numpy.ndarray):         Zenith hydrostatic delay for each observation in [m]
        latitude (numpy.ndarray):    Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):      Orthometric height for each observation in [m]

    Returns:
        numpy.ndarray:    Atmospheric pressure for each observation in [hPa]

    Example:
        >>> saastamoinen_pressure(2.2762945349647778, 0.83776, 200)
        1000.0
    """
    return zhd * saastamoinen_gravity_correction(latitude, height) / 0.002_276_8


def saastamoinen_zenith_hydrostatic_delay(pressure, latitude, height):
    r"""Calculates zenith hydrostatic delay with Saastamoinen model

    The zenith hydrostatic delay :math:`zhd` is calculated using equation (9.11) from Saastamoinen in the IERS
    Conventions :cite:`iers2010`:

    .. math::
         \Delta T_h^z = zhd = (0.0022768 \pm 0.0000005) P_0 / f_s(\phi, H)

    with the atmospheric pressure :math:`P_0` and the correction :math:`f_s(\phi, H)`
    (see equation (9.4) in :cite:`iers2010`):

    .. math::
       f_s(\phi, H) = 1 - 0.00266 \cos 2 \phi - 0.00000028 H ,

    where :math:`\phi` is the geodetic latitude of the station and :math:`H` is the geodetic height of the station.

    The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to
    100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant by using
    equation (9.11). The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.

    Args:
        pressure (numpy.ndarray):      Atmospheric pressure for each observation in [hPa]
        latitude (numpy.ndarray):      Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):        Orthometric height for each observation in [m]

    Returns:
        numpy.ndarray:    Zenith hydrostatic delay for each observation in [m]

    Example:
        >>> saastamoinen_zenith_hydrostatic_delay(1000, 0.83776, 200)
        2.2762945349647778

    TODO:
        Example is not compatible to http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w/FORTRAN/saasthyd.f
        example (zhd = 2.2695 m).
    """
    return 0.002_276_8 * pressure / saastamoinen_gravity_correction(latitude, height)


def saastamoinen_zenith_wet_delay(latitude, height, temperature, e):
    r"""Calculates zenith wet delay based Saastamoinen model

    The total tropospheric delay for a given zenith distance :math:`z` is determined after Equation (19a) in
    Saastamoinen :cite:`saastamoinen1972`:

    .. math::
       \Delta T = 0.002277 \cdot \sec z \cdot (p + (1255/T + 0.05) \cdot e) - 1.16 \cdot \tan^2 z

    The zenith tropospheric delay is determined with :math:`z = 0`:

    .. math::
       \Delta T^z = \Delta T_h^z + \Delta T_w^z = zhd + zwd

    with the zenith hydrostatic delay :math:`zhd = 0.002277 \cdot p` and zenith wet delay
    :math:`zwd = 0.002277 \cdot (1255/T + 0.05) \cdot e`.

    Args:
        latitude (numpy.ndarray):     Geodetic latitude for each observation in [rad]
        height (numpy.ndarray):       Orthometric height for each observation in [m]
        temperature (numpy.ndarray):  Temperature for each observation in [Celsius]
        e (numpy.ndarray):            Water vapor pressure for each observation in [hPa]

    Returns:
        numpy.ndarray:     Zenith wet delay for each observation in [m]
    """
    # Zenith wet delay based on Eq. (19a) in Saastamoinen :cite:`saastamoinen1972`
    return 0.002_276_8 * (1255 / Unit.celsius_to_kelvin(temperature) + 0.05) * e


def site_pressure(dset):
    """Get atmospheric pressure from local site measurements

    If local atmospheric pressure measurements on a site are not available an alternative model given in configuration
    file is used to determine atmospheric pressure.

    TODO:
        So far only gridded VMF1 model is used, if local pressure data are not available.  Which alternative is used,
        should be decided via the configuration file. How to check after an alternative model in configuration file?
        model_list = config.tech[MODEL].meteorological_data.list???

    Args:
        dset (Dataset): A Dataset containing model data.

    Returns:
        numpy.ndarray: Atmospheric pressure for each observation in [hPa]
    """
    pressure = np.zeros(dset.num_obs)

    i_given = np.zeros(dset.num_obs, dtype=bool)
    if "pressure" + (dset.default_field_suffix or "") in dset.fields:
        i_given[np.logical_not(np.isnan(dset.pressure))] = True
        pressure[i_given] = dset.pressure[i_given]

    i_missing = np.logical_not(i_given)
    if i_missing.any():
        pressure[i_missing] = vmf1_gridded_pressure(dset)[i_missing]

    return pressure


def vmf1_gridded_pressure(dset):
    """Calulates VMF1 gridded atmospheric pressure

    Equation (9.11) in IERS conventions 2010 :cite:`iers2010` is converted to 'pressure' and the VMF1 gridded
    atmospheric pressure can be determined with help of the gridded zenith hydrostatic delays from VMF1. The gridded
    VMF1 pressure values needs to be rescaled to actual station height. The method is described with Eq. (4) in Kouba
    :cite:`kouba2007`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        numpy.ndarray:  Atmospheric pressure for each observation in [hPa].
    """
    # Get gridded VMF1 data
    vmf1 = apriori.get("vmf1", time=dset.time)

    lat, lon, height = dset.site_pos.pos.llh.T
    grid_zhd = vmf1["zh"](dset.time, lon, lat)  # Interpolation in time and space in VMF1 grid
    grid_height = vmf1["ell"](lon, lat, grid=False)
    grid_pressure = pressure_zhd(grid_zhd, lat, grid_height)

    # Rescale gridded pressure to station pressure
    pressure = pressure_height_correction(grid_pressure, grid_height, height)

    return pressure


def vmf1_mapping_function(dset):
    """Calculates VMF1 hydrostatic and wet mapping functions

    This routine determines the VMF1 (Vienna Mapping Functions 1) described in Boehm et al. :cite:`boehm2006a` and
    Kouba :cite:`kouba2007` and uses the 'vmf1_ht.f' from the IERS software library :cite:`iers2010`. The mapping
    function coefficients 'ah' and 'aw' are taken from VMF1 gridded files.

    TODO: Multiplication with 1.e-8 should maybe be handled when reading gridded VMF1-files and not here.

    Args:
        dset (Dataset):    Model data.

    Returns:
        tuple of Numpy Arrays: Includes the following elements, each with entries for each observation

    ============  ===========  =======================================================
     Element       Unit         Description
    ============  ===========  =======================================================
     mh                         Hydrostatic mapping function coefficient ah
     mw                         Wet mapping function coefficient aw
    ============  ===========  =======================================================
    """
    # Get gridded VMF1 data
    vmf1 = apriori.get("vmf1", time=dset.time)

    mh = np.empty(dset.num_obs)
    mw = np.empty(dset.num_obs)

    lat, lon, height = dset.site_pos.pos.llh.T
    zd = dset.site_pos.zenith_distance

    for obs in range(dset.num_obs):

        # Interpolation in time and space in VMF1 grid
        # t = dset.time[obs]
        # TODO vlbi t_2
        # if dset._default_field_suffix == '_2':
        # baseline_gcrs = dset.site_pos_2.gcrs - dset.site_pos_1.gcrs
        # delta_t = (dset.src_dir.unit_vector[:, None, :] @ baseline_gcrs[:, :, None])[:, 0, 0] / constant.c

        # t -= timedelta(seconds=delta_t[obs])
        mh[obs], mw[obs] = iers.vmf1_ht(
            vmf1["ah"](dset.time[obs], lon[obs], lat[obs]) * 1e-8,
            vmf1["aw"](dset.time[obs], lon[obs], lat[obs]) * 1e-8,
            dset.time.utc.mjd_int[obs],
            lat[obs],
            height[obs],
            zd[obs],
        )

    return mh, mw


def vmf1_zenith_wet_delay(dset):
    """Calculates zenith wet delay based on gridded zenith wet delays from VMF1

    Uses gridded zenith wet delays from VMF1, which are rescaled from the gridded height to actual station height by
    using Equation(5) described in Kouba :cite:`kouba2007`.

    Args:
        dset (Dataset):    Model data.

    Returns:
        numpy.ndarray:     Zenith wet delay for each observation in [m]
    """
    # Get gridded VMF1 data
    vmf1 = apriori.get("vmf1", time=dset.time)

    lat, lon, height = dset.site_pos.pos.llh.T
    grid_zwd = vmf1["zw"](dset.time, lon, lat)  # Interpolation in time and space in VMF1 grid
    grid_height = vmf1["ell"](lon, lat, grid=False)

    # Zenith Wet delay. Eq. (5) in Kouba :cite:`kouba2007`
    zwd = grid_zwd * np.exp(-(height - grid_height) / 2000)

    return zwd


def _rounded_dates(datetimes):
    """Calculate days including the given datetimes

    Args:
        datetimes (numpy.float64):   np.array with datetimes as floats

    Returns:
        numpy.ndarray:  Unique days representing the datetimes.

    Example:
        >>> _rounded_dates(np.array([54001.2, 54001.4, 54001.8, 54002, 54005.5, 54006.]))
        array([ 54001.,  54002.,  54003.,  54005.,  54006.,  54007.])
    """
    return np.unique(np.hstack((np.floor(datetimes), np.floor(datetimes) + 1)))

Calculate the delay caused by the troposphere for radio techniques  Description: ------------  Different models are available to determine the troposphere delay for radio techniques. The main models are GPT/GMF, GPT2, GPT2w and VMF1 gridded.  References: ----------- .. [1] IERS Conventions (2010) Software Library.        http://maia.usno.navy.mil/conv2010/software.html    TODO:     * Maybe it could be better to implement a "troposphere" class. For example with an object the information about       used models can be retrieved, which is useful for example in routine :func:`zenith_wet_delay` (see TODO).     * Test Dataset needed to assure better unit testing. So far all models are tested against output of the GIPSY       program 'tropnominal', except GPT2w model. # External library imports # Midgard imports # Where imports # Name of model # Available troposphere models # Default relation between used mapping function model and zenith wet delay model # Default relation between used mapping function model and meteorological data # Cache for GPT2 model # Cache for GPT2W model Calculate tropospheric delay for all stations      Depending on the used technique a Dataset can include one or more stations. In the case of VLBI a baseline between     two stations are analysed. Here the relative tropospheric delay is needed between two stations and not the absolute     one.  Therefore if a Dataset has two stations then the difference between the tropospheric delay is calculated by     using ``multiplier``, which is '-1' for station 1 and '1' for station 2. For one or more than two stations in a     Dataset the ``multiplier`` is equal to '1'.      Args:         dset (Dataset):     Model data.      Returns:         numpy.ndarray: Corrections in meters for each observation. Calculate tropospheric delay for a station      The tropospheric delay :math:`\Delta T` is calculated using equation (9.12) in the IERS Conventions     :cite:`iers2010`,      .. math:: \Delta T = \Delta T_h^z \cdot m_h(e) + \Delta T_w \cdot m_w + \Delta T_{asym} ,      with :math:`m_h` and :math:`m_w` denoting the hydrostatic and wet mapping functions, :math:`\Delta T_h^z` zenith     hydrostatic delay, :math:`\Delta T_w^z` zenith wet delay and :math:`\Delta T_{asym}` asymmetric delay.      Different models for the mapping functions, the zenith hydrostatic and wet delay, and the asymmetric delay     determination are available and are defined in the configuration file.      The different terms of the tropospheric delay are stored in the Dataset in a table called ``troposphere`` in fields     prefixed by ``troposphere_``.      Args:         dset (Dataset):     Model data.      Returns:         numpy.ndarray:  Total tropospheric delay for each observation and one station in [m]. # Line of sight delay # dset.add_float(field, table="troposphere", val=locals()[term], write_level=write_level) # +DEBUG # if True: #    refdate_mjd = date.dt_to_mjd(dset.ref_date) #    delta_mjd = dset.time_unit.seconds / 86400 #    for obs, sat in enumerate(dset.satellite): #        mjd = dset.time.utc.mjd[obs] #        zd = dset.site_pos.zenith_distance[obs] #        sat_id = dset.satellite_list[sat] #        print('MURKS: '+sat_id,mjd,zd,dT[obs],zhd[obs],zwd[obs],mh[obs],mw[obs],trop_asym[obs],pressure[obs]) # -DEBUG Determine meteorological data (atmospheric pressure, ...) based on configuration file definition.      Args:         dset (Dataset): Model data.      Returns:         tuple of Numpy Arrays: Includes the following elements, each with entries for each observation      ============  ===========  =====================================================================================      Element       Unit         Description     ============  ===========  =====================================================================================      pressure      hPa          Atmospheric pressure value      temperature   Celsius      Temperature values. Marked with 'None', if unknown.      e             hPa          Water vapor pressure. Marked with 'None', if unknown.      tm            K            Mean temperature of the water vapor. Marked with 'None', if unknown.      lambd                      Water vapor decrease factor for each observation. Marked with 'None', if unknown.     ============  ===========  ===================================================================================== # Use default meteorological models, if no model is defined in configuration file Calculates asymmetric delay based on gradient model given in configuration file      Args:         dset (Dataset):       Model data.      Returns:         numpy.ndarray:  Troposphere asymmetric delay in [m] for each observation # Note: Be aware, that the apg.f function uses c = 0.0031 instead of c = 0.0032. Calculates hydrostatic and wet mapping functions based on configuration file definition      Args:         dset (Dataset):    A Dataset containing model data.      Returns:         tuple of Numpy Arrays: Includes the following elements, each with entries for each observation      ============  ===========  =====================================================================================      Element       Unit         Description     ============  ===========  =====================================================================================      mh                         Hydrostatic mapping function values      mw                         Wet mapping function values     ============  ===========  ===================================================================================== Calculates zenith hydrostatic delay based on configuration file definition      Args:         pressure:    Array with atmospheric pressure for each observation in [hPa]         latitude:    Array with geodetic latitude for each observation in [rad]         height:      Array with orthometric height for each observation in [m]      Returns:         numpy.ndarray:         Array with zenith hydrostatic delay for each observation in [m] Calculates zenith wet delay based on configuration file definition      Args:         dset (Dataset):                A Dataset containing model data         temperature (numpy.ndarray):   Temperature for each observation in [Celsius]         e (numpy.ndarray):             Water vapor pressure for each observation in [hPa]         tm (numpy.ndarray):            Mean temperature of the water vapor for each observation in [K]         lambd (numpy.ndarray):         Water vapor decrease factor for each observation       Returns:         numpy.ndarray:    Zenith wet delay values for each observation in [m]. # Use default zenith wet delay models, if no model is defined in configuration file # TODO: log.fatal(f"Meteorological model {met_model!r} does not provide input parameters for using Askne and " #                 "Nordius zenith wet delay model. Use 'gpt2w' model") # TODO: log.fatal(f"Meteorological model {met_model} does not provide input parameters for using Saastamoinen " #                 "zenith wet delay model. Use 'gpt2' or 'gpt2w' model") Calculates ECMWF gradient model based on Fortran routine 'apg.f'      The tropospheric asymmetric delay and the horizontal delay gradients ``G_N`` and ``G_E`` are computed by using the     ``apg.f`` function from the IERS software library [1]_. We use the asymmetric delay, which is based on Equation     (9.12) described in IERS conventions 2010 in Section 9.2 :cite:`iers2010`.      Args:         dset (Dataset):       A Dataset containing model data.      Returns:         tuple of Numpy Arrays: Includes the following elements, each with entries for each observation      ============  ===========  =====================================================================================      Element       Unit         Description     ============  ===========  =====================================================================================      gn            m            Horizontal delay gradient in the North direction      ge            m            Horizontal delay gradient in the East direction     ============  ===========  ===================================================================================== # Get horizontal gradients G_N and G_E Calculates zenith wet delay based on Askne and Nordius :cite:`askne1987`      The Fortran routine 'asknewet.f' of the GPT2w library is used, which is based on Equation (18) of Askne and Nordius     :cite:`askne1987`.      Args:         e (numpy.ndarray):      Water vapor pressure for each observation in [hPa]         tm (numpy.ndarray):     Mean temperature of the water vapor for each observation in [K]         lambd (numpy.ndarray):  Water vapor decrease factor for each observation      Returns:         numpy.ndarray:    Zenith wet delay for each observation in [m] # Loop over all observations Calculates zenith wet delay based on Saastamoinen/Davis model      The total tropospheric delay for a given zenith distance :math:`z` is determined after Equation (19a) in     Saastamoinen :cite:`saastamoinen1972`:      .. math::        \Delta T = 0.002277 \cdot \sec z \cdot (p + (1255/T + 0.05) \cdot e) - 1.16 \cdot \tan^2 z      The zenith tropospheric delay is determined with :math:`z = 0`:      .. math::        \Delta T^z = \Delta T_h^z + \Delta T_w^z = zhd + zwd      with the zenith hydrostatic delay :math:`zhd = 0.002277 \cdot p` and zenith wet delay :math:`zwd = 0.002277 \cdot     (1255/T + 0.05) \cdot e`.      A Fortran routine written by Davis corrects also for the gravity effect due to height :math:`H` and latitude     :math:`\phi` (see http://acc.igs.org/tropo/wetsaas.f) and uses the constant :math:`0.0022768` instead of     :math:`0.002277`, which leads to      .. math::        zwd = 0.0022768 \cdot (1255/T + 0.05) \cdot e / (1 - 0.00266 \cos 2 \phi - 0.00000028 H) .      The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to     100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant in this     equation.  The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.      Args:         latitude (numpy.ndarray):     Geodetic latitude for each observation in [rad]         height (numpy.ndarray):       Orthometric height for each observation in [m]         temperature (numpy.ndarray):  Temperature for each observation in [Celsius]         e (numpy.ndarray):            Water vapor pressure for each observation in [hPa]      Returns:         numpy.ndarray:    Zenith wet delay for each observation in [m] # Zenith wet delay based on Eq. (19a) in Saastamoinen :cite:`saastamoinen1972` with additional gravity # correction Calculates GMF hydrostatic and wet mapping functions      Use the 'gmf.f' Fortran routine from the IERS software library to calculate the Global Mapping Function (see     Section 9.2 in :cite:`iers2010`), which are described in Boehm et al. :cite:`boehm2006b`.      Args:         dset (Dataset):    Model data.      Returns:         tuple of Numpy Arrays: Includes the following elements, each with entries for each observation      ============  ===========  =======================================================      Element       Unit         Description     ============  ===========  =======================================================      mh                         Hydrostatic mapping function coefficient ah      mw                         Wet mapping function coefficient aw     ============  ===========  ======================================================= Calculates Global Pressure and Temperature (GPT)      This subroutine determines atmospheric pressure and temperature globally based on spherical harmonics up to degree     and order 9 (see Boehm et al. :cite:`boehm2007`) for a given latitude, longitude and ellipsoidal height.      Args:         dset (Dataset):    Model data.      Returns:         tuple of Numpy Arrays: Includes the following elements, each with entries for each observation      ============  ===========  =======================================================      Element       Unit         Description     ============  ===========  =======================================================      pressure      hPa          Pressure value      temperature   Celsius      Temperature values      geoid_undu    m            Geoid undulation (based on 9x9 EGM model)     ============  ===========  ======================================================= # Note: For GPT2 and GPT2w linear interpolation is done between daily solutions. Here #      it does not seem to be necessary, because the performance of GPT Fortran routine #      is not so worse than for GPT2 and GPT2w. Calculates meteorological data and mapping function coefficients based on GPT2 model      The GPT2 model is described in Lagler et al. :cite:`lagler2013`.      Args:         dset (Dataset):    Model data.      Returns:         tuple of Numpy Arrays: Includes the following elements, each with entries for each observation      ============  ===========  =======================================================      Element       Unit         Description     ============  ===========  =======================================================      pressure      hPa          Pressure value      temperature   Celsius      Temperature values      dt            degree/km    Temperature lapse rate      e             hPa          Water vapor pressure      mh                         Hydrostatic mapping function coefficient ah      mw                         Wet mapping function coefficient aw      geoid_undu    m            Geoid undulation (based on 9x9 EGM model)     ============  ===========  ======================================================= # Determine GPT2 values for each observation by interpolating between two unique # daily solutions # Start 'gpt2.f' day-by-day in folder where 'gpt2_5.grd' is placed and carry out # linear interpolation # Determine mapping function values based on coefficients 'ah' and 'aw' Calculates meteorological data and mapping function coefficients based on GPT2 model      The functions calls the IERS library routine 'gpt2.f' (see Section 9.2 in :cite:`iers2010`). The Fortran routine     ``gpt2.f`` reads the grid file ``gpt2_5.grd``, which should be available in the same folder, where the Fortran     programs runs. Therefore we change the current directory to the IERS source directory, so that 'gpt2.f' can read     the grid file.      Due to performance reasons the GPT2 values are not determined for each observation. The call of the Fortran routine     ``gpt2.f`` takes time, because the grid file ``gpt2_5.grd`` has to be read for each observation. Instead the GPT2     values are calculated only once for each unique day (modified Julian date rounded to integer) and saved in the     cache _GPT2. The final GPT2 values are computed by a linear interpolation of the daily determined GPT2 values. The     difference between the use of routine ``gpt2.f`` for each observation and the linear interpolation between daily     solution is on the submillimeter level and can therefore be neglected.      Args:         mjd (numpy.float64):  Modified Julian date.         lat (list):           Array with latitude for each station in [rad].         lon (list):           Array with longitude for each station in [rad].         hell (list):          Array with height for each station in [m].      Returns:         numpy.ndarray:  Array with following entries:      =======  ===========  =======================================================      Index    Unit         Description     =======  ===========  =======================================================      [0]      hPa          Pressure value      [1]      Celsius      Temperature values      [2]      degree/km    Temperature lapse rate      [3]      hPa          Water vapor pressure      [4]                   Hydrostatic mapping function coefficient ah      [5]                   Wet mapping function coefficient aw      [6]      m            Geoid undulation (based on 9x9 EGM model)     =======  ===========  ======================================================= # Number of stations # Use of time variations (annual and semiannual terms) # TODO: Case not handled if several stations are included in Dataset. Is that the case for VLBI? # Change directory so that gpt2.f can read the gpt2_5.grd-file in the IERS source directory # Loop over all unique dates (rounded to integer value) # Check if date is already included in cache # Linear interpolation between two daily GPT2 solutions Calculates meteorological data and mapping function coefficients based on GPT2w model      The GPT2w model is described in Boehm et al. :cite:`boehm2015`.      Args:         dset (Dataset):    Model data.      Returns:         tuple of Numpy Arrays: Includes the following elements, each with entries for each observation      ============  ===========  =======================================================      Element       Unit         Description     ============  ===========  =======================================================      pressure      hPa          Pressure value      temperature   Celsius      Temperature values      dt            degree/km    Temperature lapse rate      tm            K            Mean temperature of the water vapor      e             hPa          Water vapor pressure      mh                         Hydrostatic mapping function coefficient ah      mw                         Wet mapping function coefficient aw      la                         Water vapor decrease factor      geoid_undu    m            Geoid undulation (based on 9x9 EGM model)     ============  ===========  ======================================================= # Determine GPT2W values for each observation by interpolating between two unique # daily solutions # Start 'gpt2.f' day-by-day in folder where 'gpt2_5.grd' is placed and carry out # linear interpolation # Determine mapping function values based on coefficients 'ah' and 'aw' Calculates meteorological data and mapping function coefficients based on GPT2w model      The functions calls the GPT2w library routine ``gpt2w_1w.f`` (see     http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w). The Fortran routine ``gpt2w_1w.f`` reads the grid file     ``gpt2_1wA.grd``, which should be available in the same folder, where the Fortran programs runs. Therefore we     change the current directory to the GPT2w source directory, so that ``gpt2w_1w.f`` can read the grid file.      Due to performance reasons the GPT2w values are not determined for each observation. The call of the Fortran     routine ``gpt2w_1w.f`` takes time, because the grid file ``gpt2_1wA.grd`` has to be read for each     observation. Instead the GPT2w values are calculated only once for each unique day (modified Julian date rounded to     integer) and saved in the cache _GPT2W. The final GPT2W values are computed by a linear interpolation of the daily     determined GPT2W values.  The difference between the use of routine ``gpt2w_1w.f`` for each observation and the     linear interpolation between daily solution is on the submillimeter level and can therefore be neglected.      Args:         mjd (numpy.float64):  Modified Julian date.         lat (list):           Array with latitude for each station in [rad].         lon (list):           Array with longitude for each station in [rad].         hell (list):          Array with height for each station in [m].      Returns:         numpy.ndarray:  Array with following entries:      =======  ===========  =======================================================      Index    Unit         Description     =======  ===========  =======================================================      [0]      hPa          Pressure value      [1]      Celsius      Temperature values      [2]      degree/km    Temperature lapse rate      [3]      K            Mean temperature of the water vapor      [4]      hPa          Water vapor pressure      [5]                   Hydrostatic mapping function coefficient ah      [6]                   Wet mapping function coefficient aw      [7]                   Water vapor decrease factor      [8]      m            Geoid undulation (based on 9x9 EGM model)     =======  ===========  ======================================================= # Number of stations # Use of time variations (annual and semiannual terms) # TODO: Case not handled if several stations are included in Dateset. Is that the case for VLBI? # Change directory so that gpt2w.f can read the gpt2_5.grd-file in the GPT2w source directory # Loop over all unique dates (rounded to integer value) # Check if date is already included in cache # Linear interpolation between two daily GPT2W solutions Calculates atmospheric pressure based on zenith hydrostatic troposphere delay      Args:         zhd (numpy.ndarray):         Zenith hydrostatic delay for each observation in [m]         latitude (numpy.ndarray):    Geodetic latitude for each observation in [rad]         height (numpy.ndarray):      Orthometric height for each observation in [m]      Returns:         numpy.ndarray:    Atmospheric pressure for each observation in [hPa] Atmospheric pressure is corrected from reference height to given height      The atmospheric pressure height correction is based on a standard atmosphere model of Berg. See equation (1) in     Boehm et al. :cite:`boehm2007` or equation (4) in Kouba :cite:`kouba2007`.      The method uses equation (4) in Kouba :cite:`kouba2007` to determine the pressure for the reference height     :math:`p_r` and for the station height :math:`p_s`:      .. math::         p_r &= 1013.25 (1 - 0.0000226 h_r)^{5.225} ,         p_s &= 1013.25 (1 - 0.0000226 h_s)^{5.225} .      The pressure at the station height is determined by setting these two equations     in following relation:      .. math:: p_s = (1 - 0.0000226h_r)^{5.225} / (1 - 0.0000226h_s)^{5.225} * p_r .      TODO:         The implementation of the pressure height correction is not completely described in Kouba :cite:`kouba2007`.  I         found different solutions from Trimble, <NAME> and Bernese, whereby Bernese uses a similiar approach for the         computation (see BERN52/LIB/FOR/VMF1ELL.f90). The use of the Trimble solution leads to differences of several         millimeters for the zenith hydrostatic troposphere delay compared to our solution. Our solution compared to Per         Helge's solution shows differences on the millimeter level.  What is the correct solution?          Trimble solution::              pressure = pressure_ref + 1013.25 * ((1 - 0.0000226 * height)**5.225 - (1 - 0.0000226 * height_ref)**5.225)          Per Helge's solution::              pressure = pressure_ref * (1 - 0.0000226 * (height-height_ref))**5.225      Args:         pressure_ref (numpy.ndarray):    Atmospheric pressure given for reference height for each observation in [hPa]         height_ref (numpy.ndarray):      Reference height for each observation in [m]         height (numpy.ndarray):          Height for each observation in [m]      Returns:         numpy.ndarray:      Pressure for given height for each observation in [hPa] Calculates gravity correction due to height and latitude      The gravity correction due to height and latitude changes is shown in Davis et al. :cite:`davis1985` (see Equation     (A13)) and is based on Saastamoinen :cite:`saastamoinen1972`.      Args:         latitude (numpy.ndarray):      Geodetic latitude for each observation in [rad]         height (numpy.ndarray):        Orthometric height for each observation in [m]      Returns:         numpy.ndarray: Gravity correction for each observation in [m] # Gravity correction based on Eq. (A13) in Davis et al. :cite:`davis1985` Calculates atmospheric pressure based on Saastamoinen model      The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to     100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant by using     equation (9.11). The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.      Args:         zhd (numpy.ndarray):         Zenith hydrostatic delay for each observation in [m]         latitude (numpy.ndarray):    Geodetic latitude for each observation in [rad]         height (numpy.ndarray):      Orthometric height for each observation in [m]      Returns:         numpy.ndarray:    Atmospheric pressure for each observation in [hPa]      Example:         >>> saastamoinen_pressure(2.2762945349647778, 0.83776, 200)         1000.0 Calculates zenith hydrostatic delay with Saastamoinen model      The zenith hydrostatic delay :math:`zhd` is calculated using equation (9.11) from Saastamoinen in the IERS     Conventions :cite:`iers2010`:      .. math::          \Delta T_h^z = zhd = (0.0022768 \pm 0.0000005) P_0 / f_s(\phi, H)      with the atmospheric pressure :math:`P_0` and the correction :math:`f_s(\phi, H)`     (see equation (9.4) in :cite:`iers2010`):      .. math::        f_s(\phi, H) = 1 - 0.00266 \cos 2 \phi - 0.00000028 H ,      where :math:`\phi` is the geodetic latitude of the station and :math:`H` is the geodetic height of the station.      The difference between the orthometric height and geodetic height is called geoid undulation and can reach up to     100 m due to Boehm et al. :cite:`boehm2007`. The influence of height difference is not significant by using     equation (9.11). The geodetic (ellipsoidal) height is therefore often used instead of the orthometric height.      Args:         pressure (numpy.ndarray):      Atmospheric pressure for each observation in [hPa]         latitude (numpy.ndarray):      Geodetic latitude for each observation in [rad]         height (numpy.ndarray):        Orthometric height for each observation in [m]      Returns:         numpy.ndarray:    Zenith hydrostatic delay for each observation in [m]      Example:         >>> saastamoinen_zenith_hydrostatic_delay(1000, 0.83776, 200)         2.2762945349647778      TODO:         Example is not compatible to http://ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w/FORTRAN/saasthyd.f         example (zhd = 2.2695 m). Calculates zenith wet delay based Saastamoinen model      The total tropospheric delay for a given zenith distance :math:`z` is determined after Equation (19a) in     Saastamoinen :cite:`saastamoinen1972`:      .. math::        \Delta T = 0.002277 \cdot \sec z \cdot (p + (1255/T + 0.05) \cdot e) - 1.16 \cdot \tan^2 z      The zenith tropospheric delay is determined with :math:`z = 0`:      .. math::        \Delta T^z = \Delta T_h^z + \Delta T_w^z = zhd + zwd      with the zenith hydrostatic delay :math:`zhd = 0.002277 \cdot p` and zenith wet delay     :math:`zwd = 0.002277 \cdot (1255/T + 0.05) \cdot e`.      Args:         latitude (numpy.ndarray):     Geodetic latitude for each observation in [rad]         height (numpy.ndarray):       Orthometric height for each observation in [m]         temperature (numpy.ndarray):  Temperature for each observation in [Celsius]         e (numpy.ndarray):            Water vapor pressure for each observation in [hPa]      Returns:         numpy.ndarray:     Zenith wet delay for each observation in [m] # Zenith wet delay based on Eq. (19a) in Saastamoinen :cite:`saastamoinen1972` Get atmospheric pressure from local site measurements      If local atmospheric pressure measurements on a site are not available an alternative model given in configuration     file is used to determine atmospheric pressure.      TODO:         So far only gridded VMF1 model is used, if local pressure data are not available.  Which alternative is used,         should be decided via the configuration file. How to check after an alternative model in configuration file?         model_list = config.tech[MODEL].meteorological_data.list???      Args:         dset (Dataset): A Dataset containing model data.      Returns:         numpy.ndarray: Atmospheric pressure for each observation in [hPa] Calulates VMF1 gridded atmospheric pressure      Equation (9.11) in IERS conventions 2010 :cite:`iers2010` is converted to 'pressure' and the VMF1 gridded     atmospheric pressure can be determined with help of the gridded zenith hydrostatic delays from VMF1. The gridded     VMF1 pressure values needs to be rescaled to actual station height. The method is described with Eq. (4) in Kouba     :cite:`kouba2007`.      Args:         dset (Dataset):    Model data.      Returns:         numpy.ndarray:  Atmospheric pressure for each observation in [hPa]. # Get gridded VMF1 data # Interpolation in time and space in VMF1 grid # Rescale gridded pressure to station pressure Calculates VMF1 hydrostatic and wet mapping functions      This routine determines the VMF1 (Vienna Mapping Functions 1) described in Boehm et al. :cite:`boehm2006a` and     Kouba :cite:`kouba2007` and uses the 'vmf1_ht.f' from the IERS software library :cite:`iers2010`. The mapping     function coefficients 'ah' and 'aw' are taken from VMF1 gridded files.      TODO: Multiplication with 1.e-8 should maybe be handled when reading gridded VMF1-files and not here.      Args:         dset (Dataset):    Model data.      Returns:         tuple of Numpy Arrays: Includes the following elements, each with entries for each observation      ============  ===========  =======================================================      Element       Unit         Description     ============  ===========  =======================================================      mh                         Hydrostatic mapping function coefficient ah      mw                         Wet mapping function coefficient aw     ============  ===========  ======================================================= # Get gridded VMF1 data # Interpolation in time and space in VMF1 grid # t = dset.time[obs] # TODO vlbi t_2 # if dset._default_field_suffix == '_2': # baseline_gcrs = dset.site_pos_2.gcrs - dset.site_pos_1.gcrs # delta_t = (dset.src_dir.unit_vector[:, None, :] @ baseline_gcrs[:, :, None])[:, 0, 0] / constant.c # t -= timedelta(seconds=delta_t[obs]) Calculates zenith wet delay based on gridded zenith wet delays from VMF1      Uses gridded zenith wet delays from VMF1, which are rescaled from the gridded height to actual station height by     using Equation(5) described in Kouba :cite:`kouba2007`.      Args:         dset (Dataset):    Model data.      Returns:         numpy.ndarray:     Zenith wet delay for each observation in [m] # Get gridded VMF1 data # Interpolation in time and space in VMF1 grid # Zenith Wet delay. Eq. (5) in Kouba :cite:`kouba2007` Calculate days including the given datetimes      Args:         datetimes (numpy.float64):   np.array with datetimes as floats      Returns:         numpy.ndarray:  Unique days representing the datetimes.      Example:         >>> _rounded_dates(np.array([54001.2, 54001.4, 54001.8, 54002, 54005.5, 54006.]))         array([ 54001.,  54002.,  54003.,  54005.,  54006.,  54007.])