from datetime import datetime, timedelta
import os
import numpy as np
from glob import glob
from re import sub
from scipy.interpolate import CubicSpline
from tools import read_gnss

def dt2float(dt,dt0):
    return ((dt - dt0) / timedelta(seconds=1)).astype('float')

def cal_zenith(dt):
    Gdata, Gdt, Gsats = read_gnss(dt, prev=True)
    idx = np.all(np.isfinite(Gdata[:-1, :, :3]), axis=(0, 2))
    Gdata = Gdata[:, idx]
    Gsats = Gsats[idx]
    for leofile in glob(f's01_leoOrb/*.{dt:%Y.%j}.npz'):
        leosat = sub(r'.*_(....)\..*',r'\1',leofile)
        with np.load(leofile, allow_pickle=True) as fid:
            Ldata = fid['data']
            Ldt = fid['dts']
        Tbdry = [max(Gdt[0],Ldt[0]), min(Gdt[-1],Ldt[-1])]
        Tbdry_i = [x.replace(second=10*(x.second//10), microsecond=0) for x in Tbdry]
        if Tbdry_i[0] != Tbdry[0]:
            Tbdry_i[0] = Tbdry_i[0] + timedelta(seconds=10)
        Tnew = np.arange(Tbdry_i[0], Tbdry_i[1] + timedelta(seconds=10), timedelta(seconds=10)).astype(datetime)
        Gnew = CubicSpline(dt2float(Gdt[:-1], Tnew[0]), Gdata[:-1, :, :3])(dt2float(Tnew, Tnew[0]))
        Lnew = CubicSpline(dt2float(Ldt, Tnew[0]), Ldata[:, :3])(dt2float(Tnew, Tnew[0]))[:, None, :]
        LG = Gnew - Lnew
        zenith = np.arccos(np.sum(LG * Lnew, axis=2) /
                           np.sum(LG * LG, axis=2) ** 0.5 / np.sum(Lnew * Lnew, axis=2) ** 0.5) / np.pi * 180
        np.savez(f's02_zenith/zenith_{leosat}.{dt:%Y.%j}_.npz', zenith=zenith,Tnew=Tnew,Gsats=Gsats)