#!/usr/bin/python3
print('*** FM vremenski radar - S53MV 05.05.2022 ***')

ime=input('Vhodni zapis (prazno=nova meritev): ')	#odpri vhodni zapis
if len(ime)==0:
	import time,serial	#uporabljene funkcije
	cime='/dev/ttyUSB0'	#ime zaporednega vmesnika vrtiljaka
	print('Počakaj na lenobni USB...')
	com=serial.Serial(cime,9600,timeout=0)	#odpri comport,9600bps,8bit,noparity
	time.sleep(2)		#počakaj na lenobni USB
	s=com.read(999)		#izprazni vmesnik comporta
	com.write(bytes([83,8]))	#Dušenje AZ - postavi časovne konstante KR5600
	com.write(bytes([84,8]))	#Dušenje EL
	com.write(bytes([85,22]))	#Vztrajnost AZ
	com.write(bytes([86,22]))	#Vztrajnost EL
	el0=16			#konstante pretvorbe podatkov vrtiljaka (cela števila)
	el90=130
	azm180=17		#min 15
	az180=247		#max 249

	import sounddevice as sd	#poiskati na spletu in naložiti s pip3!
	import numpy as np		#uporaba učinkovitih funkcij numpy
	fs=96000		#frekvenca vzročenja zvočnega signala (max 192000Hz)
	pts=2048		#število točk spektra 0...fs/2 (potenca 2 za učinkovit FFT)
	avg=8			#število povprečenj spektra
	sd.default.device=0	#hardverska izbira vhoda
	Hzkm=1627		#faktor pretvorbe frekvence v dolžino [km]
	domet=fs/2/Hzkm		#izračunaj domet [km]
	fskala=np.linspace(0,domet,pts)	#izračunaj dolžinsko skalo [km]

	t=time.ctime()			#ime zapisa iz časa meritve
	ime=t[20:24]+t[4:7]+t[8:10]+'-'+t[11:13]+t[14:16]+t[17:19]
	ime=ime.replace(' ','0')
	zapis=open(ime+'.csv','w')	#ustvari izhodni zapis
	zapis.write(ime+'\n')		#ime meritve v prvo vrstico
	print('Ime izhodnega zapisa: '+ime)
	com.write(bytes([80]))		#določi smer vrtenja iz začetnega azimuta
	time.sleep(0.1)
	s=com.read(99)
	smer=1
	smerbyte=az180
	if int(s[3])>128:
		smer=-1
		smerbyte=azm180
	zapis.write(str(smer)+'\n')	#smer v drugo vrstico
	print('Smer='+str(smer))
	d='Domet '+"%.2f"%domet+' km'	#domet v tretjo vrstico
	zapis.write(d+'\n')
	print(d)

	stevec=0			#zanka meritev spektra
	dolzina=0
	while stevec<1555:
		com.write(bytes([80]))	#vprašaj za parametre vrtiljaka KR5600
		
		s=sd.rec(pts*(1+avg),samplerate=fs,channels=1)	#vzorčenje ADC zvočne kartice, 2D polje!
		sd.wait()		#počakaj do konca vzorčenja
		nz=1.e-12		#dodatek za neničelni argument logaritma
		s2=np.square(s)		#izpis jakosti/max [dB]
		jakost=' Jakost '+"%.1f"%(10*np.log10(nz+np.mean(s2)))+'dB Vrh '+"%.1f"%(10*np.log10(nz+np.amax(s2)))+'dB'
		o=1+np.cos(np.linspace(-np.pi,np.pi,2*pts))	#okno dvignjeni kosinus
		k=0
		s2=np.zeros(pts)
		while k<pts*avg:	#seštevanje moči FFT
			s2=s2+np.square(np.abs(np.fft.rfft(o*s.reshape(-1)[k:k+2*pts])[:-1]))
			k+=pts
		logsp=10*np.log10(nz+0.5*s2/pts/pts/avg)	#povprečje [dB]
		jakost+=' Max '+"%.1f"%(np.amax(logsp))+'dB @'+"%.1f"%(fskala[np.argmax(logsp)])+'km'

		azel='Vrtiljak '	#prečitaj azimut in elevacijo vrtiljaka KR5600
		azb=-1
		s=com.read(99)
		if len(s)!=20:		#napaka pri prenosu podatkov?
			azel+='ni podatka!!!'
		else:
			azb=int(s[3])	#byte azimut
			elb=int(s[4])	#byte elevacija
			azel+=str(azb)+'/'+str(elb)
			az=(azb-azm180)*360.0/(az180-azm180)-180	#preračunaj azimut
			if az<0:
				az+=360
			el=(elb-el0)*90.0/(el90-el0)			#preračunaj elevacijo
			azel+=str(' AZ=%05.1f'%az)+str('/EL=%04.1f'%el)

		zapis.write(azel+jakost+'\n')	#izpis vrtiljaka in jakosti
		print(azel+jakost)
		j=0				#zanka za zapis spektra
		ss='['
		while j<len(logsp):
			ss+="%.1f"%(logsp[j])+','
			j+=1
		zapis.write(ss[:-1]+']\n')
		com.write(bytes([81,smerbyte]))	#poženi vrtiljak KR5600
		stevec+=1			#števec meritev
		dolzina+=1
		if smer==1 and azb>=az180:	#konec vrtenja?
			stevec+=888
		if smer==-1 and azb<=azm180:
			stevec+=888

	print('Dolzina zapisa '+str(dolzina)+' azimutov')
	com.write(bytes([0]))	#ustavi vrtiljak KR5600
	zapis.close()		#zapri izhodni zapis
	com.close()		#zapri comport

