mrtk_atmos.c

Go to the documentation of this file.

/*------------------------------------------------------------------------------
 * mrtk_atmos.c : atmosphere model functions
 *
 * Copyright (C) 2026 H.SHIONO (MRTKLIB Project)
 * Copyright (C) 2023-2025 Cabinet Office, Japan
 * Copyright (C) 2024-2025 Lighthouse Technology & Consulting Co. Ltd.
 * Copyright (C) 2023-2025 Japan Aerospace Exploration Agency
 * Copyright (C) 2023-2025 TOSHIBA ELECTRONIC TECHNOLOGIES CORPORATION
 * Copyright (C) 2015- Mitsubishi Electric Corp.
 * Copyright (C) 2014 Geospatial Information Authority of Japan
 * Copyright (C) 2014 T.SUZUKI
 * Copyright (C) 2007-2023 T.TAKASU
 *
 * SPDX-License-Identifier: BSD-2-Clause
 *----------------------------------------------------------------------------*/
/**
 * @file mrtk_atmos.c
 * @brief MRTKLIB Atmosphere Module — Ionospheric and tropospheric models.
 *
 * This file contains all broadcast ionosphere and troposphere model functions
 * extracted from rtkcmn.c with zero algorithmic changes (pure cut-and-paste).
 */
#include "mrtklib/mrtk_atmos.h"
 
#include <math.h>
 
#include "mrtklib/mrtk_ionex.h"
#include "mrtklib/mrtk_sbas.h"
#include "mrtklib/mrtk_trace.h"
 
/*============================================================================
 * Constants (duplicated from rtklib.h to avoid header dependency)
 *===========================================================================*/
 
#define PI 3.1415926535897932 /* pi */
#define D2R (PI / 180.0)      /* deg to rad */
#define R2D (180.0 / PI)      /* rad to deg */
#define CLIGHT 299792458.0    /* speed of light (m/s) */
#define RE_WGS84 6378137.0    /* earth semimajor axis (WGS84) (m) */
#define HION 350000.0         /* ionosphere height (m) */
 
/*============================================================================
 * Forward declarations for legacy functions (resolved at link time)
 *===========================================================================*/
 
#ifdef IERS_MODEL
extern int gmf_(double* mjd, double* lat, double* lon, double* hgt, double* zd, double* gmfh, double* gmfw);
extern double geoidh(const double* pos);
#endif
 
/*============================================================================
 * Ionosphere Models
 *===========================================================================*/
 
/* ionosphere model ------------------------------------------------------------
 * compute ionospheric delay by broadcast ionosphere model (klobuchar model)
 * args   : gtime_t t        I   time (gpst)
 *          double *ion      I   iono model parameters {a0,a1,a2,a3,b0,b1,b2,b3}
 *          double *pos      I   receiver position {lat,lon,h} (rad,m)
 *          double *azel     I   azimuth/elevation angle {az,el} (rad)
 * return : ionospheric delay (L1) (m)
 *-----------------------------------------------------------------------------*/
extern double ionmodel(gtime_t t, const double* ion, const double* pos, const double* azel) {
    const double ion_default[] = {
        /* 2004/1/1 */
        0.1118E-07, -0.7451E-08, -0.5961E-07, 0.1192E-06, 0.1167E+06, -0.2294E+06, -0.1311E+06, 0.1049E+07};
    double tt, f, psi, phi, lam, amp, per, x;
    int week;
 
    if (pos[2] < -1E3 || azel[1] <= 0) {
        return 0.0;
    }
    if (norm(ion, 8) <= 0.0) {
        ion = ion_default;
    }
 
    /* earth centered angle (semi-circle) */
    psi = 0.0137 / (azel[1] / PI + 0.11) - 0.022;
 
    /* subionospheric latitude/longitude (semi-circle) */
    phi = pos[0] / PI + psi * cos(azel[0]);
    if (phi > 0.416) {
        phi = 0.416;
    } else if (phi < -0.416) {
        phi = -0.416;
    }
    lam = pos[1] / PI + psi * sin(azel[0]) / cos(phi * PI);
 
    /* geomagnetic latitude (semi-circle) */
    phi += 0.064 * cos((lam - 1.617) * PI);
 
    /* local time (s) */
    tt = 43200.0 * lam + time2gpst(t, &week);
    tt -= floor(tt / 86400.0) * 86400.0; /* 0<=tt<86400 */
 
    /* slant factor */
    f = 1.0 + 16.0 * pow(0.53 - azel[1] / PI, 3.0);
 
    /* ionospheric delay */
    amp = ion[0] + phi * (ion[1] + phi * (ion[2] + phi * ion[3]));
    per = ion[4] + phi * (ion[5] + phi * (ion[6] + phi * ion[7]));
    amp = amp < 0.0 ? 0.0 : amp;
    per = per < 72000.0 ? 72000.0 : per;
    x = 2.0 * PI * (tt - 50400.0) / per;
 
    return CLIGHT * f * (fabs(x) < 1.57 ? 5E-9 + amp * (1.0 + x * x * (-0.5 + x * x / 24.0)) : 5E-9);
}
/* ionosphere mapping function -------------------------------------------------
 * compute ionospheric delay mapping function by single layer model
 * args   : double *pos      I   receiver position {lat,lon,h} (rad,m)
 *          double *azel     I   azimuth/elevation angle {az,el} (rad)
 * return : ionospheric mapping function
 *-----------------------------------------------------------------------------*/
