mrtk_peph.c

Go to the documentation of this file.

/*------------------------------------------------------------------------------
 * mrtk_peph.c : precise ephemeris and earth rotation parameter 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_peph.c
 * @brief MRTKLIB Precise Ephemeris Module — Earth rotation parameter I/O,
 *        precise ephemeris/clock computation, and satellite antenna offset
 *        functions.
 *
 * This file contains the ERP functions extracted from rtkcmn.c, plus all
 * precise ephemeris functions absorbed from preceph.c with zero algorithmic
 * changes (pure cut-and-paste).
 */
#include "mrtklib/mrtk_peph.h"
 
#include <stdio.h>
#include <stdlib.h>
 
#include "mrtklib/mrtk_antenna.h"
#include "mrtklib/mrtk_astro.h"
#include "mrtklib/mrtk_mat.h"
#include "mrtklib/mrtk_nav.h"
#include "mrtklib/mrtk_sys.h"
 
/* forward declarations for functions in mrtk_nav.c */
extern int satsys_bd2(int sat, int* prn);
extern int freq_idx2freq_num(int sys, int freq_idx);
extern double freq_num2freq(int sys, int freq_num, int fcn);
extern int freq_idx2ant_idx(int sys, int freq_idx);
#include <math.h>
#include <string.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 AS2R (D2R / 3600.0)   /* arc sec to radian */
#define CLIGHT 299792458.0    /* speed of light (m/s) */
#define OMGE 7.2921151467E-5  /* earth angular velocity (IS-GPS) (rad/s) */
 
#define SQR(x) ((x) * (x))
 
/* navigation system identifiers */
#define SYS_NONE 0x00
#define SYS_GPS 0x01
#define SYS_SBS 0x02
#define SYS_GLO 0x04
#define SYS_GAL 0x08
#define SYS_QZS 0x10
#define SYS_CMP 0x20
#define SYS_IRN 0x40
#define SYS_LEO 0x80
 
/* carrier frequencies */
#define FREQ1 1.57542E9      /* L1/E1/B1C  frequency (Hz) */
#define FREQ2 1.22760E9      /* L2         frequency (Hz) */
#define FREQ5 1.17645E9      /* L5/E5a/B2a frequency (Hz) */
#define FREQ7 1.20714E9      /* E5b        frequency (Hz) */
#define FREQ9 2.492028E9     /* S          frequency (Hz) */
#define FREQ1_CMP 1.561098E9 /* BDS B1I    frequency (Hz) */
#define FREQ2_CMP 1.20714E9  /* BDS B2I/B2b frequency (Hz) */
 
/* precise ephemeris constants */
#define NMAX 10         /* order of polynomial interpolation */
#define MAXDTE 900.0    /* max time difference to ephem time (s) */
#define EXTERR_CLK 1E-3 /* extrapolation error for clock (m/s) */
#define EXTERR_EPH 5E-7 /* extrapolation error for ephem (m/s^2) */
 
/*============================================================================
 * Forward declarations for legacy functions (resolved at link time)
 *===========================================================================*/
 
/*============================================================================
 * Public Functions
 *===========================================================================*/
 
/* read earth rotation parameters ----------------------------------------------
 * read earth rotation parameters
 * args   : char   *file       I   IGS ERP file (IGS ERP ver.2)
 *          erp_t  *erp        O   earth rotation parameters
 * return : status (1:ok,0:file open error)
 *-----------------------------------------------------------------------------*/
extern int readerp(const char* file, erp_t* erp) {
    FILE* fp;
    erpd_t* erp_data;
    double v[14] = {0};
    char buff[256];
 
    trace(NULL, 3, "readerp: file=%s\n", file);
 
    if (!(fp = fopen(file, "r"))) {
        trace(NULL, 2, "erp file open error: file=%s\n", file);
        return 0;
    }
    while (fgets(buff, sizeof(buff), fp)) {
        if (sscanf(buff, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", v, v + 1, v + 2, v + 3, v + 4,
                   v + 5, v + 6, v + 7, v + 8, v + 9, v + 10, v + 11, v + 12, v + 13) < 5) {
            continue;
        }
        if (erp->n >= erp->nmax) {
            erp->nmax = erp->nmax <= 0 ? 128 : erp->nmax * 2;
            erp_data = (erpd_t*)realloc(erp->data, sizeof(erpd_t) * erp->nmax);
            if (!erp_data) {
                free(erp->data);
                erp->data = NULL;
                erp->n = erp->nmax = 0;
                fclose(fp);
                return 0;
            }
            erp->data = erp_data;
        }
        erp->data[erp->n].mjd = v[0];
        erp->data[erp->n].xp = v[1] * 1E-6 * AS2R;
        erp->data[erp->n].yp = v[2] * 1E-6 * AS2R;
        erp->data[erp->n].ut1_utc = v[3] * 1E-7;
        erp->data[erp->n].lod = v[4] * 1E-7;
        erp->data[erp->n].xpr = v[12] * 1E-6 * AS2R;
        erp->data[erp->n++].ypr = v[13] * 1E-6 * AS2R;
    }
    fclose(fp);
    return 1;
}
/* get earth rotation parameter values -----------------------------------------
 * get earth rotation parameter values
 * args   : erp_t  *erp        I   earth rotation parameters
 *          gtime_t time       I   time (gpst)
 *          double *erpv       O   erp values {xp,yp,ut1_utc,lod} (rad,rad,s,s/d)
 * return : status (1:ok,0:error)
 *-----------------------------------------------------------------------------*/
