import os
from datetime import datetime, timedelta
from re import sub

def load_dcb(dt_ref):
    import pickle
    dcbdir = '../gnss/gnssDCB/'
    dcbfile = dcbdir+f'{dt_ref:%Y-%m}_DCB.pkl'
    if os.path.exists(dcbfile):
        with open(dcbfile, 'rb') as fid:
            dcb = pickle.load(fid)
    else:
        from subprocess import run
        from re import search
        dcb = {'P1C1': {}, 'P2C2': {}, 'P1P2': {}}
        if datetime.now() - dt_ref < timedelta(days=32) and datetime.now().month == dt_ref.month:
            zfiles = ["P1C1_RINEX.DCB",
                      "P2C2_RINEX.DCB",
                      "P1P2_ALL.DCB"]
            rmtdir = "http://ftp.aiub.unibe.ch/CODE/"

        else:
            zfiles = [f"P1C1{dt_ref:%y%m}_RINEX.DCB.Z",
                      f"P2C2{dt_ref:%y%m}_RINEX.DCB.Z",
                      f"P1P2{dt_ref:%y%m}_ALL.DCB.Z"]
            rmtdir = f"http://ftp.aiub.unibe.ch/CODE/{dt_ref:%Y}/"
        for zfile in zfiles:
            run(["wget", rmtdir + zfile], cwd=dcbdir, capture_output=True)
            tfile = dcbdir + sub(r'\.Z$', r'', zfile)
            if zfile[-1]=='Z':
                run(["gunzip", zfile], cwd=dcbdir, capture_output=True)
            for line in open(tfile, 'r'):
                if len(line) > 47 and search('[GR][0-3][0-9]', line[:3]):
                    dcb[zfile[:4]][line[:3]] = [float(line[ll:ll + 9].strip()) for ll in [26, 38]]
            os.remove(tfile)
        with open(dcbfile, 'wb') as fid:
            pickle.dump(dcb, fid)

