##### Sledenje zvezd z vrtiljakom EGIS - S53MV 24.08.2016 #####

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

port='COM3'	#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=25	#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 datumura(t):		#izpisi trenutni datum in uro
	leto=t[0]
	mesec=t[1]
	dan=t[2]
	print str(leto)+'-'+str(mesec)+'-'+str(dan),		#izpisi datum
	ura=t[3]
	minuta=t[4]
	sekunda=t[5]
	print str(ura)+':'+str(minuta)+':'+str(sekunda),	#izpisi uro

def smerzvezde(z,t):		#izracunja smerni vektor zvezde 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+z[1]				#izracunaj smerni vektor zvezde
	cr=math.cos(math.radians(zdolzina))
	sr=math.sin(math.radians(zdolzina))
	cd=math.cos(math.radians(z[2]))
	sd=math.sin(math.radians(z[2]))
	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))

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

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

t=time.gmtime()			#precitaj datum in uro realnega casa racunalnika (UTC)
datumura(t)			#izpisi trenutni datum in uro
vrtiljak()			#izpisi azimut in elevacijo vrtiljaka

#seznam imen, rektascenzij (stopinje) in deklinacij (stopinje) opazovanih zvezd
zvezda=['Zvezda',0.0,0.0]
sirius=['Sirius',101.284,-16.743]
alnath=['Alnath',81.838,28.620]
rigel=['Rigel',78.835,-8.186]
betelgeuse=['Betelgeuse',89.019,7.406]
aldebaran=['Aldebaran',69.221,16.539]
capella=['Capella',79.482,46.011]
altair=['Altair',297.900,8.916]
deneb=['Deneb',310.502,45.343]
sadr=['Sadr',305.708,40.313]
scutum=['Scutum',279.029,-8.231]
polaris=['Polaris',43.272,89.333]
arcturus=['Arcturus',214.104,19.096]
spica=['Spica',201.517,-11.248]
vega=['Vega',279.375,38.800]
altair=['Altair',297.900,8.913]
regulus=['Regulus',152.313,11.886]
mizar=['Mizar',201.105,54.839]
procyon=['Procyon',115.042,5.181]
etapup=['Eta Puppis',117.500,-24.902]
mirphak=['Mirphak',51.379,49.919]
sagittariusa=['Sagittarius A',265.538,-28.833]
cygnusa=['Cygnus A',299.867,40.734]
cassiopeiaa=['Cassiopeia A',350.850,58.815]
crabnebula=['Crab Nebula',290.925,22.014]
tricstiristo=['3C400',290.925,14.509]
oriona=['Nebula Orion A',83.821,-5.391]

seznam=[zvezda,sirius,alnath,rigel,betelgeuse,aldebaran,capella,altair,deneb,sadr,scutum,polaris,arcturus,
spica,vega,altair,regulus,mizar,procyon,etapup,mirphak,sagittariusa,cygnusa,cassiopeiaa,crabnebula,tricstiristo,oriona]

print 'Seznam zvezd (0 vstavi rektascenzijo in deklinacijo)'
izbira=1	#izpisi seznam zvezd in izberi eno zvezdo
while izbira<27:
	zvezda=seznam[izbira]
	smer=smerzvezde(zvezda,t)	#izracunja smerni vektor zvezde ob dani uri in datumu
	azimut,elevacija=azel(smer,ovzhod,osever,ogor)	#izracunaj azimut in elevacijo za opazovalca
	print izbira,zvezda[0],str('%5.1f'%(azimut))+str('/%4.1f'%(elevacija))
	izbira=izbira+1
izbira=input('Izberi zvezdo: ')
zvezda=seznam[izbira]
if izbira==0:
	zvezda[1]=input('Rektascenzija: ')
	zvezda[2]=input('Deklinacija: ')
	zvezda[0]=str('Zvezda(%5.1f'%(zvezda[1]))+str('/%4.1f'%(zvezda[2])+')')

while stracun>0:			#zanka racunanja
	t=time.gmtime()			#precitaj datum in uro realnega casa racunalnika (UTC)
	datumura(t)			#izpisi trenutni datum in uro
	smer=smerzvezde(zvezda,t)	#izracunja smerni vektor zvezde ob dani uri in datumu

	azimut,elevacija=azel(smer,ovzhod,osever,ogor)	#izracunaj azimut in elevacijo za opazovalca
	print zvezda[0],str('%5.1f'%(azimut))+str('/%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)
		print 'Popravki',
		print str('%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))+'\r'
		else:			#lihi stracun za elevacijo
			s='&*EL='+str('%06.2f'%(el))+'\r'
		comport.write(s)
		print '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

#konec programa