extern int geterp(const erp_t* erp, gtime_t time, double* erpv) {
    const double ep[] = {2000, 1, 1, 12, 0, 0};
    double mjd, day, a;
    int i, j, k;
 
    trace(NULL, 4, "geterp:\n");
 
    if (erp->n <= 0) {
        return 0;
    }
 
    mjd = 51544.5 + (timediff(gpst2utc(time), epoch2time(ep))) / 86400.0;
 
    if (mjd <= erp->data[0].mjd) {
        day = mjd - erp->data[0].mjd;
        erpv[0] = erp->data[0].xp + erp->data[0].xpr * day;
        erpv[1] = erp->data[0].yp + erp->data[0].ypr * day;
        erpv[2] = erp->data[0].ut1_utc - erp->data[0].lod * day;
        erpv[3] = erp->data[0].lod;
        return 1;
    }
    if (mjd >= erp->data[erp->n - 1].mjd) {
        day = mjd - erp->data[erp->n - 1].mjd;
        erpv[0] = erp->data[erp->n - 1].xp + erp->data[erp->n - 1].xpr * day;
        erpv[1] = erp->data[erp->n - 1].yp + erp->data[erp->n - 1].ypr * day;
        erpv[2] = erp->data[erp->n - 1].ut1_utc - erp->data[erp->n - 1].lod * day;
        erpv[3] = erp->data[erp->n - 1].lod;
        return 1;
    }
    for (j = 0, k = erp->n - 1; j < k - 1;) {
        i = (j + k) / 2;
        if (mjd < erp->data[i].mjd) {
            k = i;
        } else {
            j = i;
        }
    }
    if (erp->data[j].mjd == erp->data[j + 1].mjd) {
        a = 0.5;
    } else {
        a = (mjd - erp->data[j].mjd) / (erp->data[j + 1].mjd - erp->data[j].mjd);
    }
    erpv[0] = (1.0 - a) * erp->data[j].xp + a * erp->data[j + 1].xp;
    erpv[1] = (1.0 - a) * erp->data[j].yp + a * erp->data[j + 1].yp;
    erpv[2] = (1.0 - a) * erp->data[j].ut1_utc + a * erp->data[j + 1].ut1_utc;
    erpv[3] = (1.0 - a) * erp->data[j].lod + a * erp->data[j + 1].lod;
    return 1;
}
 
/*============================================================================
 * Precise Ephemeris Functions
 * (absorbed from preceph.c — pure cut-and-paste)
 *===========================================================================*/
 
