from datetime import datetime, timedelta
# from glob import glob
import gzip
# import os
from subprocess import run
from re import sub
from gnss.tools import dt2float, lagrange
import numpy as np
from gnss.tools import read_GNSS, gps_week
from pathlib import Path
from fnmatch import fnmatch
from gnss.satkit_ext import extend_gnss, update_satkit_data
from configuration import gnss_dir, gnss_download_url1, gnss_download_url2
EXT_HOURS = 24  # prediction appended after the last epoch of each product
EXT_METHOD = 'ecom'  # 'ecom' (ECOM SRP, 48 h fit) or 'cannonball' (satkit propagator, 12 h fit)
EXT_NPROC = 16  # worker processes for the orbit fits (capped at MAX_NPROC=16 in satkit_ext)


GEO_SPEED = 0.3   # km/s: below this Earth-fixed speed (BeiDou GEO) positions use a cubic spline, not Lagrange


def near_stationary(Gdata0, Gdt0):
    """Satellites whose Earth-fixed speed stays below GEO_SPEED (geostationary)."""
    t = np.array([(x - Gdt0[0]).total_seconds() for x in Gdt0])
    sp = np.linalg.norm(np.diff(Gdata0[:, :, :3], axis=0), axis=2) / np.diff(t)[:, None]
    with np.errstate(invalid='ignore'):
        return np.nanmax(np.where(np.isfinite(sp), sp, -1), axis=0) < GEO_SPEED


def interp30(Gdata0, Gdt0, Gdt, dt):
    """5-min -> 30-s: positions/velocities by 22-point Lagrange (as before), clocks/clock rates by a
    cubic spline (Lagrange amplifies clock noise to the microsecond level at the ends of the span).
    Near-stationary (GEO) satellites also get spline positions: their Earth-fixed orbit is smooth enough
    (<= 2 cm from Lagrange mid-span) and the spline has no metre-level oscillation at the span ends.
    Returns (len(Gdt), nsat, 8): x y z [km], clk [us], vx vy vz [dm/s], clock rate [1e-4 us/s]."""
    from scipy.interpolate import CubicSpline
    x0, x = dt2float(Gdt0, dt), dt2float(Gdt, dt)
    p, v = lagrange(x0 / 86400, Gdata0[:, :, :4].reshape((Gdata0.shape[0], -1)), x / 86400, 11)
    out = np.concatenate((p.reshape(p.shape[0], -1, 4), (v * 10000 / 86400).reshape(v.shape[0], -1, 4)), axis=2)
    geo = near_stationary(Gdata0, Gdt0)
    for j in range(Gdata0.shape[1]):
        cols = [(3, 7)] + ([(0, 4), (1, 5), (2, 6)] if geo[j] else [])
        for c, cv in cols:
            m = np.isfinite(Gdata0[:, j, c])
            out[:, j, [c, cv]] = np.nan
            if m.sum() >= 4:
                cs = CubicSpline(x0[m], Gdata0[m, j, c])
                inside = (x >= x0[m][0]) & (x <= x0[m][-1])
                out[inside, j, c] = cs(x[inside])
                out[inside, j, cv] = cs(x[inside], 1) * 1e4
    return out
# from matplotlib import pyplot as plt
# from scipy.interpolate import CubicSpline
# from frame_conversions import ro_ecef2eci, ro_eci2ecef
# from astropy.time import Time
# from astropy.coordinates import ITRS, GCRS, CartesianRepresentation
# import astropy.units as u

working_dir = Path(gnss_dir)
working_dir.mkdir(parents=True, exist_ok=True)
for folder in ('gnss5min', 'gnssLrg', 'log'):
    (working_dir / folder).mkdir(exist_ok=True)

logfile = working_dir / f"log/interp_gnss_{datetime.now():%Y%m%d}.txt"
def log_message(msg):
    with open(logfile, 'a') as fid:
        fid.write(f"{datetime.now():%H:%M:%S} {msg}\n")

log_message('- Downloading...')
tar_dirs = sorted(p for p in (working_dir / 'gnssLrg').iterdir()
                  if p.is_dir() and len(p.name) == 8 and p.name[4] == '.')
