#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:\ATT\ATT_Rosetta_populated.csv' GPD_FILE = r'D:\ATT\ATT_GPD_ALL.xlsx' ANTENNAS_DIR = r'D:\ATT\AntennasD2' 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.") # Build NSMA gain lookup from antenna WC files # The gain comes from the NSMA header MDGAIN field # We store it keyed by antenna model name (same as WC filename without -WC.csv) print("Loading antenna gains from NSMA files...") nsma_gain_lookup = {} import os, glob for nsma_dir in NSMA_DIRS: try: for fpath in glob.glob(os.path.join(nsma_dir, '*.dat')) + \ glob.glob(os.path.join(nsma_dir, '*.DAT')): model = None gain = None try: with open(fpath, encoding='latin-1') as f: for line in f: line = line.strip() if line.startswith('MODNUM:'): model = line.split(',', 1)[1].strip() elif line.startswith('MDGAIN:'): try: gain = float(line.split(',', 1)[1].strip()) except ValueError: pass if model and gain is not None: if model not in nsma_gain_lookup: nsma_gain_lookup[model] = gain break except Exception: pass except Exception: pass print(f"NSMA gains loaded: {len(nsma_gain_lookup)} antenna models.") 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() print("Restart Program and Try Again") input("Press Enter to exit...") exit() # 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 (Interference) 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: pass # Loading Progress Bar # Use preloaded data if available, otherwise read from disk if preloaded_data is not None: all_rows = preloaded_data else: all_rows = read_csv_safe(interference_path_input) total_lines = len(all_rows) - 1 # subtract header with open(interference_path_output, mode='w', newline='') as output_file: csv_writer = csv.writer(output_file) # Write the header row to the 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" ]) # Initialize variables old_building = "xxx" start_flag = False input_line_count = 0 Spinner = 0 # Skip the header line data_rows = iter(all_rows[1:]) pbar = tqdm(total=total_lines, desc="Processing Rows", bar_format="{desc}: {percentage:3.0f}%|{bar}|") # Loop through each row in the input file for row in data_rows: pbar.update(1) #print("") # Increment the input line count input_line_count += 1 input_line_count_str = str(input_line_count).strip() # Strip leading and trailing whitespace from each field 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 = [field.strip() for field in row] # Convert the numeric fields to the desired format bldg_long = f"{float(bldg_long):.6f}" bldg_lat = f"{float(bldg_lat):.6f}" bldg_maxheight_ft = f"{float(bldg_maxheight_ft):.4f}" bldg_meanelev_ft = f"{float(bldg_meanelev_ft):.4f}" 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}" # Parse the building, receive site, and transmit site from the path ID building, rx_site, tx_site = parse_sites_simplified(a_rcvr_path_id) # Check if the building number is the same as the previous building if float(bldg_maxheight_ft) == 0.0: bldg_maxheight_ft = "10" if float(bldg_maxheight_ft) < 10.0: bldg_maxheight_ft = "10" ClearanceFlag = 1 if ClearanceFlag ==1: # Check if the building number is the same as the previous building if building != old_building: # Check if this is the first building evaluation if start_flag == 1: #Spinner = Spinner + 1 # Print the completion line every 10 buildings #if Spinner == 1000: # print("Inter: " + input_line_count_str + " " + bldg_number + " " + folder + " " + subfolder + " " + fileAIntVER) #file_extension + " " + the_notes) # Spinner = 0 #determine fade margin effect rfn10dB, rfn20dB, rfn30dB = determine_fade_margin_effect(interference_dbm) more = ",".join([more_custom_format(value) for value in [rfn10dB, rfn20dB, rfn30dB]]) more_parts = more.split(',') # Write the output row to the output file print_file_io = StringIO(print_file) print_file_parts = next(csv.reader(print_file_io)) csv_writer.writerow([ old_building, old_rx_site, old_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}" ]) # Initialize the print file string print_file = f'{bldg_long},{bldg_lat},{int(float(bldg_maxheight_ft))},"{bldg_desc}","{bldg_type}"' int_path_dist_miles = float(bldg_a_rcvr_distance_km) * 0.621371 # interference path distance in miles from interference source to victim antenna # Calculate victim antenna height AMSL in feet ant_ht_ft = (3.28084 * float(a_rcvr_grd_elev_ft)) + (3.28084 * float(a_rcvr_ant_raat_ft)) # Calculate interference heights AMSL in feet int_ht_ft = float(bldg_maxheight_ft) + float(bldg_meanelev_ft) - 5.0 # Convert interference heights to string and trim whitespace int_ht_ft_str = str(int_ht_ft).strip() Int_Building = bldg_number # path_height = ant_ht_ft + ((int_ht_ft - ant_ht_ft) / int_path_dist_miles) * current_distance_from_victim_antenna #2 determine interference angle and tilt relative to victim antenna #------ #GOSUB 9000 #determine bearings #calculate_distance_and_bearing(lat_a, lon_a, lat_b, lon_b) tx_bearing_deg, rx_bearing_deg = calculate_distance_and_bearing(bldg_lat, bldg_long, a_rcvr_lat, a_rcvr_long) #Interference: TxBearing_deg$ #Victim: RxBearing_deg$ #GOSUB 9500 #determine tilts 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, old_building) #Interference: TxTilt_deg$ #Victim: RxTilt_deg$ #GOSUB 9600 #find desired mw path information rcv_gain, rcv_tilt, rcv_bearing, rcv_ant_model, fr_frequency_assigned_MHzMW, rxan_gain_dBiMW, channels, bandwidth, site_name_and_direction = find_mw_path_information(gpd_lookup, rx_site, tx_site) #4 determine antenna discrimination #maybe rx_bearing_deg math got messed up # print("RCV BEARING: ", rcv_bearing) # print("RX BEARING: ", rx_bearing_deg) # print((abs(float(rcv_bearing)))) # print(float(rx_bearing_deg)) hor_off_axis = abs(float(rcv_bearing) - float(rx_bearing_deg)) if hor_off_axis > 180: hor_off_axis = 360 - hor_off_axis hor_off_axis = round(hor_off_axis, 1) hor_off_axis_str = str(hor_off_axis).strip() if building != old_building: int_hor_off_axis = hor_off_axis 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) # Calculate angles angle1 = math.cos(delta_d) * math.cos(alpha_d) * math.cos(phi_i) * math.cos(alpha_i) angle2 = math.cos(delta_d) * math.sin(alpha_d) * math.cos(phi_i) * math.sin(alpha_i) angle3 = math.sin(delta_d) * math.sin(phi_i) x = angle1 + angle2 + angle3 ARCCOS = math.acos(x) #9400, arccosine function example total_angle = math.degrees(ARCCOS) total_angle_str = str(total_angle).strip() #determine antenna discrimination, 9700 ant_disc = determine_antenna_discrimination(drive, folder, rcv_ant_model, total_angle) # Use gain from NSMA lookup if available, otherwise fall back to GPD value 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)) fsl_db = round(fsl_db, 4) fsl_db_str = str(fsl_db).strip() hor_off_axis = abs(float(rcv_bearing) - float(rx_bearing_deg)) if hor_off_axis > 180: hor_off_axis = 360 - hor_off_axis hor_off_axis = round(hor_off_axis, 1) hor_off_axis_str = str(hor_off_axis).strip() interference_dbm = 19.7 - fsl_db + effective_gain + float(ant_disc) interference_dbm = round(interference_dbm, 5) interference_dbm_str = str(interference_dbm).strip() # Check if the building number is the same as the previous building if building == old_building: #9800, knife edge loss (dB) GKE, D1, D2, D, PathHeightFt, H1, H2, H, X, F1 = 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) # Update interference dBm interference_dbm += GKE interference_dbm = round(interference_dbm, 4) interference_dbm_str = str(interference_dbm).strip() # Convert the numeric fields to the desired format D1_str = str(D1).strip() D2_str = str(D2).strip() D_str = str(D).strip() PathHeightFt_str = str(PathHeightFt).strip() H1_str = str(H1).strip() H2_str = str(H2).strip() H_str = str(H).strip() X_str = str(X).strip() GKE_str = str(GKE).strip() F1_str = str(F1).strip() old_building = building old_rx_site = rx_site old_tx_site = tx_site start_flag = 1 #not the first building evalualtion #9900 determine fade margin effect rfn10dB, rfn20dB, rfn30dB = determine_fade_margin_effect(interference_dbm) more = ",".join([more_custom_format(value) for value in [rfn10dB, rfn20dB, rfn30dB]]) more_parts = more.split(',') print_file_io = StringIO(print_file) print_file_parts = next(csv.reader(print_file_io)) #gonna have to getting the rounding right later csv_writer.writerow([building, rx_site, tx_site, site_name_and_direction, f"{float(interference_dbm):.4f}", *more_parts, *print_file_parts, channels, bandwidth, f"{float(InterferenceDistanceMiles):.5f}", int_hor_off_axis, f"{float(InterferenceTiltDeg):.6f}"]) pbar.close() #catches if file not found except FileNotFoundError: print(f"Error: One of the files '{interference_path_input}' or '{interference_path_output}' was not found.") print("Restart Program and Try Again") input("Press Enter to exit...") exit() #catches any other error except Exception as e: print(f"An error occurred: {e}") traceback.print_exc() print("Restart Program and Try Again") input("Press Enter to exit...") exit() #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()) print("Restart Program and Try Again") input("Press Enter to exit...") exit() 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) 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 PI = 3.1415926 ClearCounter = 0 ObstructedCounter = 0 d_miles_factor = 0.621371 try: pass # Loading Progress Bar # Use preloaded data if available, otherwise read from disk if preloaded_data is not None: all_rows = preloaded_data else: all_rows = read_csv_safe(terrain_path_input) total_lines = len(all_rows) - 1 # subtract header with open(terrain_path_output, mode='w', newline='') as output_file: csv_writer = csv.writer(output_file) # Write the header row to the output file csv_writer.writerow([ "Path", "Int Building", "Near Rx Site", "Far Tx Site", "Path Clearance Description", "Obstruction Loss (dB)" ]) header = all_rows[0] data_rows = iter(all_rows[1:]) pbar = tqdm(total=total_lines, desc="Processing Rows", bar_format="{desc}: {percentage:3.0f}%|{bar}|") OldPath = "xxx" StartFlag = 0 Description = "Free Space Propagation" ObstructionLoss = 0.0 GRE = 0.0 #InputLineCount = 0 #Spinner = 0 counter = 0 #add in progress bar for row in data_rows: pbar.update(1) #Unpack the input line into variables 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 = [field.strip() for field in row] # # Strip leading and trailing whitespace from each field # 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 = [field.strip() for field in row] #unique_path_id, bldg_number, bldg_long, bldg_lat, bldg_max_elev_ft, bldg_maxheight_ft = [field.strip() for field in row] # Skip the line if the incremental distance exceeds or equals the path distance if float(incr_distance_km) >= float(path_distance_km): continue #InputLineCount += 1 # unique_path_id = unique_path_id.strip() # bldg_number = bldg_number.strip() # InputLineCount_str = str(InputLineCount).strip() # Spinner += 1 # if Spinner == 10000: # #see how long it takes with 1000 or whatever the planned number is # #print(f"Terr: {operation} {InputLineCount_str} {bldg_number} {folder} {subfolder} {file_extension} {notes}") # print("did 10,000") # Spinner = 0 building, rx_site, tx_site = parse_sites_simplified(unique_path_id) if unique_path_id != OldPath: # If not the first evaluation, write the previous evaluation's results if StartFlag == 1: ObstructionLoss = round(ObstructionLoss, 1) ObstructionLoss_str = str(ObstructionLoss).strip() # Write the previous results to the output file # csv_writer.writerow([OldPath, OldBuilding, OldRxSite, OldTxSite, Description, ObstructionLoss_str]) # # Create clear and obstructed path files # create_clear_obstructed_path_files(ObstructionLoss, OldBuilding, bldg_lat, bldg_long, OldRxSite, a_rcvr_lat, a_rcvr_long, IntHtFtAGL, AntHtFtAGL, ClearCounter, ObstructedCounter) # buffer.append([OldPath, OldBuilding, OldRxSite, OldTxSite, Description, ObstructionLoss_str]) # batch_count += 1 csv_writer.writerow([OldPath, OldBuilding, OldRxSite, OldTxSite, Description, ObstructionLoss_str]) # Write the batch when buffer is full # Start a new evaluation Description = "Free Space Propagation" ObstructionLoss = 0.0 DMiles = d_miles_factor * float(path_distance_km) # Convert path distance to miles # Ensure building height is valid, default to 15 if zero, or at least 10 if float(bldg_maxheight_ft) == 0: bldg_maxheight_ft = "15" elif float(bldg_maxheight_ft) < 10: bldg_maxheight_ft = "10" # Calculate obstruction loss (calls to the relevant subroutines in the original code) ObstructionLoss, IntHtFtAGL ,AntHtFtAGL, 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) elif unique_path_id == OldPath: ObstructionLoss, IntHtFtAGL ,AntHtFtAGL, 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) # Update OldPath and related variables for the next iteration OldPath = unique_path_id OldBuilding = building OldRxSite = rx_site OldTxSite = tx_site # Set StartFlag to indicate the first building evaluation is done StartFlag = 1 #lines 227-237, ends at END, similar to interference code #add some more here # Round ObstructionLoss to 1 decimal place (equivalent to INT((ObstructionLoss * 10.) + 0.5)) ObstructionLoss = round(ObstructionLoss, 1) # Convert ObstructionLoss to a string and remove leading and trailing spaces ObstructionLoss_str = str(ObstructionLoss).strip() # Write the results to the output file (OldPath, OldBuilding, OldRxSite, OldTxSite, Description, ObstructionLoss) #output_file.write(f"{OldPath},{OldBuilding},{OldRxSite},{OldTxSite},{Description},{ObstructionLoss_str}") csv_writer.writerow([OldPath, OldBuilding, OldRxSite, OldTxSite, Description, ObstructionLoss_str]) # # Call the subroutine to create clear and obstructed path files (equivalent to GOSUB 9600) # create_clear_obstructed_path_files(ObstructionLoss, building, bldg_lat, bldg_long, rx_site, a_rcvr_lat, a_rcvr_long, IntHtFtAGL, AntHtFtAGL, ClearCounter, ObstructedCounter) # if buffer: # csv_writer.writerows(buffer) # Close the output file (if not using 'with' statement for file handling) #output_file.close() pbar.close() except FileNotFoundError: print(f"Error: One of the files '{terrain_path_input}' or '{terrain_path_output}' was not found.") print("Restart Program and Try Again") input("Press Enter to exit...") exit() #catches any other error except Exception as e: print(f"An error occurred: {e}") traceback.print_exc() print("Restart Program and Try Again") input("Press Enter to exit...") exit() #(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()