/* satellite code to satellite system ----------------------------------------*/
static int code2sys(char code) {
    if (code == 'G' || code == ' ') {
        return SYS_GPS;
    }
    if (code == 'R') {
        return SYS_GLO;
    }
    if (code == 'E') {
        return SYS_GAL; /* SP3-d */
    }
    if (code == 'J') {
        return SYS_QZS; /* SP3-d */
    }
    if (code == 'C') {
        return SYS_CMP; /* SP3-d */
    }
    if (code == 'I') {
        return SYS_IRN; /* SP3-d */
    }
    if (code == 'L') {
        return SYS_LEO; /* SP3-d */
    }
    return SYS_NONE;
}
/* read SP3 header -----------------------------------------------------------*/
static int readsp3h(FILE* fp, gtime_t* time, char* type, int* sats, double* bfact, char* tsys) {
    int i, j, k = 0, ns = 0, sys, prn;
    char buff[1024];
 
    trace(NULL, 3, "readsp3h:\n");
 
    for (i = 0;; i++) {
        if (!fgets(buff, sizeof(buff), fp)) {
            break;
        }
 
        if (i == 0) {
            *type = buff[2];
            if (str2time(buff, 3, 28, time)) {
                return 0;
            }
        } else if (!strncmp(buff, "+ ", 2)) { /* satellite id */
            if (ns == 0) {
                ns = (int)str2num(buff, 3, 3);
            }
            for (j = 0; j < 17 && k < ns; j++) {
                sys = code2sys(buff[9 + 3 * j]);
                prn = (int)str2num(buff, 10 + 3 * j, 2);
                if (k < MAXSAT) {
                    sats[k++] = satno(sys, prn);
                }
            }
        } else if (!strncmp(buff, "++", 2)) { /* orbit accuracy */
            continue;
        } else if (!strncmp(buff, "%c", 2)) { /* time system */
            sprintf(tsys, "%.3s", buff + 9);
        } else if (!strncmp(buff, "%f", 2) && bfact[0] == 0.0) { /* fp base number */
            bfact[0] = str2num(buff, 3, 10);
            bfact[1] = str2num(buff, 14, 12);
        } else if (!strncmp(buff, "%i", 2)) {
            continue;
        } else if (!strncmp(buff, "/*", 2)) { /* comment */
            continue;
        } else if (!strncmp(buff, "* ", 2)) { /* first record */
            rewind(fp);                       /* rollback to start */
            break;
        }
    }
    return ns;
}
/* add precise ephemeris -----------------------------------------------------*/
static int addpeph(nav_t* nav, peph_t* peph) {
    peph_t* nav_peph;
 
    if (nav->ne >= nav->nemax) {
        nav->nemax += 256;
        if (!(nav_peph = (peph_t*)realloc(nav->peph, sizeof(peph_t) * nav->nemax))) {
            trace(NULL, 1, "readsp3b malloc error n=%d\n", nav->nemax);
            free(nav->peph);
            nav->peph = NULL;
            nav->ne = nav->nemax = 0;
            return 0;
        }
        nav->peph = nav_peph;
    }
    nav->peph[nav->ne++] = *peph;
    return 1;
}
/* read SP3 body -------------------------------------------------------------*/
static void readsp3b(FILE* fp, char type, int* sats, int ns, double* bfact, char* tsys, int index, int opt,
                     nav_t* nav) {
    peph_t peph;
    gtime_t time;
    double val, std, base;
    int i, j, sat, sys, prn, n = ns * (type == 'P' ? 1 : 2), pred_o, pred_c, v;
    char buff[1024];
 
    trace(NULL, 3, "readsp3b: type=%c ns=%d index=%d opt=%d\n", type, ns, index, opt);
 
    while (fgets(buff, sizeof(buff), fp)) {
        if (!strncmp(buff, "EOF", 3)) {
            break;
        }
 
        /* skip header or invalid records */
        if (buff[0] != '*' || str2time(buff, 3, 28, &time)) {
            continue;
        }
        if (!strcmp(tsys, "UTC")) {
            time = utc2gpst(time); /* utc->gpst */
        }
        peph.time = time;
        peph.index = index;
 
        for (i = 0; i < MAXSAT; i++) {
            for (j = 0; j < 4; j++) {
                peph.pos[i][j] = 0.0;
                peph.std[i][j] = 0.0f;
                peph.vel[i][j] = 0.0;
                peph.vst[i][j] = 0.0f;
            }
            for (j = 0; j < 3; j++) {
                peph.cov[i][j] = 0.0f;
                peph.vco[i][j] = 0.0f;
            }
        }
        for (i = pred_o = pred_c = v = 0; i < n && fgets(buff, sizeof(buff), fp); i++) {
            if (strlen(buff) < 4 || (buff[0] != 'P' && buff[0] != 'V')) {
                continue;
            }
 
            sys = buff[1] == ' ' ? SYS_GPS : code2sys(buff[1]);
            prn = (int)str2num(buff, 2, 2);
            if (sys == SYS_SBS) {
                prn += 100;
            } else if (sys == SYS_QZS) {
                prn += 192; /* extension to sp3-c */
            }
 
            if (!(sat = satno(sys, prn))) {
                continue;
            }
 
            if (buff[0] == 'P') {
                pred_c = strlen(buff) >= 76 && buff[75] == 'P';
                pred_o = strlen(buff) >= 80 && buff[79] == 'P';
            }
            for (j = 0; j < 4; j++) {
                /* read option for predicted value */
                if (j < 3 && (opt & 1) && pred_o) {
                    continue;
                }
                if (j < 3 && (opt & 2) && !pred_o) {
                    continue;
                }
                if (j == 3 && (opt & 1) && pred_c) {
                    continue;
                }
                if (j == 3 && (opt & 2) && !pred_c) {
                    continue;
                }
 
                val = str2num(buff, 4 + j * 14, 14);
                std = str2num(buff, 61 + j * 3, j < 3 ? 2 : 3);
 
                if (buff[0] == 'P') { /* position */
                    if (val != 0.0 && fabs(val - 999999.999999) >= 1E-6) {
                        peph.pos[sat - 1][j] = val * (j < 3 ? 1000.0 : 1E-6);
                        v = 1; /* valid epoch */
                    }
                    if ((base = bfact[j < 3 ? 0 : 1]) > 0.0 && std > 0.0) {
                        peph.std[sat - 1][j] = (float)(pow(base, std) * (j < 3 ? 1E-3 : 1E-12));
                    }
                } else if (v) { /* velocity */
                    if (val != 0.0 && fabs(val - 999999.999999) >= 1E-6) {
                        peph.vel[sat - 1][j] = val * (j < 3 ? 0.1 : 1E-10);
                    }
                    if ((base = bfact[j < 3 ? 0 : 1]) > 0.0 && std > 0.0) {
                        peph.vst[sat - 1][j] = (float)(pow(base, std) * (j < 3 ? 1E-7 : 1E-16));
                    }
                }
            }
        }
        if (v) {
            if (!addpeph(nav, &peph)) {
                return;
            }
        }
    }
}
/* compare precise ephemeris -------------------------------------------------*/
static int cmppeph(const void* p1, const void* p2) {
    peph_t *q1 = (peph_t*)p1, *q2 = (peph_t*)p2;
    double tt = timediff(q1->time, q2->time);
    return tt < -1E-9 ? -1 : (tt > 1E-9 ? 1 : q1->index - q2->index);
}
/* combine precise ephemeris -------------------------------------------------*/
static void combpeph(nav_t* nav, int opt) {
    int i, j, k, m;
 
    trace(NULL, 3, "combpeph: ne=%d\n", nav->ne);
 
    qsort(nav->peph, nav->ne, sizeof(peph_t), cmppeph);
 
    if (opt & 4) {
        return;
    }
 
    for (i = 0, j = 1; j < nav->ne; j++) {
        if (fabs(timediff(nav->peph[i].time, nav->peph[j].time)) < 1E-9) {
            for (k = 0; k < MAXSAT; k++) {
                if (norm(nav->peph[j].pos[k], 4) <= 0.0) {
                    continue;
                }
                for (m = 0; m < 4; m++) {
                    nav->peph[i].pos[k][m] = nav->peph[j].pos[k][m];
                }
                for (m = 0; m < 4; m++) {
                    nav->peph[i].std[k][m] = nav->peph[j].std[k][m];
                }
                for (m = 0; m < 4; m++) {
                    nav->peph[i].vel[k][m] = nav->peph[j].vel[k][m];
                }
                for (m = 0; m < 4; m++) {
                    nav->peph[i].vst[k][m] = nav->peph[j].vst[k][m];
                }
            }
        } else if (++i < j) {
            nav->peph[i] = nav->peph[j];
        }
    }
    nav->ne = i + 1;
 
    trace(NULL, 4, "combpeph: ne=%d\n", nav->ne);
}
/* read sp3 precise ephemeris file ---------------------------------------------
 * read sp3 precise ephemeris/clock files and set them to navigation data
 * args   : char   *file       I   sp3-c precise ephemeris file
 *                                 (wind-card * is expanded)
 *          nav_t  *nav        IO  navigation data
 *          int    opt         I   options (1: only observed + 2: only predicted +
 *                                 4: not combined)
 * return : none
 * notes  : see ref [1]
 *          precise ephemeris is appended and combined
 *          nav->peph and nav->ne must by properly initialized before calling the
 *          function
 *          only files with extensions of .sp3, .SP3, .eph* and .EPH* are read
 *-----------------------------------------------------------------------------*/