if tar_dirs:
    tar_sdmin = datetime.strptime(tar_dirs[-1].name, '%Y.%j') - timedelta(days=1)
    date_len = np.ceil((datetime.now() - tar_sdmin) / timedelta(days=1)).astype(int) + 1
    date_rng = [tar_sdmin + timedelta(days=x) for x in range(date_len)]
else:
    dt_today = datetime.today()
    date_rng = [dt_today + timedelta(days=x) for x in range(-5, 2)]

# dt_proc = datetime(2025, 11, 6)
# dt_file = [dt_proc + timedelta(hours=x) for x in range(-120, 300, 6)]
# dt_range = [dt_proc + timedelta(hours=x) for x in [-1,25]]

urls = np.hstack([[f'{gnss_download_url1}latest/ultra/GRG0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz',
        f'{gnss_download_url1}latest/ultra/COD0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz',
        f'{gnss_download_url1}latest/ultra/GFZ0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz',
        f'{gnss_download_url1}{gps_week(x):d}/WUM0MGXNRT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz'] for x in date_rng])

# urls = ['ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/final/GRG0OPSFIN_20253*00_01D_30S_CLK.CLK.gz',
#         'ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/final/GRG0OPSFIN_20253*00_01D_05M_ORB.SP3.gz',
#         'ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/final/COD0OPSFIN_20253*00_01D_05M_ORB.SP3.gz']

for ftps_url in urls:
    # Retain and reuse archived SP3 downloads older than 7 days.
    # Recent products still use wget -N so updated ultra-rapid files are retrieved.
    stamp = sub(r'.*_(\d{7})\d{4}_.*', r'\1', ftps_url)
    if stamp != ftps_url:
        product_day = datetime.strptime(stamp, '%Y%j')
        if datetime.now() - product_day > timedelta(days=7):
            archive = working_dir / 'gnss5min' / f'{product_day:%Y.%j}'
            pattern = ftps_url.rsplit('/', 1)[-1]
            if archive.is_dir() and any(fnmatch(f.name, pattern) for f in archive.glob('*.gz')):
                log_message(f'- Reusing archived SP3: {pattern}')
                continue
    command = [
        "wget",
        "-N",
        "--ftp-user=anonymous",  # Specify FTP username
        "--ftp-password=passwd",  # Specify FTP password
        ftps_url
    ]
    run(command, cwd=working_dir / 'gnss5min', capture_output=True)
urls = [f'{gnss_download_url2}GNSSOrb_{x-timedelta(days=2):%Y%j}0_withCLK.SP3' for x in date_rng]
for ftps_url in urls:
    command = ['wget',ftps_url]
    run(command, cwd=working_dir / 'gnss5min', capture_output=True)
    lfile = sub(r'.*/GNSSOrb_(\d{7})0_.*', r'GBM0MGXRAP_\g<1>0000_04D_05M_ORB.SP3', ftps_url)
    command = ['mv', sub(r'.*/', r'', ftps_url), lfile]
    run(command, cwd=working_dir / 'gnss5min', capture_output=True)
    command = ['gzip', lfile]
    run(command, cwd=working_dir / 'gnss5min', capture_output=True)
log_message('- Download is completed.')
update_satkit_data(log=log_message)  # EOP predictions are needed for the ITRF<->GCRF rotations
#
# for dt_roll in dt_file:
#     # https_url = f'https://isdc-data.gfz.de/gnss/products/ultra/w{(dt_roll-datetime(1980,1,6))//timedelta(weeks=1)}/GFZ0OPSULT_{dt_roll:%Y%j%H}00_02D_05M_ORB.SP3.gz'
#     # run(["wget","-N",https_url], cwd='./gnss5min/', capture_output=True)
#     https_url = f'https://isdc-data.gfz.de/gnss/products/final/w{(dt_roll - datetime(1980, 1, 6)) // timedelta(weeks=1)}/GFZ0OPSFIN_{dt_roll:%Y%j%H}00_01D_15M_ORB.SP3.gz'
#     run(["wget", "-N", https_url], cwd='./gnss5min/', capture_output=True)
#     https_url = f'https://isdc-data.gfz.de/gnss/products/final/w{(dt_roll - datetime(1980, 1, 6)) // timedelta(weeks=1)}/GFZ0OPSFIN_{dt_roll:%Y%j%H}00_01D_30S_CLK.CLK.gz'
#     run(["wget", "-N", https_url], cwd='./gnss5min/', capture_output=True)