extern double ionmapf(const double* pos, const double* azel) {
    if (pos[2] >= HION) {
        return 1.0;
    }
    return 1.0 / cos(asin((RE_WGS84 + pos[2]) / (RE_WGS84 + HION) * sin(PI / 2.0 - azel[1])));
}
/* ionospheric pierce point position -------------------------------------------
 * compute ionospheric pierce point (ipp) position and slant factor
 * args   : double *pos      I   receiver position {lat,lon,h} (rad,m)
 *          double *azel     I   azimuth/elevation angle {az,el} (rad)
 *          double re        I   earth radius (km)
 *          double hion      I   altitude of ionosphere (km)
 *          double *posp     O   pierce point position {lat,lon,h} (rad,m)
 * return : slant factor
 * notes  : see ref [2], only valid on the earth surface
 *          fixing bug on ref [2] A.4.4.10.1 A-22,23
 *-----------------------------------------------------------------------------*/
extern double ionppp(const double* pos, const double* azel, double re, double hion, double* posp) {
    double cosaz, rp, ap, sinap, tanap;
 
    rp = re / (re + hion) * cos(azel[1]);
    ap = PI / 2.0 - azel[1] - asin(rp);
    sinap = sin(ap);
    tanap = tan(ap);
    cosaz = cos(azel[0]);
    posp[0] = asin(sin(pos[0]) * cos(ap) + cos(pos[0]) * sinap * cosaz);
 
    if ((pos[0] > 70.0 * D2R && tanap * cosaz > tan(PI / 2.0 - pos[0])) ||
        (pos[0] < -70.0 * D2R && -tanap * cosaz > tan(PI / 2.0 + pos[0]))) {
        posp[1] = pos[1] + PI - asin(sinap * sin(azel[0]) / cos(posp[0]));
    } else {
        posp[1] = pos[1] + asin(sinap * sin(azel[0]) / cos(posp[0]));
    }
    return 1.0 / sqrt(1.0 - rp * rp);
}
 
/*============================================================================
 * Troposphere Models
 *===========================================================================*/
 
/* troposphere model -----------------------------------------------------------
 * compute tropospheric delay by standard atmosphere and saastamoinen model
 * args   : gtime_t time     I   time
 *          double *pos      I   receiver position {lat,lon,h} (rad,m)
 *          double *azel     I   azimuth/elevation angle {az,el} (rad)
 *          double humi      I   relative humidity
 * return : tropospheric delay (m)
 *-----------------------------------------------------------------------------*/
extern double tropmodel(gtime_t time, const double* pos, const double* azel, double humi) {
    const double temp0 = 15.0; /* temparature at sea level */
    double hgt, pres, temp, e, z, trph, trpw;
 
    if (pos[2] < -100.0 || 1E4 < pos[2] || azel[1] <= 0) {
        return 0.0;
    }
 
    /* standard atmosphere */
    hgt = pos[2] < 0.0 ? 0.0 : pos[2];
 
    pres = 1013.25 * pow(1.0 - 2.2557E-5 * hgt, 5.2568);
    temp = temp0 - 6.5E-3 * hgt + 273.16;
    e = 6.108 * humi * exp((17.15 * temp - 4684.0) / (temp - 38.45));
 
    /* saastamoninen model */
    z = PI / 2.0 - azel[1];
    trph = 0.0022768 * pres / (1.0 - 0.00266 * cos(2.0 * pos[0]) - 0.00028 * hgt / 1E3) / cos(z);
    trpw = 0.002277 * (1255.0 / temp + 0.05) * e / cos(z);
    return trph + trpw;
}
#ifndef IERS_MODEL
 
