############################################################################
#
#   Name:           free_space_loss.py
#   Python:         3.3.2 on win32
#   Autor:          Adrian Haas
#   
############################################################################
#Free space loss
############################################################################
#
# The script calculates a free space loss map, depending of antenna type
# and frequency. The result could be mapped due html/javascript to
# swisstopo/googlemaps (KLM mapping not supported anymore). The mapping is
# disabled in this script.
#
############################################################################


import numpy as np
import warnings
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib as mpl

#from wgs84_distance import*
#from html import*

############################################################################
#free space loss equation d: distance in km, f: frequency in GHz
############################################################################

fsl=lambda d,f: 92.45+20*np.log10(d)+20*np.log10(f)

############################################################################
#antenna model
############################################################################


#create gaussian antenna pattern  (fantasy model)
pattern_deg=np.arange(-np.pi,np.pi,np.pi/160)
func0=lambda t,mean,var: np.exp(-(t-mean)**2/2/var)
pattern_att=40*(func0(pattern_deg,0,np.pi/2)-1)
 #-1: 0 as max.; 40 for lower the side att. a little

#direction of antenna
dir_gamma=3*np.pi/4
pattern_att=np.roll(pattern_att, int(len(pattern_att)*dir_gamma/2/np.pi))


############################################################################
#Create and calculate tables
############################################################################


#creating area: distance between points=res*1km
res=0.004
xlen=2000
ylen=2000

#geocoding
place=np.array([245055,679460]) #Uetliberg

#place antenna
centerx=xlen/2
centery=ylen/2


#distance matrix
dp=lambda i,j:((i-centerx)**2+(j-centery)**2)**0.5*res

#set frequency GHz
f=1.8

# ignore warning due log at center (idle problem)
warnings.filterwarnings("ignore")

fsl_table=np.fromfunction(lambda i,j:fsl(dp(i,j),f),(xlen,ylen),
                          dtype=int)

#set center from -inf to nb. value
fsl_table[centerx,centery]=fsl_table[centerx+1,centery+1]

#calculate isotropic radiation, outputpower in dBm
outputpower=45
field=outputpower-fsl_table


#degree matrix and antenna model
degree_table=np.fromfunction(lambda i,j:
                             np.angle((i-centerx)+1j*(j-centery))
                             ,(xlen,ylen),dtype=int)

#map antenna attenuation depending on degree

matrix=degree_table.copy()
matrix_ref=matrix.copy()
k=0
while k<(len(pattern_att))/2:
    mask=matrix_ref>=pattern_deg[k]
    matrix[mask]=pattern_att[k]
    k+=1

while k<(len(pattern_att))-1:
    mask=matrix_ref>=pattern_deg[k]
    matrix[mask]=pattern_att[k+1]
    k+=1


#caluclate antenna radiation; antenna gain in dBi
antennagain=20
field2=field+matrix

x=place[0]
y=place[1]
pix=400
name='out'
map_nr=1
xlen=res*xlen*1000
ylen=res*ylen*1000


############################################################################
#matplotlib output
############################################################################

plt.figure("Antenna pattern")
plt.plot(pattern_deg,pattern_att)
plt.xlabel("rad")
plt.ylabel("dB")

plt.figure("GIS Model für Funkfeld")
plt.subplot(221)
plt.title("isotropischer Kugelstrahler")
plt.imshow(field,interpolation='nearest', cmap=cm.jet,vmin=-80,
           vmax=-40,  aspect=1)
plt.colorbar().set_label("$dBm$",rotation=0, labelpad=20)
plt.subplot(222)
plt.title("Richtungsmatrix")
plt.imshow(degree_table,interpolation='nearest', cmap=cm.jet,  aspect=1)
plt.colorbar().set_label("$rad$",rotation=0, labelpad=20)
plt.subplot(223)
plt.title("Patternmatrix")
plt.imshow(matrix,interpolation='nearest', cmap=cm.jet,  aspect=1)
plt.colorbar().set_label("$dB$",rotation=0, labelpad=20)
plt.subplot(224)
plt.title("Funkfeld mit Antenne")
plt.imshow(field2,interpolation='nearest', cmap=cm.jet,vmin=-80,
           vmax=-40,  aspect=1)
plt.colorbar().set_label("$dBm$",rotation=0, labelpad=20)
plt.show()