MGR - 19.04

#!/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 = 200
mi = 0.7
centre_range = mi * R_c
edge_range = R_c - centre_range
f = 2.5 #2.5 # 2.1 * 10**9 #[GHz]
h_UE = 1.5 #[m]
h_BS = 12 # [m]
d = R_c ### ?
R_0 = 100 #[m], nominal reference distance
gamma = 4 #path loss exponent                          # !!!!!!!!!!
alfa = ...
omega = 1.4
C = 3*(10**8) #[m/s]

#noise
k = 1.38 * (10**(-23))
T_0 = 297
BW = 20 #[MHz] ?, w atykule BW = 5MHz i S = 25 sub-ch
S = 100 #number of subchannels in the whole bandwidth
N = k * T_0 * ((BW * (10**6))/S) * (10**3)                           # funkcja zależna od schematu 1 i 2 !!!!!!!

#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 = 15
G_R = 0
L_T = 5 #7 #?
L_R = 0 #1# ?
P_T = 10 #[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
    print("potega :",str(power_1))

    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): #do wywalenia
    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 = centre_range #R_0_prim = calculateFresnelZoneBreakpoint(h_BS, h_UE, C, f)
    fi_result = calculateFi((R_0_prim - 1), R_0_prim, C, f) # fi dla center
    fi_2 = calculateFi((R_0_prim + 1), R_0_prim, C, f) # fi dla edge
    print("fi_gwiazdka: ", str(fi_2))
    print("fi_center: ", str(fi_result))

    ni_result_c = calculateNi((R_0_prim - 1), R_0_prim)
    ni_result_e = calculateNi((R_0_prim + 1), R_0_prim)

    ### 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):
        x.append(r/R_c)

        if(r < R_0_prim):
            sinr = calculateSINRinCellCenter(r, mi, omega, P_tc_1, ni_result_c, fi_result, fi_2, scheme = 1)
        else:
            sinr = calculateSINRinCellEdge(r, mi, omega, P_tc_1, ni_result_e, 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)
    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)
    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 intelligenFrequencyReuse():
    file_3 = open("intelligent_FR.txt", "w")
    file_3.write("INTELLIGENT FREQUENCY REUSE - wyniki w kolejnych krokach "); file_3.write(" \n")
    #### dopisać algorytm
    file_3.close()

def calculateGainInICIC():
    file_4 = open("calculate_gain.txt", "w")
    file_4.write("Oblicz  zysk w ICIC"); file_4.write(" \n")

    R_0_prim = centre_range  #R_0_prim = calculateFresnelZoneBreakpoint(h_BS, h_UE, C, f)
    calculated_difference_1 = []
    calculated_difference_2 = []
    fi_result = calculateFi((R_0_prim-10), R_0_prim, C, f)
    fi_2 = calculateFi((R_0_prim+10), R_0_prim, C, f)
    x_ch =[]
    ni_result = 2 #calculateNi(R_0_prim, R_0_prim)
    omega_ch = 1.01

    for omega_ch in range(101, 280, 1):
        omega_ch = omega_ch/100
        P_tc = calculateP_Tc(P_T, omega_ch, scheme = 1)
        file_4.write("Dla: "); file_4.write(str(P_tc)); file_4.write(" sinr na granicy centrum = ");
        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)
        file_4.write(str(sinr_c_in_dB))
        sinr_e = calculateSINRinCellEdge((R_0_prim + 1), mi, omega_ch, P_tc, gamma, fi_result, fi_2, scheme = 1)
        sinr_e_in_dB = 10*math.log10(sinr_e)
        file_4.write(" ,na poczatku czesci granicznej = "); file_4.write(str(sinr_e_in_dB))
        gain_1 = sinr_e_in_dB - sinr_c_in_dB
        file_4.write(" , a roznica wynosila = "); file_4.write(str(gain_1)); file_4.write(" \n")

        calculated_difference_1.append(gain_1)


        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 = 2)
        sinr_c_in_dB = 10*math.log10(sinr_c)
        sinr_e = calculateSINRinCellEdge((R_0_prim + 1), mi, omega_ch, P_tc, gamma, fi_result, fi_2, scheme = 2)
        sinr_e_in_dB = 10*math.log10(sinr_e)
        gain_2 = sinr_e_in_dB - sinr_c_in_dB

        calculated_difference_2.append(gain_2)
        x_ch.append(omega_ch)
        print(omega_ch)

    ## PLOT BOTH GRAPHS ON ONE ##
    plt.figure(4)
    plt.plot(x_ch, calculated_difference_1, label='Scheme I', markersize=20)
    plt.plot(x_ch, calculated_difference_2, label='Schemat II', linestyle='dashed')
    plt.legend(loc='upper right')
    plt.xlabel('omega')
    plt.ylabel('Gain in SINR value [dB]')
    plt.title("Wzrost w wartości SINR w zależności od wzrostu wartości omegi")
    plt.savefig('SINR diff.png')

    file_4.close()

