z parametrami Natalii

#!/usr/bin/env python
import pandas as panda
import numpy as np
import matplotlib.pyplot as plt
import math


#constants
freq_reuse_factor = 3
R_c = 1200
mi = 0.6
centre_range = 981 #mi * R_c
edge_range = R_c - centre_range
f = 2.3 #* 10**9 #[GHz]
h_UE = 2 #[m]
h_BS = 16 # [m]
d = R_c ### ?
R_0 = 100 #[m], nominal reference distance
gamma = 5.056 #path loss exponent
alfa = ...
omega = 1.5
C = 3*10**8 #[m/s]

#noise
k = 1.38 * 10**(-23)
T_0 = 293
BW = 5 #[MHz] ?, w atykule BW = 5MHz i S = 25 sub-ch
S = 25 #number of subchannels in the whole bandwidth
N = k * T_0 * (BW * 10**6/S)

#values needed to calculate fi and ni
Xf = 6 * math.log10(f/2) #correction factor for frequencies above 2 GHz
Xh = -R_c * math.log10(h_UE/2) #correction factor for the receiver
s = 5 #shadowing term, normal random variable (value from normal distribution) ???
T = 16 #number of subchannels available at the edge of a cell, T=4 w artykule

#bilans
G_T = 10
G_R = 0
L_T = 5 #?
L_R = 0# ?
P_T = 23 #[dBm]


def calculateFi(R, R_0_prim, c, freq):
    power_1 = (G_T + G_R - (s + Xf + Xh + L_T + L_R))/10
    power_2 = (G_T + G_R - (s + L_T + L_R))/10

    wavelght = c / (freq * 10**9)

    if(R/R_c > R_0_prim/R_c):
        fi = (R_0/R_c)**gamma * (wavelght/(4*math.pi*R_0_prim))**2 * 10**power_1
    else:
        fi = (wavelght/(4*math.pi*R_c))**2 * 10**power_2

    return fi

def changeToLinear(x):
    linear = 10**(x/10)
    return linear

def calculateNi(R, R_0_prim):
    if(R/R_c > R_0_prim/R_c):
        ni = gamma
    else:
        ni = 2

    return ni

def calculateP_Tc(P_T, omega, scheme = 1):
    P_T_Wat = 10**(P_T/10 - 3)
    if (scheme == 1):
        P_Tc = P_T_Wat/(omega*T + S - freq_reuse_factor*T)
    else:
        P_Tc = P_T_Wat/(omega*T + S - T)

    return P_Tc

def calculateFresnelZoneBreakpoint(h_t, h_r, c, freq):
    sum = h_t + h_r
    diff = h_t - h_r
    wavelght = c / (freq * 10**9)
    R_0_prim = 1/wavelght*math.sqrt((sum**2 - diff**2)**2 -2*(sum**2 + diff**2)*(wavelght/2)**2 + (wavelght/2)**4)
    #R_0_prim = (4*h_t*h_r)/wavelght
    return R_0_prim

def calculateSINRinCellCenter(R, mi, omega, P_Tc, ni, fi, fi_2, scheme = 1):
    if (scheme == 1):
        S = (R/R_c)**(-ni)
        I_c1 = 6*(fi_2/fi) * math.sqrt(3)**(-gamma)
        I_c2 = 12*(fi_2/fi) * math.sqrt(3*freq_reuse_factor)**(-gamma)
        Noise = N / (fi *P_Tc)

        SINR = S/(I_c1 + I_c2 + Noise)
    else:
        S = (R/R_c)**(-ni)
        I_c = 3*(fi_2/fi) * ((1 + omega)*math.sqrt(3)**(-gamma) + (3 + omega)*math.sqrt(3*freq_reuse_factor)**(-gamma))
        Noise = N / (fi *P_Tc)

        SINR = S/(I_c + Noise)

    return SINR

def calculateSINRinCellEdge(R, mi, omega, P_Tc, ni, fi, fi_2, scheme = 1):
    if (scheme == 1):
        S = (R/R_c)**(-ni)
        I_c = 6*(fi_2/fi) * math.sqrt(3*freq_reuse_factor)**(-gamma)
        Noise = N / (omega*fi *P_Tc)

        SINR = S/(I_c + Noise)
    else:
        S = (R/R_c)**(-ni)
        I_c = 6*(fi_2/fi) * ((1/omega)*math.sqrt(3)**(-gamma) + (1 + 1/omega)*math.sqrt(3*freq_reuse_factor)**(-gamma))
        Noise = N / (omega*fi *P_Tc)

        SINR = S/(I_c + Noise)

    return SINR

def plotGraph(x, y, name):
    graph = panda.DataFrame(y, x)
    graph.plot(kind='line', grid=True, label='SINR [dB]')
    plt.show()
    plt.savefig(name)