def calculate_tec(reffile):
    import numpy as np
    from scipy.interpolate import CubicSpline
    from frame_conversions import ro_ecef2eci, ro_eci2ecef
    from tools import dt2float, get_sat_freq, read_gnss_v2, conv_gnss_ant, phase_correction, remove_cycleslip
    # Constants
    kappa = 40.30819294518825
    dt0_gps = datetime(1980,1,6)
    antennas = ['POD0', 'POD1', 'OCC2', 'OCC3', 'OCC4', 'OCC5', 'OCC6', 'OCC7']
    # Choose reference file
    dt_ref = [datetime.strptime(x,'%Y%m%d%H%M%S') for x in sub(r'.*(\d{14}\.\d{14}).*',r'\1',reffile).split('.')]
    dt_rr = datetime.strptime(sub(r'.*\.(\d{4}\.\d{3})\..*',r'\1',reffile),'%Y.%j')
    sec_rr = (dt_rr - dt0_gps) / timedelta(seconds=1)
    with np.load(reffile, allow_pickle=True) as fid:
        data_ref = fid['data']
    yy = data_ref[:, 3] - data_ref[:, 2]
    if max(yy)-min(yy)>1000:
        raise ValueError('L1 L2 difference out of range!')

    # Extract info from file names
    leoName = sub(r'.*occTab_(....)\..*',r'\1',reffile)
    ant_leo_ref = int(sub(r'.*\.A(\d\d)\..*',r'\1',reffile))
    gnss_ref = sub(r'.*\.(...).npz', r'\1', reffile)
    f1_ref, f2_ref, gprn_ref = get_sat_freq(gnss_ref)
    c2_ref = 1 / ((f1_ref / f2_ref) ** 2 - 1)
    c1_ref = c2_ref + 1

    # ant_gps = np.loadtxt("GPS_Antenna_Infos.Final")
    # ant_gps_start_time = [datetime(*map(int, x[1:7])) for x in ant_gps]
    # ant_gps_end_time = [datetime(*map(int, x[7:13])) for x in ant_gps]
    with np.load('../gnss/GPS_Antenna_Infos.npz', allow_pickle=True) as fid:
        antoffset_gnss_ref = fid[gnss_ref][()]
    antoffset_gnss_ref = antoffset_gnss_ref['antenna'][np.where(np.logical_and(
                         antoffset_gnss_ref['dt_rng'][:, 0] < dt_rr,
                         antoffset_gnss_ref['dt_rng'][:, 1] > dt_rr))[0][0]][:3] / 1E3
    # Read in gnss position/velocity file accordingly
    # Gdata, Gdt, Gsats = read_gnss(dt_ref[1] + timedelta(minutes=10),
    #                               (dt_ref[0] - timedelta(minutes=10)).day != (dt_ref[1] + timedelta(minutes=10)).day)
    # if gnss_ro not in Gsats:
    #     raise LookupError(f'RO satellite {gnss_ro} not in Gsats')
    # if gnss_ref not in Gsats:
    #     raise LookupError(f'REF satellite {gnss_ref} not in Gsats')
    # gidx = np.hstack([np.where(Gsats == gnss_ro), np.where(Gsats == gnss_ref)])[0]
    # Gdata = Gdata[:,gidx,:]
    Gdata, Gdt = read_gnss_v2(dt_ref, [gnss_ref])
    tidx = np.logical_and(Gdt > dt_ref[0] - timedelta(minutes=10), Gdt < dt_ref[1] + timedelta(minutes=10))
    dts_gnss = Gdt[tidx]
    sec_gnss = dt2float(dts_gnss, dt_rr)
    pos_gnss = Gdata[tidx,:,:4]
    vel_gnss = Gdata[tidx,:,4:]
    del Gdata, Gdt
    # Read in leo position/velocity/attitude files accordingly
    with np.load(f"s01_leoOrb/leoOrb_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid:
        Ldata = fid['data']
        Ldts = fid['dts']
    tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1))
    dts_leoorb = Ldts[tidx]
    pos_leo = Ldata[tidx,:4]
    vel_leo = Ldata[tidx,4:]
    del Ldata, Ldts
    with np.load(f"s01_leoAtt/leoAtt_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid:
        Ldata = fid['data']
        Ldts = fid['dts']
        Lvar = fid['vars']
    tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1))
    dts_leoatt = Ldts[tidx]
    att_leo = Ldata[tidx,:][:,['att' in x for x in Lvar]]
    aridx = np.cumsum(np.hstack(([False], np.any(np.abs(np.diff(att_leo, axis=0)) > 0.1, axis=1)))) % 2 == 1
    att_leo[aridx] *= -1
    del Ldata, Ldts, Lvar
    tidx = np.logical_and(dts_leoorb >= dts_leoatt[0], dts_leoorb <= dts_leoatt[-1])
    dts_leo = dts_leoorb[tidx]
    sec_leo = dt2float(dts_leo, dt_rr)
    pos_leo = pos_leo[tidx]
    vel_leo = vel_leo[tidx]
    att_leo = CubicSpline(dt2float(dts_leoatt,dts_leo[0]), att_leo)(dt2float(dts_leo,dts_leo[0]))
    att_leo /= np.sum(att_leo ** 2, axis=1)[:,None] ** 0.5
    del dts_leoorb, dts_leoatt
    tsmp = 1 # 30
    tspc = np.round(np.linspace(0,len(sec_leo)-1,min(max(-((len(sec_leo)-1)//-tsmp)+1,20),len(sec_leo)))).astype(int)
    dts_leo = dts_leo[tspc]
    sec_leo = sec_leo[tspc]
    pos_leo = pos_leo[tspc]
    vel_leo = vel_leo[tspc]
    att_leo = att_leo[tspc]
    # Load ro/ref files
    with np.load(reffile, allow_pickle=True) as fid:
        data_ref = fid['data']
        header_ref = fid['header']
    l_ref = np.where([x[0] == 'L' for x in header_ref])[0]
    s_ref = np.where([x[0] == 'S' for x in header_ref])[0]
    c_ref = np.where([x[0] == 'C' for x in header_ref])[0]
    sec_ref = data_ref[:, 0] - sec_rr + data_ref[:, 1]
    if sec_leo[0]>sec_ref[0] or sec_leo[-1]<sec_ref[-1]:
        raise ValueError(f'LEO time ({sec_leo[0]:.1f} >> {sec_leo[-1]:.1f}) could not cover REF event ({sec_ref[0]:.1f} >> {sec_ref[-1]:.1f})!')
    sec_ref -= CubicSpline(sec_leo, pos_leo[:,-1])(sec_ref)/1E6
    if ~np.all(np.isfinite(pos_gnss[:, 0, 3])):
        raise ValueError(f'{np.sum(~np.isfinite(pos_gnss[:, 0, 3]))}/{sec_gnss.size} clock bias value for {gnss_ref} is not finite.')

    ant_ref_ecef = conv_gnss_ant(pos_gnss[:, 0, :3], antoffset_gnss_ref, dts_gnss)
    pos_leo_eci, vel_leo_eci = ro_ecef2eci(pos_leo, vel_leo, dts_leo)
    pos_ref_eci, vel_ref_eci = ro_ecef2eci(pos_gnss[:, 0, :3] + ant_ref_ecef, vel_gnss[:, 0, :], dts_gnss)
    pos_leo_eci = np.hstack([pos_leo_eci, pos_leo[:,[-1]]])
    vel_leo_eci = np.hstack([vel_leo_eci, pos_leo[:, [-1]]])
    pos_ref_eci = np.hstack([pos_ref_eci, pos_gnss[:, 0, [3]]])
    vel_ref_eci = np.hstack([vel_ref_eci, vel_gnss[:, 0, [3]]])
    with np.load('antenna_offset.npz') as fid:
        antoffset_leo_ref = fid[f'ANT_{leoName}_{antennas[ant_leo_ref]}'][12:15]
    obj_gnss_ref, obj_leo_ref = phase_correction(data_ref[:, l_ref[0]], pos_ref_eci, vel_ref_eci, sec_gnss, pos_leo_eci,
                                         vel_leo_eci, att_leo, sec_leo, sec_ref, antoffset_leo_ref, True)
    phL1_ref = data_ref[:, l_ref[0]] - obj_gnss_ref['range_diff_T'] + obj_gnss_ref['GL_term']
    # coef = np.polyfit(sec_ref, phL1_ref, 1)
    # if leoName.startswith('YM') and (np.abs(coef[0]) > 1 or np.max(np.abs(phL1_ref - np.polyval(coef, sec_ref))) > 10):
    #     raise ValueError('L1 non-linear!')
    phL2_ref = data_ref[:, l_ref[1]] - obj_gnss_ref['range_diff_T'] + obj_gnss_ref['GL_term']
    dph_ref = phL1_ref - phL2_ref
    # dph_ref = np.hstack(
    #     [np.cumsum(dph_ref[:49])[::2] / np.arange(1, 51, 2), np.convolve(dph_ref, np.ones((51,)), 'valid') / 51,
    #      (np.cumsum(dph_ref[:-50:-1])[::2] / np.arange(1, 51, 2))[::-1]])
    termL1_ref = phL1_ref + c2_ref * dph_ref
    termL2_ref = phL2_ref + c1_ref * dph_ref
    # termL1_ro = obj_leo_ro['clk_m'] + obj_leo_ro['delt_t_rel_r'] * c0 # Alternative method
    # termL2_ro = obj_leo_ro['clk_m'] + obj_leo_ro['delt_t_rel_r'] * c0
    exL1 = phL1_ref # - termL1_ref
    exL2 = phL2_ref # - termL2_ref
    OL = obj_leo_ref['pos']
    OG = obj_gnss_ref['pos']
    vec_lg = OG - OL
    proj_lg = OL - vec_lg * (np.sum(OL * vec_lg, axis=1) / np.sum(vec_lg * vec_lg, axis=1))[:, None]
    from astropy.coordinates import EarthLocation
    altp_km = EarthLocation.from_geocentric(*proj_lg.T, unit='m').geodetic.height.value / 1E3
    altl_km = EarthLocation.from_geocentric(*OL.T, unit='m').geodetic.height.value / 1E3
    newZ = proj_lg[0,:].copy()
    newY = np.cross(OL[0,:],OG[0,:])
    newX = np.cross(newY,newZ)
    newX /= sum(newX ** 2) ** .5
    newY /= sum(newY ** 2) ** .5
    newZ /= sum(newZ ** 2) ** .5
    mrot = np.array([newX,newY,newZ])
    # cosPL = np.sum(proj_lg*OL,axis=1)/np.sum(OL*OL,axis=1)**.5/np.sum(proj_lg*proj_lg,axis=1)**.5
    # cosPG = np.sum(proj_lg*OG,axis=1)/np.sum(OG*OG,axis=1)**.5/np.sum(proj_lg*proj_lg,axis=1)**.5
    # RL = (np.sum(OL ** 2, axis=1) ** .5)[:, None] * np.array([-(1 - cosPL ** 2) ** .5, cosPL]).T / 6378131
    # RG = (np.sum(OG ** 2, axis=1) ** .5)[:, None] * np.array([ (1 - cosPG ** 2) ** .5, cosPG]).T / 6378131
    RL = (mrot@OL.T).T/6378137
    RG = (mrot@OG.T).T/6378137
    # exL1 -= exL1[0 if obj_leo_ro['ocd'] else -1]
    # exL2 -= exL2[0 if obj_leo_ro['ocd'] else -1]
    # exL1_LM = data_ro[:,m_ro[0]] - (data_ro[:,l_ro[0]] - exL1)
    # exL1 = remove_cycleslip(exL1, exL1_LM, f1_ro, gnss_ro[0] + '1', obj_leo_ro['ocd'])
    # exL2 = remove_cycleslip(exL2, exL1_LM, f2_ro, gnss_ro[0] + '2', obj_leo_ro['ocd'])
    # # rofile_tkn = sub(r'.*_(.*).npz',r'\1',rofile)
    # fortran_input = f'/dev/shm/tcliu/in_{rofile_tkn}.txt'
    # fortran_output = f'/dev/shm/tcliu/out_{rofile_tkn}.txt'
    # with open(fortran_input, 'w') as fid:
    #     fid.write(f'{len(sec_ro)} {1 if obj_leo_ro["ocd"] else -1} {int(f1_ro)} {int(f2_ro)}\n')
    #     for n in range(len(sec_ro)):
    #         fid.write(f"{sec_ro[n]:15.8f}{exL1[n]:20.8f}{exL2[n]:20.8f}{data_ro[n,m_ro[0]]:20.8f}"
    #                   f"{pos_leo_ecef[n,0]:15.5f}{pos_leo_ecef[n,1]:15.5f}{pos_leo_ecef[n,2]:15.5f}"
    #                   f"{pos_gnss_ecef[n,0]:15.5f}{pos_gnss_ecef[n,1]:15.5f}{pos_gnss_ecef[n,2]:15.5f}3{dt_rr.month:2d}\n")
    # run(['PhaseModel/cyclesliperemoval', '-i', fortran_input, '-o', fortran_output], stdout=-3, stderr=-3, timeout=5)
    # from os import remove, makedirs
    # remove(fortran_input)
    # fortran_data = np.loadtxt(fortran_output)
    # remove(fortran_output)
    # if np.max(np.abs(fortran_data[:, 0] - exL2)) < np.max(np.abs(exL1 - exL2)):
    #     exL1 = fortran_data[:, 0]
    #     exL1_LM = fortran_data[:, 2]
    idx = np.where(np.isfinite(exL1))[0]
    exL1 -= exL1[idx[0 if obj_leo_ref['ocd'] else -1]]
    exL2 -= exL2[idx[0 if obj_leo_ref['ocd'] else -1]]
    exLC = exL1 + c2_ref * (exL1 - exL2)

    from matplotlib import pyplot as plt
    plt.close('all')
    fig, axs = plt.subplots(3, 2, figsize=(12, 9))
    axs[0,0].plot(sec_ref, data_ref[:,l_ref[0]], '.', markeredgewidth=0)
    axs[0,0].plot(sec_ref, data_ref[:,l_ref[1]], '.', markeredgewidth=0)
    axs[0,0].grid(True)
    axs[0,0].legend(['L1', 'L2'])
    axs[0,0].set_ylabel('L1 L2')
    axs[1,0].plot(sec_ref, exL1, '.', markeredgewidth=0)
    axs[1,0].plot(sec_ref, exL2, '.', markeredgewidth=0)
    axs[1,0].grid(True)
    axs[1,0].legend(['exL1','exL2'])
    axs[1,0].set_ylabel('Excess Phase')
    idx = np.sum(vec_lg*OL,axis=1)>=0
    axs[2,0].plot(sec_ref, altp_km, '.', markeredgewidth=0)
    axs[2,0].plot(sec_ref[idx], altl_km[idx], '.', markeredgewidth=0)
    axs[2,0].grid(True)
    axs[2,0].set_ylabel('Direct Height')
    idxL = np.argmax(np.abs(RL[:, 1]))
    idxG = np.argmax(np.abs(RG[:, 1]))
    if abs(RL[idxL,1])>abs(RG[idxG,1]):
        x0, y0, x1, y1 = RL[0, 1], 25, RL[idxL, 1], 1 if RL[idxL, 1] > 0 else 50
    else:
        x0, y0, x1, y1 = RL[0, 1], 25, RG[idxG, 1], 1 if RG[idxG, 1] > 0 else 50
    sca = axs[0,1].scatter(RL[:, 0], RL[:, 2], s=(RL[:, 1] - x0) / (x1 - x0) * (y1 - y0) + y0, c=sec_ref, edgecolors='none')
    axs[0,1].scatter(RG[:, 0], RG[:, 2], s=(RG[:, 1] - x0) / (x1 - x0) * (y1 - y0) + y0, c=sec_ref,edgecolors='none')
    axs[0,1].add_patch(plt.Circle((0, 0), 1, color='k', fill=False))
    axs[0,1].set_aspect('equal')
    axs[0,1].grid(True)
    axs[0,1].set_title('XZ View')
    # fig.colorbar(sca, ax=axs[0, 1], location='bottom')
    x0, y0, x1, y1 = RL[0, 2], 25, np.min(RL[:, 2]), 1
    sca=axs[1,1].scatter(RL[:, 0], RL[:, 1], s=(RL[:, 2] - x0) / (x1 - x0) * (y1 - y0) + y0, c=sec_ref,edgecolors='none')
    axs[1,1].scatter(RG[:, 0], RG[:, 1], s=(RG[:, 2] - x0) / (x1 - x0) * (y1 - y0) + y0, c=sec_ref,edgecolors='none')
    axs[1,1].add_patch(plt.Circle((0, 0), 1, color='k', fill=False))
    axs[1,1].set_aspect('equal')
    axs[1,1].grid(True)
    axs[1,1].set_title('XY View')
    # fig.colorbar(sca, ax=axs[1, 1], location='bottom')
    x0, y0, x1, y1 = RL[0, 0], 25, RG[0, 0], 50
    sca=axs[2,1].scatter(RL[:, 1], RL[:, 2], s=(RL[:, 0] - x0) / (x1 - x0) * (y1 - y0) + y0, c=sec_ref,edgecolors='none')
    axs[2,1].scatter(RG[:, 1], RG[:, 2], s=(RG[:, 0] - x0) / (x1 - x0) * (y1 - y0) + y0, c=sec_ref,edgecolors='none')
    axs[2,1].add_patch(plt.Circle((0, 0), 1, color='k', fill=False))
    axs[2,1].set_aspect('equal')
    axs[2,1].grid(True)
    axs[2,1].set_title('YZ View')
    fig.colorbar(sca, ax=axs[2, 1])
    fig.suptitle(sub(r'.*/', r'', reffile))
    fig.savefig(sub(r'.*/(.*)\.npz', r'\1.png', reffile), dpi=300)
    plt.close('all')

    Pdiff = data_ref[:,l_ref[0]] - data_ref[:,l_ref[1]]
    fr_ref = 2 / ((f1_ref / f2_ref) ** 2 - 1) + 1
    S = data_ref[:, s_ref] ** 2
    S /= np.sum(S, axis=0)[None,:]
    D1 = data_ref[:, c_ref[0]] - fr_ref * data_ref[:, l_ref[0]] + (fr_ref - 1) * data_ref[:, l_ref[1]]
    D2 = data_ref[:, c_ref[1]] - (fr_ref + 1) * data_ref[:, l_ref[0]] + fr_ref * data_ref[:, l_ref[1]]
    Ddiff = np.sum(D1 * S[:, 0]) - np.sum(D2 * S[:, 1])
    TEC = (Pdiff - Ddiff) / (f2_ref ** -2 - f1_ref ** -2) / kappa
    dcb = load_dcb(dt_ref[0])

    from matplotlib import pyplot as plt
    plt.close('all')
    fig, axs = plt.subplots(3, 3, figsize=(12, 9))
    axs[0,0].plot(data_ref[:, 1], data_ref[:, c_ref[0]], '.', markeredgewidth=0)
    axs[0,0].set_title('C1')
    axs[0,1].plot(data_ref[:, 1], data_ref[:, c_ref[1]], '.', markeredgewidth=0)
    axs[0,1].set_title('C2')
    axs[0,2].plot(data_ref[:, 1], -np.diff(data_ref[:, c_ref], axis=1), '.', markeredgewidth=0)
    axs[0,2].set_title('C1 - C2')
    axs[1,0].plot(data_ref[:, 1], data_ref[:, l_ref[0]], '.', markeredgewidth=0)
    axs[1,0].set_title('L1')
    axs[1,1].plot(data_ref[:, 1], data_ref[:, l_ref[1]], '.', markeredgewidth=0)
    axs[1,1].set_title('L2')
    axs[1,2].plot(data_ref[:, 1], -np.diff(data_ref[:, l_ref], axis=1), '.', markeredgewidth=0)
    axs[1,2].set_title('L1 - L2')
    axs[2,0].plot(data_ref[:, 1], data_ref[:, c_ref[0]] - data_ref[:, l_ref[0]], '.', markeredgewidth=0)
    axs[2,0].set_title('C1 - L1')
    axs[2,1].plot(data_ref[:, 1], data_ref[:, c_ref[1]] - data_ref[:, l_ref[1]], '.', markeredgewidth=0)
    axs[2,1].set_title('C2 - L2')
    # axs[2,2].plot(data_ref[:, 1], data_ref[:, s_ref]**2/np.sum(data_ref[:, s_ref]**2,axis=0)[None,:], '.', markeredgewidth=0)
    # axs[2,2].set_title('normalized S1² S2²')
    axs[2,2].plot(data_ref[:, 1], TEC/1E16+40, '.', markeredgewidth=0)
    axs[2,2].set_title('TEC+40')
    for ax in axs.ravel():
        ax.grid(True)
    fig.suptitle(sub(r'.*/', r'', reffile))
    fig.savefig(sub(r'.*/(.*)\.npz', r'\1_L.png', reffile), dpi=300)
    plt.close('all')

    #
    # # vec_lg = pos_gnss_ecef - pos_leo_ecef
    # # proj_lg = pos_leo_ecef - vec_lg * (np.sum(pos_leo_ecef * vec_lg, axis=1) / np.sum(vec_lg * vec_lg, axis=1))[:, None]
    # # alt_km = EarthLocation.from_geocentric(*proj_lg.T, unit='m').geodetic.height.value / 1E3
    # # if np.logical_or(np.all(alt_km < 90), np.all(alt_km > 0)):
    # #     raise ValueError(f'Non-RO direct height range: {np.min(alt_km):.3f} >> {np.max(alt_km):.3f} (km)')
    # # xx = sec_ro[::1 if obj_leo_ro['ocd'] else -1]
    # # yy = exLC[::1 if obj_leo_ro['ocd'] else -1]
    # # idx = np.isfinite(yy)
    # # xx = xx[idx]
    # # yy = yy[idx]
    # # len_old = np.where(np.abs(yy)>0.1)[0][0]
    # # coef = np.polyfit(xx[:len_old], yy[:len_old], 1)
    # # len_new = np.where(np.abs(np.polyval(coef, xx) - yy) > 0.1)[0][0]
    # # while len_new > len_old:
    # #     len_old = len_new
    # #     coef = np.polyfit(xx[:len_old], yy[:len_old], 1)
    # #     len_new = np.where(np.abs(np.polyval(coef, xx) - yy) > 0.1)[0][0]
    # # trend_LC = np.polyval(coef, sec_ro)
    # # # idx = np.where(np.logical_and(alt_km > max(90, max(alt_km) - 20), np.isfinite(exLC)))[0]
    # # # idx = np.where(np.logical_and(np.abs(alt_km - 100) < 10, np.isfinite(exLC)))[0]
    # # # trend_LC = np.polyval(np.polyfit(sec_ro[idx], exLC[idx], 1), sec_ro)
    # # exL1 -= trend_LC
    # # exL2 -= trend_LC
    # # exLC -= trend_LC
    # # idx = np.where(np.isfinite(exL1))[0]
    # # exL1 -= exL1[idx[0 if obj_leo_ro['ocd'] else -1]]
    # # exL2 -= exL2[idx[0 if obj_leo_ro['ocd'] else -1]]
    # # exLC -= exLC[idx[0 if obj_leo_ro['ocd'] else -1]]
    # # out_dir = f's04_atmPhs/{dt_rr:%Y-%m-%d}/'
    # # makedirs(out_dir, exist_ok=True)
    # # filestamp = f'{leoName}.{dt_ro[0 if obj_leo_ro["ocd"] else -1]:%Y.%j.%H.%M}.{gnss_ro}'
    # # out_file = out_dir + f'atmPhs_STAR_{filestamp}_v{version}.nc'