##### Potovanje po galaksiji z vrtiljakom EGIS - S53MV 03.12.2016 #####

#izvajanje potrebuje najmanj python-2.7.12 in pyserial-2.7.win32_py3k
import time, math, serial	#uporabljene funkcije

port='COM7'	#ime zaporednega vmesnika
pomlad=0.275	#poenostavljeno pomladisce za Julijanski koledar (del kroga)
elmin=10.0	#minimalna varna elevacija vrtiljaka

stracun=9999	#stevec racunov, liho zacni z elevacijo vrtiljaka!
casracun=18	#interval racunanja (sekunde)

sirina=45.948	#zemljepisna sirina opazovalca v stopinjah (+sever)
dolzina=13.634	#zemljepisna dolzina opazovalca v stopinjah (+vzhod)

cs=math.cos(math.radians(sirina))	#smerni vektorji opazovalca ovzhod, osever in ogor
ss=math.sin(math.radians(sirina))
cd=math.cos(math.radians(dolzina))
sd=math.sin(math.radians(dolzina))
ovzhod=[-sd,cd,0]
osever=[-ss*cd,-ss*sd,cs]
ogor=[cs*cd,cs*sd,ss]

nagibew=-1.0	#smerni vektorji nagnjenega vrtiljaka vvzhod, vsever in vgor
c=math.cos(math.radians(nagibew))
s=math.sin(math.radians(nagibew))
vvzhod=[c*ovzhod[0]-s*ogor[0],c*ovzhod[1]-s*ogor[1],c*ovzhod[2]-s*ogor[2]]
vsever=osever
vgor=[s*ovzhod[0]+c*ogor[0],s*ovzhod[1]+c*ogor[1],s*ovzhod[2]+c*ogor[2]]

def skalar(a,b):		#skalarni produkt dveh vektorjev
	return (a[0]*b[0])+(a[1]*b[1])+(a[2]*b[2])

def azel(s,ov,os,og):		#izracunaj azimut in elevacijo za opazovalca
	az=math.degrees(math.atan2(skalar(s,ov),skalar(s,os)))	#azimut
	if az<0:
		az=az+360
	el=math.degrees(math.asin(skalar(s,og)))		#elevacija
	return az,el

def pw(s):			#izpis na zaslon in v izhodni zapis brez nove vrste
	print s,
	iz.write(s+' ')

def pwn(s):			#izpis na zaslon in v izhodni zapis z novo vrsto
	print s
	iz.write(s+'\r\n')

def datumura(t):		#izpisi trenutni datum in uro
	leto=t[0]
	mesec=t[1]
	dan=t[2]
	pw(str(leto)+'-'+str(mesec)+'-'+str(dan))	#izpisi datum
	ura=t[3]
	minuta=t[4]
	sekunda=t[5]
	pw(str(ura)+':'+str(minuta)+':'+str(sekunda))	#izpisi uro

def smerneba(ra,dec,t):		#izracunja smerni vektor neba ob dani uri in datumu
	ddan=(((t[5]/60.0)+t[4])/60.0+t[3])/24.0	#decimalni del dneva
	pleto=t[0]-int((t[0]/4)*4)			#dan v stiriletju
	sdan=t[2]
	dmesec=[31,28,31,30,31,30,31,31,30,31,30,31]	#dodaj stevilo dni v mesecih
	n=1	
	while n<t[1]:
		sdan=sdan+dmesec[n-1]
		if pleto==0:				#popravek za prestopno leto
			if n==2:
				sdan=sdan+1
		n=n+1
	dleto=[366,365,365,365]				#dodaj stevilo dni v letih
	n=0
	while n<pleto:
		sdan=sdan+dleto[n]
		n=n+1
	zdolzina=(ddan*(1+(1/365.25)))+(sdan/365.25)	#izracunaj Julijansko zemljepisno dolzino
	zdolzina=zdolzina+pomlad			#dodaj zacetni polozaj pomladisca
	zdolzina=(1+int(zdolzina)-zdolzina)*360.0	#rezultat v stopinje (+vzhod)
	zdolzina=zdolzina+ra				#izracunaj smerni vektor neba
	cr=math.cos(math.radians(zdolzina))
	sr=math.sin(math.radians(zdolzina))
	cd=math.cos(math.radians(dec))
	sd=math.sin(math.radians(dec))
	return [cd*cr,cd*sr,sd]

def popravek(e):		#popravek elevacije vrtiljaka
	if e<0:			#elevacija ni negativna
		e=0
	p=[0.0,0.01,2.6,6.7,10.8,14.7,18.7,22.7,26.7,30.5,34.0,37.8,41.6,45.9,50.6,55.1,60.5,65.6,71.1,72.0]
	korak=5.0		#korak tabele
	n=int(e/korak)
	return p[n]+((p[n+1]-p[n])*(e/korak-n))