gdir = working_dir / 'gnss5min'
idir = working_dir / 'gnssLrg'
gfiles = list(gdir.glob('*SP3.gz'))
gfiles.sort()
for gfile in gfiles:
    dt = datetime.strptime(sub(r'.*_(\d{7})\d{4}_.*',r'\1',gfile.name),'%Y%j')
    tdir = gdir / f'{dt:%Y.%j}'
    tdir.mkdir(exist_ok=True)
    ldir = idir / f'{dt:%Y.%j}'
    ldir.mkdir(exist_ok=True)
    tfile = tdir / gfile.with_suffix('.txt').name # sub(r'.*/(.*)\.gz',f'{tdir}/\\1.txt',gfile)
    ftoken = sub(r'(.*)0[A-Z]{6}(_\d{9}).*', r'\1\2', tfile.name)
    hr = int(ftoken[-2:])
    if hr%6 != 0:
        gfile.unlink()
        continue
    ofile = list(ldir.glob(f'{ftoken}_*.npz'))
    mdt = datetime.fromtimestamp(gfile.stat().st_mtime)
    if len(ofile) > 0:
        odt = datetime.strptime(sub(r'.*_c(\d+)\.npz',r'\1',ofile[-1].name), '%Y%j%H%M')
        if mdt - odt > timedelta(minutes=2):
            ofile[-1].unlink()
            log_message(f'- Updating {ftoken} ({odt:%m-%d %H:%M} >> {mdt:%H:%M})...')
        else:
            gfile.unlink()
            log_message(f'- Skipping {ftoken} ({odt:%m-%d %H:%M})...')
            continue
    else:
        log_message(f'- Creating {ftoken} ({mdt:%m-%d %H:%M})...')
    with gzip.open(gfile,'rb') as fin, open(tfile,'wb') as fout:
        fout.write(fin.read())
    [Gdata0, Gdt0, Gsats] = read_GNSS(tfile)
    sidx = np.all(np.isfinite(Gdata0[:,:,:4]),axis=(0,2))
    Gsats = Gsats[sidx]
    Gdata0 = Gdata0[:,sidx,:]
    tidx = np.argsort(Gdt0)
    sidx = np.argsort(Gsats)
    Gdt0 = Gdt0[tidx]
    Gsats = Gsats[sidx]
    Gdata0 = Gdata0[np.ix_(tidx, sidx)]
    g2file = gfile.parent / sub(r'ORB\.SP3',r'CLK.CLK',gfile.name)
    if g2file.is_file():
        t2file = tdir / g2file.with_suffix('.txt').name
        with gzip.open(g2file, 'rb') as fin, open(t2file, 'wb') as fout:
            fout.write(fin.read())
        for line in open(t2file,'r'):
            if line.startswith('AS'):
                tt = datetime.strptime(line[8:24],'%Y %m %d %H %M')
                ss = line[3:6]
                tidx = np.searchsorted(Gdt0, tt)
                sidx = np.searchsorted(Gsats, ss)
                if Gdt0[tidx] == tt and Gsats[sidx] == ss:
                    Gdata0[tidx,sidx,3] = float(line[40:59])*1E6
    # --- satkit: extend the product by EXT_HOURS of dynamic prediction ---
    Gdt0_orig, Gdata0_orig = Gdt0, Gdata0
    ext_info = {}
    try:
        Gdata0, Gdt0, ext_info = extend_gnss(Gdata0, Gdt0, Gsats, ext_hours=EXT_HOURS, method=EXT_METHOD,
                                             nproc=EXT_NPROC, log=log_message)
        log_message(f'  extended to {Gdt0[-1]:%m-%d %H:%M}, median fit rms '
                    f'{np.nanmedian(ext_info["fit_rms"]):.2f} m, models: '
                    + ', '.join(f'{m} {np.sum(ext_info["model"] == m)}' for m in ('ecom', 'cannonball')))
    except Exception as e:
        log_message(f'  satkit extension failed ({e}); saving original span only.')
    Gdt = np.arange(Gdt0[0],Gdt0[-1]+timedelta(seconds=30),timedelta(seconds=30)).astype(datetime)
    Gdata = interp30(Gdata0, Gdt0, Gdt, dt)
    # satellites without a valid extension: re-interpolate the original span so the NaN
    # extension does not leak into its last hour; they stay NaN after the original end
    bad = ~np.all(np.isfinite(Gdata0[:, :, :4]), axis=(0, 2))
    if np.any(bad):
        tidx = np.where(Gdt <= Gdt0_orig[-1])[0]
        Gdata[np.ix_(tidx, np.where(bad)[0])] = interp30(Gdata0_orig[:, bad], Gdt0_orig, Gdt[tidx], dt)
    # orbit: fitted ECOM orbit + spline of the residual (no end oscillation); Lagrange stays only for
    # satellites without an ECOM model
    dense = ext_info.get('dense')
    if dense is not None and np.array_equal(dense['t'], Gdt) and dense['ok'].any():
        ok = dense['ok']
        Gdata[:, ok, :3] = dense['pos'][:, ok]
        Gdata[:, ok, 4:7] = dense['vel'][:, ok]
        geo = near_stationary(Gdata0_orig, Gdt0_orig) & ~ok
        log_message(f'  30-s orbit from fitted model + residual spline for {ok.sum()} sats, '
                    f'cubic spline (GEO) for {geo.sum()}' + (f' ({" ".join(Gsats[geo])})' if geo.any() else '')
                    + f', Lagrange for {(~ok & ~geo).sum()}')
    ofile = ldir / f'{ftoken}_s{Gdt[0]:%Y%j%H%M}_e{Gdt[-1]:%Y%j%H%M}_c{mdt:%Y%j%H%M}.npz'
    np.savez(ofile, Gdata=Gdata, Gdt=Gdt, Gsats=Gsats,
             Gext_start=ext_info.get('ext_start', Gdt[-1] + timedelta(seconds=30)),
             Gext_fitrms=ext_info.get('fit_rms', np.full(len(Gsats), np.nan)),
             Gext_model=ext_info.get('model', np.full(len(Gsats), '')))
    tfile.unlink()
    gfile.replace(tfile.with_suffix('.gz'))
    if g2file.is_file():
        t2file.unlink()
        g2file.replace(t2file.with_suffix('.gz'))

