#ifndef _Additional_h_
#define _Additional_h_
#include <cmath>
#include <filesystem>
#include <iostream>
#include <map>

#include "Structures.h"
#include "Antenna.h"
#include "Antenna_correction.h"


const double RadToDeg = 180.0/3.14159265358979323846;
const double DegToRad = 3.14159265358979323846/180.0;

struct Coordinate {
	double longitude;
	double latitude;
};
static std::map<std::string, std::vector<std::vector<std::pair<double, double>>>> borderCacheAll;

/**
 * Checks if a file exists at the given path
 * 
 * @param filePath The path to check for file existence
 * @return true if file exists, false otherwise
 */
bool fileExists(const std::string& filePath) {
	return std::filesystem::exists(filePath);
}

/**
 * Replaces all occurrences of a substring in a string
 * 
 * @param str The string to modify
 * @param from The substring to replace
 * @param to The replacement substring
 */
void replaceAll(std::string& str, const std::string& from, const std::string& to) {
	size_t start_pos = 0;
	while((start_pos = str.find(from, start_pos)) != std::string::npos) {
		str.replace(start_pos, from.length(), to);
		start_pos += to.length();
	}
}

/**
 * Checks if two numbers differ by more than 10%
 * 
 * @param a First number to compare
 * @param b Second number to compare
 * @return true if difference exceeds 10%, false otherwise
 */
bool isDiff10Percent(double a, double b) { return std::abs(a - b) > 0.1 * std::max(std::abs(a), std::abs(b)); }

/**
 * Trims whitespace characters from both ends of a string
 * 
 * @param s The string to trim
 * @return Trimmed string
 */
std::string Trimmer(const std::string &s) {
	std::string ss;
	size_t start = s.find_first_not_of(" \n\r\t\f\v");
	ss = (start == std::string::npos) ? "" : s.substr(start);
	size_t end = ss.find_last_not_of(" \n\r\t\f\v");
	return (end == std::string::npos) ? "" : ss.substr(0, end + 1);
}

/**
 * Calculates great-circle distance between two geographic coordinates using Haversine formula
 * 
 * @param Longitude1 Longitude of first point in degrees
 * @param Latitude1 Latitude of first point in degrees
 * @param Longitude2 Longitude of second point in degrees
 * @param Latitude2 Latitude of second point in degrees
 * @return Distance in kilometers
 */
double CalcDistanceM(double Longitude1, double Latitude1, double Longitude2, double Latitude2) {
	const double R = 6371.29; // Earth radius in kilometers

	double delta_lat = Latitude2 - Latitude1;
	double delta_lon = Longitude2 - Longitude1;

	// Haversine formula components
	double a = pow(sin(delta_lat / 2.0 * DegToRad), 2) +
			cos(Latitude1 * DegToRad) * cos(Latitude2 * DegToRad) * pow(sin(delta_lon / 2.0 * DegToRad), 2);

	// Calculate central angle and distance
	return 2.0 * R * atan2(sqrt(a), sqrt(1.0 - a));
}

/**
 * Calculates initial bearing (azimuth) from first point to second point
 * 
 * @param Longitude1 Longitude of starting point in degrees
 * @param Latitude1 Latitude of starting point in degrees
 * @param Longitude2 Longitude of destination point in degrees
 * @param Latitude2 Latitude of destination point in degrees
 * @return Bearing angle in degrees from north (0-360)
 */
double CalcDirectionM(double Longitude1, double Latitude1, double Longitude2, double Latitude2) {
	double delta_lon = Longitude2 - Longitude1;

	// Compute numerator (y) and denominator (x) for atan2
	double y = sin(delta_lon * DegToRad) * cos(Latitude2 * DegToRad);
	double x = cos(Latitude1 * DegToRad) * sin(Latitude2 * DegToRad) - sin(Latitude1 * DegToRad) * cos(Latitude2 * DegToRad) * cos(delta_lon * DegToRad);

	// Calculate azimuth in degrees and adjust to [0, 360)
	double Azi = std::atan2(y, x) * RadToDeg;
	if (Azi < 0.0) Azi += 360.0;
	return Azi;
}

/**
 * Calculates new coordinates from starting point given direction and distance
 * 
 * @param longitude Starting longitude in degrees
 * @param latitude Starting latitude in degrees
 * @param direction Travel direction in degrees
 * @param distance Travel distance in kilometers
 * @param nLong Output parameter for new longitude
 * @param nLat Output parameter for new latitude
 */
void New_coordinatesA(double longitude, double latitude, double direction, double distance, double &nLong, double &nLat) {
	if (distance == 0.0) {
		nLong = longitude;
		nLat = latitude;
		return;
	}
	double dp = distance / 111.2;
	double t1 = cos(latitude * DegToRad) * sin(dp * DegToRad);
	double t2 = (sin(latitude * DegToRad) * cos(dp * DegToRad)) / t1;
	double t = t1 * (cos(direction * DegToRad) + t2);
	if (t > 1.0) t = 1.0;
	if (t < -1.0) t = -1.0;
	nLat = asin(t) * RadToDeg;

	if (direction == 0.0 || direction == 360.0 || direction == 180.0) {
		nLong = longitude;
		return;
	}

	t = (cos(dp * DegToRad) - sin(latitude * DegToRad) * sin(nLat * DegToRad)) / (cos(latitude * DegToRad) * cos(nLat * DegToRad));
	if (t > 1.0) t = 1.0;
	if (t < -1.0) t = -1.0;
	nLong = acos(t) * RadToDeg;
	if (direction < 180.0) nLong = longitude + nLong;
	if (direction > 180.0) nLong = longitude - nLong;
}

/**
 * Tests how many times a ray from a point in given direction crosses border polygons
 * 
 * @param dir Direction of the ray in degrees
 * @param longitude Starting longitude in degrees
 * @param latitude Starting latitude in degrees
 * @param serviceArea Length of the ray in kilometers
 * @param nCut Output parameter for number of border crossings
 * @param filename Path to border definition file
 */