extern void readsp3(const char* file, nav_t* nav, int opt) {
    FILE* fp;
    gtime_t time = {0};
    double bfact[2] = {0};
    int i, j, n, ns, sats[MAXSAT] = {0};
    char *efiles[MAXEXFILE], *ext, type = ' ', tsys[4] = "";
 
    trace(NULL, 3, "readpephs: file=%s\n", file);
 
    for (i = 0; i < MAXEXFILE; i++) {
        if (!(efiles[i] = (char*)malloc(1024))) {
            for (i--; i >= 0; i--) {
                free(efiles[i]);
            }
            return;
        }
    }
    /* expand wild card in file path */
    n = expath(file, efiles, MAXEXFILE);
 
    for (i = j = 0; i < n; i++) {
        if (!(ext = strrchr(efiles[i], '.'))) {
            continue;
        }
 
        if (!strstr(ext, ".sp3") && !strstr(ext, ".SP3") && !strstr(ext, ".eph") && !strstr(ext, ".EPH")) {
            continue;
        }
 
        if (!(fp = fopen(efiles[i], "r"))) {
            trace(NULL, 2, "sp3 file open error %s\n", efiles[i]);
            continue;
        }
        /* read sp3 header */
        ns = readsp3h(fp, &time, &type, sats, bfact, tsys);
 
        /* read sp3 body */
        readsp3b(fp, type, sats, ns, bfact, tsys, j++, opt, nav);
 
        fclose(fp);
    }
    for (i = 0; i < MAXEXFILE; i++) {
        free(efiles[i]);
    }
 
    /* combine precise ephemeris */
    if (nav->ne > 0) {
        combpeph(nav, opt);
    }
}
/* read satellite antenna parameters -------------------------------------------
 * read satellite antenna parameters
 * args   : char   *file       I   antenna parameter file
 *          gtime_t time       I   time
 *          nav_t  *nav        IO  navigation data
 * return : status (1:ok,0:error)
 * notes  : only support antex format for the antenna parameter file
 *-----------------------------------------------------------------------------*/