log_message('- Progress is completed.')
# gdir = './gnss5min/'
# idir = './gnssFin/'
# gfiles = glob(gdir + '*SP3.gz')
# for gfile in gfiles:
#     dt = datetime.strptime(sub(r'.*_(\d{7})\d{4}_.*',r'\1',gfile),'%Y%j')
#     tdir = f'{gdir}{dt:%Y.%j}'
#     os.makedirs(tdir, exist_ok=True)
#     ldir = f'{idir}{dt:%Y.%j}'
#     os.makedirs(ldir, exist_ok=True)
#     tfile = sub(r'.*/(.*)\.gz',f'{tdir}/\\1.txt',gfile)
#     with gzip.open(gfile,'rb') as fin, open(tfile,'wb') as fout:
#         fout.write(fin.read())
#     [Gdata, Gdt, Gsats] = read_GNSS(tfile)
#     ofile = sub(r'.*/(.*)\.gz',f'{ldir}/\\1.npz',gfile)
#     np.savez(ofile, Gdata=Gdata, Gdt=Gdt, Gsats=Gsats)
#     os.remove(tfile)
#     os.replace(gfile, tfile[:-3]+'gz')

# for ss in ['COD','GRG','GFZ']:
#     Gdatam = np.empty((0,0,8))
#     Gdtm = np.empty((0,))
#     Gsats = np.empty((0,))
#     Gdt = np.array([dt_proc + timedelta(seconds=x) for x in range(0,86430,30)])
#     for dt_roll in dt_file:
#         # dt_now = None
#         gfile = f'./gnss5min/{dt_roll:%Y.%j}/{ss}0OPSULT_{dt_roll:%Y%j%H}00_02D_05M_ORB.SP3.txt'
#         if not os.path.exists(gfile):
#             continue
#         [Gdata0, Gdt0, Gsats0] = read_GNSS(gfile)
#         tidx = np.logical_and(Gdt0 >= dt_range[0], Gdt0 <= dt_range[1])
#         Gdata0 = Gdata0[tidx]
#         Gdt0 = Gdt0[tidx]
#         tidx = np.argsort(Gdt0)
#         sidx = np.argsort(Gsats0)
#         Gdt0 = Gdt0[tidx]
#         Gsats0 = Gsats0[sidx]
#         Gdata0 = Gdata0[np.ix_(tidx,sidx)]
#         sidx = ~np.isin(Gsats0, Gsats)
#         Gsats = np.hstack((Gsats, Gsats0[sidx]))
#         Gdatam = np.hstack((Gdatam, np.full((Gdatam.shape[0],np.sum(sidx),Gdatam.shape[2]),np.nan)))
#         tidx = ~np.isin(Gdt0, Gdtm)
#         Gdtm = np.hstack((Gdtm, Gdt0[tidx]))
#         Gdatam = np.vstack((Gdatam, np.full((np.sum(tidx),Gdatam.shape[1],Gdatam.shape[2]),np.nan)))
#         tidx = np.argsort(Gdtm)
#         sidx = np.argsort(Gsats)
#         Gdtm = Gdtm[tidx]
#         Gsats = Gsats[sidx]
#         Gdatam = Gdatam[np.ix_(tidx,sidx)]
#         Gdatam[np.ix_(np.isin(Gdtm, Gdt0),np.isin(Gsats, Gsats0))] = Gdata0
#
#     tidx = Gdt<Gdtm[-12]
#     Gdt = Gdt[tidx]
#     p, v = lagrange(dt2float(Gdtm,dt_proc)/86400, Gdatam[:,:,:4].reshape((Gdatam.shape[0],-1)), dt2float(Gdt,dt_proc)/86400,11)
#     Gdata = np.concatenate((p.reshape(p.shape[0],-1,4),(v*10000/86400).reshape(v.shape[0],-1,4)),axis=2)
#
#     np.savez(f'gnssOrb/GNSS_{dt_proc:%Y.%j}_{ss}.npz', Gdata=Gdata, Gdt=Gdt, Gsats=Gsats)

