import math,sys,csv,scipy,numpy,Gnuplot,Gnuplot.funcutils,TRCurves150512
from scipy import integrate,special
from numpy import vectorize
#This script can analyze multi-stage thermal rectification measurements and correct for initial temperature difference, adjusts Ch2 to Ch1
CL1L2=15.5#Rod Leg 1=15.5 (found 161103 57K)#Tube Leg1=22 (found 150611 and most used 80K),Leg2=18.5 (found 150610 80K)#Thermal capacitance of leg through which heat flows as measaured J/K
#RL1L2=63.95#Leg1=64.Leg2=63.95#Thermal resistance of leg through which heat flows as measured 150611,150610
CL2L1=16#Rod Leg 2=16 (found 161104 57K)#Tube Leg1=21.5 (found 150611 and most used 80K),Leg2=18.5 (found 150610 80K)#Thermal capacitance of leg through which heat flows as measaured in J/K
#RL2L1=64#Leg1=64.Leg2=63.95#Thermal resistance of leg through which heat flows as measured 150611,150610
IntegrationTime=60#seconds spanned in dTdt calculation 
plotAll=1#1=plot all measurements k vs delT, 2=plot individual measurement k vs delT, 3=plot final set of measurements on k vs aveT
measNo=1#which measurement to plot for individual measurements
geom=4.8#geometric factor of sample equal to sample's surface area / thickness in cm
filename="161111Run6CorrectedCut.txt"
outputFileNamePrefix='161111Run6'
outputFileName=str(outputFileNamePrefix)+str(measNo)+'KVsDelT.ps'#output file name for single measurement plot
outputFileNameAll=str(outputFileNamePrefix)+'AllKVsDelT.ps'#output file name for thermal conduct vs delT plot
outputFileNameAllVT=str(outputFileNamePrefix)+'AllKVsAveT.ps'#output file name for thermal conduct. vs temperature plot
plotLabelIntro=r'161111 Si Etch Run 6 Pattern Facing Ch1\n Measurement '+str(measNo)+', A/L= '+str(geom)+'cm, dQmax=0.5W, V0=V'
plotLabelIntroAll=r'161111 Si Etch Run 6 Pattern Facing Ch1\n A/L='+str(geom)+'cm, dQmax=0.5W'#rest of labels auto made

def MeasurementTimer():#finds the starting and stopping times of each measurement
	prevFlag=0.0
	startTime=[]#must be times at which heaters still on but stable on step zeros
	stopTime=[]#must be times at which heaters still on but stable on final steps
	with open(filename,'r') as rdata:
		reader=csv.reader(rdata,delimiter='\t')
		for Ch1,Ch2,Ch3,Ch4,WattA,WattB,Time,Ch1Ch2Stable,StateFlag in reader:
			if (prevFlag==3.0)&(float(StateFlag)==0.0):#end of series of measurements
				stopTime.append(float(Time))
				break
			if ((prevFlag==0.0)&(float(StateFlag)==1.0))|((float(StateFlag)==1.0)&(len(startTime)==0)):#start of measurement
				startTime.append(float(Time))
			if (prevFlag==3.0)&(float(StateFlag)==1.0):#end of measurement
				startTime.append(float(Time))
				stopTime.append(prevTime)
			prevFlag=float(StateFlag)
			prevTime=float(Time)
		return startTime,stopTime