extern int readsap(const char* file, gtime_t time, nav_t* nav) {
    pcvs_t pcvs = {0};
    pcv_t pcv0 = {0}, *pcv;
    int i;
 
    trace(NULL, 3, "readsap : file=%s time=%s\n", file, time_str(time, 0));
 
    if (!readpcv(file, &pcvs)) {
        return 0;
    }
 
    for (i = 0; i < MAXSAT; i++) {
        pcv = searchpcv(i + 1, "", time, &pcvs);
        nav->pcvs[i] = pcv ? *pcv : pcv0;
    }
    free(pcvs.pcv);
    return 1;
}
/* read DCB parameters file --------------------------------------------------*/
static int readdcbf(const char* file, nav_t* nav, const sta_t* sta) {
    FILE* fp;
    double cbias;
    char buff[256], str1[32], str2[32] = "";
    int i, j, sat, type = 0;
 
    trace(NULL, 3, "readdcbf: file=%s\n", file);
 
    if (!(fp = fopen(file, "r"))) {
        trace(NULL, 2, "dcb parameters file open error: %s\n", file);
        return 0;
    }
    while (fgets(buff, sizeof(buff), fp)) {
        if (strstr(buff, "DIFFERENTIAL (P1-P2) CODE BIASES")) {
            type = 1;
        } else if (strstr(buff, "DIFFERENTIAL (P1-C1) CODE BIASES")) {
            type = 2;
        } else if (strstr(buff, "DIFFERENTIAL (P2-C2) CODE BIASES")) {
            type = 3;
        }
 
        if (!type || sscanf(buff, "%s %s", str1, str2) < 1) {
            continue;
        }
 
        if ((cbias = str2num(buff, 26, 9)) == 0.0) {
            continue;
        }
 
        if (sta && (!strcmp(str1, "G") || !strcmp(str1, "R"))) { /* receiver DCB */
            for (i = 0; i < MAXRCV; i++) {
                if (!strcmp(sta[i].name, str2)) {
                    break;
                }
            }
            if (i < MAXRCV) {
                j = !strcmp(str1, "G") ? 0 : 1;
                nav->rbias[i][j][type - 1] = cbias * 1E-9 * CLIGHT; /* ns -> m */
            }
        } else if ((sat = satid2no(str1))) {                       /* satellite dcb */
            nav->cbias[sat - 1][type - 1] = cbias * 1E-9 * CLIGHT; /* ns -> m */
        }
    }
    fclose(fp);
 
    return 1;
}
/* read DCB parameters ---------------------------------------------------------
 * read differential code bias (DCB) parameters
 * args   : char   *file       I   DCB parameters file (wild-card * expanded)
 *          nav_t  *nav        IO  navigation data
 *          sta_t  *sta        I   station info data to inport receiver DCB
 *                                 (NULL: no use)
 * return : status (1:ok,0:error)
 * notes  : currently only support P1-P2, P1-C1, P2-C2, bias in DCB file
 *-----------------------------------------------------------------------------*/
