# -*- coding:Utf-8 -*-
#from numpy import *
from __future__ import print_function
"""
.. currentmodule:: pylayers.antprop.loss
.. autosummary::
:members:
"""
import doctest
import os
import glob
import doctest
import logging
import numpy as np
from scipy import io
import matplotlib.pylab as plt
import pylayers.gis.gisutil as gu
from pylayers.util.project import logger
import numpy.linalg as la
import pdb
import time
from numba import jit
[docs]def PL0(fGHz,GtdB=0,GrdB=0,R=1):
""" Path Loss at frequency fGHZ @ R
Parameters
----------
fGHz: float
frequency GHz
GtdB: float
transmitting antenna gain dB (default 0 dB)
GrdB: float
receiving antenna gain dB (default 0 dB)
R : float
distance in m
Returns
-------
PL0 : float
path @ R
Notes
-----
.. math:: PL_0 = -20 log_{10}(\\frac{\\lambda}{4\\pi}) - GtdB -GrdB
Examples
--------
>>> fGHz = 2.4
>>> PL = PL0(fGHz)
>>> assert (PL<41)&(PL>40),"something wrong"
"""
if not isinstance(fGHz,np.ndarray):
fGHz=np.array([fGHz])
ld = 0.3/fGHz
PL0 = -20*np.log10(ld/(4.0*np.pi*R))-GtdB-GrdB
return PL0
[docs]def Dgrid_points(points,Px):
""" distance point to grid
Parameters
----------
points : np.array
grid Np x 2 array
Px : np.array
point 2 x 1 array
Returns
-------
D: Euclidian distance matrix
"""
Dx = points[:,0] - Px[0]
Dy = points[:,1] - Px[1]
D = np.sqrt( Dx*Dx + Dy*Dy )
return(D)
[docs]def FMetisShad2(fGHz,r,D,sign=1):
""" F Metis shadowing function
Parameters
----------
fGHz : np.array(Nf)
frequency GHz
r : np.array(Nseg,)
distance between Tx and Rx
D : np.array(Nseg,Nscreen)
indirect distance between Tx and Rx (screen effect)
sign : np.array(Nseg,Nscreen)
== 1 : Shadowing NLOS situation
==-1 : No shadowing LOS situation
Returns
-------
F : np.array(Nseg,Nscreen,Nf)
Notes
-----
Provides an implementation of formula (6.6) in D1.4 of METIS project
See Also
--------
LossMetisShadowing
The arctan is approximating the enveloppe of the Fresnel integral.
"""
lamda = 0.3/fGHz[None,None,:]
F = np.arctan(sign[:,:,None]*np.pi/2.*(np.sqrt((np.pi/lamda)*(D[:,:,None]-r[:,None,None])))) / np.pi
return(F)
[docs]def FMetisShad(fGHz,r,D,sign=1):
""" F Metis shadowing function
Parameters
----------
fGHz : float
frequency GHz
r : float
distance between Tx and Rx
D : float
indirect distance between Tx and Rx (screen effect)
sign : int
== 1 : Shadowing NLOS situation
==-1 : No shadowing LOS situation
Notes
-----
Provides an implementation of formula (6.6) in D1.4 of METIS project
See Also
--------
LossMetisShadowing
"""
lamda = 0.3/fGHz
F = np.arctan(sign*np.pi/2.*(np.sqrt((np.pi/lamda)*(D-r)))) / np.pi
return(F)
[docs]def LossMetisShadowing(fGHz,tx,rx,pg,uw,uh,w,h):
""" Calculate the Loss from
Parameters
----------
fGHz : float
tx : np.array (,3) of floats
transmiter coordinates
rx : np.array (,3) of floats
receiver coordinates
pg : np.array (,3) of floats
center of gravity of the screen
uw : np.array (,3) of floats
unitary vector along width dimension
uh : np.array (,3) of floats
unitary vector along height dimension
w : float
width in meters
h : float
height in meters
Returns
-------
Lsh : float
Loss in dB to add to the FS path Loss
Notes
-----
This function provides an implementation of formula 6.5 of D1.4 deliverable of METIS project
[Metis D1.4](Ahttps://www.metis2020.com/wp-content/uploads/METIS_D1.4_v3.pdf)
# geometry parametric issue : find M in [tx-rx] defined as M = alpha*rx + (1-alpha)tx where alpha in [0-1].
# if alpha = 0 then M = tx ; if alpha = 1 then M = rx.
# Besides, M is defined as M = pg + beta*uw + gamma*uh then alpha*rx + (1-alpha)tx = pg + beta*uw + gamma*uh
# [rx-tx , -uw, -uh]*[alpha,beta,gamma].T = pg - tx <==> Ax = b solved by la.solve ; x[0]=alpha, x[1]=beta and
TODO To be vectorized
"""
rxtx = rx - tx # LOS distance
# x[2]=gamma.
A = np.vstack((rxtx,-uw,-uh)).T
b = pg - tx
x = la.solve(A,b)
# condition of shadowing
condseg = ((x[0]>1) or (x[0]<0))
condw = ((x[1]>w/2.) or (x[1]<-w/2.))
condh = ((x[2]>h/2.) or (x[2]<-h/2.))
visi = condseg or condw or condh
if visi:
shad = -1
else:
shad = 1
r = np.dot(rxtx,rxtx)**0.5
w1 = pg + uw*w/2.
w2 = pg - uw*w/2.
h1 = pg + uh*h/2.
h2 = pg - uh*h/2.
Dtw1 = np.dot(tx-w1,tx-w1)**0.5
Drw1 = np.dot(rx-w1,rx-w1)**0.5
Dtw2 = np.dot(tx-w2,tx-w2)**0.5
Drw2 = np.dot(rx-w2,rx-w2)**0.5
Dth1 = np.dot(tx-h1,tx-h1)**0.5
Drh1 = np.dot(rx-h1,rx-h1)**0.5
Dth2 = np.dot(tx-h2,tx-h2)**0.5
Drh2 = np.dot(rx-h2,rx-h2)**0.5
D1w = Dtw1+Drw1
D1h = Dth1+Drh1
D2w = Dtw2+Drw2
D2h = Dth2+Drh2
if shad == 1:
signw1 = 1
signw2 = 1
signh1 = 1
signh2 = 1
else:
if condw:
if D1w>D2w:
signw1=1
signw2=-1
else:
signw1=-1
signw2=1
else:
signw1 = 1
signw2 = 1
if condh:
if D1h>D2h:
signh1=1
signh2=-1
else:
signh1=-1
signh2=1
else:
signh1 = 1
signh2 = 1
Fw1 = FMetisShad(fGHz,r,D1w,sign=signw1)
Fh1 = FMetisShad(fGHz,r,D1h,sign=signh1)
Fw2 = FMetisShad(fGHz,r,D2w,sign=signw2)
Fh2 = FMetisShad(fGHz,r,D2h,sign=signh2)
tmp = (Fh1+Fh2)*(Fw1+Fw2)
Lsh = -20*np.log10(1-tmp)
#return(Lsh,shad,tmp,Fw1,Fh1,Fw2,Fh2,condh,condw)
return(Lsh)
[docs]def LossMetisShadowing2(fGHz,tx,rx,pg,uw,uh,w,h):
""" Calculate the Loss from
Parameters
----------
fGHz : np.array(,Nf)
tx : np.array (3,Nseg) of floats
transmiter coordinates
rx : np.array (3,Nseg) of floats
receiver coordinates
pg : np.array (3,Nscreen) of floats
center of gravity of the screen
uw : np.array (3,Nscreen) of floats
unitary vector along width dimension
uh : np.array (3,Nscreen) of floats
unitary vector along height dimension
w : np.array (,Nscreen)
width in meters
h : np.array (,Nscreen)
height in meters
Returns
-------
Lsh : np.array (Nseg,Nscreen,Nf)
Loss in dB to add to the FS path Loss
Notes
-----
This function provides an implementation of formula 6.5 of D1.4 deliverable of METIS project
[Metis D1.4](Ahttps://www.metis2020.com/wp-content/uploads/METIS_D1.4_v3.pdf)
# geometry parametric issue : find M in [tx-rx] defined as M = alpha*rx + (1-alpha)tx where alpha in [0-1].
# if alpha = 0 then M = tx ; if alpha = 1 then M = rx.
# Besides, M is defined as M = pg + beta*uw + gamma*uh then alpha*rx + (1-alpha)tx = pg + beta*uw + gamma*uh
# [rx-tx , -uw, -uh]*[alpha,beta,gamma].T = pg - tx <==> Ax = b solved by la.solve ; x[0]=alpha, x[1]=beta and
"""
Nseg = tx.shape[1]
Nscreen = uw.shape[1]
rxtx = rx - tx # (3,Nseg) LOS distance
# A : (Nseg,Nscreen,3,3)
# b : (Nseg,Nscreen,3)
# rxtx.T (Nseg,3)
# uw.T (Nscreen, 3)
# uh.T (Nscreen, 3)
U = rxtx.T[:,None,:,None]
W = uw.T[None,:,:,None]
H = uh.T[None,:,:,None]
We = W + np.zeros(U.shape)
He = H + np.zeros(U.shape)
Ue = U + np.zeros(He.shape)
A = np.concatenate((Ue,-We,-He),axis=3)
#A = np.vstack((rxtx,-uw,-uh)).T
# pg.T Nscreen, 3
# tx.T Nseg,3
b = pg.T[None,:,:]-tx.T[:,None,:]
#b = pg - tx
x = la.solve(A,b)
# condition of shadowing
condseg = ((x[:,:,0]>1) + (x[:,:,0]<0))
condw = ((x[:,:,1]>w[None,:]/2.) + (x[:,:,1]<-w[None,:]/2.))
condh = ((x[:,:,2]>h[None,:]/2.) + (x[:,:,2]<-h[None,:]/2.))
visi = (condseg + condw + condh)%2
# if visi:
# shad = -1
# else:
# shad = 1
#shad = - visi
r = np.sum(rxtx*rxtx,axis=0)**0.5
w1 = pg + uw*w[None,:]/2.
w2 = pg - uw*w[None,:]/2.
h1 = pg + uh*h[None,:]/2.
h2 = pg - uh*h[None,:]/2.
Dtw1 = np.sum((tx[...,None]-w1[:,None,:])*(tx[...,None]-w1[:,None,:]),axis=0)**0.5
Drw1 = np.sum((rx[...,None]-w1[:,None,:])*(rx[...,None]-w1[:,None,:]),axis=0)**0.5
Dtw2 = np.sum((tx[...,None]-w2[:,None,:])*(tx[...,None]-w2[:,None,:]),axis=0)**0.5
Drw2 = np.sum((rx[...,None]-w2[:,None,:])*(rx[...,None]-w2[:,None,:]),axis=0)**0.5
Dth1 = np.sum((tx[...,None]-h1[:,None,:])*(tx[...,None]-h1[:,None,:]),axis=0)**0.5
Drh1 = np.sum((rx[...,None]-h1[:,None,:])*(rx[...,None]-h1[:,None,:]),axis=0)**0.5
Dth2 = np.sum((tx[...,None]-h2[:,None,:])*(tx[...,None]-h2[:,None,:]),axis=0)**0.5
Drh2 = np.sum((rx[...,None]-h2[:,None,:])*(rx[...,None]-h2[:,None,:]),axis=0)**0.5
# Drw1 = np.dot(rx-w1,rx-w1)**0.5
# Dtw2 = np.dot(tx-w2,tx-w2)**0.5
# Drw2 = np.dot(rx-w2,rx-w2)**0.5
# Dth1 = np.dot(tx-h1,tx-h1)**0.5
# Drh1 = np.dot(rx-h1,rx-h1)**0.5
# Dth2 = np.dot(tx-h2,tx-h2)**0.5
# Drh2 = np.dot(rx-h2,rx-h2)**0.5
D1w = Dtw1+Drw1
D1h = Dth1+Drh1
D2w = Dtw2+Drw2
D2h = Dth2+Drh2
signw1 = np.ones((Nseg,Nscreen))
signw2 = np.ones((Nseg,Nscreen))
signh1 = np.ones((Nseg,Nscreen))
signh2 = np.ones((Nseg,Nscreen))
condw1 = (visi*condw*(D1w<=D2w)).astype(bool)
condw2 = (visi*condw*(D1w>D2w)).astype(bool)
signw1[condw1]=-1
signw2[condw2]=-1
condh1 = (visi*condh*(D1h<=D2h)).astype(bool)
condh2 = (visi*condh*(D1h>D2h)).astype(bool)
signh1[condh1]=-1
signh2[condh2]=-1
Fw1 = FMetisShad2(fGHz,r,D1w,sign=signw1)
Fh1 = FMetisShad2(fGHz,r,D1h,sign=signh1)
Fw2 = FMetisShad2(fGHz,r,D2w,sign=signw2)
Fh2 = FMetisShad2(fGHz,r,D2h,sign=signh2)
tmp = (Fh1+Fh2)*(Fw1+Fw2)
Lsh = -20*np.log10(1-tmp)
#return(Lsh,shad,tmp,Fw1,Fh1,Fw2,Fh2,condh,condw)
return(Lsh)
[docs]def Dgrid_zone(zone,Px):
""" Distance point to zone
A zone is a quadrilateral zone.
Parameters
----------
zone : dictionnary
xmin xmax Nx
ymin ymax Ny
Px : np.array
point
Returns
-------
D : np.array Nx x Ny
Euclidian distance matrix
Notes
-----
Build the distance matrix between Tx and points in the zone
use broadcasting instead
"""
rx = np.linspace(zone['xmin'],zone['xmax'],zone['Nx'])
ry = np.linspace(zone['ymin'],zone['ymax'],zone['Ny'])
R_x = np.outer(np.ones(len(ry)),rx)
R_y = np.outer(ry,np.ones(len(rx)))
Dx = R_x - Px[0]
Dy = R_y - Px[1]
D = np.sqrt(Dx*Dx+Dy*Dy)
return (D)
[docs]def OneSlopeMdl(D,n,fGHz):
""" one slope model
Parameters
----------
D : np.array
distance array
n : float
path loss exponent
fGHz : np.array
frequency GHz
Returns
-------
PL : np.array
path loss as a function of distance
"""
PL = PL0(fGHz) + 10*n*np.log10(D)
return(PL)
[docs]def cost231(pBS,pMS,hroof,phir,wr,fMHz,wb=20,dB=True,city='medium'):
""" Walfish Ikegami model (COST 231)
Parameters
----------
pBS : np.array (3xNlink)
pMS : np.array (3xNlink)
hroof : np.array (1xNlink)
phir : np.array (1xNlink)
degrees
wr : np.array (1xNlink)
fMHz : np.array (1xNf)
wb : float
average building separation
dB : boolean
Returns
-------
PathLoss (Nlink,Nf)
References
----------
http://morse.colorado.edu/~tlen5510/text/classwebch3.html
Examples
--------
.. plot::
:include-source:
>>> from pylayers.antprop.loss import *
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> # Number of links and BS and MS heights
>>> Nlink = 100
>>> hBS = 300
>>> hMS = 1.5
>>> # hroof and phir are drawn uniformily at random
>>> hroof = 40*np.random.rand(Nlink)
>>> wr = 10*np.ones(Nlink)
>>> phir = 90*np.random.rand(Nlink)
>>> pMS = np.vstack((np.linspace(10,2500,Nlink),np.zeros(Nlink),hMS*np.ones(Nlink)))
>>> pBS = np.vstack((np.zeros(Nlink),np.zeros(Nlink),hBS*np.ones(Nlink)))
>>> # frequency range
>>> fMHz = np.linspace(700,1900,120)
>>> pl = cost231(pBS,pMS,hroof,phir,wr,fMHz)
>>> im = plt.imshow(pl,extent=(0,100,0.7,1.9))
>>> cb = plt.colorbar()
>>> cb.set_label('Loss (dB)')
>>> plt.axis('tight')
>>> plt.xlabel('Frequency (GHz)')
>>> plt.ylabel('Link Number')
>>> plt.title('100 WI Path Loss realizations ')
>>> plt.show()
"""
hBS = pBS[2,:][:,np.newaxis]
hMS = pMS[2,:][:,np.newaxis]
wr = wr[:,np.newaxis]
hroof = hroof[:,np.newaxis]
phir = phir[:,np.newaxis]
fMHz = fMHz[np.newaxis,:]
dm = np.sqrt(np.sum((pBS-pMS)*(pBS-pMS),axis=0))[:,np.newaxis]
dkm = dm/1000.
Nlink = len(dm)
pl0 = 32.4 + 20*np.log10(dkm) + 20*np.log10(fMHz)
delta_base = hBS-hroof
u035 = np.where((phir>=0) & (phir<35))
u3555 = np.where((phir>=35) & (phir<55))
u5590 = np.where((phir>=55) & (phir<90))
plori = np.zeros(Nlink)[:,np.newaxis]
# path loss due to orientation w.r.t road
plori[u035] = -10+0.354*phir[u035]
plori[u3555] = 2.5+0.075*phir[u3555]
plori[u5590] = 4.0-0.114*(phir[u5590]-55)
# rooftop to street
plrts = -16.9-10*np.log10(wr)+10*np.log10(fMHz)+20*np.log10(hroof-hMS)+plori
uroofsupBS = np.where(hBS>hroof)
uroofinfBS = np.where(hBS<=hroof)
udistsup500 = np.where((hBS<=hroof)&(dkm>0.5))
udistinf500 = np.where((hBS<=hroof)&(dkm<0.5))
plbsh = np.zeros((Nlink,1))
plbsh[uroofsupBS] = -18*np.log10(1+delta_base[uroofsupBS])
ka = 54*np.ones((Nlink,1))
ka[udistsup500] = ka[udistsup500]-0.8*delta_base[udistsup500]
ka[udistinf500] = ka[udistinf500]-0.8*delta_base[udistinf500]*dkm[udistinf500]/0.5
kd = 18*np.ones((Nlink,1))
kd[uroofinfBS] = kd[uroofinfBS]-15*delta_base[uroofinfBS]/hroof[uroofinfBS]
if city=='medium':
kf = -4+0.7*(fMHz/925.-1)
else:
kf = -4+1.5*(fMHz/925.-1)
plmsd = plbsh+ka+kd*np.log10(dkm)+kf*np.log10(fMHz)-9*np.log10(wb)
pl = pl0
padd = plmsd + plrts
ulosspos = np.where(padd>0)[0]
pl[ulosspos]=pl[ulosspos]+padd[ulosspos]
if not dB:
pl = 10**(-pl/20.)
return(pl)
[docs]def cost259(pMS,pBS,fMHz):
""" cost259 model
Parameters
----------
pMS : np.array (position of Mobile Station)
pBS : np.array (position of Base station)
fMHz : float
"""
dm = np.sqrt((pBS-pMS)*(pBS-pMS))
lmbd = 300/fMHz
pl = 10*2.6*np.log10(dm)+20*log10(4*np.pi/lmbd)
if not dB:
pl = 10**(-pl/20.);
return(pl)
[docs]def hata(pMS,pBS,fGHz,hMS,hBS,typ):
""" Hata Path loss model
Parameters
----------
pMS : np.array
Mobile position (meters)
pBS : np.array
Base station position (meters)
fGHz : np.array
hMS : height mobile station (m)
hBS : height base station (m)
typ : 'small' or 'big'
Returns
-------
L : Attenuation (dB)
Examples
--------
>>> d = np.linspace(100,5000,120)
>>> hBS = 30
>>> hMS = 1.5
>>> fGHz = 0.9
>>> pMS = np.array([d,0,hMS])
>>> pBS = np.array([d,0,hBS])
>>> L = hata(pMS,pBS,fGHz,hMS,hBS,'small')
Notes
-----
This model is valid until 1.5GHz, for higher frequency see
COST231-Hata model
References
----------
OKUMURA (Y.), OHMORI (E.), KAWANO (T.)
et FUKUA (K.). Field strength and its varia-
bility in UHF and VHF land-mobile radio ser-
vice. Rev. Elec. Commun. Lab., vol. 16, n o 9,
1968.
HATA (M.). Empirical formula for propaga-
tion loss in land mobile radio services. IEEE
Trans. Veh. Technol., vol. 29, pp. 317-325,
Aug. 1980
"""
dm = np.sqrt((pBS-pMS)*(pBS-pMS))
if (typ=='small'):
CH = (1.1*np.log10(fGHz*1000)-0.7)*hMS-(1.56*np.log10(fGHz*1000)-0.8)
if (typ=='big'):
if fGHz<0.2:
CH = 8.29*(np.log10(1.54*hMS)**2)-1.1
else:# valid until 1.5GHz
CH = 3.2*(np.log10(11.75*hMS)**2)-4.97
L = 69.55+26.16*np.log10(fGHz*1000)-13.82*np.log10(hBS)+(44.9-6.55*np.log10(hBS))*np.log10(dm/1000.)-CH
return(L)
[docs]def cost2100(pMS,pBS,fGHz,nfloor=1,dB=True):
""" cost 2100 model
Parameters
----------
pMS :
pBS :
fGHz : float
nfloor : int
dB : boolean
"""
# distance (meters)
dm = np.sqrt((pBS-pMS)*(pBS-pMS))
pl0 = 32.4+20*log10(dm)+20*np.log10(fGHz)
pld = nfloor*30
pl = pl0+pld
if not dB:
pl = 10**(-pl/20.)
return(pl)
[docs]def PL(fGHz,pts,p,n=2.0,dB=True,d0=1):
""" calculate Free Space Path Loss
Parameters
----------
fGHz : float
frequency (GHz)
pts : np.array (2xNp)
points
p : np.array (2x1) or (2xNp)
n : float
path loss exponent (default = 2)
dB : : boolean
return result in dB
Returns
-------
PL : np.array
path loss w.r.t distance and frequency
"""
shp = np.shape(p)
# assert(shp[0]==2)
D = np.sqrt(np.sum((pts-p)**2,axis=0))
# f x grid x ap
#PL = np.array([PL0(fGHz)])[:,np.newaxis] + 10*n*np.log10(D)[np.newaxis,:]
PL = PL0(fGHz,d0)[:,np.newaxis] + 10*n*np.log10(D/d0)[np.newaxis,:]
if not dB:
PL=10**(-PL/10)
return(PL)
[docs]def Lconcrete(fGHz):
""" 3GPP above 6GHz concrete loss
From Table 7.4.3.1 of 3GPP TR38.900
Study on channel model for frequency spectrum above 6 GHz
"""
return(5+4*fGHz)
[docs]def Lglass(fGHz):
""" 3GPP above 6GHz glass loss
From Table 7.4.3.1 of 3GPP TR38.900
Study on channel model for frequency spectrum above 6 GHz
"""
return(2 + 0.2* fGHz)
[docs]def LIRRglass(fGHz):
""" 3GPP above 6GHz IRR (Infra Red Reflecting) glass loss
From Table 7.4.3.1 of 3GPP TR38.900
Study on channel model for frequency spectrum above 6 GHz
"""
return(23 + 0.3* fGHz)
[docs]def Lwood(fGHz):
""" 3GPP above 6GHz wood loss
From Table 7.4.3.1 of 3GPP TR38.900
Study on channel model for frequency spectrum above 6 GHz
"""
return(4.85 + 0.12* fGHz)
[docs]def LossPenetration(fGHz, alpha = 0.3, typ='low'):
if typ=='low':
PLTW = 5 - 10*np.log10(alpha*10**(-Lglass(fGHz)/10)+(1-alpha)*10**(-Lconcrete(fGHz)))
if typ=='high':
PLTW = 5 - 10*np.log10((1-alpha)*10**(-LIRRglass(fGHz)/10)+alpha*10**(-Lconcrete(fGHz)))
return(PLTW)
[docs]def Losst(L,fGHz,p1,p2,dB=True,bceilfloor=False):
""" calculate Losses between links p1-p2
Parameters
----------
L : Layout object
fGHz : np.array
frequency GHz
p1 : source points
(3 x Np1) array or (3,) array
p2 : observation points
(3 x Np2) array or (3,) array
dB : boolean
bceilfloor : boolean
Examples
--------
.. plot::
:include-source:
>>> import matplotlib.pyplot as plt
>>> from pylayers.measures.mesuwb import *
>>> from pylayers.antprop.loss import *
>>> S = Simul()
>>> S.layout('WHERE1.lay')
>>> fGHz = 4
>>> Tx,Rx = ptw1()
>>> Lwo,Lwp,Edo,Edp = Losst(S.L,fGHz,Tx.T,Rx[1,0:2],dB=True)
>>> fig=plt.figure(figsize=(20,10))
>>> fig,ax = S.L.showGs(fig=fig)
>>> tit = plt.title('test Losst')
>>> sc2 = ax.scatter(Rx[1,0],Rx[1,1],s=20,marker='x',c='k')
>>> sc1 = ax.scatter(Tx[:,0],Tx[:,1],s=20,c=Lwo,linewidth=0)
>>> cb = plt.colorbar(sc1)
>>> cb.set_label('dB')
>>> plt.show()
See Also
--------
pylayers.antprop.coverage
pylayers.slab.Interface.losst
"""
if (type(fGHz)==float) | (type(fGHz)==int):
fGHz=np.array([fGHz],dtype=float)
sh1 = np.shape(p1)
sh2 = np.shape(p2)
if (len(sh1)>1) & (len(sh2)>1):
Nlink = max(sh1[1],sh2[1])
if (len(sh1)>1) & (len(sh2)<2):
Nlink = sh1[1]
if (len(sh1)<2) & (len(sh2)>1):
Nlink = sh2[1]
if (len(sh1)<2) & (len(sh2)<2):
Nlink = 1
# determine incidence angles on segment crossing p1-p2 segment
#data = L.angleonlink(p1,p2)
logger.debug('losst before angleonlink3')
data = L.angleonlink3(p1,p2)
# as many slabs as segments and subsegments
us = data['s']
slabs = np.array([ L.Gs.node[x]['name'] for x in us ])
#slabs = L.sla[us]
check = np.where(slabs=='')
#
# As segment numbering is not necessarily contiguous
# there exist void string '' in slabs
cslab = list(np.unique(slabs))
if '' in cslab:
cslab.remove('')
if 'AIR' in cslab:
cslab.remove('AIR')
if '_AIR' in cslab:
cslab.remove('_AIR')
LossWallo = np.zeros((len(fGHz),Nlink))
LossWallp = np.zeros((len(fGHz),Nlink))
EdWallo = np.zeros((len(fGHz),Nlink))
EdWallp = np.zeros((len(fGHz),Nlink))
for slname in cslab:
# u index of slabs of name slname
# data['a'][u] angle
# data['s'][u] segment number including subsegment
u = np.nonzero(np.array(slabs)==slname)[0]
#
# calculate Loss for slab slname
#
lko,lkp = L.sl[slname].losst(fGHz,data['a'][u])
#
# calculate Excess delay for slab slname
#
do , dp = L.sl[slname].excess_grdelay(theta=data['a'][u])
# data['i'][u] links number
indexu = data['i'][u]
# reduce to involved links
involved_links, indices = np.unique(indexu,return_index=True)
indicep = np.hstack((indices[1:],np.array([len(indexu)])))
# range on involved links
irange = np.arange(len(involved_links))
#
# sum contribution of slab of a same link
#
Wallo = np.array([ np.sum(lko[:,indices[x]:indicep[x]],axis=1) for x in irange ] ).T
Wallp = np.array([ np.sum(lkp[:,indices[x]:indicep[x]],axis=1) for x in irange ] ).T
Edo = np.array([np.sum(do[indices[x]:indicep[x]]) for x in irange]).T
Edp = np.array([np.sum(dp[indices[x]:indicep[x]]) for x in irange]).T
LossWallo[:,involved_links] = LossWallo[:,involved_links] + Wallo
LossWallp[:,involved_links] = LossWallp[:,involved_links] + Wallp
EdWallo[:,involved_links] = EdWallo[:,involved_links] + Edo
EdWallp[:,involved_links] = EdWallp[:,involved_links] + Edp
if bceilfloor:
# Managing Ceil / Floor transmission
# At that point there is only a single type of ceil and floor
# it should be defined ideally as a specific entity
#
# TODO fix shape error p2 is not always (3 x N)
#
if (p1[2,:]> L.zceil).any() or (p2[2,:]> L.zceil).any():
# WARNING : this test should be done individually
if (p1[2]>p2[2]).all():
v0 = p1
v1 = p2
else:
v0 = p2
v1 = p1
uu = v0 - v1
# 1 x N
nu = np.sqrt(np.sum(uu * uu, axis=0))
# 3 x N
un = uu / nu[np.newaxis, :]
dotp = np.einsum('ij,i->j',un,np.array([0,0,1]))
alphas = np.arccos(dotp)
#
# calculate Loss for slab CEIL
#
lkco,lkcp = L.sl['CEIL'].losst(fGHz,alphas)
#
# calculate Excess delay for slab CEIL
#
dco , dcp = L.sl['CEIL'].excess_grdelay(theta=alphas)
LossWallo = LossWallo + lkco
LossWallp = LossWallp + lkcp
EdWallo = EdWallo + dco
EdWallp = EdWallp + dcp
# check crossing floor
if (p1[2,:]< L.zfloor).any() or (p2[2,:]< L.zfloor).any():
# WARNING : this test should be done individually
if (p1[2]>p2[2]).all():
v0 = p1
v1 = p2
else:
v0 = p2
v1 = p1
uu = v0 - v1
# 1 x N
nu = np.sqrt(np.sum(uu * uu, axis=0))
# 3 x N
un = uu / nu[np.newaxis, :]
dotp = np.einsum('ij,i->j',un,np.array([0,0,1]))
alphas = np.arccos(dotp)
#
# calculate Loss for slab CEIL
#
lkfo,lkfp = L.sl['FLOOR'].losst(fGHz,alphas)
#
# calculate Excess delay for slab CEIL
#
dfo , dfp = L.sl['FLOOR'].excess_grdelay(theta=alphas)
LossWallo = LossWallo + lkfo
LossWallp = LossWallp + lkfp
EdWallo = EdWallo + dfo
EdWallp = EdWallp + dfp
if not dB:
LossWallo = 10**(-LossWallo/10)
LossWallp = 10**(-LossWallp/10)
print(LossWallo)
return(LossWallo,LossWallp,EdWallo,EdWallp)
[docs]def gaspl(d,fGHz,T,PhPa,wvden):
""" attenuation due to atmospheric gases
Parameters
----------
d : np.array
range (meters)
fGHz : np.array
frequency (GHz)
T : float
Temprature in degree Celcius
PhPa : float
Pressure in hPa
wvden : float
Water vapor density (g/m**3)
Examples
--------
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> T = 15
>>> PhPa = 1013
>>> wvden = 7.5
>>> d = 1000
>>> fGHz = np.linspace(1,1000,100)
>>> L = gaspl(d,fGHz,T,PhPa,wvden)
>>> plt.plot(fGHz,L)
Notes
-----
This function implements the recommandation UIT-P676-10
"""
affO2=np.array([
[50.474214,0.975,9.651,6.690,0.0,2.566,6.850],
[50.987745,2.529,8.653,7.170,0.0,2.246,6.800],
[51.503360,6.193,7.709,7.640,0.0,1.947,6.729],
[52.021429,14.320,6.819,8.110,0.0,1.667,6.640],
[52.542418,31.240,5.983,8.580,0.0,1.388,6.526],
[53.066934,64.290,5.201,9.060,0.0,1.349,6.206],
[53.595775,124.600,4.474,9.550,0.0,2.227,5.085],
[54.130025,227.300,3.800,9.960,0.0,3.170,3.750],
[54.671180,389.700,3.182,10.370,0.0,3.558,2.654],
[55.221384,627.100,2.618,10.890,0.0,2.560,2.952],
[55.783815,945.300,2.109,11.340,0.0,-1.172,6.135],
[56.264774,543.400,0.014,17.030,0.0,3.525,-0.978],
[56.363399,1331.800,1.654,11.890,0.0,-2.378,6.547],
[56.968211,1746.600,1.255,12.230,0.0,-3.545,6.451],
[57.612486,2120.100,0.910,12.620,0.0,-5.416,6.056],
[58.323877,2363.700,0.621,12.950,0.0,-1.932,0.436],
[58.446588,1442.100,0.083,14.910,0.0,6.768,-1.273],
[59.164204,2379.900,0.387,13.530,0.0,-6.561,2.309],
[59.590983,2090.700,0.207,14.080,0.0,6.957,-0.776],
[60.306056,2103.400,0.207,14.150,0.0,-6.395,0.699],
[60.434778,2438.000,0.386,13.390,0.0,6.342,-2.825],
[61.150562,2479.500,0.621,12.920,0.0,1.014,-0.584],
[61.800158,2275.900,0.910,12.630,0.0,5.014,-6.619],
[62.411220,1915.400,1.255,12.170,0.0,3.029,-6.759],
[62.486253,1503.000,0.083,15.130,0.0,-4.499,0.844],
[62.997984,1490.200,1.654,11.740,0.0,1.856,-6.675],
[63.568526,1078.000,2.108,11.340,0.0,0.658,-6.139],
[64.127775,728.700,2.617,10.880,0.0,-3.036,-2.895],
[64.678910,461.300,3.181,10.380,0.0,-3.968,-2.590],
[65.224078,274.000,3.800,9.960,0.0,-3.528,-3.680],
[65.764779,153.000,4.473,9.550,0.0,-2.548,-5.002],
[66.302096,80.400,5.200,9.060,0.0,-1.660,-6.091],
[66.836834,39.800,5.982,8.580,0.0,-1.680,-6.393],
[67.369601,18.560,6.818,8.110,0.0,-1.956,-6.475],
[67.900868,8.172,7.708,7.640,0.0,-2.216,-6.545],
[68.431006,3.397,8.652,7.170,0.0,-2.492,-6.600],
[68.960312,1.334,9.650,6.690,0.0,-2.773,-6.650],
[118.750334,940.300,0.010,16.640,0.0,-0.439,0.079],
[368.498246,67.400,0.048,16.400,0.0,0.000,0.000],
[424.763020,637.700,0.044,16.400,0.0,0.000,0.000],
[487.249273,237.400,0.049,16.000,0.0,0.000,0.000],
[715.392902,98.100,0.145,16.000,0.0,0.000,0.000],
[773.839490,572.300,0.141,16.200,0.0,0.000,0.000],
[834.145546,183.100,0.145,14.700,0.0,0.000,0.000]])
## spectroscopic data for H20
## f0 b1 b2 b3 b4 b5 b6
affH2O=np.array([
[22.235080,0.1130,2.143,28.11,0.69,4.800,1.00],
[67.803960,0.0012,8.735,28.58,0.69,4.930,0.82],
[119.995940,0.0008,8.356,29.48,0.70,4.780,0.79],
[183.310091,2.4200,0.668,30.50,0.64,5.300,0.85,],
[321.225644,0.0483,6.181,23.03,0.67,4.690,0.54],
[325.152919,1.4990,1.540,27.83,0.68,4.850,0.74],
[336.222601,0.0011,9.829,26.93,0.69,4.740,0.61],
[380.197372,11.5200,1.048,28.73,0.54,5.380,0.89,],
[390.134508,0.0046,7.350,21.52,0.63,4.810,0.55],
[437.346667,0.0650,5.050,18.45,0.60,4.230,0.48,],
[439.150812,0.9218,3.596,21.00,0.63,4.290,0.52,],
[443.018295,0.1976,5.050,18.60,0.60,4.230,0.50],
[448.001075,10.3200,1.405,26.32,0.66,4.840,0.67],
[470.888947,0.3297,3.599,21.52,0.66,4.570,0.65,],
[474.689127,1.2620,2.381,23.55,0.65,4.650,0.64],
[488.491133,0.2520,2.853,26.02,0.69,5.040,0.72],
[503.568532,0.0390,6.733,16.12,0.61,3.980,0.43],
[504.482692,0.0130,6.733,16.12,0.61,4.010,0.45],
[547.676440,9.7010,0.114,26.00,0.70,4.500,1.00],
[552.020960,14.7700,0.114,26.00,0.70,4.500,1.00],
[556.936002,487.4000,0.159,32.10,0.69,4.110,1.00],
[620.700807,5.0120,2.200,24.38,0.71,4.680,0.68],
[645.866155,0.0713,8.580,18.00,0.60,4.000,0.50],
[658.005280,0.3022,7.820,32.10,0.69,4.140,1.00],
[752.033227,239.6000,0.396,30.60,0.68,4.090,0.84],
[841.053973,0.0140,8.180,15.90,0.33,5.760,0.45],
[859.962313,0.1472,7.989,30.60,0.68,4.090,0.84],
[899.306675,0.0605,7.917,29.85,0.68,4.530,0.90],
[902.616173,0.0426,8.432,28.65,0.70,5.100,0.95],
[906.207325,0.1876,5.111,24.08,0.70,4.700,0.53],
[916.171582,8.3400,1.442,26.70,0.70,4.780,0.78],
[923.118427,0.0869,10.220,29.00,0.70,5.000,0.80],
[970.315022,8.9720,1.920,25.50,0.64,4.940,0.67],
[987.926764,132.1000,0.258,29.85,0.68,4.550,0.90],
[1780.000000,22,300.0000,0.952,176.20,0.50,30.500]])
dkm = d/1000.
TK = T + 273.15
theta = 300./TK #3
fO2 = affO2[:,0]
a1 = affO2[:,1]
a2 = affO2[:,2]
a3 = affO2[:,3]
a4 = affO2[:,4]
a5 = affO2[:,5]
a6 = affO2[:,6]
fH2O = affH2O[:,0]
b1 = affH2O[:,1]
b2 = affH2O[:,2]
b3 = affH2O[:,3]
b4 = affH2O[:,4]
b5 = affH2O[:,5]
b6 = affH2O[:,6]
e = wvden*TK/216.7 # 4
SO2 = a1*1e-7*PhPa*(theta**3)*np.exp(a2*(1-theta)) # 3
DO2 = a3*1e-4*(PhPa*(theta**(0.8-a4))+1.1*e*theta) # 6a
SH2O = b1*1e-1*e*(theta**(3.5))*np.exp(b2*(1-theta)) # 3
DH2O = b3*1e-4*(PhPa*theta**b4+b5*e*theta**b6) # 6a
DO2_m = np.sqrt(DO2**2+2.25e-6) # 6b
DH2O_m = 0.535*DH2O+np.sqrt(0.217*DH2O**2+(2.1316*1e-12*fH2O**2)/theta)
deltaO2 = (a5+a6*theta)*1e-4*(PhPa+e)*theta**(0.8)
#
# O2
#
uO2 = fO2[:,None]-fGHz[None,:]
vO2 = fO2[:,None]+fGHz[None,:]
n1O2 = DO2_m[:,None]-deltaO2[:,None]*uO2
n2O2 = DO2_m[:,None]-deltaO2[:,None]*vO2
d1O2 = uO2**2 + DO2_m[:,None]**2
d2O2 = vO2**2 + DO2_m[:,None]**2
FO2 = (fGHz[None,:]/fO2[:,None])*(n1O2/d1O2+n2O2/d2O2)
UO2 = SO2[:,None]*FO2
#
# H2O
#
uH2O = fH2O[:,None]-fGHz[None,:]
vH2O = fH2O[:,None]+fGHz[None,:]
nH2O = DH2O_m[:,None]
d1H2O = uH2O**2 + DH2O_m[:,None]**2
d2H2O = vH2O**2 + DH2O_m[:,None]**2
FH2O = (fGHz[None,:]/fH2O[:,None])*(nH2O/d1H2O+nH2O/d2H2O)
UH2O = SH2O[:,None]*FH2O
# Nsec (8)
dD = 5.6e-4*(PhPa+e)*theta**(0.8)
t1 = 6.14e-5/(dD*(1.+(fGHz/dD)**2))
t2 = 1.4e-12*PhPa*(theta**(1.5))/(1+1.9e-5*fGHz**(1.5))
Nsec = fGHz*PhPa*(theta**2)*(t1+t2) # 9
ulow = np.where(fGHz<118.750343)[0]
uhigh = np.where(fGHz>=118.750343)[0]
UO2low = UO2[:,ulow]
UO2high = UO2[:,uhigh]
SO2low = np.sum(UO2low,axis=0)
SO2high = np.sum(UO2high[38:,:],axis=0)
sSO2 = np.hstack((SO2low,SO2high))
Npp = sSO2 + np.sum(UH2O,axis=0)+Nsec
Npp = np.sum(UO2,axis=0) + np.sum(UH2O,axis=0)+Nsec
gamma = 0.1820*fGHz*Npp
LgasdB = gamma*dkm
return(LgasdB)
[docs]def Loss0(S,rx,ry,f,p):
""" calculate Loss through Layers for theta=0 deg
Parameters
----------
S : Simulation object
rx : extremity of link
ry : extremity of link
fGHz : float
frequency GHz
p :
"""
Nx = len(rx)
Ny = len(ry)
Lw = np.zeros((Nx,Ny))
print(shape(Lw))
i = 0
for x in rx:
j = 0
for y in ry:
Loss = 0
pxy = np.array([x,y])
seglist,theta = L.angleonlinkold(p,pxy)
for k in seglist:
name = L.name[k]
lk = L.sl[name].loss0(f)
Loss = Loss + lk[0]
Lw[i,j] = Loss
j = j+1
i = i+1
return(Lw)
[docs]def Loss_diff(u):
""" calculate Path Loss of the diffraction
"""
if u < -0.7:
Ld = 0
elif u > 1.5:
Ld = 13 + 20*np.log10(u)
else:
Ld = 6.9 + 20*np.log10(np.sqrt((u-0.1)**2+1)+u-0.1)
return(Ld)
[docs]def calnu(h,d1,d2,fGHz):
r""" Calculate the diffraction Fresnel parameter
Parameters
----------
h : signed height w.r.t LOS (meter)
d1 : distance 1 (meter)
d2 : distance 2 (meter)
fGHz : frequency GHz
Notes
-----
.. math:: \nu = h \sqrt{\frac{2}{\lambda} \frac{d_1+d_2}{d_1 d_2}}
"""
ld = 0.3/fGHz
nu = h*np.sqrt(2*(d1+d2)/(ld*d1*d2))
return(nu)
[docs]def route(X, Y, Z, Ha, Hb, fGHz, K, method='deygout'):
""" diffraction loss along a route
Parameters
----------
X : np.array (Nphi,Nr)
cartesian coordinate grid
Y : np.array (Nphi,Nr)
cartesian coordinate grid
Z : np.array (Nphi,Nr)
height (meters)
Ha : float
Hb : float
fGHz : np.array (,Nf)
frequency in GHz
method : 'deygout' | 'bullington'
Returns
-------
L : Losses (dB)
"""
Nphi, Nr = Z.shape
if (type(fGHz) == float):
fGHz = np.array([fGHz])
Nf = len(fGHz)
L = np.zeros((Nphi, Nf))
L0 = np.zeros(Nf)
# loop over azimut
for ip in range(Nphi):
x = X[ip, :]
y = Y[ip, :]
z = Z[ip, :]
d = np.sqrt((x-x[0])**2+(y-y[0])**2)
# effect of refraction in equivalent earth curvature
dh = d*(d[::-1])/(2*K*6375e3)
z = z + dh
LOS = 32.4 + 20*np.log10(fGHz) + 20*np.log10(d[-1])
z[0] = z[0] + Ha
z[-1] = z[-1] + Hb
if method == 'deygout':
LDiff = deygout(d, z, fGHz, L0, 0)
if method == 'bullington':
LDiff, deq, heq = bullington(d, z, fGHz)
L[ip, :] = LDiff+LOS
return(L)
[docs]def cover(X, Y, Z, Ha, Hb, fGHz, K, method='deygout'):
""" outdoor coverage on a region
Parameters
----------
X : np.array (Nphi,Nr)
cartesian coordinate grid
Y : np.array (Nphi,Nr)
cartesian coordinate grid
Z : np.array (Nphi,Nr)
height (meters)
Ha : float
Hb : float
fGHz : np.array (,Nf)
frequency in GHz
method : 'deygout' | 'bullington'
Returns
-------
L : Losses (dB)
"""
Nphi, Nr = Z.shape
if (type(fGHz) == float):
fGHz = np.array([fGHz])
Nf = len(fGHz)
L = np.zeros((Nphi, Nr, Nf))
L0 = np.zeros(Nf)
# loop over azimut
for ip in range(Nphi):
# loop over range
# il : 2 ... Nr-2
# uk : 0 ....Nr-1
for il in np.arange(2, Nr-1):
uk = np.arange(0, il+1)
z = np.empty(len(uk))
x = X[ip, uk]
y = Y[ip, uk]
z[uk] = Z[ip, uk]
d = np.sqrt((x-x[0])**2+(y-y[0])**2)
# effect of refraction in equivalent earth curvature
dh = d*(d[::-1])/(2*K*6375e3)
z = z + dh
LOS = 32.4 + 20*np.log10(fGHz) + 20*np.log10(d[-1])
z[0] = z[0] + Ha
z[-1] = z[-1] + Hb
if method == 'deygout':
LDiff = deygout(d, z, fGHz, L0, 0)
if method == 'bullington':
LDiff, deq, heq = bullington(d, z, fGHz)
L[ip, il, :] = LDiff[None, :]+LOS[None,:]
return(L)
[docs]def deygout(d, height, fGHz, L, depth):
""" Deygout attenuation
Parameters
----------
d : np.array (,N)
horizontal distance
height : np.array (,N)
height profile
fGHz : np.array (,Nf)
frequency GHz
L : np.array (,Nf)
Additional Loss
depth : recursive depth
Notes
-----
This function is recursive
"""
lmbda = 0.3/fGHz
L0 = np.zeros(len(fGHz))
depth = depth+1
N = len(height)
if depth < 3:
if N > 3:
u = np.arange(N)/(N-1.0) # float
# l : straight line between termination (LOS)
l = (height[0])*(1-u)+(height[-1])*u
# h excludes termination points
h = height[1:-1] - l[1:-1]
# Fresnel parameter (engagement)
nu = h[:, None] * np.sqrt((2/lmbda[None, :]) *
(1/d[1:-1, None]+1/(d[-1]-d[1:-1, None])))
imax = np.unique(np.nanargmax(nu, axis=0))[0]
numax = nu[imax, :]
else:
numax = -10*np.ones(len(fGHz))
if (numax > -0.78).any():
w = numax - 0.1
L = L + np.maximum(6.9 + 20*np.log10(np.sqrt(w**2+1)+w), 0)
# left link
height1 = height[0:imax+2]
d1 = d[0:imax+2]
Ll = deygout(d1, height1, fGHz, L0, depth)
# right link
height2 = height[imax+1:]
d2 = d[imax+1:]
Lr = deygout(d2, height2, fGHz, L0, depth)
# add losses
L = L + Lr + Ll
return(L)
[docs]def bullington(d, height, fGHz):
""" edges attenuation with Bullington method
Parameters
----------
d : np.array
height : np.array
antenna height is includes in height[0] and height[-1]
fGHz : np.array
Returns
-------
L : np.array
total loss
"""
def recl(d, height):
""" determine left interception point
Parameters
----------
d : np.array
height : np.array
"""
N = len(height)
u = np.arange(N)/(N-1.)
# l : straight line between termination (LOS)
l = height[0]*(1-u)+(height[-1])*u
h = height - l
# imax : index of the maximum height offset
imax = np.argmax(h)
if imax>0:
# hmax : maximum height offset
hmax = h[imax]
# parameterization from 0 to imax
ul = np.arange(imax)/(imax-1.)
# straight line
dhl = h[0]*(1-ul) + hmax*ul
# el : offset if <0 split again
el = dhl - h[0:imax]
if np.min(el) < 0:
u, v = recl(d[0:imax+1], height[0:imax+1])
else:
u = d[0:imax+1]
v = h[0:imax+1]
else:
u = d[0:1]
v = d[0:1]
return(u, v)
#if min(er)<0:
# u,v = rec(d[imax-1:-1],dhl)
#else:
#er = dhr - h[imax-1:-1]
def recr(d, height):
""" determine the right interception point
"""
N = len(height)
u = np.arange(N)/(N-1.)
l = height[0]*(1-u)+(height[-1])*u
h = height - l
imax = np.argmax(h)
hmax = h[imax]
ur = np.arange(N-imax)/(N-imax-1.)
dhr = hmax*(1-ur) + h[-1]*ur
er = dhr - h[imax:]
if np.min(er) < 0:
u, v = recr(d[imax:],h[imax:])
else:
u = d[imax:]
v = h[imax:]
return(u,v)
#if min(er)<0:
# u,v = rec(d[imax-1:-1],dhl)
#else:
#er = dhr - h[imax-1:-1]
lmbda = 0.3/fGHz
u = np.arange(len(height))/(len(height)-1.)
l = (height[0])*(1-u)+(height[-1])*u
h = height - l
if (h>0).any():
ul, vl = recl(d, height)
ur, vr = recr(d, height)
idtx = len(ul)
idrx = len(h) - len(ur)
dtx = d[idtx]
drx = d[-1]-d[idrx]
htx = h[idtx-1]
hrx = h[idrx]
deq = (dtx*hrx)*d[-1]/(drx*htx+dtx*hrx)
heq = deq*(htx/dtx)
else:
heq = -np.min(np.abs(h[1:-1]))
ieq = np.where(h==heq)[0][0]
deq = d[ieq]
nu = heq*np.sqrt((2/lmbda)*(1/deq+1/(d[-1]-deq)))
w = nu - 0.1
L = np.maximum(6.9 + 20*np.log10(np.sqrt(w**2+1)+w), 0)
return(L, deq, heq)
[docs]def two_rays_flatearth(fGHz, **kwargs):
"""
Parameters
----------
p0 : transmitter position
(3 x Np1) array or (2,) array
p1 : receiver position
(3 x Np2) array or (2,) array
OR :
d : distance between Tx and Rx
(Np1,)
ht : Tx height
hr : Rx height
(Np1)
GtdB : float (0)
Transmitter Antenna Gain (dB)
GrdB : float(0)
Receiver Antenna Gain (dB)
fGHz : float (2.4)
frequency (GHz)
gamma : complex (-1.+0.j)
Reflexion coeff
dB : boolean (True)
return result in d
Returns
-------
P :
received power
Examples
--------
.. plot::
:include-source:
>>> from pylayers.antprop.loss import *
>>> NPT=10000
>>> x=np.array([0,0,8])
>>> x=x.reshape(3,1)
>>> y = np.ones((3,NPT))
>>> y[0,:]=0
>>> y[1,:]=np.arange(NPT)
>>> y[2,:]=2
>>> g0=1
>>> g1=1
>>> fGHz=2.4
>>> PL2R=two_rays_flatearth(p0=x,p1=y,fGHz=fGHz,GtdB=g0,GrdB=g1)
>>> PL1R = PL(fGHz,x,y,2)
>>> plt.semilogx(PL2R,label='two-ray model')
>>> plt.semilogx(-PL1R[0,:],label='one slope model')
>>> plt.axis([10,NPT,-150,-50])
>>> plt.title('Loss 2-rays model vs one slope model')
>>> plt.xlabel('distance (m)')
>>> plt.ylabel('Loss Pr/Pt (dB)')
>>> plt.legend()
>>> plt.show()
>>> d=np.arange(1,1000)
>>> PL2Rd = two_rays_flatearth(d=d,ht=np.array([5]),hr=np.array([10]),fGHz=fGHz,GtdB=g0,GrdB=g1)
>>> plt.semilogx(PL2Rd,label='two-ray model')
>>> plt.semilogx(-PL1R[0,:],label='one slope model')
>>> plt.axis([10,NPT,-150,-50])
>>> plt.title('Loss 2-rays model vs one slope model')
>>> plt.xlabel('distance (m)')
>>> plt.ylabel('Loss Pr/Pt (dB)')
>>> plt.legend()
>>> plt.show()
References
----------
https://en.wikipedia.org/wiki/Two-ray_ground-reflection_model#As_a_case_of_log_distance_path_loss_model
http://morse.colorado.edu/~tlen5510/text/classwebch3.html#x15-590003.3.3
"""
defaults = { 'p0':np.array((0,0,10)),
'p1':np.array((0,10,10)),
'd':[],
'ht':10,
'hr':10,
'GtdB':0.,
'GrdB':0.,
'gamma': -1.+0.j,
'pol':'v',
'eps' :[],
'sig':0.,
'dB':True
}
for k in defaults:
if k not in kwargs:
kwargs[k]=defaults[k]
GtdB=kwargs.pop('GtdB')
GrdB=kwargs.pop('GrdB')
Gt = 10**((1.*GtdB)/10.)
Gr = 10**((1.*GrdB)/10.)
gamma=kwargs.pop('gamma')
pol=kwargs.pop('pol')
eps=kwargs.pop('eps')
sig=kwargs.pop('sig')
if kwargs['d'] == []:
p0=kwargs['p0']
p1=kwargs['p1']
assert p0.shape[0] == 3, 'p0 is not 3D'
assert p1.shape[0] == 3, 'p1 is not 3D'
if len(p0.shape) == 1:
p0=p0.reshape(p0.shape[0],1)
if len(p1.shape) == 1:
p1=p1.reshape(p1.shape[0],1)
p0=1.*p0
p1=1.*p1
ht = p0[2,:]
hr = p1[2,:]
dloss = np.sqrt(np.sum((p0-p1)**2,axis=0)) #l0
else:
dloss=kwargs['d']
ht=kwargs['ht']
hr=kwargs['hr']
Gt = 10**((1.*Gt)/10.)
Gr = 10**((1.*Gr)/10.)
d0 = np.sqrt( dloss**2 - 1.*(ht-hr)**2 ) # d0
dref = np.sqrt(d0**2+1.*(ht+hr)**2) #l0'
if eps != []:
psy = np.arcsin((ht+hr)/dref)
er = eps - 60.j*sig*0.3/fGHz
if pol == 'v':
Z = (1./er)* np.sqrt(er-np.cos(psy)**2)
elif pol == 'h':
Z = np.sqrt(er-np.cos(psy)**2)
gamma = (np.sin(psy)-Z)/((np.sin(psy)+Z))
deltad = dref-dloss
deltaphi = (2*np.pi*fGHz*deltad)/0.3
E= (0.3/(4*np.pi*fGHz)) * (np.sqrt(Gt*Gr)/dloss + gamma * np.sqrt(Gr*Gr)*(np.exp(-1.j*deltaphi))/dref)
P = abs(E)**2
# import ipdb
# ipdb.set_trace()
if kwargs['dB'] :
return 10*np.log10(P)
else:
return P
[docs]def lossref_compute(P,h0,h1,k=4/3.) :
"""
compute loss and reflection rays on curved earth
Parameters
----------
P : float |list
if len(P) == 1 => P is a distance
if len(P) == 4 => P is a list of [lon0,lat0,lon1,lat1]
where :
lat0 : float |string
latitude first point (decimal |deg min sec Direction)
lat1 : float |string
latitude second point (decimal |deg min sec Direction)
lon0 : float |string
longitude first point (decimal |deg min sec Direction)
lon1 : float |string
longitude second point (decimal |deg min sec Direction)
h0 : float:
height of 1st point
h1 : float:
height of 2nd point
k : electromagnetic earth factor
Returns
-------
dloss : float
length of direct path (meter)
dref : float
length of reflective path (meter)
psy : float
Reflection angle
References
----------
B. R. Mahafza, Radar systems analysis and design using MATLAB, Third edition. Boca Raton; London: CRC/Taylor & Francis, chapter 8, 2013.
"""
if isinstance(P,float) or isinstance(P,int) :
#P is a distance
r=P
mode = 'dist'
elif isinstance(P,np.ndarray) or isinstance(P,list):
if len(P) == 1:
#P is a distance
r=P
mode = 'dist'
elif len(P) == 4:
#P is a lonlat
lat0=P[0]
lon0=P[1]
lat1=P[2]
lon1=P[2]
mode = 'lonlat'
else :
raise AttributeError('P must be a list [lat0,lon0,lat1,lon0] or a distance')
else :
raise AttributeError('Invalid P format ( list |ndarray )')
# if h0<h1:
# h1,h0 = h0,h1
r0 = 6371.e3 # earth radius
re = k*r0 # telecom earth radius
if mode == 'lonlat':
# r = distance curvilignenp.arcsin((h1/R1)-R1/(2.*re)) entre TXetRX / geodesic
r = gu.distance_on_earth(lat0, lon0, lat1, lon1)
else :
r=P
r=1.*r
# import ipdb
# ipdb.set_trace()
p = 2/(np.sqrt(3))*np.sqrt(re*(h0+h1)+(r**2/4.)) #eq 8.45
eps = np.arcsin(2*re*r*(h1-h0)/p**3) # eq 8.46
#distance of reflection on curved earth
r1 = r/2 - p*np.sin(eps/3) #eq 8.44
r2 = r -r1
phi1 = r1/re #8.47
phi2 = r2/re # 8.48
R1 = np.sqrt(h0**2+4*re*(re+h0)*(np.sin(phi1/2))**2) # 8.51
R2 = np.sqrt(h1**2+4*re*(re+h1)*(np.sin(phi2/2))**2) #8.52
Rd = np.sqrt((h1-h0)**2+4*(re+h1)*(re+h0)*np.sin((phi1+phi2)/2.)**2) # 8.53
# tangente angle on earth
psy = np.arcsin((h1/R1)-R1/(2.*re)) #eq 8.55
deltaR = 4*R1*R2*np.sin(psy)**2/(R1+R2+Rd)
dloss = Rd
dref = R1+R2
return psy,dloss,dref
[docs]def two_rays_curvedearthold(P,h0,h1,fGHz=2.4,**kwargs):
"""
Parameters
----------
P : float |list
if len(P) == 1 => P is a distance
if len(P) == 4 => P is a list of [lon0,lat0,lon1,lat1]
where :
lat0 : float |string
latitude first point (decimal |deg min sec Direction)
lat1 : float |string
latitude second point (decimal |deg min sec Direction)
lon0 : float |string
longitude first point (decimal |deg min sec Direction)
lon1 : float |string
longitude second point (decimal |deg min sec Direction)
h0 : float:
height of 1st point
h1 : float:
height of 2nd point
fGHz : float
frequency (GHz)
k : float
electromagnetic earth factor
GtdB : float
Transmitter Antenna Gain (dB)
GrdB : float
Receiver Antenna Gain (dB)
gamma : complex (-1.+0.j)
Reflexion coeff if eps and sig are not precised
'pol': string ('v')
polarization ('v'|'h')
'eps' : float ([])
lossless relative permittivity [],
'sig': float (0.)
conductivity
dB : boolean (True)
return result in dB
Returns
-------
P :
received power
Examples
--------
.. plot::
:include-source:
>>> from pylayers.antprop.loss import *
>>> import matplotlib.pyplot as plt
>>> fGHz=2.4
>>> p0=np.array(([0,0,20]))
>>> p1=np.array(([0,1,20]))
>>> p0=p0.reshape(3,1)
>>> p1=p1.reshape(3,1)
>>> TRF = [] #Two Ray model on flat earth
>>> TRC = [] #Two Ray model on curved earth
>>> PLoss=[]
>>> for d in np.arange(1,10000,1):
>>> p1[1,:]=d
>>> TRF.append(two_rays_flatearth(p0[:,0],p1[:,0],fGHz,GtdB=0.,GrdB=0.,))
>>> TRC.append(two_rays_curvedearth(d,p0[2,:],p1[2,:],fGHz))
>>> PLoss.append(PL(fGHz, p0[:,0],p1[:,0], n=2.0, dB=True, d0=np.array([1])))
>>> PLoss=np.array(PLoss)[:,0,0]
>>> plt.semilogx(TRF,label='two-rays model flat earth')
>>> plt.semilogx(TRC,label='two-rays model curved earth')
>>> plt.semilogx(-PLoss,label='Path Loss')
>>> plt.legend()
>>> plt.show()
"""
defaults = { 'GtdB':0.,
'GrdB':0.,
'k':4/3.,
'gamma': -1.+0.j,
'pol':'v',
'eps' :[],
'sig':0.,
'mode':'PL',
'dB':True
}
for k in defaults:
if k not in kwargs:
kwargs[k]=defaults[k]
GtdB=kwargs.pop('GtdB')
GrdB=kwargs.pop('GrdB')
Gt = 10**((1.*GtdB)/10.)
Gr = 10**((1.*GrdB)/10.)
k=kwargs.pop('k')
gamma=kwargs.pop('gamma')
pol=kwargs.pop('pol')
eps=kwargs.pop('eps')
sig=kwargs.pop('sig')
h0=1.*h0
h1=1.*h1
psy,dloss,dref = lossref_compute(P,h0,h1,k)
if eps != []:
er = eps - 60.j*sig*0.3/fGHz
if pol == 'v':
Z= (1./er)* np.sqrt(er-np.cos(psy)**2)
elif pol == 'h':
Z= np.sqrt(er-np.cos(psy)**2)
gamma = (np.sin(psy)-Z)/((np.sin(psy)+Z))
deltad= dref-dloss
deltaphi = (2*np.pi*fGHz*deltad)/0.3
E= (0.3/(4*np.pi*fGHz) ) *(np.sqrt(Gt*Gr)/dloss + gamma * np.sqrt(Gr*Gr)*(np.exp(-1.j*deltaphi))/dref)
P = abs(E)**2
# import ipdb
# ipdb.set_trace()
if kwargs['dB'] :
return 10*np.log10(P)
else:
return P
[docs]def visuPts(S,nu,nd,Pts,Values,fig=[],sp=[],vmin=0,vmax=-1,label=' ',tit='',size=25,colbar=True,xticks=False):
""" visuPt : Visualization of values a given points
Parameters
----------
S : Simul
nu : useful Points
nd : Points deleted
Pts : Points coordinates
Value
"""
vx = Pts[nu,0]
vy = Pts[nu,1]
vxd = Pts[nd,0]
vyd = Pts[nd,1]
if vmax<vmin:
#vmin = min(Values)
vmin = -150
vmax = max(Values)
S.L.showGs()
if xticks:
for loc, spine in sp.spines.iteritems():
if loc in ['left','bottom']:
spine.set_position(('outward',10)) # outward by 10 points
sp.yaxis.set_ticks_position('left')
sp.xaxis.set_ticks_position('bottom')
elif loc in ['right','top']:
spine.set_color('none') # don't draw spine
else:
raise ValueError('unknown spine location: %s'%loc)
else:
for loc, spine in sp.spines.iteritems():
if loc in ['left']:
spine.set_position(('outward',10)) # outward by 10 points
sp.yaxis.set_ticks_position('left')
elif loc in ['right','top','bottom']:
spine.set_color('none') # don't draw spine
sp.xaxis.set_ticks([])
else:
raise ValueError('unknown spine location: %s'%loc)
# no xaxis ticks
#ax.xaxis.set_ticks([])
#sp.spines['left'].set_position('center')
#sp.spines['left'].set_color('none')
#sp.spines['right'].set_position('center')
#sp.spines['right'].set_color('none')
#sp.spines['bottom'].set_position('center')
#sp.xaxis.set_ticks_position('bottom')
#sp.yaxis.set_ticks_position('left')
#sp.spines['bottom'].set_color('none')
#sp.spines['top'].set_position('center')
#sp.spines['top'].set_color('none')
#
# Rooms annotation
#
annotate('R 8',xy=(-19,14.1))
annotate('R 9',xy=(-24.5,6.5))
annotate('R 14',xy=(-20,6.5))
annotate('R 10',xy=(-16.5,6.5))
annotate('R 7',xy=(-10.5,14.1))
annotate('R 11',xy=(-2.5,13.5))
annotate('R 12',xy=(-8.7,6.5))
annotate('R 13',xy=(-5.2,14.5))
annotate('R 1',xy=(3.5,8))
annotate('R 2',xy=(1.5,13.8))
annotate('R 6',xy=(-3.6,6.5))
n=scatter(vx,vy,marker='o',c=Values,s=size,vmin=vmin,vmax=vmax)
n.set_edgecolor('face')
# m=scatter(vxd,vyd,marker='o',c='k',s=size)
# m.set_edgecolor('face')
axis('scaled')
title(tit)
ylabel('meters')
if xticks:
xlabel('meters')
if colbar:
cbar=colorbar(orientation='vertical')
cbar.set_label(label)
[docs]def cdf(x,colsym="",lab="",lw=4):
""" plot the cumulative density function
Parameters
----------
x : np.array()
colsym : string
lab : string
lw : int
linewidth
Examples
--------
.. plot::
:include-source:
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.random.randn(10000)
>>> cdf(x)
>>> plt.show()
"""
rcParams['legend.fontsize']=20
rcParams['font.size']=20
x = np.sort(x)
n = len(x)
x2 = np.repeat(x, 2)
y2 = np.hstack([0.0, repeat(np.arange(1,n) / float(n), 2), 1.0])
plt.plot(x2,y2,colsym,label=lab,linewidth=lw)
plt.grid('on')
plt.legend(loc=2)
plt.xlabel('Ranging Error[m]')
plt.ylabel('Cumulative Probability')
if __name__=="__main__":
doctest.testmod()