void TestCutA(double dir, double longitude, double latitude,
			double serviceArea, int &nCut, std::string filename) {
	nCut = 0;
	double nLong, nLat;
	New_coordinatesA(longitude, latitude, dir, serviceArea, nLong, nLat);
	double AX = longitude;  
	double AY = latitude;   
	double BX = nLong;      
	double BY = nLat;       

	if (borderCacheAll.find(filename) == borderCacheAll.end()) {
		std::ifstream file(filename, std::ios::binary);
		if (!file.is_open()) {
			return;
		}

		union BorderRecord {
			char bytes[176];
			double coords[22];
		} buffer;

		std::vector<std::vector<std::pair<double, double>>> records;
		while (file.read(buffer.bytes, 176)) {
			std::vector<std::pair<double, double>> points;
			for (int i = 0; i < 22; i += 2) {
				double x = buffer.coords[i] * RadToDeg;   
				double y = buffer.coords[i+1] * RadToDeg;
				points.emplace_back(x, y);
			}
			records.push_back(points);
		}
		file.close();
		borderCacheAll[filename] = std::move(records);
	}

	const auto& records = borderCacheAll.at(filename);

	for (const auto& record : records) {
		for (size_t i = 1; i < record.size()-1; ++i) {
			double CX = record[i-1].first;  
			double CY = record[i-1].second;
			double DX = record[i].first;    
			double DY = record[i].second;

			const double RN = (BX - AX) * (DY - CY) - (BY - AY) * (DX - CX);
			if (RN != 0.0) {
				const double R = ((AY - CY) * (DX - CX) - (AX - CX) * (DY - CY)) / RN;
				const double S = ((AY - CY) * (BX - AX) - (AX - CX) * (BY - AY)) / RN;
				if (R >= 0.0 && R <= 1.0 && S >= 0.0 && S <= 1.0) {
					nCut++;
				}
			}
		}
	}
}

/**
 * Finds the nearest border point within service area from a reference point
 * 
 * @param longitude Reference longitude in degrees
 * @param latitude Reference latitude in degrees
 * @param nLong Output parameter for nearest point longitude
 * @param nLat Output parameter for nearest point latitude
 * @param longX Search center longitude in degrees
 * @param latX Search center latitude in degrees
 * @param serviceArea Search radius in kilometers
 * @param filename Path to border definition file
 */
void NearestLinePointA(double longitude, double latitude, double &nLong, double &nLat,
					double longX, double latX, double serviceArea, std::string filename) {
	double maxdi = 1.0e7;

	if (borderCacheAll.find(filename) == borderCacheAll.end()) {
		std::ifstream file(filename, std::ios::binary);
		if (!file.is_open()) {
			return;
		}

		union BorderRecord {
			char bytes[176];
			double coords[22];
		} buffer;

		std::vector<std::vector<std::pair<double, double>>> records;
		while (file.read(buffer.bytes, 176)) {
			std::vector<std::pair<double, double>> points;
			for (int i = 0; i < 22; i += 2) { 
				points.emplace_back(buffer.coords[i] * RadToDeg, buffer.coords[i+1] * RadToDeg);
			}
			records.push_back(points);
		}
		file.close();
		borderCacheAll[filename] = std::move(records);
	}


	const auto& records = borderCacheAll.at(filename);


	double dis, di;
	for (const auto& record : records) {
		for (size_t i = 0; i < record.size() - 1; ++i) {
			const auto& [x, y] = record[i];
			di = CalcDistanceM(longX, latX, x, y);
			if (di <= serviceArea) {
				dis = CalcDistanceM(longitude, latitude, x, y);
				if (dis < maxdi) {
					maxdi = dis;
					nLong = x;
					nLat = y;
				}
			}
		}
	}
}


/**
 * Calculates transmitter position considering border constraints and service area
 * 
 * @param longTx Transmitter longitude in degrees
 * @param latTx Transmitter latitude in degrees
 * @param longRx Receiver longitude in degrees
 * @param latRx Receiver latitude in degrees
 * @param Dir_Tx_Rx Direction from transmitter to receiver in degrees
 * @param Distance Distance between stations in kilometers
 * @param newLongTx Output parameter for calculated transmitter longitude
 * @param newLatTx Output parameter for calculated transmitter latitude
 * @param Tx_serv_area Transmitter service area radius in kilometers
 * @param filename Path to transmitter border definition file
 */
void Calc_Tx_posA(double longTx, double latTx, double longRx, double latRx, double Dir_Tx_Rx, double Distance,
				double &newLongTx, double &newLatTx, double Tx_serv_area, std::string filename) {

	double dp1 = Distance;
	// NO further than Tx_serv_area
	if (dp1 > Tx_serv_area) dp1 = Tx_serv_area;
	int nCut;

	// how many times line crosses border
	TestCutA(Dir_Tx_Rx, longTx, latTx, dp1, nCut, filename);

	if (nCut % 2 == 0) {
		// if in and out
		New_coordinatesA(longTx, latTx, Dir_Tx_Rx, dp1, newLongTx, newLatTx);
	} else {
		// if crosses
		NearestLinePointA(longRx, latRx, newLongTx, newLatTx, longTx, latTx, Tx_serv_area, filename);
	}
}

/**
 * Calculates receiver position considering border constraints and service area
 * 
 * @param longTx Transmitter longitude in degrees
 * @param latTx Transmitter latitude in degrees
 * @param longRx Receiver longitude in degrees
 * @param latRx Receiver latitude in degrees
 * @param Dir_Rx_Tx Direction from receiver to transmitter in degrees
 * @param Distance Distance between stations in kilometers
 * @param newLongRx Output parameter for calculated receiver longitude
 * @param newLatRx Output parameter for calculated receiver latitude
 * @param Rx_serv_area Receiver service area radius in kilometers
 * @param filename Path to receiver border definition file
 */
void Calc_Rx_posA(double longTx, double latTx, double longRx, double latRx, double Dir_Rx_Tx, double Distance,
				double &newLongRx, double &newLatRx, double Rx_serv_area, std::string filename) {
	double dp1 = Distance;
	// NO further than Rx_serv_area
	if (dp1 > Rx_serv_area) dp1 = Rx_serv_area;

	int nCut;
	// how many times line crosses border
	TestCutA(Dir_Rx_Tx, longRx, latRx, dp1, nCut, filename);
	if (nCut % 2 == 0) {
		// if in and out
		New_coordinatesA(longRx, latRx, Dir_Rx_Tx, dp1, newLongRx, newLatRx);
	} else {
		// if crosses
		NearestLinePointA(longTx, latTx, newLongRx, newLatRx, longRx, latRx, Rx_serv_area, filename);
	}
}



/**
 * Calculates final positions for mobile stations considering service areas and border constraints
 * 
 * @param input Input parameters containing station coordinates and service areas
 * @param output Output parameters with calculated positions and directions
 */