def StageTimer(stageStart,stageStop):#finds the start/stop times of each step in a stage and calculates thermal drift adjustment
	Ch1I=Ch2I=0
	Ch1F=Ch2F=0
	minTimeA=[]#list of starting times of steps with power supply A being high
	maxTimeA=[]#list of ending times of steps with power supply A being high
	minTimeB=[]#list of starting times of steps with power supply B being high
	maxTimeB=[]#list of ending times of steps with power supply B being high
	Ch12StableTimeI=0
	Ch12StableTimeF=0
	Ch2Offset=0#offset correction factor for thermal drift over stage
	Ch2Slope=0#slope correction factor for thermal drift over stage
	prevStable=0
	prevState=0
	prevWattA=0
	prevWattB=0
	prevTime=0
	with open(filename,'r') as rdata:#reads measurement and finds time of each new step
		reader=csv.reader(rdata,delimiter='\t')
		for Ch1,Ch2,Ch3,Ch4,WattA,WattB,Time,Ch1Ch2Stable,StateFlag in reader:
			if (float(Time)>=stageStart)&(float(Time)<stageStop):
				if (float(StateFlag)==1.0):#initial heating step? Will overwrite self, keeping only the last values prior new state
					Ch1I=float(Ch1)
					Ch2I=float(Ch2)
					Ch12StableTimeI=float(Time)
				if float(StateFlag)==3.0:#final heating step? Will overwrite self, keeping only the last values prior new state
					Ch1F=float(Ch1)
					Ch2F=float(Ch2)
					Ch12StableTimeF=float(Time)
				if (float(Ch1Ch2Stable)==1.0)&(prevStable==0.0)&(float(StateFlag)==2.0):#did ch1&ch2 just become stable?
					if float(WattA)>float(WattB):
						minTimeA.append(float(Time))
					if float(WattB)>float(WattA):
						minTimeB.append(float(Time))
				if (float(Ch1Ch2Stable)==0.0)&(prevStable==1.0)&(prevState==2.0):#did ch1&ch2 just stop being stable?
					if prevWattA>prevWattB:
						maxTimeA.append(prevTime)
					if prevWattB>prevWattA:
						maxTimeB.append(prevTime)
			prevWattA=float(WattA)
			prevWattB=float(WattB)
			prevTime=float(Time)
			prevStable=float(Ch1Ch2Stable)
			prevState=float(StateFlag)			

	Ch2Offset=Ch1I-Ch2I#correction trend based on linear change in average temperature over time
	Ch2Slope=(Ch1I-Ch2I-(Ch1F-Ch2F))/(Ch12StableTimeF-Ch12StableTimeI)
	return minTimeA,maxTimeA,minTimeB,maxTimeB,Ch12StableTimeI,Ch12StableTimeF,Ch2Offset,Ch2Slope

def ConductAnalyze(minT, maxT, Ch12StableTimeI, Ch12StableTimeF, Ch2Offset, Ch2Slope, orient):#Analyzes step for resistances. True orient indicates power supply A/Ch1 is high
	timeOut=[]#Time of resistance calculation
	rOut=[]#Thermal resistance of sample
	delTOut=[]#Temperature difference across sample
	aveTOut=[]#Average temperature across the sample
	oldTempHigh=0
	oldTempLow=0
	IntIndex=1
	Ch1Ave=0
	Ch2Ave=0
	Ch3Ave=0
	oldTime=0
	deltdTdt=0
	avedTdt=0
	highTemp=0
	lowTemp=0
	oldTemp=0
	capac=0
	with open(filename,'r') as rdata:
		reader=csv.reader(rdata,delimiter='\t')
		for Ch1,Ch2,Ch3,Ch4,WattA,WattB,Time,Ch1Ch2Stable,StateFlag in reader:
			if (float(Time)>=minT)&(float(Time)<=maxT):
				if (float(Time)>=Ch12StableTimeI)&(float(Time)<=Ch12StableTimeF):#correcting initial temp difference
					Ch2=float(Ch2)+(Ch2Slope*(float(Time)-Ch12StableTimeI)+Ch2Offset)
				if orient:#Ch1 high step
					highTemp=float(Ch1)
					lowTemp=float(Ch2)
					capac=CL1L2
				if not orient:#Ch2 high step
					highTemp=float(Ch2)
					lowTemp=float(Ch1)
					capac=CL2L1
				if IntIndex==1:
					oldTempHigh=highTemp
					oldTempLow=lowTemp
					oldTime=float(Time)
					Ch3Ave=float(Ch3)
				if IntIndex>1:
					Ch3Ave+=float(Ch3)
				if (float(Time)-oldTime)>IntegrationTime:
					deldTdt=(highTemp-oldTempHigh-lowTemp+oldTempLow)/(float(Time)-oldTime)
					avedTdt=0.5*(highTemp-oldTempHigh+lowTemp-oldTempLow)/(float(Time)-oldTime)
					Ch3Ave=Ch3Ave/IntIndex
					delQ=abs(float(WattA)-float(WattB))
					aveQ=0.5*(float(WattA)+float(WattB))
					delT=highTemp-lowTemp
					aveT=0.5*(highTemp+lowTemp)
					conductLow=(aveQ-capac*avedTdt)/(aveT-Ch3Ave)
					conductSample=0.5*((delQ-capac*deldTdt)/delT-conductLow)
					#print 2*conductSample-conductLow,conductLow
					rOut.append(1.0/conductSample)
					timeOut.append(float(Time))
					delTOut.append(delT)
					aveTOut.append(aveT)
					#rL2L1.append(1.0/(0.5*(float(WattB)-float(WattA)-CL2L1*deldTdt)/(Ch2Ave-Ch1Ave)-(0.5*float(WattA)+float(WattB)-CL2L1*avedTdt)/(0.5*(Ch1Ave+Ch2Ave)-Ch3Ave)))
					IntIndex=0
				IntIndex+=1
	#print timeOut
	return rOut, timeOut, delTOut, aveTOut

