##### Slabljenje klinaste ovire - S53MV 18.03.2018 #####
##### Prirejeno za navodila vaj: Aljaz Blatnik 01.12.2018 #####

import math,cmath,numpy	#uporabljene funkcije

kx=0.001	#korak po zaslonu [ro1]
xmin=-5		#zacetek zaslona [ro1]
xmax=5		#konec zaslona [ro1]

def integral(xx):	#izracun integrala exp(-ju**2) v obmocju (0..xx)
	ku=0.01
	u=ku/2
	s=0
	while u<math.sqrt(math.pi)*xx:
		s=s+cmath.exp(complex(0,-u*u))
		u=u+ku
	return	s*ku

def priblizek(xx):	#priblizek slabljenja v dB (max +6dB)
	if xx>0.001:
		ww=-16-20*math.log10(abs(xx))
	else:
		ww=6
	if ww>6:
		ww=6
	return ww

x=[]		#polja za risanje grafa
y=[]
z=[]
w=[]

nn=complex(math.sqrt(math.pi/8),-math.sqrt(math.pi/8))
xx=xmin
while xx<xmax:			#zanka cez zaslon
	ss=integral(abs(xx))
	if xx<0:
		ss=-ss-nn
	else:
		ss=-ss+nn
	y=y+[abs(ss)/math.sqrt(math.pi)]
	z=z+[20*math.log10(abs(ss)/math.sqrt(math.pi))]
	w=w+[priblizek(xx)]
	x=x+[xx]
	xx=xx+kx

print(numpy.amax(y),r' @ $h$/$\rho_1$',x[numpy.argmax(y)])
print(numpy.amax(z),r'dB @ $h$/$\rho_1$',x[numpy.argmax(z)])

import matplotlib.pyplot as plt		#risanje grafa funkcije
fig=plt.figure(figsize=(8,5))
plt.plot(x,y,'k-')
plt.axis([xmin,xmax,0,1.2])
plt.xticks([-5,-4,-3,-2,-1,0,1,2,3,4,5])
#plt.title('Slabljenje klinaste ovire v linearnem merilu')
plt.xlabel(r'Velikost ovire $H$/$\rho_1$')
plt.ylabel(r'Slabljenje ovire |$E$/$E_\infty$|')
plt.grid(True)
fig.savefig('slabljenje.png', transparent=True)

fig=plt.figure(figsize=(8,5))
plt.plot(x,z,'k-')
plt.axis([xmin,xmax,-30,5])
plt.xticks([-5,-3.3333333,-1.6666666,0,1.6666666,3.3333333,5],[-30,-20,-10,0,10,20,30])
plt.xlabel(r'Višina ovire H[cm]')
plt.ylabel('Jakost sprejetega polja 20log|$E$/$E_\infty$| [dB]')
plt.grid(False)
fig.savefig('slabljenje2.png', transparent=True)
plt.show()