Mgr

#!/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 = 600
mi = 0.6
centre_range = mi * R_c
edge_range = R_c - centre_range
f = 2.1 #* 10**9 #[GHz]
h_UE = 1.5 #[m]
h_BS = 10 # [m]
d = 1 ### ?
R_0 = 100 #[m], nominal reference distance
gamma = 4 #path loss exponent
alfa = ...
omega = 1.25
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)

#values needed to calculate fi and ni
Xf = 6 * math.log10(f/2) #correction factor for frequencies above 2 GHz
Xh = -d * 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 = 7 #?
L_R = 1# ?
P_T = 1 #[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)

    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):
    graph = panda.DataFrame(y, x)
    graph.plot(kind='line', grid=True, label='SINR [dB]')
    plt.show()
    #plt.savefig('wykres SINR.png')


def main():
    R_0_prim = calculateFresnelZoneBreakpoint(h_BS, h_UE, C, f)
    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)
        #x.append(r)

        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")

    plotGraph(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")

    plotGraph(x, calculated_SINR_sch_2)
    plt.savefig('wykres SINR (sch_2).png')
    file_2.close()

    ## PLOT BOTH GRAPHS ON ONE ##
    #graph.plot(kind='line', grid=True, label='SINR [dB]')
    plt.plot(x, calculated_SINR_sch_1, label='Scheme I')
    plt.plot(x, calculated_SINR_sch_2, label='Scheme II')
    plt.legend(loc='upper right')
    plt.xlabel('d/R')
    plt.ylabel('SINR [dB]')
    plt.savefig('Both schemes.png')

if __name__ == "__main__":
    main()