void PositionOfMobileA(PointToPoint_Input_Parameters &input, PointToPoint_Output_Parameters &output){
	std::string filenameTx = std::string(input.Border_path) + std::string(input.Land_from) + ".all";
	std::string filenameRx = std::string(input.Border_path) + std::string(input.Land_to) + ".all";
	if (input.Tx_Serv_area > 0.0){
		if (!fileExists(filenameTx)){
			output.HCM_error = 1036;// "+The 'xxx.ALL' borderline file for Tx is missing"
			return;
		}
	}
	if (input.Rx_Serv_area > 0.0){
		if (!fileExists(filenameRx)){
			output.HCM_error = 1037;// "+The 'xxx.ALL' borderline file for Rx is missing"
			return;
		}
	}

	output.Distance = CalcDistanceM(input.Tx_Longitude, input.Tx_Latitude, input.Rx_Longitude, input.Rx_Latitude);
	output.Dir_Rx_Tx = CalcDirectionM(input.Rx_Longitude, input.Rx_Latitude, input.Tx_Longitude, input.Tx_Latitude);
	output.Dir_Tx_Rx = CalcDirectionM(input.Tx_Longitude, input.Tx_Latitude, input.Rx_Longitude, input.Rx_Latitude);

	int nCutTx = 0;
	int nCutRx = 0;
	// Check if Serv_area crosses border line
	if (input.Tx_Serv_area > 0.0) {
		double dp1 = output.Distance;
		if (dp1 > input.Tx_Serv_area) dp1 = input.Tx_Serv_area;

		TestCutA(output.Dir_Tx_Rx, input.Tx_Longitude, input.Tx_Latitude, dp1, nCutTx, filenameTx);
	}
	if (input.Rx_Serv_area > 0.0) {
		double dp1 = output.Distance;
		if (dp1 > input.Rx_Serv_area) dp1 = input.Rx_Serv_area;
		TestCutA(output.Dir_Rx_Tx, input.Rx_Longitude, input.Rx_Latitude, dp1, nCutRx, filenameRx);
	}
	if (nCutTx == 0 && nCutRx == 0) {
		if (output.Distance <= input.Tx_Serv_area + input.Rx_Serv_area) {
			output.Distance = 0.0;
			output.Info[7] = true;
			output.HCM_error = 1000;//"Distance between Tx and Rx = 0. Calculations not possible"
			output.Perm_FS = HCM4A_ERROR_VALUE_DOUBLE;
		} else {
			if (input.Tx_Serv_area > 0.0) {
				New_coordinatesA(input.Tx_Longitude, input.Tx_Latitude, output.Dir_Tx_Rx, input.Tx_Serv_area,
								output.Tx_Longitude_Calc, output.Tx_Latitude_Calc);
			}
			if (input.Rx_Serv_area > 0.0) {
				New_coordinatesA(input.Rx_Longitude, input.Rx_Latitude, output.Dir_Rx_Tx, input.Rx_Serv_area,
								output.Rx_Longitude_Calc, output.Rx_Latitude_Calc);
			}
		}
		return;
	}


	if (input.Tx_Serv_area > 0.0 && input.Rx_Serv_area > 0.0) {
		if (input.Tx_Serv_area > input.Rx_Serv_area) {

			//Calculate new coor output.Tx_Longitude_Calc, Tx_Latitude_Calc
			Calc_Tx_posA(input.Tx_Longitude, input.Tx_Latitude,
						input.Rx_Longitude, input.Rx_Latitude,
						output.Dir_Tx_Rx, output.Distance,
						output.Tx_Longitude_Calc, output.Tx_Latitude_Calc, input.Tx_Serv_area, filenameTx );

			output.Distance = CalcDistanceM(output.Tx_Longitude_Calc, output.Tx_Latitude_Calc,
											input.Rx_Longitude, input.Rx_Latitude);
			output.Dir_Tx_Rx = CalcDirectionM(output.Tx_Longitude_Calc, output.Tx_Latitude_Calc,
											input.Rx_Longitude, input.Rx_Latitude);
			output.Dir_Rx_Tx = CalcDirectionM(input.Rx_Longitude, input.Rx_Latitude,
											output.Tx_Longitude_Calc, output.Tx_Latitude_Calc);
			Calc_Rx_posA(output.Tx_Longitude_Calc, output.Tx_Latitude_Calc,
						input.Rx_Longitude, input.Rx_Latitude,
						output.Dir_Rx_Tx, output.Distance,
						output.Rx_Longitude_Calc, output.Rx_Latitude_Calc,
						input.Rx_Serv_area, filenameRx );
		} else {
			//Calculate new coor output.Rx_Longitude_Calc, Rx_Latitude_Calc
			Calc_Rx_posA(input.Tx_Longitude, input.Tx_Latitude,
						input.Rx_Longitude, input.Rx_Latitude,
						output.Dir_Rx_Tx, output.Distance,
						output.Rx_Longitude_Calc, output.Rx_Latitude_Calc,
						input.Rx_Serv_area, filenameRx );

			output.Distance = CalcDistanceM(output.Rx_Longitude_Calc, output.Rx_Latitude_Calc,
											input.Tx_Longitude, input.Tx_Latitude);
			output.Dir_Tx_Rx = CalcDirectionM(input.Tx_Longitude, input.Tx_Latitude,
											output.Rx_Longitude_Calc, output.Rx_Latitude_Calc);
			output.Dir_Rx_Tx = CalcDirectionM(output.Rx_Longitude_Calc, output.Rx_Latitude_Calc,
											input.Tx_Longitude, input.Tx_Latitude);
			Calc_Tx_posA(input.Tx_Longitude, input.Tx_Latitude,
						output.Rx_Longitude_Calc, output.Rx_Latitude_Calc,
						output.Dir_Tx_Rx, output.Distance,
						output.Tx_Longitude_Calc, output.Tx_Latitude_Calc,
						input.Tx_Serv_area, filenameTx );
		}
	} else if (input.Tx_Serv_area > 0.0) {
		Calc_Tx_posA(input.Tx_Longitude, input.Tx_Latitude,
					input.Rx_Longitude, input.Rx_Latitude,
					output.Dir_Tx_Rx, output.Distance,
					output.Tx_Longitude_Calc, output.Tx_Latitude_Calc,
					input.Tx_Serv_area, filenameTx );
	} else if (input.Rx_Serv_area > 0.0) {
		Calc_Rx_posA(input.Tx_Longitude, input.Tx_Latitude,
					input.Rx_Longitude, input.Rx_Latitude,
					output.Dir_Rx_Tx, output.Distance,
					output.Rx_Longitude_Calc, output.Rx_Latitude_Calc,
					input.Rx_Serv_area, filenameRx );
	}
};

/**
 * Converts coordinate string (DDMMSS) to decimal degrees
 * 
 * @param source Coordinate string in "DDDMMSSXDDDMMSSY" format
 * @param lo Output parameter for longitude in decimal degrees
 * @param la Output parameter for latitude in decimal degrees
 */
void StrtoCoor(std::string source, double& lo, double& la) {
	lo = 0;
	la = 0;

	lo = std::stod(source.substr(0, 3)) +
				std::stod(source.substr(4, 2)) / 60.0 +
				std::stod(source.substr(6, 2)) / 3600.0;
	if (source[3] == 'W') lo = -lo;
	la = std::stod(source.substr(8, 2)) +
				std::stod(source.substr(11, 2)) / 60.0 +
				std::stod(source.substr(13, 2)) / 3600.0;
	if (source[10] == 'S') la = -la;
}

/**
 * Converts decimal degrees coordinate to string format (DDMMSS)
 * 
 * @param lo Longitude in decimal degrees
 * @param la Latitude in decimal degrees
 * @return Coordinate string in "DDDMMSSXDDDMMSSY" format
 */