def plotSINRforDifferentOmegas():
    R_0_prim = centre_range
    fi_result = calculateFi((R_0_prim - 1), R_0_prim, C, f) # fi dla center
    fi_2 = calculateFi((R_0_prim + 1), R_0_prim, C, f) # fi dla edge

    ni_result_c = calculateNi((R_0_prim - 1), R_0_prim)
    ni_result_e = calculateNi((R_0_prim + 1), R_0_prim)

    # omegas_for_plotting = [1.0, 1.2, 1.5, 1.8, 2.0]
    omegas_for_plotting = [1.0, 1.5, 2.0, 2.5, 3.5]

    dict_scheme_1 = {}

    ### SCHEME 1 ###
    for omega_1 in omegas_for_plotting:
        calculated_SINR_sch_1 = []
        x = []
        P_tc_1 = calculateP_Tc(P_T, omega_1, scheme = 1)

        for r in range(1, R_c):
            x.append(r/R_c)

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

        dict_scheme_1[omega_1] = calculated_SINR_sch_1


    ## GRAPHS FOR SCHEME 1 ##
    plt.figure(5)
    for omega_1 in omegas_for_plotting:
        plt.plot(x, dict_scheme_1[omega_1], label='Schemat I - omega: '+str(omega_1))

    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('SINR_for_few_omegas_SCH_1.png')


    ### SCHEME 2 ###
    dict_scheme_2 = {}

    for omega_1 in omegas_for_plotting:
        P_tc_1 = calculateP_Tc(P_T, omega_1, scheme = 2)
        calculated_SINR_sch_2 = []

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

            if(r < R_0_prim):
                sinr = calculateSINRinCellCenter(r, mi, omega_1, P_tc_1, ni_result, fi_result, fi_2, scheme = 2)
            else:
                sinr = calculateSINRinCellEdge(r, mi, omega_1, 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)

        dict_scheme_2[omega_1] = calculated_SINR_sch_2

    ## GRAPHS FOR SCHEME 2 ##
    plt.figure(6)
    for omega_1 in omegas_for_plotting:
        plt.plot(x, dict_scheme_2[omega_1], label='Schemat I - omega: '+str(omega_1))

    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('SINR_for_few_omegas_SCH_2.png')


def intelligentFrequencyReuse():
    R_0_prim = centre_range
    fi_result = calculateFi((R_0_prim - 1), R_0_prim, C, f) # fi dla center
    fi_2 = calculateFi((R_0_prim + 1), R_0_prim, C, f) # fi dla edge

    ni_result_c = calculateNi((R_0_prim - 1), R_0_prim)
    ni_result_e = calculateNi((R_0_prim + 1), R_0_prim)

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

    for r in range(1, R_c):
        x.append(r/R_c)

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

    plt.figure(1)
    plt.plot(x, calculated_SINR_sch_1)
    plt.savefig('IntelligentFR.png')


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


if __name__ == "__main__":
    main()