#Combined Interference and Terrain Calculation (kizer_krunch_code) # Note: the Munge cuts the A and B buildings and the INT/TER queries use them. As a result, this code below does not change. # ============================================================================ # CLIENT CONFIGURATION — change these to switch between ATT and Oncor # ============================================================================ ROSETTA_FILE = r'D:\KinderMorgan\KinderMorgan_Rosetta_populated.csv' GPD_FILE = r'D:\KinderMorgan\KinderMorgan_GPD_ALL.xlsx' ANTENNAS_DIR = r'D:\antennas\AntennasD2' WARNING_LOG = r'D:\KinderMorgan\warnings.log' NSMA_DIRS = [r'D:\antennas\Andrew', r'D:\antennas\Radiowaves', r'D:\antennas\RFS', r'D:\antennas\Gabriel'] EM_DAT_FILE = r'D:\ATT\fcc_data\EM.dat' # ============================================================================ import csv from io import StringIO from tqdm import tqdm # This is a progress bar library import math import time import traceback #used for error handling from geopy.distance import geodesic #used for a commented out line import pandas as pd import numpy as np # ============================================================================ # HEADER NAME CONFIGURATION # ---------------------------------------------------------------------------- # The two CSV files below are now read by COLUMN HEADER NAME instead of # positional index. If your spreadsheets use different header text than # what's listed here, update the strings on the right-hand side to match # the exact column headers in row 1 of each file. # ============================================================================ # --- SCE_Path_File_Naming_Convention_v8.csv (item 1) --- # Original positional indices -> header names used by this script FNC_COLUMNS = { "AReceiver": "A_Rcvr_Call_Sign", "BReceiver": "B_Rcvr_Call_Sign", "fileATerr": "A_Receiver_Terrain_File_Name", "fileBTerr": "B_Receiver_Terrain_File_Name", "fileAInt": "A_Receiver_Interference_File_Name", "fileBInt": "B_Receiver_Interference_File_Name", "amergefilename": "A_Rcvr_Merged_File", "bmergefilename": "B_Rcvr_Merged_File", "AIntFunction": "IntFunction", "ATerrFunction": "TerrFunction", "drive": "Drive", "folder": "Folder", "subfolder": "Subfolder", } # --- GeneralPathDataALL.csv (item 2) --- # Original positional order -> header names used by this script. # This file must now include a header row with these exact column names # (order no longer matters, since columns are looked up by name). GENERALPATHDATA_COLUMNS = { "tx_call_sign": "tx_call_sign", "rx_call_sign": "rx_call_sign", "txan_location_name": "txan_location_name", "rxan_location_name": "rxan_location_name", "channels": "channels", "bandwidth": "bandwidth", "fr_frequency_assigned_MHz": "fr_frequency_assigned_mhz", "rxan_gain_dBi": "rxan_gain_dBi", "RxTilt_deg": "RxTilt_deg", "RxBearing_deg": "RxBearing_deg", "AltRcvAntenna": "AltRcvAntenna", } def read_csv_safe(filepath): """Read CSV file handling encoding issues.""" for encoding in ('utf-8', 'latin-1', 'cp1252', 'utf-8-sig'): try: with open(filepath, mode='r', encoding=encoding) as f: return list(csv.reader(f)) except (UnicodeDecodeError, LookupError): continue # Last resort: read with errors ignored with open(filepath, mode='r', encoding='utf-8', errors='ignore') as f: return list(csv.reader(f)) def run_program(): while True: main() # Call the main function to execute the program # After the program finishes, ask the user if they want to run it again run_again = input("\nDo you want to run the program again? (Y/N): ").strip().lower() if run_again == 'n': input("Exiting the program. Press Enter to Exit!") break # Exit the loop and end the program elif run_again != 'y': input("Invalid input. Exiting the program. Press Enter to Exit") break #Evaluate potential interference from RLANs in buildings def main(): print("Interference/Terrain Calculation") print("Note VERSION 8 of the File naming convention is being used.") print("This program will produce A and B for both Interference and Terrain.") # Pre-load GPD into a dict keyed by (tx_call_sign, rx_call_sign) for fast lookup import re print("Loading General Path Data...") gpd_df = pd.read_excel(GPD_FILE, dtype=str) gpd_lookup = {} for _, row in gpd_df.iterrows(): tx = str(row[GENERALPATHDATA_COLUMNS["tx_call_sign"]]).strip() rx = str(row[GENERALPATHDATA_COLUMNS["rx_call_sign"]]).strip() gpd_lookup[(tx, rx)] = { "channels": row[GENERALPATHDATA_COLUMNS["channels"]], "bandwidth": row[GENERALPATHDATA_COLUMNS["bandwidth"]], "frequency": row[GENERALPATHDATA_COLUMNS["fr_frequency_assigned_MHz"]], "gain": row[GENERALPATHDATA_COLUMNS["rxan_gain_dBi"]], "tilt": row[GENERALPATHDATA_COLUMNS["RxTilt_deg"]], "bearing": row[GENERALPATHDATA_COLUMNS["RxBearing_deg"]], "ant_model": str(row[GENERALPATHDATA_COLUMNS["AltRcvAntenna"]]).strip(), "tx_location_name": str(row[GENERALPATHDATA_COLUMNS["txan_location_name"]]).strip(), "rx_location_name": str(row[GENERALPATHDATA_COLUMNS["rxan_location_name"]]).strip(), } # Count actual number of channels per call sign from EM.dat print("Loading channel counts from EM.dat...") em_channels = {} try: with open(EM_DAT_FILE, encoding='latin-1') as f: for line in f: c = line.rstrip('\r\n').split('|') if len(c) < 13: continue call = c[4].strip() ch_num = c[12].strip() if call and ch_num: try: ch_int = int(ch_num) if call not in em_channels or ch_int > em_channels[call]: em_channels[call] = ch_int except ValueError: pass # Update gpd_lookup with correct channel counts for key, rec in gpd_lookup.items(): tx = key[0] if tx in em_channels: rec["channels"] = str(em_channels[tx]) except FileNotFoundError: print("Warning: EM.dat not found, using channel values from GPD.") print(f"GPD loaded: {len(gpd_lookup)} paths.") # Load manufacturer midband gains from AntennasD2_MidBandGain.csv print("Loading antenna midband gains...") nsma_gain_lookup = {} MIDBAND_GAIN_FILE = r'D:\antennas\AntennasD2_MidBandGain.csv' try: mg_df = pd.read_csv(MIDBAND_GAIN_FILE, dtype=str) for _, row in mg_df.iterrows(): gain = float(row['Mid_Band_Gain_dBi']) wc = str(row['WC_Filename']).strip() import re as _re key = _re.sub(r'-WC\.csv$', '', wc, flags=_re.IGNORECASE).strip() nsma_gain_lookup[key] = gain model = str(row['Model']).strip() if model and model != 'nan': nsma_gain_lookup[model] = gain print(f"Midband gains loaded: {len(nsma_gain_lookup)} entries.") except FileNotFoundError: print(f"Warning: {MIDBAND_GAIN_FILE} not found. Using FCC-filed gains.") try: # Load Rosetta (naming convention) filenaming = pd.read_csv(ROSETTA_FILE) #get the path number, and watch the row number to make sure they don't need an offset match_flag = False valid_path_nums = filenaming['path_num'].astype(str).str.strip().tolist() while not match_flag: pathnumber = input("Enter the PATH number from the NAMING CONVENTIONS SPREADSHEET: ").strip() if pathnumber not in valid_path_nums: print(f"Path number not found. Valid values: {valid_path_nums[0]} to {valid_path_nums[-1]}") print(f"Try again.") else: match_flag = True # look up row by path_num value rownumber = filenaming.index[filenaming['path_num'].astype(str).str.strip() == str(pathnumber)][0] AReceiver = str(filenaming.loc[rownumber, FNC_COLUMNS["AReceiver"]]) BReceiver = str(filenaming.loc[rownumber, FNC_COLUMNS["BReceiver"]]) amergefilename = str(filenaming.loc[rownumber, FNC_COLUMNS["amergefilename"]]) bmergefilename = str(filenaming.loc[rownumber, FNC_COLUMNS["bmergefilename"]]) print() # ----------------------------------------------------------------------------------- # these four are used to calculate composite interference. AIntFunction = str(filenaming.loc[rownumber, FNC_COLUMNS["AIntFunction"]]) BIntFunction = str(filenaming.loc[rownumber, FNC_COLUMNS["AIntFunction"]]) ATerrFunction = str(filenaming.loc[rownumber, FNC_COLUMNS["ATerrFunction"]]) BTerrFunction = str(filenaming.loc[rownumber, FNC_COLUMNS["ATerrFunction"]]) drive = str(filenaming.loc[rownumber, FNC_COLUMNS["drive"]]).strip() folder = str(filenaming.loc[rownumber, FNC_COLUMNS["folder"]]).strip().strip('\\') subfolder = str(filenaming.loc[rownumber, FNC_COLUMNS["subfolder"]]).strip() #this gets the file names (full paths from Rosetta) fileAInt = str(filenaming.loc[rownumber, FNC_COLUMNS["fileAInt"]]) fileBInt = str(filenaming.loc[rownumber, FNC_COLUMNS["fileBInt"]]) fileATerr = str(filenaming.loc[rownumber, FNC_COLUMNS["fileATerr"]]) fileBTerr = str(filenaming.loc[rownumber, FNC_COLUMNS["fileBTerr"]]) # input files - full paths come directly from Rosetta AIntInputfile = fileAInt BIntInputfile = fileBInt ATerrInputfile = fileATerr BTerrInputfile = fileBTerr # output files - same folder, A and B distinguished by receiver call sign base_path = f"{drive}\\{folder}\\{subfolder}" AIntOutputfile = f"{base_path}\\{subfolder}_A_{AIntFunction}.csv" BIntOutputfile = f"{base_path}\\{subfolder}_B_{BIntFunction}.csv" ATerrOutputfile = f"{base_path}\\{subfolder}_A_{ATerrFunction}.csv" BTerrOutputfile = f"{base_path}\\{subfolder}_B_{BTerrFunction}.csv" # configuration files. genpathdata = f"{drive}:\\{folder}\\ + 'AntennasD2csv'" print(f"Here's what about to run...") print(f"AReceiver: ", AReceiver) print(f"BReceiver: ", BReceiver) print(f"A Interference Input File: ", AIntInputfile) print(f"B Interference Input File: ", BIntInputfile) print(f"A Interference Output File: ", AIntOutputfile) print(f"B Interference Output File: ", BIntOutputfile) print(f"A Terrain Input File: ", ATerrInputfile) print(f"B Terrain Input File: ", BTerrInputfile) print(f"A Terrain Output File: ", ATerrOutputfile) print(f"B Terrain Output File: ", BTerrOutputfile) print(f"A merge Output file: ", amergefilename) print(f"B merge Output file: ", bmergefilename) print() print("Configuration Files are here: ") print(f"Rosetta Stone: {ROSETTA_FILE}") print(f"General Path Data: {GPD_FILE}") print(f"Antenna Patterns: {ANTENNAS_DIR}") print() print(f"If you need to change file names or versions, do it in the Rosetta file.") #print(f"The output files will be in a folder called MERGED in the PATH directory") # Open files for interference paths of interest and evaluations interference_path_inputA = AIntInputfile interference_path_inputB = BIntInputfile #interference_path_input = "C:\SCEProject\DataSCE14\WPNN436_WPOS895\WPNN436_WPOS895_A_INT_1-005.csv" interference_path_outputA = AIntOutputfile interference_path_outputB = BIntOutputfile # Open files for interference paths of interest and evaluations terrain_path_inputA = ATerrInputfile terrain_path_inputB = BTerrInputfile #terrain_path_input = "C:\SCEProject\DataSCE14\WPNN436_WPOS895\WPNN436_WPOS895_A_INT_1-005.csv" terrain_path_outputA = ATerrOutputfile terrain_path_outputB = BTerrOutputfile except Exception as e: print(f"An error occurred: {e}") traceback.print_exc() raise # Ask the user if they want to proceed with the operation while True: proceed = input("Proceed? Enter Y or N and then press Enter: ") if proceed.lower() == "y": print() print("Beginning Interference/Terrain Calculation") break elif proceed.lower() == "n": print("Please enter the path number again.") # If they enter 'N', it will re-prompt for the path number by restarting the `main()` logic main() return # Make sure to return after restarting to avoid executing the rest of the function. else: print("Invalid input, please enter 'Y' or 'N'.") # Ask the user what they want to run: Interference, Terrain, or Both while True: choice = input("Do you want to run Interference, Terrain, or Both? (Enter I, T, or B): ") if choice.lower() == 'i': print("You have chosen to run Interference calculations.") run_interference = True run_terrain = False break elif choice.lower() == 't': print("You have chosen to run Terrain calculations.") run_interference = False run_terrain = True break elif choice.lower() == 'b': print("You have chosen to run both Interference and Terrain calculations.") run_interference = True run_terrain = True break else: print("Invalid input. Please enter 'I' for Interference, 'T' for Terrain, or 'B' for Both.") start_time = time.time() # Start the timer # Only run interference calculations if the user selects Interference or Both if run_interference: #Runs A Interference Calculation print("Running A Interference Calculation:") print("Interference Path Input: ", interference_path_inputA) print("Interference Path Output: ", interference_path_outputA) run_interference_calculation(interference_path_inputA, interference_path_outputA, drive, folder, gpd_lookup, nsma_gain_lookup=nsma_gain_lookup) #Runs B Interference Calculation print("Running B Interference Calculation") print("Interference Path Input: ", interference_path_inputB) print("Interference Path Output: ", interference_path_outputB) run_interference_calculation(interference_path_inputB, interference_path_outputB, drive, folder, gpd_lookup, nsma_gain_lookup=nsma_gain_lookup) # Only run terrain calculations if the user selects Terrain or Both if run_terrain: print() print("Running A Terrain Calculation:") print("Terrain Path Input: ", terrain_path_inputA) print("Terrain Path Output: ", terrain_path_outputA) run_terrain_calculation(terrain_path_inputA, terrain_path_outputA) #Runs B Terrain Calculation print("Running B Terrain Calculation:") print("Terrain Path Input: ", terrain_path_inputB) print("Terrain Path Output: ", terrain_path_outputB) run_terrain_calculation(terrain_path_inputB, terrain_path_outputB) # Stop the timer and print the elapsed time since Main function is now ended # At the end of your script, after all processing is done: endtime = time.time() elapsed_time = endtime - start_time # Convert elapsed time to hours, minutes, and seconds hours, rem = divmod(elapsed_time, 3600) minutes, seconds = divmod(rem, 60) # Print the formatted elapsed time print(f"Elapsed time: {int(hours)} hours, {int(minutes)} minutes, {seconds:.2f} seconds") print("Program Completed") #Interference Calculation def run_interference_calculation(interference_path_input, interference_path_output, drive, folder, gpd_lookup, preloaded_data=None, nsma_gain_lookup=None): if nsma_gain_lookup is None: nsma_gain_lookup = {} PI = 3.1415926 try: all_rows = preloaded_data if preloaded_data is not None else read_csv_safe(interference_path_input) if not all_rows or len(all_rows) < 2: # Empty file - write header only with open(interference_path_output, mode='w', newline='') as f: csv.writer(f).writerow([ "Building", "Receiver", "Transmitter", "Site Name and Direction", "Raw Interference (dBm)", "RFM 10 dB BEL", "RFM 20 dB BEL", "RFM 30 dB BEL", "Longitude (deg)", "Latitude (deg)", "Height (ft)", "Interference Address", "Interference Building Type", "Channels", "Ch Bandwidth", "Int Dist (miles)", "Int Bearing", "Int Tilt" ]) msg = f"WARNING: Empty interference file: {interference_path_input}" print(msg) with open(WARNING_LOG, 'a') as wf: wf.write(msg + '\n') return data_rows = all_rows[1:] total_lines = len(data_rows) with open(interference_path_output, mode='w', newline='') as output_file: csv_writer = csv.writer(output_file) csv_writer.writerow([ "Building", "Receiver", "Transmitter", "Site Name and Direction", "Raw Interference (dBm)", "RFM 10 dB BEL", "RFM 20 dB BEL", "RFM 30 dB BEL", "Longitude (deg)", "Latitude (deg)", "Height (ft)", "Interference Address", "Interference Building Type", "Channels", "Ch Bandwidth", "Int Dist (miles)", "Int Bearing", "Int Tilt" ]) pbar = tqdm(total=total_lines, desc="Processing Rows", bar_format="{desc}: {percentage:3.0f}%|{bar}|") # State for current building group cur_building = None cur_rx_site = None cur_tx_site = None interference_dbm = 0.0 int_hor_off_axis = 0.0 InterferenceDistanceMiles = 0.0 InterferenceTiltDeg = 0.0 print_file = '' site_name_and_direction = '' channels = '' bandwidth = '' # GPD lookup results cached per path gpd_cache = {} def flush_building(): """Write the current building's result to output.""" if cur_building is None: return rfn10, rfn20, rfn30 = determine_fade_margin_effect(interference_dbm) more_parts = [more_custom_format(v) for v in [rfn10, rfn20, rfn30]] print_file_parts = next(csv.reader(StringIO(print_file))) csv_writer.writerow([ cur_building, cur_rx_site, cur_tx_site, site_name_and_direction, f"{interference_dbm:.4f}", *more_parts, *print_file_parts, channels, bandwidth, f"{float(InterferenceDistanceMiles):.5f}", int_hor_off_axis, f"{float(InterferenceTiltDeg):.6f}" ]) for row in data_rows: pbar.update(1) try: fields = [f.strip() for f in row] if len(fields) < 15: continue (a_rcvr_path_id, bldg_number, bldg_long, bldg_lat, bldg_maxheight_ft, bldg_meanelev_ft, bldg_desc, bldg_type, a_rcvr_call_sign, a_rcvr_desc, a_rcvr_long, a_rcvr_lat, a_rcvr_grd_elev_ft, a_rcvr_ant_raat_ft, bldg_a_rcvr_distance_km) = fields[:15] # Sanitize numerics bldg_maxheight_ft = str(max(10.0, float(bldg_maxheight_ft) or 10.0)) bldg_meanelev_ft = f"{float(bldg_meanelev_ft):.4f}" bldg_long = f"{float(bldg_long):.6f}" bldg_lat = f"{float(bldg_lat):.6f}" a_rcvr_long = f"{float(a_rcvr_long):.6f}" a_rcvr_lat = f"{float(a_rcvr_lat):.6f}" a_rcvr_grd_elev_ft = f"{float(a_rcvr_grd_elev_ft):.4f}" a_rcvr_ant_raat_ft = f"{float(a_rcvr_ant_raat_ft):.4f}" bldg_a_rcvr_distance_km = f"{float(bldg_a_rcvr_distance_km):.4f}" building, rx_site, tx_site = parse_sites_simplified(a_rcvr_path_id) if building != cur_building: # Write previous building flush_building() # Reset for new building int_path_dist_miles = float(bldg_a_rcvr_distance_km) * 0.621371 ant_ht_ft = (3.28084 * float(a_rcvr_grd_elev_ft)) + (3.28084 * float(a_rcvr_ant_raat_ft)) int_ht_ft = float(bldg_maxheight_ft) + float(bldg_meanelev_ft) - 5.0 print_file = f'{bldg_long},{bldg_lat},{int(float(bldg_maxheight_ft))},"{bldg_desc}","{bldg_type}"' # Bearings and tilts _, rx_bearing_deg = calculate_distance_and_bearing(bldg_lat, bldg_long, a_rcvr_lat, a_rcvr_long) TxTilt_deg, rx_tilt_deg, InterferenceDistanceMiles, InterferenceTiltDeg = calculate_antenna_tilt( bldg_meanelev_ft, bldg_maxheight_ft, a_rcvr_grd_elev_ft, a_rcvr_ant_raat_ft, bldg_a_rcvr_distance_km, a_rcvr_path_id, building) # GPD lookup (cached per rx/tx pair) cache_key = (rx_site, tx_site) if cache_key not in gpd_cache: result = find_mw_path_information(gpd_lookup, rx_site, tx_site) gpd_cache[cache_key] = result gpd_result = gpd_cache[cache_key] if gpd_result is None: cur_building = building cur_rx_site = rx_site cur_tx_site = tx_site continue (rcv_gain, rcv_tilt, rcv_bearing, rcv_ant_model, fr_frequency_assigned_MHzMW, rxan_gain_dBiMW, channels, bandwidth, site_name_and_direction) = gpd_result # Off-axis angle hor_off_axis = abs(float(rcv_bearing) - float(rx_bearing_deg)) if hor_off_axis > 180: hor_off_axis = 360 - hor_off_axis int_hor_off_axis = round(hor_off_axis, 1) # 3D total angle alpha_i = (PI/180) * float(rx_bearing_deg) alpha_d = (PI/180) * float(rcv_bearing) delta_d = (PI/180) * float(rcv_tilt) phi_i = (PI/180) * float(rx_tilt_deg) x = (math.cos(delta_d)*math.cos(alpha_d)*math.cos(phi_i)*math.cos(alpha_i) + math.cos(delta_d)*math.sin(alpha_d)*math.cos(phi_i)*math.sin(alpha_i) + math.sin(delta_d)*math.sin(phi_i)) x = max(-1.0, min(1.0, x)) total_angle = math.degrees(math.acos(x)) # Antenna discrimination and interference ant_disc = determine_antenna_discrimination(drive, folder, rcv_ant_model, total_angle) effective_gain = nsma_gain_lookup.get(rcv_ant_model, float(rxan_gain_dBiMW)) fsl_db = 96.58 + 20*math.log10(float(fr_frequency_assigned_MHzMW)/1000) + 20*math.log10(int_path_dist_miles) interference_dbm = round(19.7 - fsl_db + effective_gain + float(ant_disc), 5) cur_building = building cur_rx_site = rx_site cur_tx_site = tx_site else: # Same building — apply knife-edge loss int_path_dist_miles = float(bldg_a_rcvr_distance_km) * 0.621371 ant_ht_ft = (3.28084 * float(a_rcvr_grd_elev_ft)) + (3.28084 * float(a_rcvr_ant_raat_ft)) int_ht_ft = float(bldg_maxheight_ft) + float(bldg_meanelev_ft) - 5.0 GKE, *_ = knife_edge_loss(int_path_dist_miles, bldg_a_rcvr_distance_km, fr_frequency_assigned_MHzMW, bldg_maxheight_ft, bldg_meanelev_ft, ant_ht_ft, int_ht_ft) interference_dbm = round(interference_dbm + GKE, 4) except Exception as row_err: print(f" ERROR in row {fields[0] if fields else "?"} of {interference_path_input}: {row_err}") import traceback; traceback.print_exc() raise # Write the last building flush_building() pbar.close() except FileNotFoundError: print(f"Error: File not found: {interference_path_input}") raise except Exception as e: print(f"An error occurred in interference calculation: {e}") import traceback; traceback.print_exc() raise #formats RFM 10 dB BEL RFM 20 dB BEL RFM 30 dB BEL to 1 decimal place #(Interference/Terrain) def more_custom_format(value): if value == int(value): return str(int(value)) else: return str(value).lstrip('0') #1100, Parse building, receive site and trasmit site (Interference/Terrain) def parse_sites_simplified(a_rcvr_path_id): # Split the input string by the period ('.') delimiter parts = a_rcvr_path_id.split('.') # Check that the input string is in the correct format with exactly three parts if len(parts) != 3: raise ValueError("Input string must be in the format 'Building.RxSite.TxSite'") # Assign each part to its respective variable building = parts[0] # First part is the Building rx_site = parts[1] # Second part is the RxSite (Receive Site) tx_site = parts[2] # Third part is the TxSite (Transmit Site) # Return the parsed components as a tuple return building, rx_site, tx_site #9000, path distance and site bearing (Interference) def calculate_distance_and_bearing(bldg_lat, bldg_long, a_rcvr_lat, a_rcvr_long): #print the lat and long lat_a = float(bldg_lat) lon_a = float(bldg_long) lat_b = float(a_rcvr_lat) lon_b = float(a_rcvr_long) # lat_a = np.float32(bldg_lat) # lon_a = np.float32(bldg_long) # lat_b = np.float32(a_rcvr_lat) # lon_b = np.float32(a_rcvr_long) # Check if the two sites are the same if abs(lat_a - lat_b) < 0.000001 and abs(lon_a - lon_b) < 0.000001: bearing_a = "0.0" bearing_b = "0.0" distance = "0.0" TxBearing_deg = bearing_a RxBearing_deg = bearing_b return TxBearing_deg, RxBearing_deg # Constants PI = 3.1415926 oneeightyDivPI = 180.0/PI # Calculate intermediate terms for bearing calculations a0 = math.sin(PI * lat_a / 180) * math.sin(PI * lat_b / 180) + \ math.cos(PI * lat_a / 180) * math.cos(PI * lat_b / 180) * \ math.cos(PI * (lon_a - lon_b) / 180) # a0 = math.sin(np.float32(np.float32(PI) * np.float32(lat_a) / np.float32(180))) * \ # math.sin(np.float32(np.float32(PI) * np.float32(lat_b) / np.float32(180))) + \ # math.cos(np.float32(np.float32(PI) * np.float32(lat_a) / np.float32(180))) * \ # math.cos(np.float32(np.float32(PI) * np.float32(lat_b) / np.float32(180))) * \ # math.cos(np.float32(np.float32(PI) * (np.float32(lon_a) - np.float32(lon_b)) / np.float32(180))) y=math.acos(a0) #need to fix y # y = math.acos(np.float32(a0)) # print(y) #print(a0, y) #qbasic ex. .0010335983 #python ex. .000935331 #math acos, arccos, precise arccos a1 = (180.0 / PI) * y # a1 = np.float32((np.float32(180.0) / np.float32(PI)) * np.float32(y)) #a1 = round(a1, 7) # Round a1 to match precision # Calculate bearing from A to B denominator = math.sin(PI * a1 / 180) * math.cos(PI * lat_a / 180) # Check if the denominator is zero if denominator == 0: #raise ValueError("Denominator is zero, cannot perform division. Check your input values.") a2 = 0 else: a2 = (math.sin(PI * lat_b / 180) - math.sin(PI * lat_a / 180) * math.cos(PI * a1 / 180)) / \ (math.sin(PI * a1 / 180) * math.cos(PI * lat_a / 180)) # a2 = (math.sin(np.float32(PI) * np.float32(lat_b) / np.float32(180)) - # math.sin(np.float32(PI) * np.float32(lat_a) / np.float32(180)) * # math.cos(np.float32(PI) * np.float32(a1) / np.float32(180))) / \ # (math.sin(np.float32(PI) * np.float32(a1) / np.float32(180)) * # math.cos(np.float32(PI) * np.float32(lat_a) / np.float32(180))) #a2 = round(a2, 6) #y = math.acos(a2) if -1 <= a2 <= 1: y = math.acos(a2) # y = math.acos(np.float32(a2)) else: if a2 < -1: a2 = -1 if a2 > 1: a2 = 1 y = math.acos(a2) # y = math.acos(np.float32(a2)) bearing_a = (180.0 / PI) * y # bearing_a = np.float32(np.float32(180.0) / np.float32(PI)) * np.float32(y) denominator = (math.sin(PI * a1 / 180) * math.cos(PI * lat_b / 180)) # Check if the denominator is zero if denominator == 0: #raise ValueError("Denominator is zero, cannot perform division. Check your input values.") a4 = 0 else: # Calculate bearing from B to A a4 = (math.sin(PI * lat_a / 180) - math.sin(PI * lat_b / 180) * math.cos(PI * a1 / 180)) / (math.sin(PI * a1 / 180) * math.cos(PI * lat_b / 180)) # a4 = (math.sin(np.float32(PI) * np.float32(lat_a) / np.float32(180)) - # math.sin(np.float32(PI) * np.float32(lat_b) / np.float32(180)) * # math.cos(np.float32(PI) * np.float32(a1) / np.float32(180))) / \ # (math.sin(np.float32(PI) * np.float32(a1) / np.float32(180)) * # math.cos(np.float32(PI) * np.float32(lat_b) / np.float32(180))) #y = math.acos(a4) if -1 <= a4 <= 1: y = math.acos(a4) # y = math.acos(np.float32(a4)) else: if a4 < -1: a4 = -1 if a4 > 1: a4 = 1 y = math.acos(a4) # y = math.acos(np.float32(a4)) bearing_b = (180.0 / PI) * y # bearing_b = np.float32(np.float32(180.0) / np.float32(PI)) * np.float32(y) if math.sin(PI * (lon_b - lon_a) / 180) < 0: #< isntead of >= but works bearing_a = 360 - bearing_a else: bearing_b = 360 - bearing_b # if math.sin(np.float32(PI) * (np.float32(lon_b) - np.float32(lon_a)) / np.float32(180)) < 0: # bearing_a = np.float32(360) - np.float32(bearing_a) # else: # bearing_b = np.float32(360) - np.float32(bearing_b) # Calculate the great-circle distance #distance doesnt even get used so it doesnt matter #distance = geodesic((lat_a, lon_a), (lat_b, lon_b)).miles distance = calculate_distance(lat_a, lat_b, lon_a, lon_b, unit='M') #print(distance) distance = int((distance * 1000) + .5) distance = distance / 1000 distance = str(distance) bearing_a = int((bearing_a * 10) + .5) bearing_a = bearing_a / 10 bearing_a = str(bearing_a) bearing_b = int((bearing_b * 10) + .5) bearing_b = bearing_b / 10 bearing_b = str(bearing_b) TxBearing_deg = bearing_a RxBearing_deg = bearing_b # print(f"TxBearing_deg: {TxBearing_deg} RxBearing_deg: {RxBearing_deg}") return TxBearing_deg, RxBearing_deg #9100 #should run fine with the other function, delete this later # def calculate_distance(lat_a, lat_b, lon_a, lon_b, unit='M'): # PI = 3.1415926 # # High accuracy formula # z = (math.sin((PI * (lat_a - lat_b) / 180.) / 2)) ** 2 # z += math.cos(PI * lat_a / 180) * math.cos(PI * lat_b / 180) * (math.sin((PI * (lon_a - lon_b) / 180) / 2)) ** 2 # z = math.sqrt(z) # # Calculate arcsin # z_short = 2 * (180 / PI) * math.asin(z) # # Distance in kilometers and miles # z_km = 111.1 * z_short # z_miles = 69.06 * z_short # if unit == "M": # distance = z_miles # elif unit == "K": # distance = z_km # else: # raise ValueError("Unit must be 'M' for miles or 'K' for kilometers") # return distance #9200/9400, arccosine function example #test if this works the same for both functions #maybe i try out both python functions that were translated, then test if arccos works #result = math.acos(x) #(Interference) def custom_copysign(x): if x >= 0: return 1 else: return -1 #9400, arccosine function example (Interference) def ARCCOS(X): PI = 3.1415926 if abs(X) > 1: XMOD = X / abs(X) else: XMOD = X if XMOD == 1: ACOS = 0 elif XMOD == -1: ACOS = PI elif abs(XMOD) < 1: ACOS = PI / 2 - math.atan(XMOD / math.sqrt(1 - XMOD * XMOD)) return ACOS #9300, arcsin function example #result = math.asin(1, 3.1415926) #9500, calculate antenna tilt (Interference) def calculate_antenna_tilt(bldg_meanelev_ft, bldg_maxheight_ft, a_rcvr_grd_elev_ft, a_rcvr_ant_raat_ft, bldg_a_rcvr_distance_km, a_rcvr_path_id, old_building): # Constants PI = 3.1415926 K = 4.0 / 3.0 # Convert heights to meters and kilometers to miles Htx = (float(bldg_meanelev_ft) + float(bldg_maxheight_ft) - 5.0) * 0.3048 # meters Hrx = (float(a_rcvr_grd_elev_ft)) + (float(a_rcvr_ant_raat_ft)) # meters PathDistance = float(bldg_a_rcvr_distance_km) # kilometers # Calculate interference distance in miles if conditions match if a_rcvr_path_id != old_building: InterferenceDistanceMiles = str(0.621371 * PathDistance).strip() else: InterferenceDistanceMiles = "" # Calculate antenna tilt for transmitter H = Hrx - Htx Phi1 = 0.001 * H / PathDistance Phi2 = (math.sqrt((PathDistance * PathDistance) + (0.000001 * H * H))) / (12735. * K) Phi1 = math.atan(Phi1) # arctangent Phi2 = math.atan(Phi2 / math.sqrt(1. - (Phi2 * Phi2))) # arcsine AntTiltTx = (Phi1 - Phi2) * (180. / PI) TxTilt_deg = str(AntTiltTx).strip() # Update interference tilt if conditions match if a_rcvr_path_id != old_building: InterferenceTiltDeg = TxTilt_deg else: InterferenceTiltDeg = "" # Calculate antenna tilt for receiver H = Htx - Hrx Phi1 = 0.001 * H / PathDistance Phi2 = (math.sqrt((PathDistance * PathDistance) + (0.000001 * H * H))) / (12735. * K) Phi1 = math.atan(Phi1) # arctangent Phi2 = math.atan(Phi2 / math.sqrt(1. - (Phi2 * Phi2))) # arcsine AntTiltRx = (Phi1 - Phi2) * (180. / PI) RxTilt_deg = str(AntTiltRx).strip() return f"{float(TxTilt_deg):.6f}", f"{float(RxTilt_deg):.6f}", f"{float(InterferenceDistanceMiles):.6f}", f"{float(InterferenceTiltDeg):.6f}" #9600, find desired mw path information (Interference) def find_mw_path_information(gpd_lookup, rx_site, tx_site): rec = gpd_lookup.get((tx_site, rx_site)) if rec is None: print(f"Subroutine 9600 failed to find matching mw path <{tx_site}> <{rx_site}>.") print("Program terminated.") return None rcv_gain = rec["gain"] rcv_tilt = rec["tilt"] rcv_bearing = rec["bearing"] rcv_ant_model = rec["ant_model"] fr_frequency_assigned_MHzMW = rec["frequency"] rxan_gain_dBiMW = rec["gain"] channels = rec["channels"] bandwidth = rec["bandwidth"] tx_location_name = rec["tx_location_name"] rx_location_name = rec["rx_location_name"] site_name_and_direction = f"{tx_location_name} TO {rx_location_name}" return rcv_gain, rcv_tilt, rcv_bearing, rcv_ant_model, fr_frequency_assigned_MHzMW, rxan_gain_dBiMW, channels, bandwidth, site_name_and_direction #9700, determine antenna discrimination #(Interference) #takes the value from the spreedsheet, then does some calculations #to see which one gets picked, arc coseine function used for total_angle def determine_antenna_discrimination(drive, folder, rcv_ant_model, total_angle): total_angle = round(total_angle, 6) antenna_file = ANTENNAS_DIR + '\\' + rcv_ant_model + '-WC.csv' ant_disc = "0" old_delta = 999.0 try: with open(antenna_file, mode='r') as file: csv_reader = csv.reader(file) for row in csv_reader: row = [field.strip() for field in row if field.strip() != ''] if len(row) < 2: continue an_angle, a_value = row[0], row[1] if not an_angle or not a_value: continue new_delta = abs(float(an_angle) - total_angle) if new_delta > old_delta: break if new_delta < old_delta: old_delta = new_delta ant_disc = a_value #catches if file not found except FileNotFoundError: print(f"Error: The file '{antenna_file}' was not found.") input("Press Enter to exit...") return None #catches any other error except Exception as e: print(f"An error occurred: {e}") import traceback print(traceback.format_exc()) raise return ant_disc #9800, knife edge loss (dB) (Interference/Terrain) #needs to be tested def knife_edge_loss(int_path_dist_miles, bldg_a_rcvr_distance_km, fr_frequency_assigned_mhz, bldg_maxheight_ft, bldg_meanelev_ft, ant_ht_ft, int_ht_ft): pi = 3.1415926 # Path distance in miles d = int_path_dist_miles # Distance from one end to point on path in miles d1 = 0.621371 * float(bldg_a_rcvr_distance_km) if d1 < 0.00947: d1 = 0.00947 # Distance from other end to point on path in miles d2 = d - d1 if d2 < 0.00947: d2 = 0.00947 # Operating frequency in GHz freq = float(fr_frequency_assigned_mhz) / 1000.0 # First Fresnel zone in feet f1 = 72.1 * math.sqrt(d1 * d2 / (freq * d)) # Building obstruction clearance in feet h1 = float(bldg_maxheight_ft) + float(bldg_meanelev_ft) # Propagation height factor h2 = (d1 * d2) / (1.5 * (4.0 / 3.0)) # Current distance from victim antenna in miles current_distance_from_victim_antenna_miles = 0.621371 * float(bldg_a_rcvr_distance_km) # Path height in feet path_height_ft = ant_ht_ft + (((int_ht_ft - ant_ht_ft) / int_path_dist_miles) * current_distance_from_victim_antenna_miles) # Calculate H #print(f"h1: {h1}, h2: {h2}, path_height_ft: {path_height_ft}") h = path_height_ft - (float(h1) + h2) # Calculate X x = h / f1 # Calculate knife edge loss (GKE) if x > 0.6: gke = 0.0 elif x < -600.0: gke = -73.0 elif x <= 0.6 and x >= -1.0: a = -6.013050244594124 b = 12.56373978136271 c = -4.638348100303959 d = -14.27304778829325 e = 4.624584658556855 f = 26.08257628627929 g = 14.16222106310102 gke = a + b * x + c * x**2 + d * x**3 + e * x**4 + f * x**5 + g * x**6 elif x <= -1.0 and x >= -600.0: a = -10.13843752145488 b = -0.2211687451291319 c = 10.31071865966655 d = 0.0009104403873283171 e = -0.0762342416178572 gke = (a + c * x + e * x**2) / (1.0 + b * x + d * x**2) else: gke = 0.0 return gke, d1, d2, d, path_height_ft, h1, h2, h, x, f1 #9900, determine effect of interference on fade margin (Interference/Terrain) def determine_fade_margin_effect(interference_dbm): rcv_noise = -94.0 interference_10db = interference_dbm - 10.0 interference_20db = interference_dbm - 20.0 interference_30db = interference_dbm - 30.0 log10 = math.log(10) rfn_10db = (10.0 * math.log10((10.0 ** (rcv_noise / 10.0)) + (10.0 ** (interference_10db / 10.0))) ) - rcv_noise rfn_20db = (10.0 * math.log10((10.0 ** (rcv_noise / 10.0)) + (10.0 ** (interference_20db / 10.0))) ) - rcv_noise rfn_30db = (10.0 * math.log10((10.0 ** (rcv_noise / 10.0)) + (10.0 ** (interference_30db / 10.0))) ) - rcv_noise # if rfn_10db < 0.1: # rfn_10db = 0.1 # if rfn_20db < 0.1: # rfn_20db = 0.1 # if rfn_30db < 0.1: # rfn_30db = 0.1 # Apply the one significant decimal conversion # rfn_10db = convert_to_one_significant_decimal(rfn_10db) # rfn_20db = convert_to_one_significant_decimal(rfn_20db) # rfn_30db = convert_to_one_significant_decimal(rfn_30db) rfn_10db = round(rfn_10db, 1) rfn_20db = round(rfn_20db, 1) rfn_30db = round(rfn_30db, 1) return rfn_10db, rfn_20db, rfn_30db #9950, convert to one significant decimal (Interference) def convert_to_one_significant_decimal(variable): variable = int(0.5 + (variable * 10.0)) variable = variable / 10.0 return variable #9000, path distance and site bearing (Interference) def calculate_distance_and_bearing(bldg_lat, bldg_long, a_rcvr_lat, a_rcvr_long): #print the lat and long lat_a = float(bldg_lat) lon_a = float(bldg_long) lat_b = float(a_rcvr_lat) lon_b = float(a_rcvr_long) # lat_a = np.float32(bldg_lat) # lon_a = np.float32(bldg_long) # lat_b = np.float32(a_rcvr_lat) # lon_b = np.float32(a_rcvr_long) # Check if the two sites are the same if abs(lat_a - lat_b) < 0.000001 and abs(lon_a - lon_b) < 0.000001: bearing_a = "0.0" bearing_b = "0.0" distance = "0.0" TxBearing_deg = bearing_a RxBearing_deg = bearing_b return TxBearing_deg, RxBearing_deg # Constants PI = 3.1415926 oneeightyDivPI = 180.0/PI # Calculate intermediate terms for bearing calculations a0 = math.sin(PI * lat_a / 180) * math.sin(PI * lat_b / 180) + \ math.cos(PI * lat_a / 180) * math.cos(PI * lat_b / 180) * \ math.cos(PI * (lon_a - lon_b) / 180) # a0 = math.sin(np.float32(np.float32(PI) * np.float32(lat_a) / np.float32(180))) * \ # math.sin(np.float32(np.float32(PI) * np.float32(lat_b) / np.float32(180))) + \ # math.cos(np.float32(np.float32(PI) * np.float32(lat_a) / np.float32(180))) * \ # math.cos(np.float32(np.float32(PI) * np.float32(lat_b) / np.float32(180))) * \ # math.cos(np.float32(np.float32(PI) * (np.float32(lon_a) - np.float32(lon_b)) / np.float32(180))) y=math.acos(a0) #need to fix y # y = math.acos(np.float32(a0)) # print(y) #print(a0, y) #qbasic ex. .0010335983 #python ex. .000935331 #math acos, arccos, precise arccos a1 = (180.0 / PI) * y # a1 = np.float32((np.float32(180.0) / np.float32(PI)) * np.float32(y)) #a1 = round(a1, 7) # Round a1 to match precision # Calculate bearing from A to B denominator = math.sin(PI * a1 / 180) * math.cos(PI * lat_a / 180) # Check if the denominator is zero if denominator == 0: #raise ValueError("Denominator is zero, cannot perform division. Check your input values.") a2 = 0 else: a2 = (math.sin(PI * lat_b / 180) - math.sin(PI * lat_a / 180) * math.cos(PI * a1 / 180)) / \ (math.sin(PI * a1 / 180) * math.cos(PI * lat_a / 180)) # a2 = (math.sin(np.float32(PI) * np.float32(lat_b) / np.float32(180)) - # math.sin(np.float32(PI) * np.float32(lat_a) / np.float32(180)) * # math.cos(np.float32(PI) * np.float32(a1) / np.float32(180))) / \ # (math.sin(np.float32(PI) * np.float32(a1) / np.float32(180)) * # math.cos(np.float32(PI) * np.float32(lat_a) / np.float32(180))) #a2 = round(a2, 6) #y = math.acos(a2) if -1 <= a2 <= 1: y = math.acos(a2) # y = math.acos(np.float32(a2)) else: if a2 < -1: a2 = -1 if a2 > 1: a2 = 1 y = math.acos(a2) # y = math.acos(np.float32(a2)) bearing_a = (180.0 / PI) * y # bearing_a = np.float32(np.float32(180.0) / np.float32(PI)) * np.float32(y) denominator = (math.sin(PI * a1 / 180) * math.cos(PI * lat_b / 180)) # Check if the denominator is zero if denominator == 0: #raise ValueError("Denominator is zero, cannot perform division. Check your input values.") a4 = 0 else: # Calculate bearing from B to A a4 = (math.sin(PI * lat_a / 180) - math.sin(PI * lat_b / 180) * math.cos(PI * a1 / 180)) / (math.sin(PI * a1 / 180) * math.cos(PI * lat_b / 180)) # a4 = (math.sin(np.float32(PI) * np.float32(lat_a) / np.float32(180)) - # math.sin(np.float32(PI) * np.float32(lat_b) / np.float32(180)) * # math.cos(np.float32(PI) * np.float32(a1) / np.float32(180))) / \ # (math.sin(np.float32(PI) * np.float32(a1) / np.float32(180)) * # math.cos(np.float32(PI) * np.float32(lat_b) / np.float32(180))) #y = math.acos(a4) if -1 <= a4 <= 1: y = math.acos(a4) # y = math.acos(np.float32(a4)) else: if a4 < -1: a4 = -1 if a4 > 1: a4 = 1 y = math.acos(a4) # y = math.acos(np.float32(a4)) bearing_b = (180.0 / PI) * y # bearing_b = np.float32(np.float32(180.0) / np.float32(PI)) * np.float32(y) if math.sin(PI * (lon_b - lon_a) / 180) < 0: #< isntead of >= but works bearing_a = 360 - bearing_a else: bearing_b = 360 - bearing_b # if math.sin(np.float32(PI) * (np.float32(lon_b) - np.float32(lon_a)) / np.float32(180)) < 0: # bearing_a = np.float32(360) - np.float32(bearing_a) # else: # bearing_b = np.float32(360) - np.float32(bearing_b) # Calculate the great-circle distance #distance doesnt even get used so it doesnt matter #distance = geodesic((lat_a, lon_a), (lat_b, lon_b)).miles distance = calculate_distance(lat_a, lat_b, lon_a, lon_b, unit='M') #print(distance) distance = int((distance * 1000) + .5) distance = distance / 1000 distance = str(distance) bearing_a = int((bearing_a * 10) + .5) bearing_a = bearing_a / 10 bearing_a = str(bearing_a) bearing_b = int((bearing_b * 10) + .5) bearing_b = bearing_b / 10 bearing_b = str(bearing_b) TxBearing_deg = bearing_a RxBearing_deg = bearing_b # print(f"TxBearing_deg: {TxBearing_deg} RxBearing_deg: {RxBearing_deg}") return TxBearing_deg, RxBearing_deg #9100 #should run fine with the other function, delete this later # def calculate_distance(lat_a, lat_b, lon_a, lon_b, unit='M'): # PI = 3.1415926 # # High accuracy formula # z = (math.sin((PI * (lat_a - lat_b) / 180.) / 2)) ** 2 # z += math.cos(PI * lat_a / 180) * math.cos(PI * lat_b / 180) * (math.sin((PI * (lon_a - lon_b) / 180) / 2)) ** 2 # z = math.sqrt(z) # # Calculate arcsin # z_short = 2 * (180 / PI) * math.asin(z) # # Distance in kilometers and miles # z_km = 111.1 * z_short # z_miles = 69.06 * z_short # if unit == "M": # distance = z_miles # elif unit == "K": # distance = z_km # else: # raise ValueError("Unit must be 'M' for miles or 'K' for kilometers") # return distance #9200/9400, arccosine function example #test if this works the same for both functions #maybe i try out both python functions that were translated, then test if arccos works #result = math.acos(x) #(Interference) def custom_copysign(x): if x >= 0: return 1 else: return -1 #9400, arccosine function example (Interference) def ARCCOS(X): PI = 3.1415926 if abs(X) > 1: XMOD = X / abs(X) else: XMOD = X if XMOD == 1: ACOS = 0 elif XMOD == -1: ACOS = PI elif abs(XMOD) < 1: ACOS = PI / 2 - math.atan(XMOD / math.sqrt(1 - XMOD * XMOD)) return ACOS #9300, arcsin function example #result = math.asin(1, 3.1415926) #9500, calculate antenna tilt (Interference) def calculate_antenna_tilt(bldg_meanelev_ft, bldg_maxheight_ft, a_rcvr_grd_elev_ft, a_rcvr_ant_raat_ft, bldg_a_rcvr_distance_km, a_rcvr_path_id, old_building): # Constants PI = 3.1415926 K = 4.0 / 3.0 # Convert heights to meters and kilometers to miles Htx = (float(bldg_meanelev_ft) + float(bldg_maxheight_ft) - 5.0) * 0.3048 # meters Hrx = (float(a_rcvr_grd_elev_ft)) + (float(a_rcvr_ant_raat_ft)) # meters PathDistance = float(bldg_a_rcvr_distance_km) # kilometers # Calculate interference distance in miles if conditions match if a_rcvr_path_id != old_building: InterferenceDistanceMiles = str(0.621371 * PathDistance).strip() else: InterferenceDistanceMiles = "" # Calculate antenna tilt for transmitter H = Hrx - Htx Phi1 = 0.001 * H / PathDistance Phi2 = (math.sqrt((PathDistance * PathDistance) + (0.000001 * H * H))) / (12735. * K) Phi1 = math.atan(Phi1) # arctangent Phi2 = math.atan(Phi2 / math.sqrt(1. - (Phi2 * Phi2))) # arcsine AntTiltTx = (Phi1 - Phi2) * (180. / PI) TxTilt_deg = str(AntTiltTx).strip() # Update interference tilt if conditions match if a_rcvr_path_id != old_building: InterferenceTiltDeg = TxTilt_deg else: InterferenceTiltDeg = "" # Calculate antenna tilt for receiver H = Htx - Hrx Phi1 = 0.001 * H / PathDistance Phi2 = (math.sqrt((PathDistance * PathDistance) + (0.000001 * H * H))) / (12735. * K) Phi1 = math.atan(Phi1) # arctangent Phi2 = math.atan(Phi2 / math.sqrt(1. - (Phi2 * Phi2))) # arcsine AntTiltRx = (Phi1 - Phi2) * (180. / PI) RxTilt_deg = str(AntTiltRx).strip() return f"{float(TxTilt_deg):.6f}", f"{float(RxTilt_deg):.6f}", f"{float(InterferenceDistanceMiles):.6f}", f"{float(InterferenceTiltDeg):.6f}" #9600, find desired mw path information (Interference) def find_mw_path_information(gpd_lookup, rx_site, tx_site): rec = gpd_lookup.get((tx_site, rx_site)) if rec is None: print(f"Subroutine 9600 failed to find matching mw path <{tx_site}> <{rx_site}>.") print("Program terminated.") return None rcv_gain = rec["gain"] rcv_tilt = rec["tilt"] rcv_bearing = rec["bearing"] rcv_ant_model = rec["ant_model"] fr_frequency_assigned_MHzMW = rec["frequency"] rxan_gain_dBiMW = rec["gain"] channels = rec["channels"] bandwidth = rec["bandwidth"] tx_location_name = rec["tx_location_name"] rx_location_name = rec["rx_location_name"] site_name_and_direction = f"{tx_location_name} TO {rx_location_name}" return rcv_gain, rcv_tilt, rcv_bearing, rcv_ant_model, fr_frequency_assigned_MHzMW, rxan_gain_dBiMW, channels, bandwidth, site_name_and_direction #9700, determine antenna discrimination #(Interference) #takes the value from the spreedsheet, then does some calculations #to see which one gets picked, arc coseine function used for total_angle def determine_antenna_discrimination(drive, folder, rcv_ant_model, total_angle): total_angle = round(total_angle, 6) antenna_file = ANTENNAS_DIR + '\\' + rcv_ant_model + '-WC.csv' ant_disc = "0" old_delta = 999.0 try: with open(antenna_file, mode='r') as file: csv_reader = csv.reader(file) for row in csv_reader: row = [field.strip() for field in row if field.strip() != ''] if len(row) < 2: continue an_angle, a_value = row[0], row[1] if not an_angle or not a_value: continue new_delta = abs(float(an_angle) - total_angle) if new_delta > old_delta: break if new_delta < old_delta: old_delta = new_delta ant_disc = a_value #catches if file not found except FileNotFoundError: print(f"Error: The file '{antenna_file}' was not found.") input("Press Enter to exit...") return None #catches any other error except Exception as e: print(f"An error occurred: {e}") import traceback print(traceback.format_exc()) raise return ant_disc #9800, knife edge loss (dB) (Interference/Terrain) #needs to be tested def knife_edge_loss(int_path_dist_miles, bldg_a_rcvr_distance_km, fr_frequency_assigned_mhz, bldg_maxheight_ft, bldg_meanelev_ft, ant_ht_ft, int_ht_ft): pi = 3.1415926 # Path distance in miles d = int_path_dist_miles # Distance from one end to point on path in miles d1 = 0.621371 * float(bldg_a_rcvr_distance_km) if d1 < 0.00947: d1 = 0.00947 # Distance from other end to point on path in miles d2 = d - d1 if d2 < 0.00947: d2 = 0.00947 # Operating frequency in GHz freq = float(fr_frequency_assigned_mhz) / 1000.0 # First Fresnel zone in feet f1 = 72.1 * math.sqrt(d1 * d2 / (freq * d)) # Building obstruction clearance in feet h1 = float(bldg_maxheight_ft) + float(bldg_meanelev_ft) # Propagation height factor h2 = (d1 * d2) / (1.5 * (4.0 / 3.0)) # Current distance from victim antenna in miles current_distance_from_victim_antenna_miles = 0.621371 * float(bldg_a_rcvr_distance_km) # Path height in feet path_height_ft = ant_ht_ft + (((int_ht_ft - ant_ht_ft) / int_path_dist_miles) * current_distance_from_victim_antenna_miles) # Calculate H #print(f"h1: {h1}, h2: {h2}, path_height_ft: {path_height_ft}") h = path_height_ft - (float(h1) + h2) # Calculate X x = h / f1 # Calculate knife edge loss (GKE) if x > 0.6: gke = 0.0 elif x < -600.0: gke = -73.0 elif x <= 0.6 and x >= -1.0: a = -6.013050244594124 b = 12.56373978136271 c = -4.638348100303959 d = -14.27304778829325 e = 4.624584658556855 f = 26.08257628627929 g = 14.16222106310102 gke = a + b * x + c * x**2 + d * x**3 + e * x**4 + f * x**5 + g * x**6 elif x <= -1.0 and x >= -600.0: a = -10.13843752145488 b = -0.2211687451291319 c = 10.31071865966655 d = 0.0009104403873283171 e = -0.0762342416178572 gke = (a + c * x + e * x**2) / (1.0 + b * x + d * x**2) else: gke = 0.0 return gke, d1, d2, d, path_height_ft, h1, h2, h, x, f1 #9900, determine effect of interference on fade margin (Interference/Terrain) def determine_fade_margin_effect(interference_dbm): rcv_noise = -94.0 interference_10db = interference_dbm - 10.0 interference_20db = interference_dbm - 20.0 interference_30db = interference_dbm - 30.0 log10 = math.log(10) rfn_10db = (10.0 * math.log10((10.0 ** (rcv_noise / 10.0)) + (10.0 ** (interference_10db / 10.0))) ) - rcv_noise rfn_20db = (10.0 * math.log10((10.0 ** (rcv_noise / 10.0)) + (10.0 ** (interference_20db / 10.0))) ) - rcv_noise rfn_30db = (10.0 * math.log10((10.0 ** (rcv_noise / 10.0)) + (10.0 ** (interference_30db / 10.0))) ) - rcv_noise # if rfn_10db < 0.1: # rfn_10db = 0.1 # if rfn_20db < 0.1: # rfn_20db = 0.1 # if rfn_30db < 0.1: # rfn_30db = 0.1 # Apply the one significant decimal conversion # rfn_10db = convert_to_one_significant_decimal(rfn_10db) # rfn_20db = convert_to_one_significant_decimal(rfn_20db) # rfn_30db = convert_to_one_significant_decimal(rfn_30db) rfn_10db = round(rfn_10db, 1) rfn_20db = round(rfn_20db, 1) rfn_30db = round(rfn_30db, 1) return rfn_10db, rfn_20db, rfn_30db #9950, convert to one significant decimal (Interference) def convert_to_one_significant_decimal(variable): variable = int(0.5 + (variable * 10.0)) variable = variable / 10.0 return variable #Terrain Calculation def run_terrain_calculation(terrain_path_input, terrain_path_output, preloaded_data=None): K = 4.0 / 3.0 FreqGHz = (6.175 + 6.7) / 2.0 d_miles_factor = 0.621371 try: all_rows = preloaded_data if preloaded_data is not None else read_csv_safe(terrain_path_input) data_rows = all_rows[1:] if all_rows and len(all_rows) > 1 else [] with open(terrain_path_output, mode='w', newline='') as output_file: csv_writer = csv.writer(output_file) csv_writer.writerow([ "Path", "Int Building", "Near Rx Site", "Far Tx Site", "Path Clearance Description", "Obstruction Loss (dB)" ]) pbar = tqdm(total=len(data_rows), desc="Processing Rows", bar_format="{desc}: {percentage:3.0f}%|{bar}|") OldPath = None OldBuilding = None OldRxSite = None OldTxSite = None StartFlag = False Description = "Free Space Propagation" ObstructionLoss = 0.0 DMiles = 0.0 def flush(): if StartFlag: csv_writer.writerow([ OldPath, OldBuilding, OldRxSite, OldTxSite, Description, str(round(ObstructionLoss, 1)) ]) for row in data_rows: pbar.update(1) try: fields = [f.strip() for f in row] if len(fields) < 17: continue (unique_path_id, bldg_number, bldg_long, bldg_lat, bldg_max_elev_ft, bldg_maxheight_ft, bldg_address, a_rcvr_call_sign, a_rcvr_long, a_rcvr_lat, a_rcvr_grd_elev_me, a_rcvr_raat_me, incr_long, incr_lat, incr_elev_me, incr_distance_km, path_distance_km) = fields[:17] # Skip points at or beyond path end if float(incr_distance_km) >= float(path_distance_km): continue building, rx_site, tx_site = parse_sites_simplified(unique_path_id) if unique_path_id != OldPath: # Write previous building flush() # Reset for new building Description = "Free Space Propagation" ObstructionLoss = 0.0 DMiles = d_miles_factor * float(path_distance_km) if float(bldg_maxheight_ft) <= 0: bldg_maxheight_ft = "15" elif float(bldg_maxheight_ft) < 10: bldg_maxheight_ft = "10" ObstructionLoss, _, _, Description = calculate_obstruction_loss( a_rcvr_grd_elev_me, a_rcvr_raat_me, bldg_maxheight_ft, bldg_max_elev_ft, incr_distance_km, incr_elev_me, DMiles, FreqGHz, K, ObstructionLoss, Description) OldPath = unique_path_id OldBuilding = building OldRxSite = rx_site OldTxSite = tx_site StartFlag = True except Exception as row_err: print(f" ERROR in terrain row of {terrain_path_input}: {row_err}") import traceback; traceback.print_exc() raise # Write last building flush() if not StartFlag: if len(data_rows) == 0: msg = f"WARNING: Terrain file is empty (header only): {terrain_path_input}" else: msg = f"WARNING: No terrain rows processed for {terrain_path_input} - path is too short for any sample point to fall within range" print(msg) with open(WARNING_LOG, 'a') as wf: wf.write(msg + '\n') pbar.close() except FileNotFoundError: print(f"Error: File not found: {terrain_path_input}") raise except Exception as e: print(f"An error occurred in terrain calculation: {e}") import traceback; traceback.print_exc() raise #(Interference/Terrain) def calculate_distance(lat_a, lat_b, lon_a, lon_b, unit='M'): PI = 3.1415926 # High accuracy formula z = (math.sin((PI * (lat_a - lat_b) / 180.) / 2)) ** 2 z += math.cos(PI * lat_a / 180) * math.cos(PI * lat_b / 180) * (math.sin((PI * (lon_a - lon_b) / 180) / 2)) ** 2 z = math.sqrt(z) # Calculate arcsin z_short = 2 * (180 / PI) * math.asin(z) # Distance in kilometers and miles z_km = 111.1 * z_short z_miles = 69.06 * z_short if unit == "M": distance = z_miles elif unit == "K": distance = z_km else: raise ValueError("Unit must be 'M' for miles or 'K' for kilometers") return distance #9500, rounded earth loss (Terrain) def rounded_earth_loss(H_F1): # Check for H_F1 > 0.6 if H_F1 > 0.6: GRE = 0.0 # Check for 0 <= H_F1 <= 0.6 elif 0.0 <= H_F1 <= 0.6: GRE = 10.0 - (20.0 * H_F1) + (5.556 * H_F1 * H_F1) # For H_F1 < 0 else: GRE = 10.0 - (20.0 * H_F1) return GRE #(Terrain) def calculate_obstruction_loss(a_rcvr_grd_elev_me, a_rcvr_raat_me, bldg_maxheight_ft, bldg_max_elev_ft, incr_distance_km, incr_elev_me, DMiles, FreqGHz, K, ObstructionLoss, Description): # Convert victim antenna height AMSL in feet AntHtFt = (float(a_rcvr_grd_elev_me) * 3.28084) + (float(a_rcvr_raat_me) * 3.28084) # Convert interference height AMSL in feet IntHtFt = float(bldg_maxheight_ft) + float(bldg_max_elev_ft) - 5.0 # Convert victim antenna height AGL in feet AntHtFtAGL = float(a_rcvr_raat_me) * 3.28084 # Convert interference height AGL in feet IntHtFtAGL = float(bldg_max_elev_ft) - 5.0 # Convert distance to miles D1Miles = float(incr_distance_km) * 0.621371 D2Miles = DMiles - D1Miles # Apply minimum distance thresholds D1Miles = max(D1Miles, 0.01) D2Miles = max(D2Miles, 0.01) # Calculate additional height AddlHt = (D1Miles * D2Miles) / (1.5 * K) # Calculate terrain height in feet TerrainHeight = (float(incr_elev_me) * 3.28084) + AddlHt # Calculate terrain distance from receiver antenna in miles TerrainDistance = float(incr_distance_km) * 0.621371 # Calculate path height PathHeight = AntHtFt + (((IntHtFt - AntHtFt) / DMiles) * TerrainDistance) # Calculate height difference Difference = PathHeight - TerrainHeight # First Fresnel zone in feet F1 = 72.1 * math.sqrt(D1Miles * D2Miles / (FreqGHz * DMiles)) # Calculate H/F1 ratio H_F1 = Difference / F1 # Call the rounded earth loss function (GOSUB 9500 in the original) GRE = rounded_earth_loss(H_F1) # Initialize obstruction description and check for obstruction if GRE > ObstructionLoss: ObstructionLoss = GRE Description = "Mildly Obstructed Path" if GRE > 10.0: Description = "Noticeably Obstructed Path" if GRE > 20.0: Description = "Strongly Obstructed Path" # Convert values to strings for output # DMiles_str = f"{DMiles:.2f}" # D1Miles_str = f"{D1Miles:.2f}" # D2Miles_str = f"{D2Miles:.2f}" # TerrainDistance_str = f"{TerrainDistance:.2f}" # PathHeight_str = f"{PathHeight:.2f}" # TerrainHeight_str = f"{TerrainHeight:.2f}" # F1_str = f"{F1:.2f}" # H_F1_str = f"{H_F1:.2f}" # GRE_str = f"{GRE:.2f}" return ObstructionLoss, IntHtFtAGL ,AntHtFtAGL, Description #9600, create clear and obstructed path files (terrain) def create_clear_obstructed_path_files(ObstructionLoss, Building, bldg_lat, bldg_long, RxSite, a_rcvr_lat, a_rcvr_long, IntHtFtAGL, AntHtFtAGL, ClearCounter, ObstructedCounter): PRINTFILE = "" # If no obstruction loss (clear path) if ObstructionLoss == 0.0: ClearCounter += 1 # Construct the line for the clear path PRINTFILE = f"{ClearCounter},{Building},{bldg_lat},{bldg_long},{RxSite},{a_rcvr_lat},{a_rcvr_long},{IntHtFtAGL},{AntHtFtAGL}" # If obstruction loss is greater than 20 (obstructed path) if ObstructionLoss > 20.0: ObstructedCounter += 1 # Construct the line for the obstructed path PRINTFILE = f"{ObstructedCounter},{Building},{bldg_lat},{bldg_long},{RxSite},{a_rcvr_lat},{a_rcvr_long},{IntHtFtAGL},{AntHtFtAGL}" return PRINTFILE, ClearCounter, ObstructedCounter #start of program if __name__ == "__main__": run_program()