std::string CoorToString(double lo, double la) {
	std::string coo = "               ";

	if (lo >= 0.0) {
		coo[3] = 'E';
	} else {
		coo[3] = 'W';
		lo = -lo;
	}

	if (la >= 0.0) {
		coo[10] = 'N';
	} else {
		coo[10] = 'S';
		la = -la;
	}

	int i = static_cast<int>(lo);
	lo = (lo - i) * 60.0;
	int j = static_cast<int>(lo);
	lo = (lo - j) * 60.0;
	int k = std::round(lo);
	if (k >= 60) {
		k -= 60;
		j += 1;
		if (j >= 60) {
			j -= 60;
			i += 1;
		}
	}

	// Write the longitude part
	std::ostringstream oss;
	oss << std::setw(3) << std::setfill('0') << i;
	coo.replace(0, 3, oss.str());
	oss.str("");
	oss.clear();
	oss << std::setw(2) << std::setfill('0') << j;
	coo.replace(4, 2, oss.str());
	oss.str("");
	oss.clear();
	oss << std::setw(2) << std::setfill('0') << k;
	coo.replace(6, 2, oss.str());

	i = static_cast<int>(la);
	la = (la - i) * 60.0;
	j = static_cast<int>(la);
	la = (la - j) * 60.0;
	k = std::round(la);
	if (k >= 60) {
		k -= 60;
		j += 1;
		if (j >= 60) {
			j -= 60;
			i += 1;
		}
	}

	// Write the latitude part
	oss.str("");
	oss.clear();
	oss << std::setw(2) << std::setfill('0') << i;
	coo.replace(8, 2, oss.str());
	oss.str("");
	oss.clear();
	oss << std::setw(2) << std::setfill('0') << j;
	coo.replace(11, 2, oss.str());
	oss.str("");
	oss.clear();
	oss << std::setw(2) << std::setfill('0') << k;
	coo.replace(13, 2, oss.str());

	return coo;
}

/**
 * Calculates terrain height difference (Dh) from terrain profile data
 * 
 * @param dist Distance between stations in kilometers
 * @param terra Terrain profile data structure
 * @return Height difference in meters after 10% trimming
 */
double DhCalculationM(double dist, HcmTerrain& terra) {
	std::vector<int> HI(1670);
	int PStart, PStop;
	int I, J, K, L, N, S1, S2, S3;

	// Prepare useful values
	N = 0;
	I = static_cast<int>(std::round(4.5 / terra.stepProfile));
	L = static_cast<int>(std::round(25.0 / terra.stepProfile));

	// Starting number of profile point 'PStart' and ending number 'PStop'
	if (dist <= 50.0) {
		// Between 10 km and 50 km: Normal calculation minus two times 4.5 km (transmitter and receiver).
		PStart = I + 1;
		PStop = terra.countSteps - I;
		for (J = PStart; J <= PStop; ++J) {
			HI[N++] = terra.T_Prof[J];
		}
	} else {
		// Distance > 50 km
		// Calculation of two parts: 4.5 to 25 km and distance - 25 km to distance - 4.5 km.
		// 1st part:
		PStart = I + 1;
		PStop = L + 1;
		for (J = PStart; J <= PStop; ++J) {
			HI[N++] = terra.T_Prof[J];
		}
		// 2nd part:
		PStart = terra.countSteps - L;
		PStop = terra.countSteps - I;
		for (J = PStart; J <= PStop; ++J) {
			HI[N++] = terra.T_Prof[J];
		}
	}

	// Sort all heights:
	S1 = N / 2;
	while (S1 != 0) {
		S2 = N - S1;
		for (S3 = 1; S3 <= S2; ++S3) {
			for (K = S3; K >= 1; K -= S1) {
				L = K + S1;
				if (HI[K] < HI[L]) break;
				std::swap(HI[K], HI[L]);
			}
		}
		S1 /= 2;
	}

	// Subtract 10% at the top and at the bottom:
	// Remaining difference between the lowest and the highest value is the value of 'Dh':
	if (N < 5) {
		return static_cast<double>(HI[N - 1] - HI[0]);
	} else {
		return static_cast<double>(HI[N - static_cast<int>(std::round(static_cast<double>(N) / 10.0))] - HI[static_cast<int>(std::round(static_cast<double>(N) / 10.0)) - 1]);
	}
}

/**
 * Applies correction factor based on terrain height difference for propagation calculations
 * 
 * @param dh Terrain height difference in meters
 * @param distance Path distance in kilometers
 * @param frequency Operating frequency in MHz
 * @param dh_corr Output parameter for calculated correction value
 */