static double interpc(const double coef[], double lat) {
    int i = (int)(lat / 15.0);
    if (i < 1) {
        return coef[0];
    } else if (i > 4) {
        return coef[4];
    }
    return coef[i - 1] * (1.0 - lat / 15.0 + i) + coef[i] * (lat / 15.0 - i);
}
static double mapf(double el, double a, double b, double c) {
    double sinel = sin(el);
    return (1.0 + a / (1.0 + b / (1.0 + c))) / (sinel + (a / (sinel + b / (sinel + c))));
}
static double nmf(gtime_t time, const double pos[], const double azel[], double* mapfw) {
    /* ref [5] table 3 */
    /* hydro-ave-a,b,c, hydro-amp-a,b,c, wet-a,b,c at latitude 15,30,45,60,75 */
    const double coef[][5] = {{1.2769934E-3, 1.2683230E-3, 1.2465397E-3, 1.2196049E-3, 1.2045996E-3},
                              {2.9153695E-3, 2.9152299E-3, 2.9288445E-3, 2.9022565E-3, 2.9024912E-3},
                              {62.610505E-3, 62.837393E-3, 63.721774E-3, 63.824265E-3, 64.258455E-3},
 
                              {0.0000000E-0, 1.2709626E-5, 2.6523662E-5, 3.4000452E-5, 4.1202191E-5},
                              {0.0000000E-0, 2.1414979E-5, 3.0160779E-5, 7.2562722E-5, 11.723375E-5},
                              {0.0000000E-0, 9.0128400E-5, 4.3497037E-5, 84.795348E-5, 170.37206E-5},
 
                              {5.8021897E-4, 5.6794847E-4, 5.8118019E-4, 5.9727542E-4, 6.1641693E-4},
                              {1.4275268E-3, 1.5138625E-3, 1.4572752E-3, 1.5007428E-3, 1.7599082E-3},
                              {4.3472961E-2, 4.6729510E-2, 4.3908931E-2, 4.4626982E-2, 5.4736038E-2}};
    const double aht[] = {2.53E-5, 5.49E-3, 1.14E-3}; /* height correction */
 
    double y, cosy, ah[3], aw[3], dm, el = azel[1], lat = pos[0] * R2D, hgt = pos[2];
    int i;
 
    if (el <= 0.0) {
        if (mapfw) {
            *mapfw = 0.0;
        }
        return 0.0;
    }
    /* year from doy 28, added half a year for southern latitudes */
    y = (time2doy(time) - 28.0) / 365.25 + (lat < 0.0 ? 0.5 : 0.0);
 
    cosy = cos(2.0 * PI * y);
    lat = fabs(lat);
 
    for (i = 0; i < 3; i++) {
        ah[i] = interpc(coef[i], lat) - interpc(coef[i + 3], lat) * cosy;
        aw[i] = interpc(coef[i + 6], lat);
    }
    /* ellipsoidal height is used instead of height above sea level */
    dm = (1.0 / sin(el) - mapf(el, aht[0], aht[1], aht[2])) * hgt / 1E3;
 
    if (mapfw) {
        *mapfw = mapf(el, aw[0], aw[1], aw[2]);
    }
 
    return mapf(el, ah[0], ah[1], ah[2]) + dm;
}
#endif /* !IERS_MODEL */
 
/* troposphere mapping function ------------------------------------------------
 * compute tropospheric mapping function by NMF
 * args   : gtime_t t        I   time
 *          double *pos      I   receiver position {lat,lon,h} (rad,m)
 *          double *azel     I   azimuth/elevation angle {az,el} (rad)
 *          double *mapfw    IO  wet mapping function (NULL: not output)
 * return : dry mapping function
 * note   : see ref [5] (NMF) and [9] (GMF)
 *          original JGR paper of [5] has bugs in eq.(4) and (5). the corrected
 *          paper is obtained from:
 *          ftp://web.haystack.edu/pub/aen/nmf/NMF_JGR.pdf
 *-----------------------------------------------------------------------------*/
extern double tropmapf(gtime_t time, const double pos[], const double azel[], double* mapfw) {
#ifdef IERS_MODEL
    const double ep[] = {2000, 1, 1, 12, 0, 0};
    double mjd, lat, lon, hgt, zd, gmfh, gmfw;
#endif
    trace(NULL, 4, "tropmapf: pos=%10.6f %11.6f %6.1f azel=%5.1f %4.1f\n", pos[0] * R2D, pos[1] * R2D, pos[2],
          azel[0] * R2D, azel[1] * R2D);
 
    if (pos[2] < -1000.0 || pos[2] > 20000.0) {
        if (mapfw) {
            *mapfw = 0.0;
        }
        return 0.0;
    }
#ifdef IERS_MODEL
    mjd = 51544.5 + (timediff(time, epoch2time(ep))) / 86400.0;
    lat = pos[0];
    lon = pos[1];
    hgt = pos[2] - geoidh(pos); /* height in m (mean sea level) */
    zd = PI / 2.0 - azel[1];
 
    /* call GMF */
    gmf_(&mjd, &lat, &lon, &hgt, &zd, &gmfh, &gmfw);
 
    if (mapfw) *mapfw = gmfw;
    return gmfh;
#else
    return nmf(time, pos, azel, mapfw); /* NMF */
#endif
}
 
