from datetime import datetime, timedelta import numpy as np from pathlib import Path working_dir = Path(__file__).absolute().parent def to_datetime(dtlist): return datetime(int(dtlist[0]), int(dtlist[1]), int(dtlist[2]),int(dtlist[3]),int(dtlist[4]),int(float(dtlist[5])),round((float(dtlist[5])%1)*1E6)) def dt2float(dt,dt0): return ((dt - dt0) / timedelta(seconds=1)).astype('float') def get_sat_freq(sat): GPRN = int(sat[1:]) if sat[0] == 'G': #GPS FL1 = 1.57542E9 FL2 = 1.22760E9 elif sat[0] == 'R': #GLONASS Kfrq = [1, -4, 5, 6, 1, -4, 5, 6, -2, -7, 0, -1, -2, -7, 0, -1, 4, -3, 3, 2, 4, -3, 3, 2] baseFL1 = 1.602E9 baseFL2 = 1.246E9 DFRQ1_GLO = 5.625E5 DFRQ2_GLO = 4.375E5 FL1 = baseFL1 + Kfrq[GPRN - 1] * DFRQ1_GLO FL2 = baseFL2 + Kfrq[GPRN - 1] * DFRQ2_GLO GPRN = GPRN + 100 elif sat[0] == 'E': #GALILEO FL1 = 1.57542E9 FL2 = 1.20714E9 GPRN = GPRN + 200 elif sat[0] == 'C': # BEIDOU FL1 = 1.57542E9 FL2 = 1.17645E9 GPRN = GPRN + 400 else: # ERROR!! FL1 = 0 FL2 = 0 GPRN = 0 return FL1, FL2, GPRN def load_gnss(file): with np.load(file, allow_pickle=True) as fid: Gdata = fid['Gdata'] Gdt = fid['Gdt'] Gsats = fid['Gsats'] return Gdata, Gdt, Gsats def read_gnss(dt,prev=False): import os Gdata1 = np.array([]) if prev: for x in ['NRT','RT']: gnssfile = f'../gnss/gnssOrb/GNSS_{dt - timedelta(days=1):%Y.%j}_{x}.npz' if os.path.isfile(gnssfile): if len(Gdata1) == 0: Gdata1, Gdt1, Gsats1 = load_gnss(gnssfile) else: Gdata0, Gdt0, Gsats0 = load_gnss(gnssfile) if Gdt1.shape == Gdt0.shape and np.all(Gdt1 == Gdt0): idx = ~np.isin(Gsats0, Gsats1) Gsats1 = np.hstack([Gsats1, Gsats0[idx]]) Gdata1 = np.hstack([Gdata1, Gdata0[:, idx, :]]) del Gdata0, Gdt0, Gsats0 Gdata2 = np.array([]) for x in ['NRT','RT']: gnssfile = f'../gnss/gnssOrb/GNSS_{dt:%Y.%j}_{x}.npz' if os.path.isfile(gnssfile): if len(Gdata2) == 0: Gdata2, Gdt2, Gsats2 = load_gnss(gnssfile) else: Gdata0, Gdt0, Gsats0 = load_gnss(gnssfile) if Gdt2.shape == Gdt0.shape and np.all(Gdt2 == Gdt0): idx = ~np.isin(Gsats0, Gsats2) Gsats2 = np.hstack([Gsats2, Gsats0[idx]]) Gdata2 = np.hstack([Gdata2, Gdata0[:, idx, :]]) del Gdata0, Gdt0, Gsats0 Gsats = Gsats2 if prev: idx = Gdt1dt_rng[1]) if np.any(tidx): gfiles = gfiles[tidx] tstart = tstart[tidx] tidx = tstart == max(tstart) gfiles = gfiles[tidx] tstart = tstart[tidx] # tend = tend[tidx] Gdt = None Gsats = [] Gdata = [] for ss in src_types: sidx = np.where([ss in gf for gf in gfiles])[0] if sidx.size == 0: continue gfile = gfiles[sidx[0]] Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0 dt_rng[0]) gfiles = gfiles[tidx] tstart = tstart[tidx] Gdt = None Gsats = [] Gdata = [] for gfile in gfiles: Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0dt_rng[1]) gfiles = gfiles[tidx] tstart = tstart[tidx] # tend = tend[tidx] Gdt = None for gsat in gsats: gfile = None for ss in src_types: sidx = np.where([ss in gf for gf in gfiles])[0] if sidx.size==0: continue sfile = gfiles[sidx[np.argmin(np.abs(dt_rng[0] - tstart[sidx] - timedelta(hours=6)))]] with np.load(sfile, allow_pickle=True) as fid: if gsat in fid['Gsats']: gfile = sfile break if gfile is None: raise FileNotFoundError(f'Could not find a GNSS Orb file containing {gsat}.') Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0 proj_lg2[1] exph -= exph[0 if ocd else -1] else: exph = np.zeros_like(l1) range_leo = np.sum(pos_leo_intp ** 2, axis=1) ** 0.5 range_gnss = np.sum(pos_gnss_intp ** 2, axis=1) ** 0.5 delta_rho = GMoCC * np.log((range_gnss + range_leo + range_lg) / (range_gnss + range_leo - range_lg)) tao_old = tao_new.copy() tao_new = (range_lg + exph + delta_rho) / C pos_gnss_intp = lagrange(sec_gnss - pos_gnss[:, 3] / 1E6, pos_gnss[:, :3], sec_ro - tao_new, 9) if not np.any(np.abs(tao_new - tao_old) > 1E-12): break vel_gnss_intp = lagrange(sec_gnss - pos_gnss[:, 3] / 1E6, vel_gnss[:, :3], sec_ro - tao_new, 9) range_dT = np.sum((pos_leo_intp - pos_gnss_intp) ** 2, axis=1) ** 0.5 # from scipy.interpolate import make_interp_spline # speed_dT = make_interp_spline(sec_ro, range_dT, k=2).derivative()(sec_ro) obj_gnss['pos'] = pos_gnss_intp obj_gnss['vel'] = vel_gnss_intp coef_leo = 1 - (GMoCC / range_leo - np.sum(vel_leo_intp ** 2, axis=1) / C / C) / 2 - const propertimedelay_gnss = (coef_leo - 1) * obj_gnss['clk_s'] * C delt_t_rel_gnss = -2 * np.sum(pos_gnss_intp * vel_gnss_intp, axis=1) / C / C delt_t_rel_leo = -2 * np.sum(pos_leo_intp * vel_leo_intp, axis=1) / C / C obj_gnss['range_diff_T'] = range_dT # obj_gnss['range'] = range_lg obj_gnss['GL_term'] = -delta_rho + C * delt_t_rel_gnss + propertimedelay_gnss + obj_gnss['clk_s'] * C obj_leo['delt_t_rel_r'] = delt_t_rel_leo obj_gnss['delt_t_rel_s'] = delt_t_rel_gnss obj_leo['ocd'] = ocd return obj_gnss, obj_leo def remove_cycleslip(ph, ph_m, freq, sat_l, ocd): c0 = 299792458 if not ocd: ph = ph[::-1] ph_m = ph_m[::-1] ph_new = (ph - ph_m) * (2 * freq / c0) if sat_l == 'G1': idx = np.unique(np.where(np.round(np.diff(ph_new[:100])) % 2 == 1)[0] % 2) if idx.size != 1: ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new))))) else: for ii in range(idx[0] + 1, len(ph_new), 2): ph_new[ii:ii + 2] -= np.round(ph_new[ii] - ph_new[ii - 1]) ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new) / 2) * 2))) elif sat_l[0] == 'R': ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new))))) else: ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new) / 2) * 2))) ph_new = ph_new * (c0 / (2 * freq)) + ph_m if not ocd: ph_new = ph_new[::-1] return ph_new # def load_dcb(dt): # import os # from re import sub # import pickle # dcbdir = '../gnss/gnssDCB/' # dcbfile = dcbdir+f'{dt:%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 < timedelta(days=32) and datetime.now().month == dt.month: # zfiles = ["P1C1_RINEX.DCB", # "P2C2_RINEX.DCB", # "P1P2_ALL.DCB"] # rmtdir = "http://ftp.aiub.unibe.ch/CODE/" # else: # zfiles = [f"P1C1{dt:%y%m}_RINEX.DCB.Z", # f"P2C2{dt:%y%m}_RINEX.DCB.Z", # f"P1P2{dt:%y%m}_ALL.DCB.Z"] # rmtdir = f"http://ftp.aiub.unibe.ch/CODE/{dt:%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) # return dcb def calculate_tec(reffile, tecfile=None, dumpnc=True): import os from re import sub from netCDF4 import Dataset from astropy.coordinates import EarthLocation from zoneinfo import ZoneInfo # if dumpnc: # out_file = sub(r'(.*)/p01_podPair/(.*)/occTab(.*)\.npz',r'\1/p05_podTec/\2/podTec\3.nc',reffile) # if os.path.exists(out_file) and datetime.fromtimestamp(os.path.getmtime(out_file))>datetime(2026,4,1): # print(f'File {out_file} already exists, skipping.') # return 0,1,1,1 if 'tecQ' in np.load(reffile): with np.load(reffile, allow_pickle=True) as fid: sec_ref = fid['sec_ref'] tec = fid['tec'] coszn = fid['coszn'] meps = fid['meps'] else: 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 # Constants c0 = 299792458 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.name).split('.')] dt_rr = datetime.strptime(sub(r'.*\.(\d{4}\.\d{3})\..*',r'\1',reffile.name),'%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.name) ant_leo_ref = int(sub(r'.*\.A(\d\d)\..*',r'\1',reffile.name)) gnss_ref = sub(r'.*\.(...).npz', r'\1', reffile.name) 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,:] try: Gdata, Gdt = read_gnss_v2(dt_ref, [gnss_ref]) except: raise LookupError(f'GNSS not found: {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"p01_leoOrb/leoOrb_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid: Ldata = fid['data'] Ldts = fid['dts'] tmsk = ~np.isin(np.arange(len(Ldts)), np.unique(np.where(np.diff(Ldts) > timedelta(seconds=5))[0] + np.arange(-30, 90)[:, None])) if np.all(Ldts[tmsk] < dt_ref[1]) or np.all(Ldts[tmsk] > dt_ref[0]): raise ValueError( f'LEO time ({Ldts[tmsk][0]:%d-%H:%M} >> {Ldts[tmsk][-1]:%d-%H:%M}) could not cover REF event ({dt_ref[0]:%H:%M} >> {dt_ref[1]:%H:%M})!') tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1)) if np.sum(tidx)<20: raise ValueError( f'LEO with large gap, no enough point in REF event ({dt_ref[0]:%H:%M} >> {dt_ref[1]:%H:%M})!') if np.all(tmsk[tidx]): dts_leoorb = Ldts[tidx] pos_leo = Ldata[tidx,:4] vel_leo = Ldata[tidx,4:] else: tidx0 = np.where(np.logical_and(Ldts < dt_ref[0], tmsk))[0] tidx1 = np.where(np.logical_and(Ldts > dt_ref[1], tmsk))[0] tidx = [tidx0[-60 if len(tidx0)>60 else 0],tidx1[min(60, len(tidx1)-1)]] dts_leoorb = Ldts[tidx[0]:tidx[1]] pos_leo = Ldata[tidx[0]:tidx[1], :4] vel_leo = Ldata[tidx[0]:tidx[1], 4:] tmsk = tmsk[tidx[0]:tidx[1]] sec_leoorb = dt2float(dts_leoorb, dt_rr) sec_mid = (sec_leoorb[0] + sec_leoorb[-1]) / 2 sec_rng = (sec_leoorb[-1] - sec_leoorb[0]) / 2 dts_leoorb = np.arange(dts_leoorb[0],dts_leoorb[-1],timedelta(seconds=1)).astype(datetime) pos_leo = np.array([np.polyval(np.polyfit((sec_leoorb[tmsk] - sec_mid) / sec_rng, x, 10), (np.arange(sec_leoorb[0], sec_leoorb[-1]) - sec_mid) / sec_rng) for x in pos_leo[tmsk, :3].T]+ [np.interp(np.arange(sec_leoorb[0], sec_leoorb[-1]), sec_leoorb, pos_leo[:, 3])]).T vel_leo = np.array([np.polyval(np.polyfit((sec_leoorb[tmsk] - sec_mid) / sec_rng, x, 10), (np.arange(sec_leoorb[0], sec_leoorb[-1]) - sec_mid) / sec_rng) for x in vel_leo[tmsk, :3].T]+ [np.interp(np.arange(sec_leoorb[0], sec_leoorb[-1]), sec_leoorb, vel_leo[:, 3])]).T del Ldata, Ldts with np.load(f"p01_leoAtt/leoAtt_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid: Ldata = fid['data'] Ldts = fid['dts'] Lvar = fid['vars'] Ldts, tidx = np.unique(Ldts, return_index=True) Ldata = Ldata[tidx] tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1)) if ~np.any(tidx): raise ValueError('No matching leoOrb/leoAtt') dts_leoatt = Ldts[tidx] att_leo = Ldata[tidx, :][:, ['att' in x for x in Lvar]] sad_leo = Ldata[tidx, :][:, ['sad_a' 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 sad_leo = CubicSpline(dt2float(dts_leoatt,dts_leo[0]), sad_leo)(dt2float(dts_leo,dts_leo[0])) 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 = np.hstack([fid['header'],['C1C','C2X']]) 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_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] antbsight_leo = fid[f'ANT_{leoName}_{antennas[ant_leo_ref]}'][15:18] 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 exL1 -= exL1[0 if obj_leo_ref['ocd'] else -1] exL2 -= exL2[0 if obj_leo_ref['ocd'] else -1] OL = obj_leo_ref['pos'] OG = obj_gnss_ref['pos'] if np.all(np.isnan(OG)) or np.all(np.isnan(OL)): raise ValueError('No finite coordinate values') LG = OG - OL antbsight_eci = inv_quat_rotate(att_leo, antbsight_leo) ex = lagrange(sec_leo / 86400, antbsight_eci, np.round(sec_ref) / 86400, 11) ex /= np.sum(ex ** 2, axis=1)[:, None] ** 0.5 ez = -OL # * [1, 1, 1.00669437999014] ey = np.cross(ez, ex) ey /= np.sum(ey ** 2, axis=1)[:, None] ** 0.5 ez = np.cross(ex, ey) LG_bsight = np.array([np.sum(LG * x, axis=1) for x in [ex, ey, ez]]).T el_bsight = -np.arcsin(LG_bsight[:, 2] / np.sum(LG_bsight ** 2, axis=1) ** 0.5) / np.pi * 180 az_bsight = np.arctan2(LG_bsight[:, 1], LG_bsight[:, 0]) / np.pi * 180 sad_ref = lagrange(sec_leo / 86400, sad_leo, np.round(sec_ref) / 86400, 11)[:, 0] bslut = f'./p03_bslut/bslut_{leoName}.2025-09-18.D125.A{ant_leo_ref:02d}.{gnss_ref[0]}--.npz' with np.load(bslut, allow_pickle=True) as fid: mpe_grid = fid['mpe_grid'] sad_bins = fid['sad_bins'] mpe_grid = mpe_grid[np.floor(el_bsight + 90).astype('int'), np.floor(az_bsight + 180).astype('int')] idx = np.any(np.isfinite(mpe_grid), axis=(0, 2)) mpe = np.array( [[np.interp(sad_ref[x], sad_bins[idx], mpe_grid[x, idx, y]) for y in range(2)] for x in range(len(el_bsight))]) proj_lg = OL - LG * (np.sum(OL * LG, axis=1) / np.sum(LG * LG, axis=1))[:, None] 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 coszn = np.sum(OL * LG, axis=1) / np.sum(LG * LG, axis=1)**0.5 / np.sum(OL * OL, axis=1)**0.5 # if np.all(coszn<0): # raise ValueError('All observation below horizon!') alti_km = 4000 # km rad_ratio = (alti_km+6378.137)/(altl_km+6378.137) meps = (coszn + (rad_ratio**2+coszn**2-1)**0.5)/(1+rad_ratio) # 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) tt = np.arange(mpe.shape[0]) idx = np.all(np.isfinite(mpe), axis=1) mpe = np.array([np.interp(tt,tt[idx],mpe[idx,x]) for x in [0,1]]).T 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[coszn<0] = 0 D1 = data_ref[:, c_ref[0]] - fr_ref * data_ref[:, l_ref[0]] + (fr_ref - 1) * data_ref[:, l_ref[1]] - mpe[:,0] D2 = data_ref[:, c_ref[1]] - (fr_ref + 1) * data_ref[:, l_ref[0]] + fr_ref * data_ref[:, l_ref[1]] - mpe[:,1] S /= np.sum(S, axis=0)[None, :] Ddiff = np.nansum(D1 * S[:, 0]) - np.nansum(D2 * S[:, 1]) Cdiff = (data_ref[:, c_ref[1]] - data_ref[:, c_ref[0]]) / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 # tec = (Pdiff - Ddiff) / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 tec = (Pdiff - Ddiff) / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 # tec0 = tec.copy() tec -= np.nanmean(tec-Cdiff) # dcb = load_dcb(dt_ref[0]) # dcb_factor = c0 / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 # dcb_value = (dcb['P1P2'][gnss_ref][0] - dcb['P1C1'][gnss_ref][0] + dcb['P2C2'][gnss_ref][0]) / 1E9 * dcb_factor # tec += dcb_value dts_ref = [dt_rr + timedelta(seconds=x) for x in sec_ref] pos_leo_ecef, vel_leo_ecef, pos_leo_tod, vel_leo_tod = ro_eci2ecef(obj_leo_ref['pos'], obj_leo_ref['vel'], dts_ref, tod=True) pos_gnss_ecef, vel_gnss_ecef, pos_gnss_tod, vel_gnss_tod = ro_eci2ecef(obj_gnss_ref['pos'], obj_gnss_ref['vel'], dts_ref, tod=True) with np.load(reffile, allow_pickle=True) as fid: data = fid['data'] header = fid['header'] np.savez(reffile, data=data, header=header, sec_ref=sec_ref, tec=tec, coszn=coszn, meps=meps) if dumpnc: out_dir = working_dir / f'p05_podTec/{dt_rr:%Y-%m-%d}' filestamp = f'{leoName}_{dt_ref[0]:%Y%m%d%H%M%S}_{dt_ref[-1]:%Y%m%d%H%M%S}_{gnss_ref}_A{ant_leo_ref:02d}_c{datetime.now().astimezone(ZoneInfo("UTC")):%Y%m%d%H%M%S}' out_file = out_dir / f'ROSW_TEC_{filestamp}_nes_ops.nc' # with np.load('dcb_default.npz', allow_pickle=True) as fid: # allgnsses = fid['allgnsses'][0 if leoName == 'GN04' else 1, ant_leo_ref] # alldcbs = fid['alldcbs'][0 if leoName == 'GN04' else 1, ant_leo_ref] # with np.load('dcb_daily.npz', allow_pickle=True) as fid: # allgnsses = fid['allgnsses'] # alldcbs = fid['alldcbs'][0 if leoName == 'GN04' else 1, ant_leo_ref] # alldts = fid['alldts'] # alldcbs = alldcbs[dt_rr==alldts][0] dcb = 0.0 for dayprev in range(366): dt_now = dt_rr - timedelta(days=dayprev) dcbfile = Path(f'p04_teczdcb/teczdcb_{leoName}.{dt_now:%Y-%m-%d}.A{ant_leo_ref:02d}.npz') if dcbfile.exists(): with np.load(dcbfile, allow_pickle=True) as fid: allgnsses = fid['gnsses'] alldcbs = fid['dcb_sol'] if gnss_ref in allgnsses: dcb = alldcbs[allgnsses==gnss_ref][0] break dt_now = datetime.now() with Dataset(out_file, 'w') as fid: fid.createDimension('time',len(tec)) fid.createDimension('attribute_dimension',1) fid.createVariable('time', 'f8', ('time', )) fid['time'].setncatts( {'_FillValue': -999.0, 'calendar': 'standard', 'long_name': f"Seconds since start time {sec_rr + sec_ref[0]:.1f}", 'units': f"Seconds since start time {sec_rr + sec_ref[0]:.1f}", 'add_offset': sec_rr + sec_ref[0]} ) fid['time'][...] = sec_ref+sec_rr fid.createVariable('TEC', 'f8', ('time',)) fid['TEC'].setncatts( {'_FillValue': -999.0, 'units': 'TECU', 'long_name': 'Total Electron Content along LEO-GPS link'} ) fid['TEC'][...] = tec+dcb fid.createVariable('antid', 'int16', ('time',)) fid['antid'].setncatts( {'_FillValue': np.int16(-999), 'units': 'Unitless', 'long_name': 'Antenna ID for each time step'} ) fid['antid'][...] = np.full((len(tec),),ant_leo_ref).astype('int16') fid.createVariable('elevation', 'f8', ('time',)) fid['elevation'].setncatts( {'_FillValue': -999.0, 'units': 'deg', 'long_name': 'Elevation angle of LEO-GPS link'} ) fid['elevation'][...] = np.arcsin(coszn)/np.pi*180 fid.createVariable('occheight', 'f8', ('time',)) fid['occheight'].setncatts( {'_FillValue': -999.0, 'units': 'km', 'long_name': 'Occultation perigee point height'} ) fid['occheight'][...] = altp_km fid.createVariable('caL1_SNR', 'f4', ('time',)) fid['caL1_SNR'].setncatts( {'_FillValue': np.float32(-999.0), 'long_name': 'First frequency signal strength'} ) fid['caL1_SNR'][...] = data_ref[:, s_ref[0]].astype('float32') fid.createVariable('pL2_SNR', 'f4', ('time',)) fid['pL2_SNR'].setncatts( {'_FillValue': np.float32(-999.0), 'long_name': 'Second frequency signal strength'} ) fid['pL2_SNR'][...] = data_ref[:, s_ref[1]].astype('float32') fid.createVariable('x_LEO', 'f8', ('time',)) fid['x_LEO'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF LEO satellite X position', 'valid_range': [-8000.0, 8000.0]} ) fid['x_LEO'][...] = pos_leo_ecef[:, 0] / 1E3 fid.createVariable('y_LEO', 'f8', ('time',)) fid['y_LEO'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF LEO satellite Y position', 'valid_range': [-8000.0, 8000.0]} ) fid['y_LEO'][...] = pos_leo_ecef[:, 1] / 1E3 fid.createVariable('z_LEO', 'f8', ('time',)) fid['z_LEO'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF LEO satellite Z position', 'valid_range': [-8000.0, 8000.0]} ) fid['z_LEO'][...] = pos_leo_ecef[:, 2] / 1E3 fid.createVariable('x_GNSS', 'f8', ('time',)) fid['x_GNSS'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF GNSS satellite X position', 'valid_range': [-30000.0, 30000.0]} ) fid['x_GNSS'][...] = pos_gnss_ecef[:, 0] / 1E3 fid.createVariable('y_GNSS', 'f8', ('time',)) fid['y_GNSS'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF GNSS satellite Y position', 'valid_range': [-30000.0, 30000.0]} ) fid['y_GNSS'][...] = pos_gnss_ecef[:, 1] / 1E3 fid.createVariable('z_GNSS', 'f8', ('time',)) fid['z_GNSS'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF GNSS satellite Z position', 'valid_range': [-30000.0, 30000.0]} ) fid['z_GNSS'][...] = pos_gnss_ecef[:, 2] / 1E3 fid.createVariable('exph_L1', 'f8', ('time',)) fid['exph_L1'].setncatts( {'units': 'm', 'long_name': 'Excess phase L1'} ) fid['exph_L1'][...] = exL1 fid.createVariable('exph_L2', 'f8', ('time',)) fid['exph_L2'].setncatts( {'units': 'm', 'long_name': 'Excess phase L2'} ) fid['exph_L2'][...] = exL2 fid.setncatts({ 'start_time': np.uint32(np.round(sec_ref[0]+sec_rr)), 'stop_time': np.uint32(np.round(sec_ref[-1]+sec_rr)), 'gobs_start': np.uint32(np.round(sec_ref[0]+sec_rr)), 'gobs_stop': np.uint32(np.round(sec_ref[-1]+sec_rr)), 'History': f'Created by atec.py {dt_now:%d/%m/%Y %H:%M:%S}', 'Title': 'Absolute TEC measurements from radio occultation', 'Conventions': 'CF-1.6', 'year': dts_ref[0].year, 'month': dts_ref[0].month, 'day': dts_ref[0].day, 'hour': dts_ref[0].hour, 'minute': dts_ref[0].minute, 'second': dts_ref[0].second+dts_ref[0].microsecond/1E6, 'tecmin': np.min(tec+dcb), 'tecmax': np.max(tec+dcb), 'elevmin': np.min(np.arcsin(coszn)/np.pi*180), 'elevmax': np.max(np.arcsin(coszn)/np.pi*180), 'occaltmin': np.min(altp_km), 'occaltmax': np.max(altp_km), 'processing_center': 'NESDIS', 'L2C': 1, 'creation_time': f'{dt_now:%Y-%m-%d %H:%M:%S.%f}', 'mission': 'gnomestec', 'dump_id': f'{dt_ref[0]:%Y.%j}.0{leoName[2:]}.{ant_leo_ref:02d}', 'leo_id': int(leoName[3]), 'prn_id': int(gnss_ref[1:]), 'duration': int(sec_ref[-1]-sec_ref[0]), 'leodcb_flag': 1, 'leodcb': dcb, 'gpsdcb_flag': 0, 'gpsdcb': -999, 'leveling_offset': -Ddiff / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16, 'leveling_err': np.sum((tec-Cdiff)**2)**0.5/len(tec), 'conid': gnss_ref[0], 'dcb_units': 'TECU', }) return sec_ref, tec, coszn, meps def combine_dcb(): from glob import glob from re import sub gnsses = np.array( ['C06', 'C07', 'C08', 'C09', 'C10', 'C11', 'C12', 'C13', 'C14', 'C16', 'C19', 'C20', 'C21', 'C22', 'C23', 'C24', 'C25', 'C26', 'C27', 'C28', 'C29', 'C30', 'C32', 'C33', 'C34', 'C35', 'C36', 'C37', 'C38', 'C39', 'C40', 'C41', 'C42', 'C43', 'C44', 'C45', 'C47', 'C48', 'C49', 'C50', 'E02', 'E03', 'E04', 'E05', 'E06', 'E07', 'E08', 'E09', 'E10', 'E11', 'E12', 'E13', 'E14', 'E15', 'E16', 'E18', 'E19', 'E21', 'E23', 'E24', 'E25', 'E26', 'E27', 'E29', 'E30', 'E31', 'E33', 'E34', 'E36', 'G01', 'G02', 'G03', 'G04', 'G05', 'G06', 'G07', 'G08', 'G09', 'G10', 'G11', 'G12', 'G13', 'G14', 'G15', 'G16', 'G17', 'G18', 'G19', 'G20', 'G21', 'G22', 'G23', 'G24', 'G25', 'G26', 'G27', 'G28', 'G29', 'G30', 'G31', 'G32', 'J02', 'J03', 'J04', 'R02', 'R03', 'R04', 'R05', 'R06', 'R07', 'R08', 'R09', 'R11', 'R12', 'R14', 'R15', 'R16', 'R17', 'R18', 'R19', 'R20', 'R21', 'R22', 'R24', 'R25', 'R26', 'R27', 'R28']) ants = [0,1] leos = ['GN04','GN05'] # gnsstypes = ['C','E'] alldcbs = [[[],[]],[[],[]]] allgnsses = [[[],[]],[[],[]]] alldts = [[[],[]],[[],[]]] kernal = np.array([1]) # np.arange(10, 0, -1) / 55 for aa,ant in enumerate(ants): for ll,leo in enumerate(leos): dcbfiles = glob(f'p04_teczdcb/teczdcb_{leo}.*.A{ant:02d}.npz') dcbfiles.sort() dts = np.array([datetime.strptime(sub(r'.*\.(\d{4}-\d\d-\d\d)\..*',r'\1',x),'%Y-%m-%d') for x in dcbfiles]) idx = dts>datetime(2025,12,15) dcbfiles = np.array(dcbfiles)[idx] dts = dts[idx] alldcb = np.full((len(dts),len(gnsses)),np.nan) for dd, dt, dcbfile in zip(range(len(dts)),dts,dcbfiles): with np.load(dcbfile, allow_pickle=True) as fid: # sec_tec = fid['sec_tec'] # tecz = fid['tecz'] dcb_gnss = fid['gnsses'] dcb_sol = fid['dcb_sol'] alldcb[dd,np.isin(gnsses,dcb_gnss)] = dcb_sol idx = np.any(np.isfinite(alldcb), axis=0) allgnss = gnsses[idx] alldcb = alldcb[:, idx] idx = np.any(np.abs(alldcb-np.nanmedian(alldcb,axis=0))>3.5,axis=1) alldcb[idx] = np.nan for ii in range(alldcb.shape[1]): idx = np.isfinite(alldcb[:,ii]) alldcb[:,ii] = np.interp(dt2float(dts,dts[0]),dt2float(dts[idx],dts[0]),alldcb[idx,ii]) alldcbs[ll][aa] = np.array([np.convolve(kernal,alldcb[:,x],'valid') for x in range(alldcb.shape[1])]).T # alldcbs[ll][aa] = np.nanmedian(alldcb, axis=0) allgnsses[ll][aa] = allgnss.copy() alldts[ll][aa] = dts[len(kernal)-1:] cmb_gnss = [] cmb_dts = [] for aa,ant in enumerate(ants): for ll,leo in enumerate(leos): cmb_gnss = np.union1d(cmb_gnss,allgnsses[aa][ll]) cmb_dts = np.union1d(cmb_dts,alldts[aa][ll]) cmb_dcbs = np.full((len(ants),len(leos),len(cmb_dts),len(cmb_gnss)),np.nan) for aa,ant in enumerate(ants): for ll,leo in enumerate(leos): tidx = np.isin(cmb_dts,alldts[aa][ll]) gidx = np.isin(cmb_gnss,allgnsses[aa][ll]) idx = np.ix_([aa],[ll],tidx,gidx) cmb_dcbs[idx] = alldcbs[aa][ll] idx = np.isfinite(cmb_dcbs) if ~np.all(idx): for tt,dts in enumerate(cmb_dts): for gg,gnss in enumerate(cmb_gnss): if np.any(~np.isfinite(cmb_dcbs[:,:,tt,gg])): for aa,ant in enumerate(ants): idx = np.where(np.isfinite(cmb_dcbs[aa, :, tt, gg]))[0] if len(idx)>0: cmb_dcbs[aa, :, tt, gg] = np.interp(np.arange(len(leos)), idx, cmb_dcbs[aa, idx, tt, gg]) for ll,leo in enumerate(leos): idx = np.where(np.isfinite(cmb_dcbs[:, ll, tt, gg]))[0] if len(idx)>0: cmb_dcbs[:, ll, tt, gg] = np.interp(np.arange(len(ants)), idx, cmb_dcbs[idx, ll, tt, gg]) np.savez('dcb_daily.npz',alldts=cmb_dts,alldcbs=cmb_dcbs,allgnsses=cmb_gnss) return def calc_bore_sight(args): from re import sub import os from glob import glob from scipy.interpolate import CubicSpline from frame_conversions import ro_ecef2eci from traceback import extract_tb # c0 = 299792458 # kappa = 40.30819294518825 leoName = args[0] ant_ref = args[1] gnss_ref = args[2] antennas = ['POD0', 'POD1', 'OCC2', 'OCC3', 'OCC4', 'OCC5', 'OCC6', 'OCC7'] dt0_gps = datetime(1980, 1, 6) ds_ref = [datetime.strptime(sub(r'.*/',r'',x),'%Y-%m-%d') for x in glob('./p01_podPair/*')] ds_ref.sort() # plt.close('all') # fig1, axs1 = plt.subplots(3, 3, figsize=(12, 9)) # fig2, axs2 = plt.subplots(3, 3, figsize=(12, 9)) for d_ref in ds_ref: out_dir = f'./p02_bsight/{d_ref:%Y-%m-%d}/' os.makedirs(out_dir, exist_ok=True) # for leoName in leos_ref: # for gnss_ref in gnsses_ref: # for ant_ref in ants_ref: out_file = out_dir + f'bsight_{leoName}_{d_ref:%Y.%j}.A{ant_ref:02d}.{gnss_ref}.npz' if os.path.exists(out_file): continue reffiles = glob(f'./p01_podPair/{d_ref:%Y-%m-%d}/occTab_{leoName}*A{ant_ref:02d}.{gnss_ref}.npz') if len(reffiles) == 0: continue reffiles.sort() mpe_0 = [] sec_ref_0 = [] sad_ref_0 = [] el_bsight_0 = [] az_bsight_0 = [] # print(f'{datetime.now():%H:%M:%S} - {d_ref:%Y-%m-%d} ({leoName}-{gnss_ref}-A{ant_ref:02d}), {len(reffiles):2d} files:', end='', flush=True) for reffile in reffiles: try: # print('.', end='', flush=True) 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) gnss_ref = sub(r'.*\.(...).npz', r'\1', reffile) ant_ref = int(sub(r'.*\.A(\d\d)\..*', r'\1', reffile)) with np.load('antenna_offset.npz') as fid: # antoffset_leo = fid[f'ANT_{leoName}_{antennas[ant_ref]}'][12:15] antbsight_leo = fid[f'ANT_{leoName}_{antennas[ant_ref]}'][15:18] # tecfiles = np.array(glob(f'/data3/xshao/planetiq/SWx/UCAR/L1B/20251220/podTec_GN05.*.{gnss_ref}.00*_nc')) # dts_ucar = np.array([datetime.strptime(sub(r'.*_GN0.\.(.*)\.\d{4}\.G.*',r'\1',tecfile),'%Y.%j.%H.%M') for tecfile in tecfiles]) # dur_ucar = np.array([timedelta(minutes=int(sub(r'.*\.(\d{4})\.G.*',r'\1',tecfile))) for tecfile in tecfiles]) # fidx = np.logical_and(dts_ucardt_ref[0]) # tecfiles = tecfiles[fidx] 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]).astype('int') f1_ref, f2_ref, _ = get_sat_freq(gnss_ref) try: Gdata, Gdt = read_gnss_v2(dt_ref, [gnss_ref]) except: continue 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).astype('int') pos_gnss = Gdata[tidx, 0, :4] vel_gnss = Gdata[tidx, 0, 4:] # OG = lagrange(sec_gnss/86400,Gdata[tidx, 0, :3],sec_ref/86400,11) # vel_gnss = Gdata[tidx, :, 4:] del Gdata, Gdt # Read in leo position/velocity/attitude files accordingly with np.load(f"p01_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:] # sec_leo = dt2float(dts_leo, dt_rr).astype('int') # OL = lagrange(sec_leo/86400,Ldata[tidx, :3], sec_ref/86400, 11) # vel_leo = Ldata[tidx, 4:] del Ldata, Ldts with np.load(f"p01_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)) if ~np.any(tidx): continue dts_leoatt = Ldts[tidx] sad_leo = Ldata[tidx, :][:, ['sad_a' in x for x in Lvar]] att_leo = Ldata[tidx, :][:, ['att' in x for x in Lvar]] aridx = np.cumsum(np.hstack(([False], np.sum(att_leo[1:] * att_leo[:-1], axis=1) < 0))) % 2 == 1 att_leo[aridx] *= -1 del Ldata, Ldts, Lvar tidx = np.logical_and(dts_leoorb >= dts_leoatt[0], dts_leoorb <= dts_leoatt[-1]) if ~np.any(tidx): continue dts_leo = dts_leoorb[tidx] sec_leo = dt2float(dts_leo, dt_rr) pos_leo = pos_leo[tidx] vel_leo = vel_leo[tidx] dts_leoatt, tidx = np.unique(dts_leoatt, return_index=True) sad_leo = CubicSpline(dt2float(dts_leoatt, dts_leo[0]), sad_leo[tidx])(dt2float(dts_leo, dts_leo[0])) att_leo = CubicSpline(dt2float(dts_leoatt, dts_leo[0]), att_leo[tidx])(dt2float(dts_leo, dts_leo[0])) att_leo /= np.sum(att_leo ** 2, axis=1)[:, None] ** 0.5 del dts_leoorb, dts_leoatt 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})!') continue pos_leo_eci, vel_leo_eci = ro_ecef2eci(pos_leo, vel_leo, dts_leo) pos_ref_eci, vel_ref_eci = ro_ecef2eci(pos_gnss[:, :3], vel_gnss, dts_gnss) OL = lagrange(sec_leo / 86400, pos_leo_eci, sec_ref / 86400, 11) OG = lagrange(sec_gnss / 86400, pos_ref_eci, sec_ref / 86400, 11) antbsight_eci = inv_quat_rotate(att_leo, antbsight_leo) ex = lagrange(sec_leo / 86400, antbsight_eci, sec_ref / 86400, 11) ex /= np.sum(ex ** 2, axis=1)[:, None] ** 0.5 ez = -OL # * [1, 1, 1.00669437999014] ey = np.cross(ez, ex) ey /= np.sum(ey ** 2, axis=1)[:, None] ** 0.5 ez = np.cross(ex, ey) LG = OG - OL LG_bsight = np.array([np.sum(LG * x, axis=1) for x in [ex, ey, ez]]).T el_bsight = -np.arcsin(LG_bsight[:, 2] / np.sum(LG_bsight ** 2, axis=1) ** 0.5) / np.pi * 180 az_bsight = np.arctan2(LG_bsight[:, 1], LG_bsight[:, 0]) / np.pi * 180 sad_ref = lagrange(sec_leo / 86400, sad_leo, sec_ref / 86400, 11)[:, 0] # idx_sad = np.isin(np.round(sad_ref * 10), np.arange(-1200, 1500, 300) - 123) proj_lg = OL - LG * (np.sum(OL * LG, axis=1) / np.sum(LG * LG, axis=1))[:, None] # altp_km = EarthLocation.from_geocentric(*proj_lg.T, unit='km').geodetic.height.value # altl_km = EarthLocation.from_geocentric(*OL.T, unit='km').geodetic.height.value # coszn = np.sum(OL * LG, axis=1) / np.sum(LG * LG, axis=1) ** 0.5 / np.sum(OL * OL, axis=1) ** 0.5 # 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] S2 = S ** 2 # cumS = np.cumsum(np.vstack(([[0, 0]], S)), axis=0) # cumS2 = np.cumsum(np.vstack(([[0, 0]], S2)), axis=0) # clen = 31 # S4 = ((cumS2[clen:] - cumS2[:-clen]) * clen / (cumS[clen:] - cumS[:-clen]) ** 2 - 1) ** .5 # L = data_ref[:, l_ref] - data_ref[:, c_ref] # L2 = L ** 2 # cumL = np.cumsum(np.vstack(([[0, 0]], L)), axis=0) # cumL2 = np.cumsum(np.vstack(([[0, 0]], L2)), axis=0) # sigma = ((cumL2[clen:] - cumL2[:-clen]) * clen / (cumL[clen:] - cumL[:-clen]) ** 2 - 1) ** .5 S2 /= np.sum(S2, 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]] sig1 = np.sum(D1 * S2[:, 0]) sig2 = np.sum(D2 * S2[:, 1]) sad_ref = np.round(sad_ref * 10) idx_sad = np.isin(sad_ref, np.arange(-1200, 1500, 300) - 123) sad_ref /= 10 mpe = np.array([D1[idx_sad] - sig1, D2[idx_sad] - sig2]).T sec_ref = sec_ref[idx_sad] sad_ref = sad_ref[idx_sad] el_bsight = el_bsight[idx_sad] az_bsight = az_bsight[idx_sad] except Exception as e: with open("bore_sight_error.txt", 'a') as fid: fid.write(f"{datetime.now().strftime('%H:%M:%S')} -- Error detected for {reffile}:\n") fid.write(f' --- {type(e).__name__} >> {str(e)}\n') for tb in extract_tb(e.__traceback__): fid.write(f' --- line {tb.lineno} of {tb.filename}\n') mpe_0.append(mpe) sec_ref_0.append(sec_ref) sad_ref_0.append(sad_ref) el_bsight_0.append(el_bsight) az_bsight_0.append(az_bsight) if len(mpe_0) > 1: mpe_0 = np.vstack(mpe_0) sec_ref_0 = np.hstack(sec_ref_0) sad_ref_0 = np.hstack(sad_ref_0) el_bsight_0 = np.hstack(el_bsight_0) az_bsight_0 = np.hstack(az_bsight_0) elif len(mpe_0) == 1: mpe_0 = mpe_0[0] sec_ref_0 = sec_ref_0[0] sad_ref_0 = sad_ref_0[0] el_bsight_0 = el_bsight_0[0] az_bsight_0 = az_bsight_0[0] else: # print('', flush=True) continue sec_ref_0, idx_sec = np.unique(sec_ref_0, return_index=True) mpe_0 = mpe_0[idx_sec] sad_ref_0 = sad_ref_0[idx_sec] el_bsight_0 = el_bsight_0[idx_sec] az_bsight_0 = az_bsight_0[idx_sec] np.savez(out_file, mpe=mpe_0, sec_ref=sec_ref_0, sad_ref=sad_ref_0, el_bsight=el_bsight_0, az_bsight=az_bsight_0) # print('', flush=True) return 0 def gen_bore_sight_lut(): from glob import glob from re import sub degsize = 1 el_rng = [-90, 90] az_rng = [-180, 180] el_size = int((el_rng[1] - el_rng[0]) // degsize) az_size = int((az_rng[1] - az_rng[0]) // degsize) el_bins = np.linspace(el_rng[0], el_rng[1], el_size + 1) az_bins = np.linspace(az_rng[0], az_rng[1], az_size + 1) for leoName in ['GN04','GN05']: for ant_ref in range(2): for gnss_ref in ['C*','E*','G*']: bsfiles = glob(f'./p02_bsight/*/bsight_{leoName}*.A{ant_ref:02d}.{gnss_ref}.npz') bsfiles.sort() bsfiles = [x for x in bsfiles if datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\..*', r'\1', x), '%Y.%j') <= datetime(2026,1,20)] date_rng = [datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\..*', r'\1', bsfiles[x]), '%Y.%j') for x in [0,-1]] gnss_fname = sub(r"\*", r"--", gnss_ref) outfile = f'./p03_bslut/bslut_{leoName}.{date_rng[0]:%Y-%m-%d}.D{int((date_rng[1]-date_rng[0])/timedelta(days=1)+1):03d}.A{ant_ref:02d}.{gnss_fname}.npz' print(f'{datetime.now():%H:%M:%S} - {outfile}') mpe = [] sad_ref = [] el_bsight = [] az_bsight = [] for bsfile in bsfiles: with np.load(bsfile, allow_pickle=True) as fid: mpe.append(fid['mpe'][...]) sad_ref.append(fid['sad_ref'][...]) el_bsight.append(fid['el_bsight'][...]) az_bsight.append(fid['az_bsight'][...]) mpe = np.vstack(mpe) sad_ref = np.hstack(sad_ref) el_bsight = np.hstack(el_bsight) az_bsight = np.hstack(az_bsight) # idx = np.all(np.abs(mpe) < 3, axis=1) # mpe = mpe[idx] # sad_ref = sad_ref[idx] # el_bsight = el_bsight[idx] # az_bsight = az_bsight[idx] sad_bins = np.unique(sad_ref) mpe_grid = np.full((el_size, az_size, 9, 2), np.nan) el_bin = ((el_bsight - el_rng[0]) // degsize).astype('int') el_bin[el_bin==180] = 179 az_bin = ((az_bsight - az_rng[0]) // degsize).astype('int') az_bin[az_bin==360] = 359 ea_bin = el_bin * 360 + az_bin sad_bin = np.round((sad_ref - sad_bins[0]) / 30).astype('int') idx = np.lexsort((sad_bin,ea_bin)) mpe = mpe[idx] # sad_ref = sad_ref[idx] # el_bsight = el_bsight[idx] # az_bsight = az_bsight[idx] ea_bin = ea_bin[idx] sad_bin = sad_bin[idx] div = np.hstack(([0],np.where(np.diff(ea_bin)>0)[0]+1,[len(ea_bin)])) for dd in range(len(div)-1): if dd%1000==0: print(dd) el_idx = int(ea_bin[div[dd]] // 360) az_idx = int(ea_bin[div[dd]] % 360) sad_slice = sad_bin[div[dd]:div[dd+1]] mpe_slice = mpe[div[dd]:div[dd+1]] sdiv = np.hstack(([0],np.where(np.diff(sad_slice)>0)[0]+1,[len(sad_slice)])) for ss in range(len(sdiv)-1): if sdiv[ss + 1] - sdiv[ss] < 3: continue sad_idx = sad_slice[sdiv[ss]] for ll in range(2): mpe_bin = mpe_slice[sdiv[ss]:sdiv[ss + 1], ll] mpe_bin = mpe_bin[np.abs(mpe_bin)<3] mpe_med = np.median(mpe_bin) mpe_std = np.std(mpe_bin) while np.any(np.abs(mpe_bin - mpe_med) > 3 * mpe_std): mpe_bin = mpe_bin[np.abs(mpe_bin - mpe_med) < 3 * mpe_std] mpe_med = np.median(mpe_bin) mpe_std = np.std(mpe_bin) mpe_grid[el_idx, az_idx, sad_idx, ll] = np.mean(mpe_bin) for ll in range(2): sad_idx = np.where(np.isfinite(mpe_grid[el_idx, az_idx, :, ll]))[0] if len(sad_idx) > 0: mpe_grid[el_idx, az_idx, :, ll] = np.interp(np.arange(9), sad_idx, mpe_grid[el_idx, az_idx, sad_idx, ll], left=np.nan, right=np.nan) np.savez(outfile, el_bins=el_bins, az_bins=az_bins, sad_bins=sad_bins, mpe_grid=mpe_grid) return