extern int readdcb(const char* file, nav_t* nav, const sta_t* sta) {
    int i, j, n;
    char* efiles[MAXEXFILE] = {0};
 
    trace(NULL, 3, "readdcb : file=%s\n", file);
 
    for (i = 0; i < MAXSAT; i++) {
        for (j = 0; j < 3; j++) {
            nav->cbias[i][j] = 0.0;
        }
    }
    for (i = 0; i < MAXEXFILE; i++) {
        if (!(efiles[i] = (char*)malloc(1024))) {
            for (i--; i >= 0; i--) {
                free(efiles[i]);
            }
            return 0;
        }
    }
    n = expath(file, efiles, MAXEXFILE);
 
    for (i = 0; i < n; i++) {
        readdcbf(efiles[i], nav, sta);
    }
    for (i = 0; i < MAXEXFILE; i++) {
        free(efiles[i]);
    }
 
    return 1;
}
/* polynomial interpolation by Neville's algorithm ---------------------------*/
static double interppol(const double* x, double* y, int n) {
    int i, j;
 
    for (j = 1; j < n; j++) {
        for (i = 0; i < n - j; i++) {
            y[i] = (x[i + j] * y[i] - x[i] * y[i + 1]) / (x[i + j] - x[i]);
        }
    }
    return y[0];
}
/* satellite position by precise ephemeris -----------------------------------*/
static int pephpos(gtime_t time, int sat, const nav_t* nav, double* rs, double* dts, double* vare, double* varc) {
    double t[NMAX + 1], p[3][NMAX + 1], c[2], *pos, std = 0.0, s[3], sinl, cosl;
    int i, j, k, index;
 
    trace(NULL, 4, "pephpos : time=%s sat=%2d\n", time_str(time, 3), sat);
 
    rs[0] = rs[1] = rs[2] = dts[0] = 0.0;
 
    if (nav->ne < NMAX + 1 || timediff(time, nav->peph[0].time) < -MAXDTE ||
        timediff(time, nav->peph[nav->ne - 1].time) > MAXDTE) {
        trace(NULL, 3, "no prec ephem %s sat=%2d\n", time_str(time, 0), sat);
        return 0;
    }
    /* binary search */
    for (i = 0, j = nav->ne - 1; i < j;) {
        k = (i + j) / 2;
        if (timediff(nav->peph[k].time, time) < 0.0) {
            i = k + 1;
        } else {
            j = k;
        }
    }
    index = i <= 0 ? 0 : i - 1;
 
    /* polynomial interpolation for orbit */
    i = index - (NMAX + 1) / 2;
    if (i < 0) {
        i = 0;
    } else if (i + NMAX >= nav->ne) {
        i = nav->ne - NMAX - 1;
    }
 
    for (j = 0; j <= NMAX; j++) {
        t[j] = timediff(nav->peph[i + j].time, time);
        if (norm(nav->peph[i + j].pos[sat - 1], 3) <= 0.0) {
            trace(NULL, 3, "prec ephem outage %s sat=%2d\n", time_str(time, 0), sat);
            return 0;
        }
    }
    for (j = 0; j <= NMAX; j++) {
        pos = nav->peph[i + j].pos[sat - 1];
        /* correciton for earh rotation ver.2.4.0 */
        sinl = sin(OMGE * t[j]);
        cosl = cos(OMGE * t[j]);
        p[0][j] = cosl * pos[0] - sinl * pos[1];
        p[1][j] = sinl * pos[0] + cosl * pos[1];
        p[2][j] = pos[2];
    }
    for (i = 0; i < 3; i++) {
        rs[i] = interppol(t, p[i], NMAX + 1);
    }
    if (vare) {
        for (i = 0; i < 3; i++) {
            s[i] = nav->peph[index].std[sat - 1][i];
        }
        std = norm(s, 3);
 
        /* extrapolation error for orbit */
        if (t[0] > 0.0) {
            std += EXTERR_EPH * SQR(t[0]) / 2.0;
        } else if (t[NMAX] < 0.0) {
            std += EXTERR_EPH * SQR(t[NMAX]) / 2.0;
        }
        *vare = SQR(std);
    }
    /* linear interpolation for clock */
    t[0] = timediff(time, nav->peph[index].time);
    t[1] = timediff(time, nav->peph[index + 1].time);
    c[0] = nav->peph[index].pos[sat - 1][3];
    c[1] = nav->peph[index + 1].pos[sat - 1][3];
 
    if (t[0] <= 0.0) {
        if ((dts[0] = c[0]) != 0.0) {
            std = nav->peph[index].std[sat - 1][3] * CLIGHT - EXTERR_CLK * t[0];
        }
    } else if (t[1] >= 0.0) {
        if ((dts[0] = c[1]) != 0.0) {
            std = nav->peph[index + 1].std[sat - 1][3] * CLIGHT + EXTERR_CLK * t[1];
        }
    } else if (c[0] != 0.0 && c[1] != 0.0) {
        dts[0] = (c[1] * t[0] - c[0] * t[1]) / (t[0] - t[1]);
        i = t[0] < -t[1] ? 0 : 1;
        std = nav->peph[index + i].std[sat - 1][3] + EXTERR_CLK * fabs(t[i]);
    } else {
        dts[0] = 0.0;
    }
    if (varc) {
        *varc = SQR(std);
    }
    return 1;
}
/* satellite clock by precise clock ------------------------------------------*/
static int pephclk(gtime_t time, int sat, const nav_t* nav, double* dts, double* varc) {
    double t[2], c[2], std;
    int i, j, k, index;
 
    trace(NULL, 4, "pephclk : time=%s sat=%2d\n", time_str(time, 3), sat);
 
    if (nav->nc < 2 || timediff(time, nav->pclk[0].time) < -MAXDTE ||
        timediff(time, nav->pclk[nav->nc - 1].time) > MAXDTE) {
        trace(NULL, 3, "no prec clock %s sat=%2d\n", time_str(time, 0), sat);
        return 1;
    }
    /* binary search */
    for (i = 0, j = nav->nc - 1; i < j;) {
        k = (i + j) / 2;
        if (timediff(nav->pclk[k].time, time) < 0.0) {
            i = k + 1;
        } else {
            j = k;
        }
    }
    index = i <= 0 ? 0 : i - 1;
 
    /* linear interpolation for clock */
    t[0] = timediff(time, nav->pclk[index].time);
    t[1] = timediff(time, nav->pclk[index + 1].time);
    c[0] = nav->pclk[index].clk[sat - 1][0];
    c[1] = nav->pclk[index + 1].clk[sat - 1][0];
 
    if (t[0] <= 0.0) {
        if ((dts[0] = c[0]) == 0.0) {
            return 0;
        }
        std = nav->pclk[index].std[sat - 1][0] * CLIGHT - EXTERR_CLK * t[0];
    } else if (t[1] >= 0.0) {
        if ((dts[0] = c[1]) == 0.0) {
            return 0;
        }
        std = nav->pclk[index + 1].std[sat - 1][0] * CLIGHT + EXTERR_CLK * t[1];
    } else if (c[0] != 0.0 && c[1] != 0.0) {
        dts[0] = (c[1] * t[0] - c[0] * t[1]) / (t[0] - t[1]);
        i = t[0] < -t[1] ? 0 : 1;
        std = nav->pclk[index + i].std[sat - 1][0] * CLIGHT + EXTERR_CLK * fabs(t[i]);
    } else {
        trace(NULL, 3, "prec clock outage %s sat=%2d\n", time_str(time, 0), sat);
        return 0;
    }
    if (varc) {
        *varc = SQR(std);
    }
    return 1;
}
/* satellite antenna phase center offset ---------------------------------------
 * compute satellite antenna phase center offset in ecef
 * args   : gtime_t time       I   time (gpst)
 *          double *rs         I   satellite position and velocity (ecef)
 *                                 {x,y,z,vx,vy,vz} (m|m/s)
 *          int    sat         I   satellite number
 *          nav_t  *nav        I   navigation data
 *          double *dant       I   satellite antenna phase center offset (ecef)
 *                                 {dx,dy,dz} (m) (iono-free LC value)
 * return : none
 * notes  : iono-free LC frequencies defined as follows:
 *            GPS/QZSS : L1-L2
 *            GLONASS  : G1-G2
 *            Galileo  : E1-E5b
 *            BDS      : B1I-B2I
 *            NavIC    : L5-S
 *-----------------------------------------------------------------------------*/