zapis=open(ime+'.csv','r')	#odpri vhodni zapis
s=zapis.readline()		#določi ime meritve iz prve vrstice
ime=s[:-1]
print('Ime meritve: '+ime)
s=zapis.readline()		#določi smer iz druge vrstice
smer=int(s[0:2])
print('Smer vrtenja '+str(smer))
s=zapis.readline()		#določi domet iz tretje vrstice
ss=s.split(' ')
domet=float(ss[1])
print('Domet '+str(domet)+' km')
s=zapis.readline()		#določi elevacijo iz četrte vrstice
ss=s.split('/')
s=ss[2]
ss=s.split(' ')
eleva=ss[0]
print(str(eleva)+'°')

def radij():			#vstavi polmer risanja
	s=input('Polmer[km] risbe: ')
	p=domet
	if len(s)>0:
		if float(s)<domet:
			p=float(s)
	print('Polmer '+str(p)+' km')	
	return p

polmer=radij()			#zanka po različnih polmerih (polmer<0=konec)
minodboj=-72			#minimum vidnega odboja
maxodboj=-16			#maximum vidnega odboja
while polmer>0:
	s=zapis.readline()		#zanka čitanja azimutov
	polje=[]			#prazno začetno polje podatkov
	meja=int(len(s[1:-2].split(','))*polmer/domet)	#meja za izbran polmer
	while len(s)>0:
		ss=s[1:-2].split(',')	#dodatek azimuta v polje float
		a=[float(n) for n in ss]
		polje+=[a[:meja]]
		s=zapis.readline()	#čitanje vrstice v prazno
		s=zapis.readline()	#čitanje vrstice
	print('Prečitano '+str(len(polje))+' azimutov')

	#https://stackoverflow.com/questions/31051882/matplotlib-density-plot-in-polar-coordinates
	import numpy as np
	import matplotlib.pyplot as plt
	import matplotlib.cm as cm
	aa=np.linspace(0,2*np.pi,len(polje))	#pripravi podatke #vrstic,#stoplcev,transponiraj
	bb=np.linspace(0,polmer,len(polje[0]))
	cc=np.transpose(np.clip(polje,minodboj,maxodboj))	#skala z rezanjem šuma
	fig=plt.figure(figsize=(16,12))		#risanje polarnega diagrama
	ax=plt.subplot(111,polar=True)
	ax.set_xticklabels([])			#izklop xticklabel
	ctf=ax.contourf(aa,bb,cc,cmap=cm.gist_ncar)	#https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html#colorcet
	ax.set_theta_zero_location("S")		#'south' location for theta=0
	ax.set_theta_direction(-smer)		#angles increase clockwise=-1
	plt.colorbar(ctf)
	plt.title(ime+'  '+eleva+'°  Smer='+str(smer)+'  '+str(polmer)+'/'+str(domet)+'km', fontsize=20)
	plt.xlabel('S', fontsize=20)
	plt.grid(color='r',linestyle='-',linewidth=0.3)
	fig.savefig(ime+'-'+str(polmer)+'.png')
	plt.show()

	zapis.seek(0)			#pojdi na začetek zapisa
	s=zapis.readline()		#preskoči glavo štirih vrstic
	s=zapis.readline()
	s=zapis.readline()
	s=zapis.readline()
	polmer=radij()			#nov polmer za risanje

zapis.close()			#zapri vhodni zapis
print('*** Konec programa ***')