/*============================================================================
 * Ionosphere / Troposphere Correction (multi-model dispatch)
 *===========================================================================*/
 
#define SQR(x) ((x) * (x))
#define ERR_ION 5.0   /* ionospheric delay Std (m) */
#define ERR_TROP 3.0  /* tropspheric delay Std (m) */
#define ERR_SAAS 0.3  /* Saastamoinen model error Std (m) */
#define ERR_BRDCI 0.5 /* broadcast ionosphere model error factor */
#define REL_HUMI 0.7  /* relative humidity for Saastamoinen model */
 
/* ionospheric correction ------------------------------------------------------
 * compute ionospheric correction
 * args   : gtime_t time     I   time
 *          nav_t  *nav      I   navigation data
 *          int    sat       I   satellite number
 *          double *pos      I   receiver position {lat,lon,h} (rad|m)
 *          double *azel     I   azimuth/elevation angle {az,el} (rad)
 *          int    ionoopt   I   ionospheric correction option (IONOOPT_???)
 *          double *ion      O   ionospheric delay (L1) (m)
 *          double *var      O   ionospheric delay (L1) variance (m^2)
 * return : status(1:ok,0:error)
 *-----------------------------------------------------------------------------*/
int ionocorr(gtime_t time, const nav_t* nav, int sat, const double* pos, const double* azel, int ionoopt, double* ion,
             double* var) {
    trace(NULL, 4, "ionocorr: time=%s opt=%d sat=%2d pos=%.3f %.3f azel=%.3f %.3f\n", time_str(time, 3), ionoopt, sat,
          pos[0] * R2D, pos[1] * R2D, azel[0] * R2D, azel[1] * R2D);
 
    /* GPS broadcast ionosphere model */
    if (ionoopt == IONOOPT_BRDC) {
        *ion = ionmodel(time, nav->ion_gps, pos, azel);
        *var = SQR(*ion * ERR_BRDCI);
        return 1;
    }
    /* SBAS ionosphere model */
    if (ionoopt == IONOOPT_SBAS) {
        return sbsioncorr(time, nav, pos, azel, ion, var);
    }
    /* IONEX TEC model */
    if (ionoopt == IONOOPT_TEC) {
        return iontec(time, nav, pos, azel, 1, ion, var);
    }
    /* QZSS broadcast ionosphere model */
    if (ionoopt == IONOOPT_QZS && norm(nav->ion_qzs, 8) > 0.0) {
        *ion = ionmodel(time, nav->ion_qzs, pos, azel);
        *var = SQR(*ion * ERR_BRDCI);
        return 1;
    }
    *ion = 0.0;
    *var = ionoopt == IONOOPT_OFF ? SQR(ERR_ION) : 0.0;
    return 1;
}
/* tropospheric correction -----------------------------------------------------
 * compute tropospheric correction
 * args   : gtime_t time     I   time
 *          nav_t  *nav      I   navigation data
 *          double *pos      I   receiver position {lat,lon,h} (rad|m)
 *          double *azel     I   azimuth/elevation angle {az,el} (rad)
 *          int    tropopt   I   tropospheric correction option (TROPOPT_???)
 *          double *trp      O   tropospheric delay (m)
 *          double *var      O   tropospheric delay variance (m^2)
 * return : status(1:ok,0:error)
 *-----------------------------------------------------------------------------*/
int tropcorr(gtime_t time, const nav_t* nav, const double* pos, const double* azel, int tropopt, double* trp,
             double* var) {
    trace(NULL, 4, "tropcorr: time=%s opt=%d pos=%.3f %.3f azel=%.3f %.3f\n", time_str(time, 3), tropopt, pos[0] * R2D,
          pos[1] * R2D, azel[0] * R2D, azel[1] * R2D);
 
    /* Saastamoinen model */
    if (tropopt == TROPOPT_SAAS || tropopt == TROPOPT_EST || tropopt == TROPOPT_ESTG) {
        *trp = tropmodel(time, pos, azel, REL_HUMI);
        *var = SQR(ERR_SAAS / (sin(azel[1]) + 0.1));
        return 1;
    }
    /* SBAS (MOPS) troposphere model */
    if (tropopt == TROPOPT_SBAS) {
        *trp = sbstropcorr(time, pos, azel, var);
        return 1;
    }
    /* no correction */
    *trp = 0.0;
    *var = tropopt == TROPOPT_OFF ? SQR(ERR_TROP) : 0.0;
    return 1;
}