"""
Simple demo of a scatter plot.
"""
import numpy as np
import math
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
from textwrap import wrap
tTable = 100
###average gap (seconds) between trams
vWalk = 0.1
###speed "walking"
newV = 0.5
###speed on tram
###Paris metro: 0.5km/minute
#######################################
###Set tTable=0 and newV=intinity for overall result
#######################################
tPerH=6
def eucDist(p1,p2):
	return ((p1[0]-p2[0])**2+(p1[1]-p2[1])**2)**0.5
#def ntDist(p1,p2,minic):
#	return eucDist(p1,p2) - 3.0/minic
	#minic trams per hour: so tram each 60/minic minutes
	#so average waiting is 60/(2*minic) minutes. Can
	#walk 0.1 km/minute. So distance is 60*(0.1)/(2*minic)
	# = 3/minic km.
def tProj(p,q1,q2):
	#projection of point p on to tram line q1---q2
	if p[0] < q1[0]:
		newP = q1
	elif p[0] > q2[0]:
		newP = q2
	else:
		newP = [p[0],q1[1]]
	return newP	

def tDist(p1,p2,e1,e2,relV):
	#pos = x coordinate where tram starts
	#freq = num. of trams per hour
	#relV = velocity relative to non-tram velocity
	np1 = tProj(p1,e1,e2)
	np2 = tProj(p2,e1,e2)
	a = eucDist(p1,np1)
	b = abs(np1[0]-np2[0])/relV
	c = eucDist(np2,p2)
	return(a+b+c)
def infeas(a,b):
	#return % of area that is infeasible
	#######################################
	#shape is array: 
	#circle: [centerX, centerY, radius]
	#rectangle: [bottomLeftX, bottomLeftY, width, height]
	########################################
	#tram: [leftX,leftY,rightX,rightY]
	pt1 = [a,b]
	rate = 2
	radius = 10
	end1=[15,20]
	end2=[25,20]
	totalC=0
	infeasC=0
	for i in range(10,31):
		dummy = int(round((100 - (i-20)**2)**0.5))
		#print dummy
		for j in range(20-dummy,20+dummy+1):
			totalC = totalC + 1
			arbP = [i,j]
			if eucDist(pt1,arbP) < tDist(pt1,arbP,end1,end2,newV/vWalk) + vWalk*tTable/2:
				infeasC = infeasC + 1
	perC = infeasC*100/totalC
	return (perC)
#################################################
u = range(40)[0:]
v = range(40)[0:]
W = [[0 for i in xrange(40)] for j in xrange(40)]
U, V = np.meshgrid(u, v)
count = 0
wSum = 0
for i in range(10,31):
	dummy = int(round((100 - (i-20)**2)**0.5))
	for j in range(20-dummy,21+dummy):
		temp9 = infeas(i,j)
		W[i][j] = temp9
		count = count + 1
		wSum = wSum + temp9
print wSum/count
fig = plt.figure()
ax = plt.axes(projection='3d')
pic = ax.plot_surface(U, V, W, rstride=1, cstride=1, linewidth=0, cmap=plt.cm.spectral_r, edgecolor='none',alpha=0.9)
title1 = ax.set_title("\n".join(wrap(str(wSum/count)+'% of journeys are infeasible. Minutes between trams:'+ str(tTable)+'. Non-tram speed (km/min): '+ str(vWalk)+'. Tram speed: '+str(newV),48)), fontweight = 'bold', fontsize = 12);
fig.colorbar(pic, shrink=0.9, aspect=2)
cset = ax.contour(U, V, W, zdir='z', offset=100, cmap=plt.cm.spectral_r)
plt.show()
#import numpy as np
#import matplotlib.pyplot as plt
#########################################
#############################################
contours = plt.contour(U, V, W, 10, colors='black')
plt.contourf(U, V, W, 10, cmap='gist_rainbow');
plt.clabel(contours, inline=True, fontsize=12)
plt.colorbar();


#fig = plt.figure()
#ax = plt.axes(projection='3d')
#ax.contour3D(X, Y, Z, 50, cmap='binary')
#ax.set_xlabel('x')
#ax.set_ylabel('y')
#ax.set_zlabel('z');
#####################################################
pt1=[20,20]
testP1=[17,13]
rate=2
fig, ax = plt.subplots()
#x = np.linspace(1, 40, 40)
end1=[15,20]
end2=[25,20]
dEnd=eucDist(pt1,end1)
print(dEnd)
#print(x)
for i in range(40)[1:]:
	#y=[i]*40
	#print(y)
	for j in range(40)[1:]:
		arbP = [i,j]
		#print arbP
		if eucDist(pt1,arbP) > tDist(pt1,arbP,end1,end2,10,3.0):
			txt='red'
		else:
			txt='green'
#		txt='orange'
		plt.plot([arbP[0]],[arbP[1]],color=txt,marker='o',linestyle='')
#	plt.plot(x,y,color='green',marker='s',linestyle=' ')
plt.plot([end1[0],end2[0]],[end1[1],end2[1]],color='red',linestyle='-',linewidth=7.0)
c=plt.Circle((20,20), radius= 10, alpha=0.5)
rect = plt.Rectangle((5, 15), 31, 10, color='k', alpha=0.5)
ax.add_artist(c)
ax.add_patch(rect)
plt.plot([pt1[0]],[pt1[1]],color='black',marker='s',linestyle='')
print('point: ',testP1)
print(eucDist(pt1,testP1))
print(tDist(pt1,testP1,end1,end2,10,3.0))
plt.show()