def gllvrd(glon,glat):				#galakticne koordinate v rektascenzijo,deklinacijo
	glon=glon-32.9				#galakticna dolzina dviznega vozla galaksije 32.9
	glon=math.radians(glon)			#pretvorba v radiane
	glat=math.radians(glat)
	x=math.cos(glon)*math.cos(glat)		#pretvorba v kartezicne koordinate
	y=math.sin(glon)*math.cos(glat)
	z=math.sin(glat)
	nagib=math.radians(62.8)		#nagib ravnine galaksije 62.8
	yy=y*math.cos(nagib)-z*math.sin(nagib)
	z=y*math.sin(nagib)+z*math.cos(nagib)
	ra=math.degrees(math.atan2(yy,x))	#pretvorba v krogelne koordinate
	dec=math.degrees(math.asin(z))
	ra=ra+282.9				#rektascenzija dviznega vozla galaksije 282.9
	if ra>360:				#omeji obmocje rektascenzije 0...360
		ra=ra-360
	return ra,dec

def vrtiljak():			#izpisi azimut in elevacijo vrtiljaka
	s=comport.read(99)		#izprazni vmesnik comporta
	comport.write('&*AZ=?\r')	#vprasaj za azimut vrtiljaka
	sa=comport.read(99)
	comport.write('&*EL=?\r')	#vprasaj za elevacijo vrtiljaka
	se=comport.read(99)
	pwn('Vrtiljak '+sa[0:9]+' '+se[0:9])	#izpisi azimut in elevacijo

#odpri comport port, 9600bps, 8bit, brez paritete
comport=serial.Serial(port,9600,timeout=0.1)

t=time.ctime()			#ime izhodnega zapisa
f=t[20:24]+t[4:7]+t[8:10]+t[11:13]+t[14:16]+t[17:19]+'.txt'
iz=open(f,'w')			#odpri izhodni zapis

pwn('*** Potovanje po galaksiji z vrtiljakom EGIS ***')
t=time.gmtime()			#precitaj datum in uro realnega casa racunalnika (UTC)
datumura(t)			#izpisi trenutni datum in uro
vrtiljak()			#izpisi azimut in elevacijo vrtiljaka

glon=input('Vstavi zacetno galakticno dolzino: ')
glat=input('Vstavi zacetno galakticno sirino: ')
klon=input('Vstavi korak galakticne dolzine: ')
klat=input('Vstavi korak galakticne sirine: ')

pwn('*** Nastavitve za korak '+str(casracun)+' sekund ***')
pwn('Zacetna galakticna dolzina '+str(glon))
pwn('Zacetna galakticna sirina '+str(glat))
pwn('Korak galakticne dolzine '+str(klon))
pwn('Korak galakticne sirine '+str(klat))

while stracun>0:			#zanka racunanja
	t=time.gmtime()			#precitaj datum in uro realnega casa racunalnika (UTC)
	datumura(t)			#izpisi trenutni datum in uro

	pw(str('GLON=%5.1f'%(glon))+str('/GLAT=%4.1f'%(glat)))

	ra,dec=gllvrd(glon,glat)	#galakticne koordinate v rektascenzijo,deklinacijo
	pwn(str('RA=%5.1f'%(ra))+str('/DEC=%4.1f'%(dec)))

	smer=smerneba(ra,dec,t)		#izracunja smerni vektor neba ob dani uri in datumu
	azimut,elevacija=azel(smer,ovzhod,osever,ogor)	#izracunaj azimut in elevacijo za opazovalca
	pw(str('AZ=%5.1f'%(azimut))+str('/EL=%4.1f'%(elevacija)))

	vrtiljak()		#izpisi azimut in elevacijo vrtiljaka

	azimut,elevacija=azel(smer,vvzhod,vsever,vgor)	#izracunaj azimut in elevacijo za vrtiljak
	if elevacija>elmin:	#postavi vrtiljak na azimut/elevacijo
		el=popravek(elevacija)
		pw(str('Popravki %5.1f'%(azimut))+str('/%4.1f'%(elevacija))+str('>%4.1f'%(el)))
		if stracun%2==0:	#sodi stracun za azimut
			s='&*AZ='+str('%06.2f'%(azimut))
		else:			#lihi stracun za elevacijo
			s='&*EL='+str('%06.2f'%(el))
		comport.write(s+'\r')
		pwn('Ukaz vrtiljaku: '+s)
	else:		#zacni vedno z lihim stevilom >>> elevacija
		stracun=(stracun//2)*2
	stracun=stracun-1		#zakljucek zanke racunanja
	time.sleep(casracun)		#cakanje do naslednjega racuna
	glon=glon+klon			#povecaj glon in glat za korak
	glat=glat+klat

iz.close()		#zapri izhodni zapis

#konec programa