# with (np.load(f'gnssOrb/GNSS_{dt_proc:%Y.%j}_RT.npz', allow_pickle=True) as fid):
#     GdataR = fid['Gdata']
#     GdtR = fid['Gdt']
#     GsatsR = fid['Gsats']
# sidx = np.argsort(GsatsR)
# GsatsR = GsatsR[sidx]
# GdataR = GdataR[:,sidx,:]
# with (np.load(f'gnssOrb/GNSS_{dt_proc:%Y.%j}_NRT.npz', allow_pickle=True) as fid):
#     GdataN = fid['Gdata']
#     GdtN = fid['Gdt']
#     GsatsN = fid['Gsats']
# sidx = np.argsort(GsatsN)
# GsatsN = GsatsN[sidx]
# GdataN = GdataN[:,sidx,:]

# with (np.load(f'gnssOrb/GNSS_{dt_proc:%Y.%j}_GFZ.npz', allow_pickle=True) as fid):
#     GdataR = fid['Gdata']
#     GdtR = fid['Gdt']
#     GsatsR = fid['Gsats']
# with (np.load(f'gnssOrb/GNSS_{dt_proc:%Y.%j}_COD.npz', allow_pickle=True) as fid):
#     GdataC = fid['Gdata']
#     GdtC = fid['Gdt']
#     GsatsC = fid['Gsats']
# with (np.load(f'gnssOrb/GNSS_{dt_proc:%Y.%j}_GRG.npz', allow_pickle=True) as fid):
#     GdataG = fid['Gdata']
#     GdtG = fid['Gdt']
#     GsatsG = fid['Gsats']

# sidx = np.all(np.isfinite(GdataG),axis=(0,2))
# Gdata = GdataG[:,sidx,:]
# Gdt = GdtG
# Gsats = GsatsG[sidx]
# sidx = ~np.isin(GsatsR,Gsats)
# Gsats = np.hstack((Gsats, GsatsR[sidx]))
# Gdata = np.hstack((Gdata, GdataR[:,sidx]))
# sidx = ~np.isin(GsatsC,Gsats)
# Gsats = np.hstack((Gsats, GsatsC[sidx]))
# Gdata = np.hstack((Gdata, GdataC[:,sidx]))
# sidx = ~np.isin(GsatsG,Gsats)
# Gsats = np.hstack((Gsats, GsatsG[sidx]))
# Gdata = np.hstack((Gdata, GdataG[:,sidx]))
# np.savez(f'gnssOrb/GNSS_{dt_proc:%Y.%j}_CMB.npz', Gdata=Gdata, Gdt=Gdt, Gsats=Gsats)