extern void satantoff(gtime_t time, const double* rs, int sat, const nav_t* nav, double* dant) {
    const pcv_t* pcv = nav->pcvs + sat - 1;
    double ex[3], ey[3], ez[3], es[3], r[3], rsun[3], gmst, erpv[5] = {0}, freq[2];
    double C1, C2, dant1, dant2;
    int i, sys;
 
    trace(NULL, 4, "satantoff: time=%s sat=%2d\n", time_str(time, 3), sat);
 
    dant[0] = dant[1] = dant[2] = 0.0;
 
    /* sun position in ecef */
    sunmoonpos(gpst2utc(time), erpv, rsun, NULL, &gmst);
 
    /* unit vectors of satellite fixed coordinates */
    for (i = 0; i < 3; i++) {
        r[i] = -rs[i];
    }
    if (!normv3(r, ez)) {
        return;
    }
    for (i = 0; i < 3; i++) {
        r[i] = rsun[i] - rs[i];
    }
    if (!normv3(r, es)) {
        return;
    }
    cross3(ez, es, r);
    if (!normv3(r, ey)) {
        return;
    }
    cross3(ey, ez, ex);
 
    /* iono-free LC coefficients */
    sys = satsys_bd2(sat, NULL);
    freq[0] = freq_num2freq(sys, freq_idx2freq_num(sys, 0), 0);
    freq[1] = freq_num2freq(sys, freq_idx2freq_num(sys, 1), 0);
 
    if (sys == SYS_GLO) { /* GLONASS: use sat-dependent freq */
        freq[0] = sat2freq(sat, CODE_L1C, nav);
        freq[1] = sat2freq(sat, CODE_L2C, nav);
    }
    if (freq[0] == 0.0 || freq[1] == 0.0) {
        return;
    }
 
    C1 = SQR(freq[0]) / (SQR(freq[0]) - SQR(freq[1]));
    C2 = -SQR(freq[1]) / (SQR(freq[0]) - SQR(freq[1]));
 
    /* antenna indices for each frequency */
    {
        int a1 = freq_idx2ant_idx(sys, 0);
        int a2 = freq_idx2ant_idx(sys, 1);
 
        /* iono-free LC */
        for (i = 0; i < 3; i++) {
            dant1 = pcv->off[a1][0] * ex[i] + pcv->off[a1][1] * ey[i] + pcv->off[a1][2] * ez[i];
            dant2 = pcv->off[a2][0] * ex[i] + pcv->off[a2][1] * ey[i] + pcv->off[a2][2] * ez[i];
            dant[i] = C1 * dant1 + C2 * dant2;
        }
    }
}
/* satellite position/clock by precise ephemeris/clock -------------------------
 * compute satellite position/clock with precise ephemeris/clock
 * args   : gtime_t time       I   time (gpst)
 *          int    sat         I   satellite number
 *          nav_t  *nav        I   navigation data
 *          int    opt         I   sat position option
 *                                 (0: center of mass, 1: antenna phase center)
 *          double *rs         O   sat position and velocity (ecef)
 *                                 {x,y,z,vx,vy,vz} (m|m/s)
 *          double *dts        O   sat clock {bias,drift} (s|s/s)
 *          double *var        IO  sat position and clock error variance (m)
 *                                 (NULL: no output)
 * return : status (1:ok,0:error or data outage)
 * notes  : clock includes relativistic correction but does not contain code bias
 *          before calling the function, nav->peph, nav->ne, nav->pclk and
 *          nav->nc must be set by calling readsp3(), readrnx() or readrnxt()
 *          if precise clocks are not set, clocks in sp3 are used instead
 *-----------------------------------------------------------------------------*/