void DhCorrectionM(double dh, double distance, double frequency, double& dh_corr) {
	if ((distance < 10.0) || (dh == 50.0)) {
		dh_corr = 0.0;
		return;
	}

	double a = dh;
	if (dh <= 10.0) a = 10.0;
	if (dh >= 500.0) a = 500.0;

	int f_inf, f_sup;
	if (frequency <= 600.0) {
		f_inf = 100;
		f_sup = 600;
	} else {
		f_inf = 600;
		f_sup = 2000;
	}

	double a1_100 = 0.0, a2_100 = 0.0;
	double a1_600 = 0.0, a2_600 = 0.0;
	double a1_2000 = 0.0, a2_2000 = 0.0;

	if (f_inf == 100) {
		if (distance < 200.0) {
			if (a <= 20.0) a1_100 = (a - 10.0) * 0.3 - 7.0;
			if ((a > 20.0) && (a <= 30.0)) a1_100 = (a - 30.0) * 0.15 - 2.5;
			if ((a > 30.0) && (a <= 50.0)) a1_100 = (a - 50.0) * 0.125;
			if ((a > 50.0) && (a <= 80.0)) a1_100 = (a - 80.0) * 0.1 + 3.0;
			if ((a > 80.0) && (a <= 100.0)) a1_100 = (a - 100.0) * 0.1 + 5.0;
			if ((a > 100.0) && (a <= 150.0)) a1_100 = (a - 150.0) * 0.06 + 8.0;
			if ((a > 150.0) && (a <= 300.0)) a1_100 = (a - 300.0) * 0.04 + 14.0;
			if (a > 300.0) a1_100 = (a - 500.0) * 0.025 + 19.0;
		}
		if (distance >= 200.0) {
			if (a <= 20.0) a2_100 = (a - 20.0) * 0.1 - 2.0;
			if ((a > 20.0) && (a <= 30.0)) a2_100 = (a - 30.0) * 0.05 - 1.5;
			if ((a > 30.0) && (a <= 50.0)) a2_100 = (a - 50.0) * 0.075;
			if ((a > 50.0) && (a <= 80.0)) a2_100 = (a - 80.0) * 0.067 + 2.0;
			if ((a > 80.0) && (a <= 100.0)) a2_100 = (a - 100.0) * 0.05 + 3.0;
			if ((a > 100.0) && (a <= 150.0)) a2_100 = (a - 150.0) * 0.03 + 4.5;
			if ((a > 150.0) && (a <= 300.0)) a2_100 = (a - 300.0) * 0.0167 + 7.0;
			if (a > 300.0) a2_100 = (a - 500.0) * 0.0125 + 9.5;
		}
	}

	if ((f_inf == 600) || (f_sup == 600)) {
		if (distance < 200.0) {
			if (a <= 20.0) a1_600 = (a - 10.0) * 0.4 - 10.0;
			if ((a > 20.0) && (a <= 30.0)) a1_600 = (a - 30.0) * 0.3 - 3.0;
			if ((a > 30.0) && (a <= 50.0)) a1_600 = (a - 50.0) * 0.15;
			if ((a > 50.0) && (a <= 80.0)) a1_600 = (a - 80.0) * 0.133 + 4.0;
			if ((a > 80.0) && (a <= 100.0)) a1_600 = (a - 100.0) * 0.15 + 7.0;
			if ((a > 100.0) && (a <= 150.0)) a1_600 = (a - 150.0) * 0.06 + 10.0;
			if ((a > 150.0) && (a <= 300.0)) a1_600 = (a - 300.0) * 0.067 + 20.0;
			if (a > 300.0) a1_600 = (a - 500.0) * 0.04 + 28.0;
		}
		if (distance >= 200.0) {
			if (a <= 20.0) a2_600 = (a - 20.0) * 0.2 - 3.0;
			if ((a > 20.0) && (a <= 30.0)) a2_600 = (a - 30.0) * 0.1 - 2.0;
			if ((a > 30.0) && (a <= 50.0)) a2_600 = (a - 50.0) * 0.1;
			if ((a > 50.0) && (a <= 80.0)) a2_600 = (a - 80.0) * 0.067 + 2.0;
			if ((a > 80.0) && (a <= 100.0)) a2_600 = (a - 100.0) * 0.075 + 3.5;
			if ((a > 100.0) && (a <= 150.0)) a2_600 = (a - 150.0) * 0.03 + 5.0;
			if ((a > 150.0) && (a <= 300.0)) a2_600 = (a - 300.0) * 0.033 + 10.0;
			if (a > 300.0) a2_600 = (a - 500.0) * 0.015 + 13.0;
		}
	}

	if (f_sup == 2000) {
		if (distance < 200.0) {
			if (a <= 20.0) a1_2000 = (a - 10.0) * 0.4 - 10.0;
			if ((a > 20.0) && (a <= 30.0)) a1_2000 = (a - 30.0) * 0.3 - 3.0;
			if ((a > 30.0) && (a <= 50.0)) a1_2000 = (a - 50.0) * 0.15;
			if ((a > 50.0) && (a <= 80.0)) a1_2000 = (a - 80.0) * 0.1667 + 5.0;
			if ((a > 80.0) && (a <= 100.0)) a1_2000 = (a - 100.0) * 0.185 + 8.7;
			if ((a > 100.0) && (a <= 150.0)) a1_2000 = (a - 150.0) * 0.074 + 12.4;
			if ((a > 150.0) && (a <= 300.0)) a1_2000 = (a - 300.0) * 0.082667 + 24.8;
			if (a > 300.0) a1_2000 = (a - 500.0) * 0.0495 + 34.7;
		}
		if (distance >= 200.0) {
			if (a <= 20.0) a2_2000 = (a - 20.0) * 0.2 - 3.0;
			if ((a > 20.0) && (a <= 30.0)) a2_2000 = (a - 30.0) * 0.1 - 2.0;
			if ((a > 30.0) && (a <= 50.0)) a2_2000 = (a - 50.0) * 0.1;
			if ((a > 50.0) && (a <= 80.0)) a2_2000 = (a - 80.0) * 0.08333 + 2.5;
			if ((a > 80.0) && (a <= 100.0)) a2_2000 = (a - 100.0) * 0.09 + 4.3;
			if ((a > 100.0) && (a <= 150.0)) a2_2000 = (a - 150.0) * 0.038 + 6.2;
			if ((a > 150.0) && (a <= 300.0)) a2_2000 = (a - 300.0) * 0.041333 + 12.4;
			if (a > 300.0) a2_2000 = (a - 500.0) * 0.0185 + 16.1;
		}
	}

	double c100, c600, c2000;
	if (distance <= 50.0) {
		c100 = a1_100 * (distance - 10.0) / 40.0;
		c600 = a1_600 * (distance - 10.0) / 40.0;
		c2000 = a1_2000 * (distance - 10.0) / 40.0;
	} else if ((distance > 50.0) && (distance <= 100.0)) {
		c100 = a1_100;
		c600 = a1_600;
		c2000 = a1_2000;
	} else if ((distance > 100.0) && (distance <= 200.0)) {
		c100 = a2_100 - (distance - 200.0) * (a1_100 - a2_100) / 100.0;
		c600 = a2_600 - (distance - 200.0) * (a1_600 - a2_600) / 100.0;
		c2000 = a2_2000 - (distance - 200.0) * (a1_2000 - a2_2000) / 100.0;
	} else {
		c100 = a2_100;
		c600 = a2_600;
		c2000 = a2_2000;
	}

	if ((f_inf == 100) && (f_sup == 600)) {
		dh_corr = c100 + (c600 - c100) * std::log10(frequency / 100.0) / std::log10(600.0 / 100.0);
	} else if ((f_inf == 600) && (f_sup == 2000)) {
		dh_corr = c600 + (c2000 - c600) * std::log10(frequency / 600.0) / std::log10(2000.0 / 600.0);
	}
}

/**
 * Calculates permissible field strength for point-to-point or line calculations
 * 
 * @param input Input parameters for stations and propagation
 * @param output Output parameters with calculated field strength and corrections
 * @param p2p Flag indicating point-to-point (true) or line (false) calculation
 */