# for ctype in ['CMB']:
#     with (np.load(f'gnssOrb/GNSS_{dt_proc:%Y.%j}_{ctype}.npz', allow_pickle=True) as fid):
#         GdataC = fid['Gdata']
#         GdtC = fid['Gdt']
#         GsatsC = fid['Gsats']
#     sidx = np.argsort(GsatsC)
#     GsatsC = GsatsC[sidx]
#     GdataC = GdataC[:,sidx,:]
#
#     ncidx = np.isin(GsatsN,GsatsC)
#     cnidx = np.isin(GsatsC,GsatsN)
#     # ngidx = np.isin(GsatsN,GsatsG)
#     # gnidx = np.isin(GsatsG,GsatsN)
#     # nridx = np.isin(GsatsN,GsatsR)
#     # rnidx = np.isin(GsatsR,GsatsN)
#
#     plt.close('all')
#     xyz = 'rpy'
#     clrs = np.array([[int(y['color'][x:x+2],16)/255 for x in range(1,7,2)] for y in plt.rcParams['axes.prop_cycle']])
#     fig, axs = plt.subplots(3, 1, figsize=(12, 9))
#     xc = np.where(ncidx)[0]
#     yc = GdataC[:, cnidx, :] - GdataN[:, ncidx, :]
#     ey = GdataN[:, ncidx, :3]
#     ey /= np.sum(ey**2, axis=2)[:,:,None]**0.5
#     ep = np.cross(ey, GdataN[:,ncidx,4:7])
#     ep /= np.sum(ep**2, axis=2)[:,:,None]**0.5
#     er = np.cross(ep, ey)
#     yc = np.concatenate((np.sum((yc[:,:,:3]*er,yc[:,:,:3]*ep,yc[:,:,:3]*ey),axis=3),
#                          yc[None,:,:,3],
#                          np.sum((yc[:,:,4:7]*er,yc[:,:,4:7]*ep,yc[:,:,4:7]*ey),axis=3)),axis=0)
#     # xg = np.where(ngidx)[0]
#     # yg = GdataG[:, gnidx, :] - GdataN[:, ngidx, :]
#     # xr = np.where(nridx)[0]
#     # yr = GdataR[:, rnidx, :] - GdataN[:, nridx, :]
#
#     # axs[0].plot(np.tile(xc, (2, 1)), [np.nanmin(yc[0], axis=0), np.nanmax(yc[0], axis=0)],
#     #                 color=(clrs[0] + 2) / 3, linewidth=3)
#     # axs[1].plot(np.tile(xc, (2, 1)), [np.nanmin(yc[1], axis=0), np.nanmax(yc[1], axis=0)],
#     #                 color=(clrs[0] + 2) / 3, linewidth=3)
#     # axs[2].plot(np.tile(xc, (2, 1)), [np.nanmin(yc[2], axis=0), np.nanmax(yc[2], axis=0)],
#     #                 color=(clrs[0] + 2) / 3, linewidth=3)
#     # axs[0].plot(xc, np.nanmean(yc[0], axis=0), '.-', color=clrs[0], label='CMB')
#     # axs[1].plot(xc, np.nanmean(yc[1], axis=0), '.-', color=clrs[0])
#     # axs[2].plot(xc, np.nanmean(yc[2], axis=0), '.-', color=clrs[0])
#
#     for ii in range(3):
#         axs[0].plot(np.tile(xc, (2, 1)), [np.min(yc[ii], axis=0), np.max(yc[ii], axis=0)],
#                     color=(clrs[ii] + 2) / 3, linewidth=3)
#         # axs[0].plot(np.tile(xg, (2, 1)), [np.min(yg[ii], axis=0), np.max(yg[ii], axis=0)],
#         #             color=(clrs[ii + 3] + 2) / 3, linewidth=3)
#         # axs[0].plot(np.tile(xr, (2, 1)), [np.min(yr[ii], axis=0), np.max(yr[ii], axis=0)],
#         #             color=(clrs[ii + 6] + 2) / 3, linewidth=3)
#         axs[2].plot(np.tile(xc, (2, 1)), [np.min(yc[ii + 4], axis=0), np.max(yc[ii + 4], axis=0)],
#                     color=(clrs[ii] + 2) / 3, linewidth=3)
#         # axs[2].plot(np.tile(xg, (2, 1)), [np.min(yg[ii + 4], axis=0), np.max(yg[ii + 4], axis=0)],
#         #             color=(clrs[ii + 3] + 2) / 3, linewidth=3)
#         # axs[2].plot(np.tile(xr, (2, 1)), [np.min(yr[ii + 4], axis=0), np.max(yr[ii + 4], axis=0)],
#         #             color=(clrs[ii + 6] + 2) / 3, linewidth=3)
#     axs[1].plot(np.tile(xc, (2, 1)), [np.nanmin(yc[3], axis=0), np.nanmax(yc[3], axis=0)],
#                     color=(clrs[0] + 2) / 3, linewidth=3)
#     # axs[1].plot(np.tile(xg, (2, 1)), [np.nanmin(yg[3], axis=0), np.nanmax(yg[3], axis=0)],
#     #                 color=(clrs[3] + 2) / 3, linewidth=3)
#     # axs[1].plot(np.tile(xr, (2, 1)), [np.nanmin(yr[3], axis=0), np.nanmax(yr[3], axis=0)],
#     #                 color=(clrs[6] + 2) / 3, linewidth=3)
#     for ii in range(3):
#         axs[0].plot(xc, np.mean(yc[ii], axis=0), color=clrs[ii], label=f'{ctype} {xyz[ii]}')
#         # axs[0].plot(xg, np.mean(yg[ii], axis=0), color=clrs[ii + 3], label=f'GRG {xyz[ii]}')
#         # axs[0].plot(xr, np.mean(yr[ii], axis=0), color=clrs[ii + 6], label=f'GFZ {xyz[ii]}')
#         axs[2].plot(xc, np.mean(yc[ii + 4], axis=0), color=clrs[ii])
#         # axs[2].plot(xg, np.mean(yg[i + 4], axis=0), color=clrs[ii + 3])
#         # axs[2].plot(xr, np.mean(yr[ii + 4], axis=0), color=clrs[ii + 6])
#     axs[1].plot(xc, np.nanmean(yc[3], axis=0), '.-', color=clrs[0])
#     # axs[1].plot(xg, np.nanmean(yg[3], axis=0), '.-', color=clrs[3])
#     # axs[1].plot(xr, np.nanmean(yr[3], axis=0), '.-', color=clrs[6])
#
#     axs[0].set_xlim((-0.5, len(ncidx) - 0.5))
#     axs[0].set_xticks(np.arange(len(ncidx)))
#     axs[0].set_xticklabels(GsatsN)
#     axs[0].set_ylim((-0.0008,0.0008))
#     axs[0].tick_params(axis='x', labelrotation=90, labelsize=8)
#     axs[0].set_ylabel('Position Distance (km)')
#     axs[0].legend(ncol=3, loc='upper right')
#     axs[0].set_title(f'COD/GRG/GFZ-NRT difference in {dt_proc:%Y-%m-%d} ({dt_proc:%Y.%j})')
#     axs[1].set_xlim((-0.5, len(ncidx) - 0.5))
#     axs[1].set_xticks(np.arange(len(ncidx)))
#     axs[1].set_xticklabels(GsatsN)
#     axs[1].set_ylim((-0.02,0.02))
#     axs[1].tick_params(axis='x', labelrotation=90, labelsize=8)
#     axs[1].set_ylabel('Clock Bias Difference (us)')
#     axs[2].set_xlim((-0.5, len(ncidx) - 0.5))
#     axs[2].set_xticks(np.arange(len(ncidx)))
#     axs[2].set_xticklabels(GsatsN)
#     axs[2].set_ylim((-0.008,0.008))
#     axs[2].tick_params(axis='x', labelrotation=90, labelsize=8)
#     axs[2].set_ylabel('Velocity Distance (dm/s)')
#     plt.savefig(f'{dt_proc:%Y.%j}_{ctype}.png',dpi=300)