#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <string>
#include <algorithm>
#include <fstream>
#include <cctype>
#include <vector>
#include <sys/stat.h>

#include "HcmTerrain.h"
#include "Structures.h"
#include "Additional.h"


void PointToPointForLine(PointToPoint_Input_Parameters& input,
						PointToPoint_Output_Parameters& output,
						HcmTerrain& terra){
	output.HCM_error = 0;

	output.Prot_margin = 0;
	
	output.Tx_Longitude_Calc = input.Tx_Longitude;
	output.Tx_Latitude_Calc = input.Tx_Latitude;
	output.Rx_Longitude_Calc = input.Rx_Longitude;
	output.Rx_Latitude_Calc = input.Rx_Latitude;

	
	// for backward compatibility
	std::string tmpS = CoorToString(output.Tx_Longitude_Calc, output.Tx_Latitude_Calc);
	StrtoCoor(tmpS, output.Tx_Longitude_Calc, output.Tx_Latitude_Calc);
	tmpS = CoorToString(output.Rx_Longitude_Calc, output.Rx_Latitude_Calc);
	StrtoCoor(tmpS, output.Rx_Longitude_Calc, output.Rx_Latitude_Calc);
	// for backward compatibility
	
	
	output.Distance = CalcDistanceM(output.Tx_Longitude_Calc, output.Tx_Latitude_Calc,
									output.Rx_Longitude_Calc, output.Rx_Latitude_Calc);
	output.Dir_Rx_Tx = CalcDirectionM(output.Rx_Longitude_Calc, output.Rx_Latitude_Calc,
									output.Tx_Longitude_Calc, output.Tx_Latitude_Calc);
	output.Dir_Tx_Rx = CalcDirectionM(output.Tx_Longitude_Calc, output.Tx_Latitude_Calc,
									output.Rx_Longitude_Calc, output.Rx_Latitude_Calc);

	if (output.Distance < input.PointDistance) {
		output.Info[7] = true;
		output.HCM_error = 1000;
		return;
	}
	if (output.Distance >= input.MaxCalcDistanceKm) {
		output.Calculated_FS = -999.9;
		output.HCM_error = 1028;
		return;
	}

	output.Rx_H_Site_Ant = std::trunc(terra.GetHeigth(output.Rx_Longitude_Calc, output.Rx_Latitude_Calc));
	output.Info[8] = true;
	if ( terra.isErrorHeigth(output.Rx_H_Site_Ant) ){
		output.HCM_error = 36;//"Error opening terrain- or morphological data file (data not available)
		return;
	}

	
	output.Tx_H_Site_Ant = input.Tx_H_Site_Ant;
	if (input.Tx_H_Site_Ant == HCM4A_NO_INPUT_INT || input.Tx_Serv_area > 0.0){
		output.Info[1] = true; //The height is taken from terrain database
		
		output.Tx_H_Site_Ant = std::trunc(terra.GetHeigth(output.Tx_Longitude_Calc, output.Tx_Latitude_Calc));
		if ( terra.isErrorHeigth(output.Tx_H_Site_Ant) ){
			output.HCM_error = 36;//"Error opening terrain- or morphological data file (data not available)
			return;
		}
	} else {
		output.Tx_H_Site_Ant = std::trunc(terra.GetHeigth(output.Tx_Longitude_Calc, output.Tx_Latitude_Calc));
		if ( terra.isErrorHeigth(output.Tx_H_Site_Ant) ){
			output.HCM_error = 36;//"Error opening terrain- or morphological data file (data not available)
			return;
		}
		if (std::abs(output.Tx_H_Site_Ant - input.Tx_H_Site_Ant) <= input.Tx_H_Site_Ant / 10) {
			output.Info[2] = true; //Height of Tx site differs from height of terrain data
		} else {
			output.Info[3] = true; // Tx site height differs more than 10%, calculated values may be (extremely) wrong!
		}
		output.Tx_H_Site_Ant = input.Tx_H_Site_Ant;
	}
	

	output.V_angle_Tx_Rx = std::atan2((output.Rx_H_Site_Ant + input.Rx_H_Ant) - (output.Tx_H_Site_Ant + input.Tx_H_Ant), (1000* output.Distance));
	output.V_angle_Tx_Rx = output.V_angle_Tx_Rx * RadToDeg;
	
	// Calculation of power in direction of the receiver 'Power_to_Rx'
	if (std::strcmp(input.Tx_Ant_typ_H, "000ND00") == 0 &&
	    std::strcmp(input.Tx_Ant_typ_V, "000ND00") == 0) {
		output.Tx_ant_corr = 0.0;
		output.H_diff_angle_Tx_Rx = 0.0;
		output.V_diff_angle_Tx_Rx = 0.0;
	} else {
		if (input.Tx_Elev <-90 || input.Tx_Elev >90){
			output.HCM_error = 1031; //Error in Tx elevation
			return;
		}
		if (input.Tx_Azi <0 || input.Tx_Azi >360){
			output.HCM_error = 1032; //Error in Tx azimuth
			return;
		}
		
		// NOT MODIFIED
		Antenna_correction(output.Dir_Tx_Rx, input.Tx_Azi, output.V_angle_Tx_Rx, input.Tx_Elev,
							input.Tx_Ant_typ_H, input.Tx_Ant_typ_V,
							//out
							output.H_diff_angle_Tx_Rx, output.V_diff_angle_Tx_Rx,
							output.Tx_ant_corr, output.HCM_error);
		if (output.HCM_error != 0) return;
	}
	
	input.Tx_Type_of_Ant == 0 ? output.Tx_ant_type_corr = 0.0 :	output.Tx_ant_type_corr = 2.1;
	output.Power_to_Rx = input.Tx_Power - output.Tx_ant_corr - output.Tx_ant_type_corr;
	
	output.Free_space_FS = 77.0 - 20.0 * std::log10(output.Distance) + output.Power_to_Rx;

	// If distance is less than 1 km, use free space field strength
	if (output.Distance < 1.0) {
		output.Calculated_FS = output.Free_space_FS;
		output.Info[11] = true; //Free space field strength used because distance < 1 km
		return;
	}
	
	terra.GetProfile( output.Tx_Longitude_Calc, output.Tx_Latitude_Calc,
					output.Rx_Longitude_Calc, output.Rx_Latitude_Calc,
					false, output.Tx_H_Site_Ant, output.Rx_H_Site_Ant);
					
	output.ProfilePoints = terra.countSteps;
	
	// Calculate Terrain Clearance Angles and TCA correction factors
	double xTCA = 0.0; //temp value
	
	output.Tx_TCA_corr = 0;
	if (input.Tx_Serv_area == 0.0) {
		output.Tx_TCA = -90.0;
		int jtx = (output.Distance >= 16.0) ? std::round(16.0 / input.PointDistance) : output.ProfilePoints - 2;
		for (int i = 1; i <= jtx; ++i) {
			xTCA = (terra.T_Prof[1 + i] - (input.Tx_H_Ant + output.Tx_H_Site_Ant)) / (input.PointDistance * i * 1000.0);
			xTCA = std::atan(xTCA) * RadToDeg; // Convert to degrees
			if (xTCA > output.Tx_TCA) output.Tx_TCA = xTCA;
		}
		
		output.Tx_TCA_corr = TCACorrectionCalculation(output.Tx_TCA, input.Tx_Freq, output.Distance);
	}
	
	output.Rx_TCA_corr = 0;

	// Effective antenna heights
	int d1,d2;
	if (output.Distance < 15.0) {
		d1 = std::round(output.ProfilePoints / 15.0);
		d2 = output.ProfilePoints - 1;
	}else{
		d1 = std::round(1.0 / input.PointDistance);
		d2 = std::round(15 / input.PointDistance);
	}
	
    // Transmitter:
    double hsum = 0.0; //temp value;
    
    if (input.Tx_Serv_area > 0.0) {
        // Height of a mobile 'hmTx'
        output.Heff_Tx = static_cast<double>(input.Tx_H_Ant);
    } else {
        hsum = 0.0;
        for (int I = d1; I <= d2; ++I) {
            hsum += terra.T_Prof[I];
        }
        output.Heff_Tx = input.Tx_H_Ant + output.Tx_H_Site_Ant - static_cast<double>(hsum) / static_cast<double>(d2 - d1 + 1);
        
    }
    if (output.Heff_Tx < 3.0) output.Heff_Tx = 3.0;

	output.Heff_Rx = input.Rx_H_Ant;
	if (output.Heff_Rx < 3.0) output.Heff_Rx = 3.0;

	output.Heff = (output.Heff_Tx * output.Heff_Rx) / 10.0;

	// Calculate distance over sea
	output.D_sea_calculated = 0.0;
	bool null = false;
	double ds1 = 0.0;
	int j1 = 0;
	if (output.Distance > 1.0) {
		for (int i = 1; i < output.ProfilePoints; ++i) {
			if (terra.M_Prof[i] == 1) { // Assuming 1 represents 'sea'
				if (null) {
					++j1;
				} else {
					j1 = 0;
					null = true;
				}
			} else {
				if (null) {
					ds1 = static_cast<double>(j1) * input.PointDistance;
					output.D_sea_calculated += ds1;
					ds1 = 0.0;
					j1 = 0;
					null = false;
				}
			}
		}
		if (null) {
			ds1 = static_cast<double>(j1) * input.PointDistance;
			output.D_sea_calculated += ds1;
			null = false;
		}
	}
	
	if (output.D_sea_calculated > output.Distance) {
		output.Info[13] = true;
		output.D_sea_calculated = output.Distance;
	}

	// Calculate Delta h and Dh correction factor
	if ((output.Distance <= 10.0) || (output.D_sea_calculated >= output.Distance) || (input.Tx_Serv_area > 0.0)) {
		output.Dh = 50.0;
		output.Dh_corr = 0.0;
	} else {
		output.Dh = DhCalculationM(output.Distance, terra);
		DhCorrectionM(output.Dh, output.Distance, input.Tx_Freq, output.Dh_corr);
	}

	// Calculate the field strength from the figures
	FS_figures fsf;
	fsf.Tx_frequency = input.Tx_Freq;
	fsf.Distance = output.Distance;     // km
	fsf.Heff = output.Heff;          // m
	fsf.Time_percentage = output.Time_percentage; // %
	fsf.LatM = 0;//(latTx + latRx)/2; // Middle Latitide of path
	fsf.Sea_temperature = input.Sea_temperature[0]; // Sea_temperature
	fsf.TxLongitude = input.Tx_Longitude;
	fsf.TxLatitude = input.Tx_Latitude;
	fsf.Get_FS_from_figures_Old();
	
	output.HCM_error = fsf.HCM_error;
	if (output.HCM_error != 0) return;
	
	output.Land_FS_old  = fsf.Land_FS_1kW - 30.0 + output.Power_to_Rx;
	
	fsf.Get_FS_from_figures();
	
	output.HCM_error = fsf.HCM_error;
	if (output.HCM_error != 0) return;
	
	// Corrections according to the power
	output.Land_FS = fsf.Land_FS_1kW - 30.0 + output.Power_to_Rx;
	output.Sea_FS = fsf.Sea_FS_1kW - 30.0 + output.Power_to_Rx;

	// Clearance angle corrections and delta h correction
	output.Land_FS += output.Tx_TCA_corr + output.Rx_TCA_corr - output.Dh_corr;

	if (output.D_sea_calculated == 0.0) {
		output.Calculated_FS = output.Land_FS;
	} else {
		// Mixed path calculation
		if (output.Time_percentage < 10) {
			std::vector<double> a = {0.0, 0.07, 0.13, 0.2, 0.29, 0.37, 0.46, 0.56, 0.67, 0.8, 1.0};
			double factorOfPathOverSea = output.D_sea_calculated / output.Distance;
			double x = factorOfPathOverSea * 10.0;
			double a1 = a[static_cast<int>(x) + 1];
			double a2, ax = 0;
			if (x < 10.0) {
				a2 = a[static_cast<int>(x) + 2];
			} else {
				a2 = a1;
			}
			x -= static_cast<int>(x);
			ax = a1 + x * (a2 - a1);
			output.Calculated_FS = output.Land_FS + ax * (output.Sea_FS - output.Land_FS);
		} else {
			output.Calculated_FS = output.Land_FS * (output.Distance - output.D_sea_calculated) / 
									output.Distance + output.Sea_FS * output.D_sea_calculated / output.Distance;
		}
		
	}
	
	if (output.HCM_error > 0) return;
	
	Permissible_FS_calculationM( input, output, false);
	if (!output.AnnexInFreqRange) {
		output.Perm_FS += input.Perm_FS_input;
	} else{
		output.Perm_FS += output.Perm_FS_from_table;
	}

	output.Prot_margin = output.Perm_FS - output.Calculated_FS;
}