# -*- coding: utf-8 -*- """ This Python program provides the data-reading routines used throughout the analysis. It imports the relevant N2 and N3d binary data, together with auxiliary ephemeris and beta-angle data, and converts them into structured NumPy arrays or dictionaries for further processing. As such, it forms the input stage of the workflow and supplies the data required by the subsequent programs. @author: Dorian Jost """ import numpy as np import datetime def read_n3d(filepath): """ Reads a N3D binary file and returns a structured NumPy array. Parameters: filepath (str): Path to the N3D binary file. Returns: data (np.ndarray): Structured array with fields: 'ydh', 'num', 's', 'q', 'u', 'v', 'th', 'ph', 'sn' """ # Define the structure of one record (40 bytes) record_dtype = np.dtype([ ('ydh', np.int32), # yyyydddhh of current file ('num', np.int32), # index of record in current file ('s', np.float32), # Intensity (V^2m^-2Hz^-1) ('q', np.float32), # Normalized linear polarization ('u', np.float32), # parameters Q & U ('v', np.float32), # Normalized circular polar. V ('th', np.float32), # Source colatitude (S/C frame) ('ph', np.float32), # Source azimuth (S/C frame) ('sn', np.float32, (2,)) # S/N ratio on autocorrel. values (2-element float array) ]) # Read the binary file into a structured array with open(filepath, 'rb') as f: data = np.fromfile(f, dtype=record_dtype) return data def read_n2(filepath): """ Reads a N2 binary file and returns a structured NumPy array. Parameters: filepath (str): Path to the N2 binary file. Returns: data (np.ndarray): Structured array with fields: 'ydh', 'num', 't97', 'f', 'dt', 'df', 'autoX', 'autoZ', 'crossR', 'crossI', 'ant' """ # Define the structure of one N2 record (45 bytes) record_dtype = np.dtype([ ('ydh', np.int32), # Timestamp ('num', np.int32), # Record index ('t97', np.float64), # Decimal time since 1997.0 ('f', np.float32), # Frequency (KHz) ('dt', np.float32), # Integration time (ms) ('df', np.float32), # Bandwidth (KHz) ('autoX', np.float32), # Autocorrelation X ('autoZ', np.float32), # Autocorrelation Z ('crossR', np.float32), # Cross-correlation (real) ('crossI', np.float32), # Cross-correlation (imag) ('ant', np.uint8) # Antenna selection ]) # Read the binary file with open(filepath, 'rb') as f: data = np.fromfile(f, dtype=record_dtype) return data def read_ephemeris_data(filename,year1,doy1,hour1,year2,doy2,hour2,coord_sys = "ecliptic"): # Lists to hold each column of data timestamps = [] x_rs, y_rs, z_rs, r_rs = [], [], [], [] vx_kms, vy_kms, vz_kms, vmag_kms = [], [], [], [] with open(filename, 'r') as f: lines = f.readlines() # Skip to data section: Find the line that starts with year (e.g., '2016') for i, line in enumerate(lines): if line.strip().startswith("20"): # assumes years start with "20" data_start = i break # Parse data lines for line in lines[data_start:]: parts = line.split() # Combine timestamp fields year = int(parts[0]) doy = int(parts[1]) hour = int(parts[2]) timestamp = ' '.join(parts[0:5]) # Create a comparable time interval for filtering start_time = year1 * 366*24 + doy1 * 24 + hour1 end_time = year2 * 366*24 + doy2 * 24 + hour2 current_time = year * 366*24 + doy * 24 + hour # Check if the current time is within the specified range if start_time <= current_time <= end_time: timestamps.append(timestamp) # Convert values to floats and append x_rs.append(float(parts[5])) y_rs.append(float(parts[6])) z_rs.append(float(parts[7])) r_rs.append(float(parts[8])) vx_kms.append(float(parts[9])) vy_kms.append(float(parts[10])) vz_kms.append(float(parts[11])) vmag_kms.append(float(parts[12])) # Convert timestamp list into datetime objects timestamps = [datetime.strptime(date_str, "%Y %j %H %M %S.%f") for date_str in timestamps] # Convert lists to NumPy arrays if coord_sys == "ecliptic": return { 'timestamps': timestamps, 'x_rs': np.array(x_rs), 'y_rs': np.array(y_rs), 'z_rs': np.array(z_rs), 'r_rs': np.array(r_rs), 'vx_kms': np.array(vx_kms), 'vy_kms': np.array(vy_kms), 'vz_kms': np.array(vz_kms), 'vmag_kms': np.array(vmag_kms), } elif coord_sys == "geographic": return { 'timestamps': timestamps, 'wl_iau': np.array(x_rs), 'wl_sls2': np.array(y_rs), 'wl_sls3': np.array(z_rs), 'lat': np.array(r_rs), 'L': np.array(vx_kms), 'local_time': np.array(vy_kms), 'alt': np.array(vz_kms), 'r_rs': np.array(vmag_kms), } else : raise ValueError("var: coord_sys must have value ecliptic or geographic") def read_beta(file_path, doy1, hour1, doy2, hour2): # Initialize empty lists for each of the data columns timestamps = [] beta_xw = [] beta_uw = [] beta_vw = [] with open(file_path, 'r') as file: lines = file.readlines() # Skip header and comment lines data_lines = lines[4:] # Starting from line 5 # Loop through the data lines and check if they fall within the specified time range for line in data_lines: # Split each line by whitespace and filter out empty strings values = list(filter(None, line.split(' '))) if len(values) == 7: # Extract the data (ignore year, DOY, hour, minute) year = int(values[0]) doy = int(values[1]) hour = int(values[2]) minute = int(values[3]) beta_xw_value = float(values[4]) beta_uw_value = float(values[5]) beta_vw_value = float(values[6]) # Create a comparable time interval for filtering start_time = doy1 * 24 + hour1 end_time = doy2 * 24 + hour2 current_time = doy * 24 + hour # Check if the current time is within the specified range if start_time <= current_time <= end_time: timestamp = f"{year} {doy} {hour} {minute}" timestamps.append(timestamp) beta_xw.append(beta_xw_value) beta_uw.append(beta_uw_value) beta_vw.append(beta_vw_value) # Convert lists to numpy arrays timestamps_array = [datetime.strptime(date_str, "%Y %j %H %M") for date_str in timestamps] beta_xw_array = np.array(beta_xw) beta_uw_array = np.array(beta_uw) beta_vw_array = np.array(beta_vw) return { 'timestamps': timestamps_array, 'beta_xw': beta_xw_array, 'beta_uw': beta_uw_array, 'beta_vw': beta_vw_array }