def calculateSINRAndPlotGraphs():
    #R_0_prim = calculateFresnelZoneBreakpoint(h_BS, h_UE, C, f)
    R_0_prim = 1201 #centre_range
    fi_result = calculateFi(200, R_0_prim, C, f)
    fi_2 = calculateFi(500, R_0_prim, C, f)

    ### SCHEME 1 ###
    calculated_SINR_sch_1 = []
    x = []
    P_tc_1 = calculateP_Tc(P_T, omega, scheme = 1)

    file = open("calculated_SINR.txt", "w")
    file.write("P_tc_1: "); file.write(str(P_tc_1)); file.write(" \n")
    file.write("R_0_prim: "); file.write(str(R_0_prim)); file.write(" \n")

    for r in range(1, R_c):
        ni_result = calculateNi(r, R_0_prim)
        x.append(r/R_c)

        if(r < R_0_prim):
            sinr = calculateSINRinCellCenter(r, mi, omega, P_tc_1, ni_result, fi_result, fi_2, scheme = 1)
        else:
            sinr = calculateSINRinCellEdge(r, mi, omega, P_tc_1, ni_result, fi_result, fi_2, scheme = 1)
        sinr_in_dB = 10*math.log10(sinr)
        calculated_SINR_sch_1.append(sinr_in_dB)

        file.write(str(sinr_in_dB))
        file.write(" \n")

    plt.figure(1)
    # plotGraph(x, calculated_SINR_sch_1, name = 'wykres SINR (sch_1).png')
    plt.plot(x, calculated_SINR_sch_1)
    plt.savefig('wykres SINR (sch_1).png')
    file.close()

    ### SCHEME 2 ###
    calculated_SINR_sch_2 = []
    P_tc_1 = calculateP_Tc(P_T, omega, scheme = 2)

    file_2 = open("calculated_SINR (sch_2).txt", "w")
    file_2.write("P_tc_1: "); file_2.write(str(P_tc_1)); file_2.write(" \n")
    file_2.write("R_0_prim: "); file_2.write(str(R_0_prim)); file_2.write(" \n")

    for r in range(1, R_c):
        ni_result = calculateNi(r, R_0_prim)

        if(r < R_0_prim):
            sinr = calculateSINRinCellCenter(r, mi, omega, P_tc_1, ni_result, fi_result, fi_2, scheme = 2)
        else:
            sinr = calculateSINRinCellEdge(r, mi, omega, P_tc_1, ni_result, fi_result, fi_2, scheme = 2)
        sinr_in_dB = 10*math.log10(sinr)
        calculated_SINR_sch_2.append(sinr_in_dB)

        file_2.write(str(sinr_in_dB))
        file_2.write(" \n")

    plt.figure(2)
    # plotGraph(x, calculated_SINR_sch_2, name = 'wykres SINR (sch_2).png')
    plt.plot(x, calculated_SINR_sch_2)
    plt.savefig('wykres SINR (sch_2).png')

    file_2.close()

    ## PLOT BOTH GRAPHS ON ONE ##
    plt.figure(3)
    plt.plot(x, calculated_SINR_sch_1, label='Schemat I')
    plt.plot(x, calculated_SINR_sch_2, label='Schemat II')
    plt.legend(loc='upper right')
    plt.title("Wartość SINR w funkcji odległości od stacji bazowej przy użyciu FFR")
    plt.xlabel('d/R')
    plt.ylabel('SINR [dB]')
    plt.savefig('Both schemes.png')

def calculateGainInICIC():
    #R_0_prim = calculateFresnelZoneBreakpoint(h_BS, h_UE, C, f)
    R_0_prim = centre_range
    calculated_difference = []
    fi_result = calculateFi(R_0_prim - 5, R_0_prim, C, f)
    fi_2 = calculateFi(R_0_prim + 5, R_0_prim, C, f)
    x_ch =[]
    ni_result = 2 #calculateNi(R_0_prim, R_0_prim)
    omega_ch = 1.01

    if(omega_ch < 1.9):
        P_tc = calculateP_Tc(P_T, omega_ch, scheme = 2)

        sinr_c = calculateSINRinCellCenter((R_0_prim - 1), mi, omega_ch, P_tc, ni_result, fi_result, fi_2, scheme = 1)
        sinr_c_in_dB = 10*math.log10(sinr_c)
        sinr_e = calculateSINRinCellEdge((R_0_prim + 1), mi, omega_ch, P_tc, ni_result, fi_result, fi_2, scheme = 1)
        sinr_e_in_dB = 10*math.log10(sinr_e)
        gain = sinr_e_in_dB - sinr_c_in_dB

        calculated_difference.append(gain)
        x_ch.append(omega_ch)
        omega_ch += 0.01

    ## PLOT BOTH GRAPHS ON ONE ##
    plt.figure(4)
    plt.plot(x_ch, calculated_difference, label='Scheme I')
    plt.legend(loc='upper right')
    plt.xlabel('omega')
    plt.ylabel('Gain in SINR value [dB]')
    plt.title("Schemat przedstawiający wzrost w wartości SINR w zależności od wzrostu wartości omegi")
    plt.savefig('SINR diff.png')

def main():
    calculateSINRAndPlotGraphs()
    #calculateGainInICIC()


if __name__ == "__main__":
    main()