import sys,csv,math

filename=sys.argv[1]
Freq=[]
Trust=[]
Time=[]
R=[]
Phase=[]
Aux2=[]
Power=[]
index=-1
TR=0.0
TRSigma=0.0
TTSigma=0.0
TT=0.0
TCount=0
Move=0

#Parameters of fits taken from PumpDirect Control measurement on 140820PumpDirect25PMT9Filter40mW, used to correct for electronics
#log10 values used to prevent issues with floating point precision
TP0=0.988297 
TP1=-0.00103723
TP2=-7.98119 #log10 of true value 1.04426e-08 
TP3=-13.42448 #log10 of true value 3.76286e-14 (negative)
TP4=-18.88566 #log10 of true value 1.30118e-19

output = open(filename[:-4]+"Analyzed.txt",'w')
#opening raw data file and reading values
with open(filename,'r') as rdata:
	reader=csv.reader(rdata,delimiter='\t')
	for blank,FReq,TRust,TIme,RR,PHase,AUx2,POwer in reader:
		index += 1
		FReq=float(FReq)
		TRust=int(float(TRust))
		TIme=float(TIme)
		RR=float(RR)
		
		PHase=float(PHase)
		AUx2=float(AUx2)
		POwer=float(POwer)
		Freq.append(FReq)
		Trust.append(TRust)
		Time.append(TIme)
		R.append(RR)
		Phase.append(PHase)
		Aux2.append(AUx2)
		Power.append(POwer)
		if filename[len(filename)-5]=="z":#amplitude modulation (AM) analyzer 
			if TRust == 1:
				TCount+=1
				TR+=RR
				TT+=PHase
			if ((TRust==0)and(TCount!=0)):
				TPower=POwer
				TR=TR/float(TCount)
				TT=TT/float(TCount)
				for x in range(0,TCount):
					TRSigma+=(R[index-TCount+x]-TR)**2
					TTSigma+=(Phase[index-TCount+x]-TT)**2
				TRSigma=(TRSigma/(float(TCount-1)))**0.5#R standard deviation
				TRSigma=TRSigma/((float(TCount))**0.5)#R statistical error
				TTSigma=(TTSigma/(float(TCount-1)))**0.5#Phase standard deviation
				TTSigma=TTSigma/((float(TCount))**0.5)#Phase statistical error
				output.write(str(TPower)+'\t'+str(TR)+'\t'+str(TRSigma)+'\t'+str(TT)+'\t'+str(TTSigma)+'\n')#Writing output file
				TR=0
				TRSigma=0
				TT=0
				TTSigma=0
				TCount=0
		if filename[len(filename)-5]!="z":#frequency modulation (FM) analyzer
			if TRust == 1:
				TCount+=1
				TR+=RR
				TT+=PHase
				TFreq=FReq
			if (((TRust==0)or(TRust==2))and(TCount!=0)):
				TR=TR/float(TCount)
				TT=TT/float(TCount)
				for x in range(0,TCount):
					TRSigma+=(R[index-TCount+x]-TR)**2
					TTSigma+=(Phase[index-TCount+x]-TT)**2
				TRSigma=(TRSigma/(float(TCount-1)))**0.5#R standard deviation
				TRSigma=TRSigma/((float(TCount))**0.5)#R statistical error
				TTSigma=(TTSigma/(float(TCount-1)))**0.5#Phase standard deviation
				TTSigma=TTSigma/((float(TCount))**0.5)#Phase statistical error
				if FReq>1000:
					PP2=2*math.log10(TFreq)+TP2
					PP3=3*math.log10(TFreq)+TP3
					PP4=4*math.log10(TFreq)+TP4
					TT=TT-(TP0+TFreq*TP1+10**(PP2)-10**(PP3)+10**(PP4))#phase correction (offset inherent to electronics)
				if TT<-180:
					TT=TT+360
				if TT>180:
					TT=TT-360
				output.write(str(TFreq)+'\t'+str(TR)+'\t'+str(TRSigma)+'\t'+str(TT)+'\t'+str(TTSigma)+'\t'+str(Move)+'\n')#Writing output file
				TR=0
				TRSigma=0
				TT=0
				TTSigma=0
				TCount=0
				if TRust == 2:
					Move+=1				
				