def ErrAveFinder(resistL,timeL,delTL,aveTL):#cacluation of ave and standard dev (statistical error) over a step
	index=0
	kAve=0
	kErr=0
	delTAve=0
	delTErr=0
	aveTAve=0#average of average of ch1, ch2
	aveTErr=0
	for h,i,j,k in zip(resistL,timeL,delTL,aveTL):#calculate average of resistance and temperature difference
		kAve+=100.0/(geom*h)#calculating average resistance and converting to m from cm
		delTAve+=j#calculating average temperature difference
		aveTAve+=k
		index+=1.0
		#print i, index, timeL[-1]
	kAve=kAve/index
	delTAve=delTAve/index
	aveTAve=aveTAve/index
	for a,b,c,d in zip(resistL,timeL,delTL,aveTL):#calculating standard dev of resistance and temperature difference
		kErr+=(100.0/(geom*a)-kAve)**2
		delTErr+=(c-delTAve)**2
		aveTErr+=(d-aveTAve)**2

	kErr=(kErr/(index-1.0))**0.5
	delTErr=(delTErr/(index-1.0))**0.5
	aveTErr=(aveTErr/(index-1.0))**0.5
	return kAve,kErr,delTAve,delTErr,aveTAve,aveTErr

kA=[]
kAErr=[]
delTA=[]
delTAErr=[]
aveTA=[]
aveTAErr=[]
kB=[]
kBErr=[]
delTB=[]
delTBErr=[]
aveTB=[]
aveTBErr=[]#list of lists, each element in greater list contains data from a single stage
measTimes=MeasurementTimer()
for a,b in zip(measTimes[0],measTimes[1]):#main run sequence
	times=StageTimer(a,b)#find the start/stop times of each step and thermal drift correction
	kAt=[]#lists which temporarily store data from each stage
	kAErrt=[]
	delTAt=[]
	delTAErrt=[]
	aveTAt=[]
	aveTAErrt=[]
	kBt=[]
	kBErrt=[]
	delTBt=[]
	delTBErrt=[]
	aveTBt=[]
	aveTBErrt=[]
	for x,y in zip(times[0],times[1]):#calculate thermal conductivity values over each step power supply A high oreintation
		rData=ConductAnalyze(x,y,times[4],times[5],times[6],times[7],True)
		statData=ErrAveFinder(rData[0],rData[1],rData[2],rData[3])
		kAt.append(statData[0])
		kAErrt.append(statData[1])
		delTAt.append(statData[2])
		delTAErrt.append(statData[3])
		aveTAt.append(statData[4])
		aveTAErrt.append(statData[5])
	kA.append(kAt)#making list of each stage's data list
	kAErr.append(kAErrt)
	delTA.append(delTAt)
	delTAErr.append(delTAErrt)
	aveTA.append(aveTAt)
	aveTAErr.append(aveTAErrt)
	for v,w in zip(times[2],times[3]):#calculate thermal conductivity values over each step power supply B high oreintation
		rData=ConductAnalyze(v,w,times[4],times[5],times[6],times[7],False)
		statData=ErrAveFinder(rData[0],rData[1],rData[2],rData[3])
		kBt.append(statData[0])
		kBErrt.append(statData[1])
		delTBt.append(statData[2])
		delTBErrt.append(statData[3])
		aveTBt.append(statData[4])
		aveTBErrt.append(statData[5])
	kB.append(kBt)
	kBErr.append(kBErrt)
	delTB.append(delTBt)
	delTBErr.append(delTBErrt)
	aveTB.append(aveTBt)
	aveTBErr.append(aveTBErrt)