void Permissible_FS_calculationM(PointToPoint_Input_Parameters &input, PointToPoint_Output_Parameters &output, bool p2p){
	bool TX_TETRA, RX_TETRA, TX_DIG, RX_DIG;
	double CDFN, CDFW, GANT, acorrB1, acorrsinus, OMEGA, B1, B2, X1, FACTOR ;
	int I, CSXR, CSXT;
	std::string DRX, DTX;

	output.Perm_FS = 0;
	output.Corr_delta_f = 0;
	output.Rx_ant_corr = 0;

	output.Delta_frequency = std::round(std::abs(input.Tx_Freq - input.Rx_Freq) * 1e6);

	// Input value for correction factor according delta frequency ?
	if (input.Cor_fact_frequ_diff != HCM4A_ERROR_VALUE_DOUBLE) {
		output.Info[14] = true;
		// Input value of correction factor according frequency difference is used
		output.Corr_delta_f = input.Cor_fact_frequ_diff;
	} else {
		std::string txDesig(input.Tx_Desig_of_emis);
		std::string rxDesig(input.Rx_Desig_of_emis);
		
		TX_DIG = (txDesig.substr(5, 1).find_first_of("1279") != std::string::npos);
		RX_DIG = (rxDesig.substr(5, 1).find_first_of("1279") != std::string::npos);
		TX_TETRA = (txDesig.substr(0, 7) == "25K0G7W");
		RX_TETRA = (rxDesig.substr(0, 7) == "25K0G7W");

		// Bandwidth of TX:
		if (TX_TETRA && !RX_TETRA) {
			CSXT = 16000;
		} else {
			DTX = std::string(input.Tx_Desig_of_emis).substr(0, 4);
			if (DTX[3] == 'H') {
				FACTOR = 1.0;
				I = 4;
			} else {
				I = DTX.find('K');
				if (I == std::string::npos) {
					I = DTX.find('M');
					if (I == std::string::npos) {
						output.HCM_error = 1041;
						// Channel spacing outside definition range (Tx)!
						return;
					} else {
						FACTOR = 1000000.0;
					}
				} else {
					FACTOR = 1000.0;
				}
			}
			// Replace 'K' or 'M' with '.':
			DTX[I] = '.';
			try {
				X1 = std::stod(DTX);
			} catch (const std::invalid_argument&) {
				output.HCM_error = 1041;
				// Channel spacing outside definition range (Tx)!
				return;
			}
			CSXT = std::round(X1 * FACTOR);
		}

		// bail out for line calculations
		if (!p2p) {
			output.Delta_frequency = 0;
			// WB/NB correction
			if (TX_DIG && (CSXT > 16000) && (input.Tx_Freq <= 470.0) && (input.Perm_FS_input == HCM4A_ERROR_VALUE_DOUBLE)) {
				output.Perm_FS += 6 * std::log10(CSXT / 25000.0);
			}
			return;
		}

		// shortcut for TETRA <> TETRA
		if (TX_TETRA && RX_TETRA) {
			// Annex3.4.3
			if (output.Delta_frequency <= 24999) {
				output.Corr_delta_f = 0.0;
			} else if (output.Delta_frequency <= 50000) {
				output.Corr_delta_f = 45.0;
			} else {
				output.Corr_delta_f = 70.0;
			}
		} else {
			// Bandwidth of Rx:
			if (RX_TETRA && !TX_TETRA) {
				CSXR = 16000;
			} else {
				DRX = std::string(input.Rx_Desig_of_emis).substr(0, 4);
				if (DRX[3] == 'H') {
					FACTOR = 1.0;
					I = 4;
				} else {
					I = DRX.find('K');
					if (I == std::string::npos) {
						I = DRX.find('M');
						if (I == std::string::npos) {
							output.HCM_error = 1040;
							// Channel spacing outside definition range (Rx)!
							return;
						} else {
							FACTOR = 1000000.0;
						}
					} else {
						FACTOR = 1000.0;
					}
				}
				// Replace 'K' or 'M' with '.':
				DRX[I] = '.';
				try {
					X1 = std::stod(DRX);
				} catch (const std::invalid_argument&) {
					output.HCM_error = 1040;
					// Channel spacing outside definition range (Rx)!
					return;
				}
				CSXR = std::round(X1 * FACTOR);
			}

			// Annex3.2
			if (CSXT > CSXR) {
				B1 = CSXT;
				B2 = CSXR;
			} else {
				B2 = CSXT;
				B1 = CSXR;
			}

			OMEGA = output.Delta_frequency / B1;

			// Annex3.4.1 TETRA Tx vs. narrowband Rx
			if (TX_TETRA && (CSXR <= 25000)) {
				// acorrB1:
				if (OMEGA < 0.5) {
					acorrB1 = 0.0;
				} else if (OMEGA <= 1.0) {
					acorrB1 = 32.0 * OMEGA - 16.0;
				} else if (OMEGA <= 1.4) {
					acorrB1 = 112.0 * OMEGA - 96.0;
				} else {
					acorrB1 = 41.0 * OMEGA;
				}

				// acorrsinus:
				if (OMEGA < 0.4) {
					acorrsinus = 0.0;
				} else if (OMEGA <= 0.7) {
					acorrsinus = 50.0 * OMEGA - 21.0;
				} else if (OMEGA <= 1.0) {
					acorrsinus = 225.0 * OMEGA - 145.0;
				} else {
					acorrsinus = -20.0 * OMEGA + 100.0;
				}
				output.Corr_delta_f = std::min(acorrsinus - (acorrsinus - acorrB1) * (B2 / B1), 70.0);
			} else if (RX_TETRA && (CSXT <= 25000)) {
				// Annex3.4.2 narrowband Tx vs. TETRA Rx
				// acorrB1:
				if (OMEGA < 0.45) {
					acorrB1 = 0.0;
				} else if (OMEGA <= 0.63) {
					acorrB1 = 55.0 * OMEGA - 23.0;
				} else if (OMEGA <= 0.93) {
					acorrB1 = 180.0 * OMEGA - 100.0;
				} else {
					acorrB1 = 12.5 * OMEGA + 57.0;
				}

				// acorrsinus:
				if (OMEGA < 0.45) {
					acorrsinus = 0.0;
				} else if (OMEGA <= 0.7) {
					acorrsinus = 225.0 * OMEGA - 101.0;
				} else {
					acorrsinus = 13.0 * OMEGA + 58.0;
				}
				output.Corr_delta_f = std::min(acorrsinus - (acorrsinus - acorrB1) * (B2 / B1), 70.0);
			} else {
				// Annex3.3 narrowband w/o TETRA
				// acorrB1:
				if (OMEGA < 0.5) {
					acorrB1 = 0.0;
				} else if (OMEGA <= 1.0) {
					acorrB1 = 47.0 * OMEGA - 24.0;
				} else if (OMEGA <= 1.3) {
					acorrB1 = 80.0 * OMEGA - 55.0;
				} else {
					acorrB1 = 38.0 * OMEGA;
				}

				// acorrsinus:
				if (OMEGA < 0.5) {
					acorrsinus = 0.0;
				} else if (OMEGA <= 1.3) {
					acorrsinus = 88.0 * OMEGA - 44.0;
				} else {
					acorrsinus = 12.0 * OMEGA + 55.0;
				}
				output.Corr_delta_f = std::min(acorrsinus - (acorrsinus - acorrB1) * (B2 / B1), 70.0);
				if ((CSXT <= 25000) && (CSXR <= 25000)) {
					CDFN = output.Corr_delta_f;
				} else {
					// Annex3.5 wideband
					// acorrB1:
					if (OMEGA < 0.5) {
						acorrB1 = 0.0;
					} else if (OMEGA <= 2.0) {
						acorrB1 = 33.3 * OMEGA - 16.7;
					} else {
						acorrB1 = 10.0 * OMEGA + 30.0;
					}

					// acorrsinus:
					if (OMEGA < 0.5) {
						acorrsinus = 0.0;
					} else if (OMEGA <= 1.25) {
						acorrsinus = 66.7 * OMEGA - 33.3;
					} else if (OMEGA <= 1.75) {
						acorrsinus = 20.0 * OMEGA + 25.0;
					} else {
						acorrsinus = 4.8 * OMEGA + 51.6;
					}
					output.Corr_delta_f = std::min(acorrsinus - (acorrsinus - acorrB1) * (B2 / B1), 70.0);
					if ((CSXR >= 200000) || (CSXT >= 200000)) {
						CDFW = output.Corr_delta_f;
					} else {
						// Annex3.6
						output.Corr_delta_f = CDFN + ((CDFW - CDFN) / 175.0) * (B1 / 1000.0 - 25.0);
					}
				}
			}
		}
	}

	output.Perm_FS += output.Corr_delta_f;


	// Calculation of antenna correction factors "Rx_ant_corr" and "Rx_ant_type_corr":
	output.V_angle_Rx_Tx = std::atan2((output.Tx_H_Site_Ant + input.Tx_H_Ant) - (output.Rx_H_Site_Ant + input.Rx_H_Ant), (1000* output.Distance));
	output.V_angle_Rx_Tx = output.V_angle_Rx_Tx * RadToDeg;

	if ((std::strcmp(input.Rx_Ant_typ_V, "000ND00") != 0) || (std::strcmp(input.Rx_Ant_typ_H, "000ND00") != 0)){

		Antenna_correction(output.Dir_Rx_Tx, input.Rx_Azi, output.V_angle_Rx_Tx, input.Rx_Elev,
							input.Rx_Ant_typ_H, input.Rx_Ant_typ_V,
							output.H_diff_angle_Rx_Tx, output.V_diff_angle_Rx_Tx,
							output.Rx_ant_corr, output.HCM_error);
		if (output.HCM_error != 0) return;
	}

	output.Perm_FS += output.Rx_ant_corr;

	if (input.Rx_Type_of_Ant == 0) {
		output.Rx_ant_type_corr = 0.0;
	} else if (input.Rx_Type_of_Ant == 1) {
		output.Rx_ant_type_corr = 2.1;
	} else {
		output.HCM_error = 1044;
		// Error in typ of Rx antenna (E/I)
		return;
	}

	output.Perm_FS += output.Rx_ant_type_corr;
	output.Perm_FS -= input.Rx_Ant_gain;

	// Depolarisation loss:
	if (input.Depol_loss == HCM4A_ERROR_VALUE_DOUBLE){
		output.HCM_error = 1046;		// Error in input data of depolarization loss
		return;
	}
	double DPN = input.Depol_loss;

	// automatic DPN calculation
	if (input.Depol_loss <= -9.8) {
		double X1 = output.Free_space_FS - output.Calculated_FS;
		if ((X1 < 50.0) && (input.Rx_H_Ant > 10) && (input.Tx_H_Ant > 10) &&
			(input.Tx_Serv_area + input.Rx_Serv_area == 0.0) &&
			((output.Rx_ant_corr <= 10.0) || (output.Tx_ant_corr <= 10.0))) {
			DPN = 25.0 - 0.5 * (X1);
		} else {
			DPN = 0.0;
		}
	}
	output.Depol_loss_calc = DPN;

	output.Perm_FS += output.Depol_loss_calc;
}


