#! /usr/bin/python

###############################
# Author : Alexandre Mouradian
# year : 2013
###############################

import random
import sys
import math
import scipy.special
import numpy

import matplotlib.pyplot as plt
H=5
NR=1/math.sqrt(math.pi)

x=[]
y=[]
graph=[]
res=[]
Nx=[]
points=20

Np=1000
R=0.07
Pt=0.02*10**-3

class Node:
	ID=0
	rank=100
	coordinate=0
	mts=1
	x=0.0
	y=0.0

def distance(x1,y1,x2,y2):
	return math.sqrt((x1-x2)**2+(y1-y2)**2)

def produceGraph():
	graph.append(Node())
	graph[0].x=0.0
	graph[0].y=0.0
	x.append(0.0)
	y.append(0.0)
	graph[0].ID=0
	for i in range(1,N):
		graph.append(Node())
		xg=random.uniform(-NR,NR)
		yg=random.uniform(-NR,NR)
		while distance(0.0,0.0,xg,yg)>=NR:
			xg=random.uniform(-NR,NR)
			yg=random.uniform(-NR,NR)
		graph[i].x=xg
		graph[i].y=yg
		x.append(xg)
		y.append(yg)
		graph[i].ID=i
	graph.append(Node())
	print len(graph)
	graph[N].rank=0
	graph[N].ID=N
	graph[N].x=0.0
	graph[N].y=0.0

def distance(x1,y1,x2,y2):
	return math.sqrt((x1-x2)**2+(y1-y2)**2)

def setRank(R):
	for i in range(len(graph)):
		graph[i].rank=math.ceil(distance(0.0,0.0,graph[i].x,graph[i].y)/R)
		if math.ceil(distance(0.0,0.0,graph[i].x,graph[i].y)/R)==0 and graph[i].ID!=0 and graph[i].ID!=N:
			graph[i].rank=1

def prob(d):
	Gt=1.0
	Gr=1.0
	ple=2.0
	Nb=800.0
	L=1
	f=0.868*10**9
	c=3*10**8
	lambd=c/f
	N0=10**(-15.4)*10**-3
	B=0.5*10**6

	sigma=6

	K=Gt*Gr*lambd**2/(((4*math.pi)**2)*L)
	SNR=K*Pt*(d**-ple)/(N0*B)
	SNRdB=10*scipy.special.log10(SNR)

	threshold=scipy.special.erfcinv(1.0/Nb)**2

	thresholddB=10*math.log10(threshold)
	pr=0.5*(1-scipy.special.erf((thresholddB-SNRdB)/(sigma*math.sqrt(2))))
	return pr

def main():
	global N
	global graph
	for i in range(40):
		res.append(points*[0.0])
	for m in range(1,41):
		N=m*50
		Nx.append(N)
		for o in range(points):
			graph=[]
			print N, o
			NbH=0
			produceGraph()
			setRank(R)
			for i in range(len(graph)):
				if graph[i].rank==H:
					NbH+=1
			sender=0
			rcvd=0
			for i in range(Np):
				if NbH>0:
					indH=random.randint(1,NbH)
					nb=0
					for j in range(len(graph)):
						if graph[j].rank==H:
							nb+=1
							if nb==indH:
								sender=j
					oldSender=sender
					s=0
					while sender!=0 and sender!=-1 and sender!=N:
						s+=1
						for l in range(int(graph[sender].rank)):
							for k in range(len(graph)):
								if graph[k].rank==l:
									d=distance(graph[sender].x,graph[sender].y,graph[k].x,graph[k].y)
									d2=(graph[sender].rank-l)*R
									p2=random.random()
									p=prob(d*1000)
									oldSender=sender
									if p2<p:
										oldSender=sender
										sender=k
										break;
							if sender==k:
								break
						if oldSender==sender:
							sender=-1

					if sender==0 or sender==N:
						rcvd+=1
				else:
					print "no node at this H", N, H

			print "DR", rcvd/float(Np)
			res[m-1][o]=rcvd/float(Np)


	maxRes=[]
	minRes=[]
	avgRes=[]
	for i in range(40):
		maxRes.append(max(res[i]))
		minRes.append(min(res[i]))
		avgRes.append(numpy.mean(res[i]))
	legends=[]
	fig=plt.figure()
	ax=fig.add_subplot(111)
	leg1=ax.plot(Nx,maxRes,'*-')
	legends.append("maximum")
	leg2=ax.plot(Nx,minRes,'p-')
	legends.append("minimum")
	leg3=ax.plot(Nx,avgRes,'o-')
	legends.append("average")

	ax.legend((leg1,leg2,leg3),legends,'lower right')
	ax.set_xlabel('Number of nodes $N$', fontsize=15)
	ax.set_ylabel('$P_{e2e\_2b}$', fontsize=15)
	plt.show()


if __name__ == '__main__':
    main()