output = open(filename[:-4]+"TCvsT.txt",'w')
for t,y,u,i,o,p in zip(kA,aveTA,delTA,kB,aveTB,delTB):
	for a,b,c,d,e,f in zip(t,y,u,i,o,p):
		output.write(str(a)+'\t'+str(b)+'\t'+str(c)+'\t'+str(d)+'\t'+str(e)+'\t'+str(f)+'\n')
g=Gnuplot.Gnuplot(persist=1)
if plotAll==1:
	kL1L2Trend1=Gnuplot.Data(delTA[0],kA[0],delTAErr[0],kAErr[0],with_="xyerrorbars ls 2 lc 0",title="Ch1 High Steps, Ave. "+str(round(aveTA[0][-1],1))+"K")
	kL2L1Trend1=Gnuplot.Data(delTB[0],kB[0],delTBErr[0],kBErr[0],with_="xyerrorbars ls 2 lc 7",title="Ch2 High Steps, Ave. "+str(round(aveTB[0][-1],1))+"K")
	kL1L2Trend2=Gnuplot.Data(delTA[1],kA[1],delTAErr[1],kAErr[1],with_="xyerrorbars ls 4 lc 0",title="Ch1 High Steps, Ave. "+str(round(aveTA[1][-1],1))+"K")
	kL2L1Trend2=Gnuplot.Data(delTB[1],kB[1],delTBErr[1],kBErr[1],with_="xyerrorbars ls 4 lc 7",title="Ch2 High Steps, Ave. "+str(round(aveTB[1][-1],1))+"K")
	kL1L2Trend3=Gnuplot.Data(delTA[2],kA[2],delTAErr[2],kAErr[2],with_="xyerrorbars ls 5 lc 0",title="Ch1 High Steps, Ave. "+str(round(aveTA[2][-1],1))+"K")
	kL2L1Trend3=Gnuplot.Data(delTB[2],kB[2],delTBErr[2],kBErr[2],with_="xyerrorbars ls 5 lc 7",title="Ch2 High Steps, Ave. "+str(round(aveTB[2][-1],1))+"K")
	kL1L2Trend4=Gnuplot.Data(delTA[3],kA[3],delTAErr[3],kAErr[3],with_="xyerrorbars ls 6 lc 0",title="Ch1 High Steps, Ave. "+str(round(aveTA[3][-1],1))+"K")
	kL2L1Trend4=Gnuplot.Data(delTB[3],kB[3],delTBErr[3],kBErr[3],with_="xyerrorbars ls 6 lc 7",title="Ch2 High Steps, Ave. "+str(round(aveTB[3][-1],1))+"K")
	kL1L2Trend5=Gnuplot.Data(delTA[4],kA[4],delTAErr[4],kAErr[4],with_="xyerrorbars ls 7 lc 0",title="Ch1 High Steps, Ave. "+str(round(aveTA[4][-1],1))+"K")
	kL2L1Trend5=Gnuplot.Data(delTB[4],kB[4],delTBErr[4],kBErr[4],with_="xyerrorbars ls 7 lc 7",title="Ch2 High Steps, Ave. "+str(round(aveTB[4][-1],1))+"K")
	#kL1L2Trend6=Gnuplot.Data(delTA[5],kA[5],delTAErr[5],kAErr[5],with_="xyerrorbars ls 8 lc 0",title="Ch1 High Steps, Ave. "+str(round(aveTA[5][-1],1))+"K C="+str(CL1L2)+"J/K")
	#kL2L1Trend6=Gnuplot.Data(delTB[5],kB[5],delTBErr[5],kBErr[5],with_="xyerrorbars ls 8 lc 7",title="Ch2 High Steps, Ave. "+str(round(aveTB[5][-1],1))+"K C="+str(CL2L1)+"J/K")
	g('set terminal postscript color eps enhanced')
	g('set output "'+str(outputFileNameAll)+'"')
	g('set multiplot layout 1,1 title "'+str(plotLabelIntroAll)+', All Meas. dt='+str(IntegrationTime)+'s"')
	g('set grid')
	#g('set yrange [0:6]')
	g('set ylabel "Thermal Conductivity (W/m-K)"')
	g('set xlabel "Temperature Difference Across Sample (K)"')
	g('set key top right')
	g.plot(kL1L2Trend1,kL2L1Trend1,kL1L2Trend2,kL2L1Trend2,kL1L2Trend3,kL2L1Trend3,kL1L2Trend4,kL2L1Trend4,kL1L2Trend5,kL2L1Trend5)#,kL1L2Trend5,kL2L1Trend5, kL1L2Trend6,kL2L1Trend6)