/**
 * Adjusts coordinates for coordination boundary region (CBR) calculations
 * 
 * @param LongLi Line point longitude in degrees
 * @param LatLi Line point latitude in degrees
 * @param LongTx Transmitter longitude in degrees
 * @param LatTx Transmitter latitude in degrees
 * @param d2b Distance to border in kilometers
 * @param Tx_serv_area Transmitter service area in kilometers
 * @param filename Path to border definition file
 */
void CBR_CoordinatesM(double& LongLi, double& LatLi, double& LongTx, double& LatTx,
						double d2b, double Tx_serv_area, std::string filename) {

	// Calculate the direction from Tx_Site to Rx in linepoint
	double Dir_Tx_Rx = CalcDirectionM(LongTx, LatTx, LongLi, LatLi);

	// In case of Tx is a mobile, calculate new Tx coordinates:
	if (Tx_serv_area > 0.0) {
		double Lo, La, Dist;
		Dist = CalcDistanceM(LongTx, LatTx, LongLi, LatLi);
		Calc_Tx_posA(LongTx, LatTx, LongLi, LatLi, Dir_Tx_Rx, Dist, Lo, La, Tx_serv_area, filename);
		LongTx = Lo;
		LatTx = La;
	}
	if (d2b < 0.0) {
		New_coordinatesA(LongLi, LatLi, Dir_Tx_Rx, std::abs(d2b), LongLi, LatLi);
	}
}

/**
 * Tests if line segment between two points crosses any border polygon segments
 * 
 * @param LONG Start point longitude in degrees
 * @param LAT Start point latitude in degrees
 * @param N_LONG End point longitude in degrees
 * @param N_LAT End point latitude in degrees
 * @param Co_cp Vector of border polygon coordinates
 * @param N_cp Number of coordinate pairs in polygon
 * @return true if no crossings detected, false if crossing exists
 */
bool Test_cutM(double LONG, double LAT, double N_LONG, double N_LAT, std::vector<Coordinate>& Co_cp, int N_cp) {
	double AX = LONG;
	double AY = LAT;
	double BX = N_LONG;
	double BY = N_LAT;

	int I = 2;
	// Take first two points:
	double CX = Co_cp[0].longitude;
	double CY = Co_cp[0].latitude;
	double DX = Co_cp[1].longitude;
	double DY = Co_cp[1].latitude;

	// Determine intersection:
	while (true) {
		double RN = (BX - AX) * (DY - CY) - (BY - AY) * (DX - CX);
		if (RN == 0.0) {
			CX = DX;
			CY = DY;
			I += 1;
			if (I > N_cp)
				return true;
			DX = Co_cp[I].longitude;
			DY = Co_cp[I].latitude;
			continue;
		}
		double R = ((AY - CY) * (DX - CX) - (AX - CX) * (DY - CY)) / RN;
		double S = ((AY - CY) * (BX - AX) - (AX - CX) * (BY - AY)) / RN;
		if ((R >= 0.0) && (R <= 1.0) && (S >= 0.0) && (S <= 1.0)) {
			return false;
		}

		CX = DX;
		CY = DY;
		I += 1;
		if (I > N_cp)
			return true;
		DX = Co_cp[I].longitude;
		DY = Co_cp[I].latitude;
	}
}