extern int peph2pos(gtime_t time, int sat, const nav_t* nav, int opt, double* rs, double* dts, double* var) {
    gtime_t time_tt;
    double rss[3], rst[3], dtss[1], dtst[1], dant[3] = {0}, vare = 0.0, varc = 0.0, tt = 1E-3;
    int i;
 
    trace(NULL, 4, "peph2pos: time=%s sat=%2d opt=%d\n", time_str(time, 3), sat, opt);
 
    if (sat <= 0 || MAXSAT < sat) {
        return 0;
    }
 
    /* satellite position and clock bias */
    if (!pephpos(time, sat, nav, rss, dtss, &vare, &varc) || !pephclk(time, sat, nav, dtss, &varc)) {
        return 0;
    }
 
    time_tt = timeadd(time, tt);
    if (!pephpos(time_tt, sat, nav, rst, dtst, NULL, NULL) || !pephclk(time_tt, sat, nav, dtst, NULL)) {
        return 0;
    }
 
    /* satellite antenna offset correction */
    if (opt) {
        satantoff(time, rss, sat, nav, dant);
    }
    for (i = 0; i < 3; i++) {
        rs[i] = rss[i] + dant[i];
        rs[i + 3] = (rst[i] - rss[i]) / tt;
    }
    /* relativistic effect correction */
    if (dtss[0] != 0.0) {
        dts[0] = dtss[0] - 2.0 * dot(rs, rs + 3, 3) / CLIGHT / CLIGHT;
        dts[1] = (dtst[0] - dtss[0]) / tt;
    } else { /* no precise clock */
        dts[0] = dts[1] = 0.0;
    }
    if (var) {
        *var = vare + varc;
    }
 
    return 1;
}