if plotAll==2:
	kL1L2Trend1=Gnuplot.Data(delTA[measNo-1],kA[measNo-1],delTAErr[measNo-1],kAErr[measNo-1],with_="xyerrorbars ls 1 lc 0",title="Ch1 High Steps, Ave. "+str(round(aveTA[measNo-1][-1],1))+"K C="+str(CL1L2)+"J/K")
	kL2L1Trend1=Gnuplot.Data(delTB[measNo-1],kB[measNo-1],delTBErr[measNo-1],kBErr[measNo-1],with_="xyerrorbars ls 1 lc 7",title="Ch2 High Steps, Ave. "+str(round(aveTB[measNo-1][-1],1))+"K C="+str(CL2L1)+"J/K")
	g('set terminal postscript color eps enhanced')
	g('set output "'+str(outputFileName)+'"')
	g('set multiplot layout 1,1 title "'+str(plotLabelIntro)+' '+str(int(measTimes[0][measNo-1]))+'s to '+str(int(measTimes[1][measNo-1]))+'s, dt='+str(IntegrationTime)+'s"')
	g('set grid')
	g('set ylabel "Thermal Conductivity (W/m-K)"')
	g('set xlabel "Temperature Difference Across Sample (K)"')
	g('set key top right')
	g.plot(kL1L2Trend1,kL2L1Trend1)
#g('set yrange [0.054:0.062]')
#for i,v in zip(maxTA,maxTB):
#	if v<i:
#		g('set arrow from '+str(i/1000.0-startTime/1000.0)+',graph(0,0) to '+str(i/1000.0-startTime/1000.0)+',graph(1,1) nohead lc rgb "blue" lw 0.5')
#	if v>i:	
#		g('set arrow from '+str(v/1000.0-startTime/1000.0)+',graph(0,0) to '+str(v/1000.0-startTime/1000.0)+',graph(1,1) nohead lc rgb "blue" lw 0.5')
if plotAll==3:
	kLastA=[]
	kLastAErr=[]
	kLastB=[]
	kLastBErr=[]
	aveTemperatureA=[]
	aveTemperatureB=[]
	for x,y,z in zip(kA,kAErr,aveTA):#making list of final thermal conductivity measurements for each V0 setting Ch1 High
		kLastA.append(x[-1])
		kLastAErr.append(y[-1])
		aveTemperatureA.append(z[-1])
	for x,y,z in zip(kB,kBErr,aveTB):#making list of final thermal conductivity measurements for each V0 setting Ch2 High
		kLastB.append(x[-1])
		kLastBErr.append(y[-1])
		aveTemperatureB.append(z[-1])
	kL1L2Last=Gnuplot.Data(aveTemperatureA,kLastA,kLastAErr,with_="yerrorbars ls 2 lc 0",title="Ch1 High Step")
	kL2L1Last=Gnuplot.Data(aveTemperatureB,kLastB,kLastBErr,with_="yerrorbars ls 2 lc 7",title="Ch2 High Step")
	g('set terminal postscript color eps enhanced')
	g('set output "'+str(outputFileNameAllVT)+'"')
	g('set multiplot layout 1,1 title "'+str(plotLabelIntroAll)+', Final Steps, dt='+str(IntegrationTime)+'s"')
	g('set grid')
	g('set ylabel "Thermal Conductivity (W/m-K)"')
	g('set xlabel "Average Sample Temperature (K)"')
	g('set key bottom right')
	g.plot(kL1L2Last,kL2L1Last)
#g('set grid')
#print "Final Temperature Percentage Difference of Each Step" 
#print "Ch1 High Step, Ch2 High Step"
#for x,y in zip(rCh1OverCh2,rCh2OverCh1):
#	print str(round(x,2)),"%, ",str(round(y,2)),"%"