/**
 * Counts number of times line segment crosses border polygon segments
 * 
 * @param LONG Start point longitude in degrees
 * @param LAT Start point latitude in degrees
 * @param N_LONG End point longitude in degrees
 * @param N_LAT End point latitude in degrees
 * @param Co_cp Vector of border polygon coordinates
 * @param N_cp Number of coordinate pairs in polygon
 * @return Number of border crossings
 */
int Test_cutCount(double LONG, double LAT, double N_LONG, double N_LAT, std::vector<Coordinate>& Co_cp, int N_cp) {
	double AX = LONG;
	double AY = LAT;
	double BX = N_LONG;
	double BY = N_LAT;

	int I = 2;
	int Count = 0;
	// Take first two points:
	double CX = Co_cp[0].longitude;
	double CY = Co_cp[0].latitude;
	double DX = Co_cp[1].longitude;
	double DY = Co_cp[1].latitude;

	// Determine intersection:
	while (true) {
		double RN = (BX - AX) * (DY - CY) - (BY - AY) * (DX - CX);
		if (RN == 0.0) {
			CX = DX;
			CY = DY;
			I += 1;
			if (I > N_cp)
				return Count;
			DX = Co_cp[I].longitude;
			DY = Co_cp[I].latitude;
			continue;
		}
		double R = ((AY - CY) * (DX - CX) - (AX - CX) * (DY - CY)) / RN;
		double S = ((AY - CY) * (BX - AX) - (AX - CX) * (BY - AY)) / RN;
		if ((R >= 0.0) && (R <= 1.0) && (S >= 0.0) && (S <= 1.0)) {
			Count++;
		}

		CX = DX;
		CY = DY;
		I += 1;
		if (I > N_cp)
			return Count;
		DX = Co_cp[I].longitude;
		DY = Co_cp[I].latitude;
	}
	return Count;
}

/**
 * Manages a list of top 3 field strength values and their corresponding record indexes
 * 
 * @param N_rec Current record index
 * @param Calculated_FS Current calculated field strength value
 * @param ListIndexes Vector of 3 record indexes (modified in-place)
 * @param ListMaxFS Vector of 3 field strength values (modified in-place)
 */

void Manage_ListM(int N_rec, double Calculated_FS, std::vector<int>& ListIndexes, std::vector<double>& ListMaxFS){

	int minIndex = 0;
	double minFS = ListMaxFS[0];

	for (int i = 0; i < 3; ++i) {
		if (ListMaxFS[i] >  minFS) {
			minFS = ListMaxFS[i];
			minIndex = i;
		}
	}

	if (Calculated_FS > minFS) {
		ListMaxFS[minIndex] = Calculated_FS;
		ListIndexes[minIndex] = N_rec;
	}
}

/**
 * Calculates terrain clearance angle (TCA) correction for propagation models
 * 
 * @param TCAI Terrain clearance angle in degrees
 * @param Frequency Operating frequency in MHz
 * @param Distance Path distance in kilometers
 * @return TCA correction value in dB
 */
double TCACorrectionCalculation(double TCAI, double Frequency, double Distance) {
	int f_inf, f_sup;
	double TCA_c_100, TCA_c_600, TCA_c_2000;
	double v100, v600, v2000;
	double TCAI_corr;

	// Calculate v values
	v100 = 0.649 * TCAI;
	v600 = 1.592 * TCAI;
	v2000 = 2.915 * TCAI;

	// Compute TCA corrections
	TCA_c_100 = 9.1 - (6.9 + 20.0 * std::log10(std::sqrt(std::pow(v100 - 0.1, 2) + 1.0) + (v100 - 0.1)));
	TCA_c_600 = 13.1 - (6.9 + 20.0 * std::log10(std::sqrt(std::pow(v600 - 0.1, 2) + 1.0) + (v600 - 0.1)));
	TCA_c_2000 = 17.3 - (6.9 + 20.0 * std::log10(std::sqrt(std::pow(v2000 - 0.1, 2) + 1.0) + (v2000 - 0.1)));

	// Apply limits to corrections
	if (TCA_c_100 >= 0.0) TCA_c_100 = 0.0;
	if (TCA_c_100 < -32.0) TCA_c_100 = -32.0;

	if (TCA_c_600 >= 0.0) TCA_c_600 = 0.0;
	if (TCA_c_600 < -35.0) TCA_c_600 = -35.0;

	if (TCA_c_2000 >= 0.0) TCA_c_2000 = 0.0;
	if (TCA_c_2000 < -36.0) TCA_c_2000 = -36.0;

	// Determine frequency range
	if (Frequency <= 600.0) {
		f_inf = 100;
		f_sup = 600;
	} else {
		f_inf = 600;
		f_sup = 2000;
	}

	// Interpolate correction based on frequency
	if (f_inf == 100 && f_sup == 600) {
		TCAI_corr = TCA_c_100 + (TCA_c_600 - TCA_c_100) * std::log10(Frequency / 100.0) * std::log10(600.0 / 100.0);
	} else if (f_inf == 600 && f_sup == 2000) {
		TCAI_corr = TCA_c_600 + (TCA_c_2000 - TCA_c_600) * std::log10(Frequency / 600.0) * std::log10(2000.0 / 600.0);
	} else {
		TCAI_corr = 0.0; // Default case, adjust if needed
	}

	// Apply distance correction if necessary
	if (Distance < 16.0) {
		TCAI_corr *= Distance / 16.0;
	}

	return TCAI_corr;
}


/**
 * Parses emission class string to extract bandwidth value in Hertz
 * 
 * @param emissionClassIn Emission designation string (e.g., "25K0G7W")
 * @return Bandwidth value in Hertz, or 0 if parsing fails
 */
int64_t parseEmissionBand(const std::string& emissionClassIn) {
	std::string emissionClass = Trimmer(emissionClassIn);
	
	if (emissionClass.empty()) return 0;

	size_t numEnd = 0;
	while (numEnd < emissionClass.size() && isdigit(emissionClass[numEnd])) {
		numEnd++;
	}

	if (numEnd == 0) return 0;

	int64_t baseValue = std::stoull(emissionClass.substr(0, numEnd));

	if (numEnd >= emissionClass.size()) {
		return baseValue;
	}

	char multiplierChar = emissionClass[numEnd];
	int64_t multiplier;

	switch (multiplierChar) {
		case 'H': multiplier = 1; break;
		case 'K': multiplier = 1'000; break;
		case 'M': multiplier = 1'000'000; break;
		case 'G': multiplier = 1'000'000'000; break;
		default : multiplier = 1;
	}

	double preciseValue = baseValue;
	if (numEnd + 1 < emissionClass.size() && isdigit(emissionClass[numEnd + 1])) {
		preciseValue += (emissionClass[numEnd + 1] - '0') * 0.1;
	}

	return static_cast<int64_t>(preciseValue * multiplier);
}

#endif
