#!/usr/bin/python
##################################################################################
# Average C1 values - Outlier v1
# Rel. 10/2015
# Dr. Rene Cejas Bolecek
# Low Temperatures Laboratory, CAB, Argentine.
# email: reneczechdev@gmail.com
##################################################################################
#Libraries
from numpy import zeros, sqrt, where, pi, average, arange, histogram,log,log10 # log is ln
import numpy as np
from PyAstronomy import pyasl
#http://www.hs.uni-hamburg.de/DE/Ins/Per/Czesla/PyA/PyA/pyaCDoc/installingPyA.html
#User settings
filename = '14_10T_norm_f1'
ColT = 0 #En los archivos de J.G. es 0
ColC1 = 6 #En los archivos de J.G. es 6
#Av. settings
Tcutoff = 0.01
#Outlier settings
polyFit = 1
stdlim = 0.05
flagDEB = False #True or False
###################################################################
#Output Files
writeIt1 = open(filename+"_Av_V1.dat", "w")
writeIt2 = open(filename+"_Outlier_V1.dat", "w")
###################################################################
#Gnuplot check
#P3(x)=a0+a1*x+a2*x**2+a3*x**3
#a0 = 1.99699e-06-8E-6
#a1 = 5.93816e-07
#a2 = 3.49326e-09
#a3= -3.94088e-11
#plot "14_10T_norm.dat" u 1:7, P3(x), "14_10T_norm_f1.dat" u 1:7, "14_10T_norm_f1_AvV1.dat" u 1:3, "14_10T_norm_f1_OutlierV1.dat" u 1:2
#Var init
Tav = 0
C1av = 0
i0 = 0
flag = True
data = np.loadtxt(filename+'.dat')
T = data[:,ColT]
C1 = data[:,ColC1]
#Calculation init
for i in range(len(data)-1):
TCalc = []
TCalc.append(T[i])
C1Calc = []
C1Calc.append(C1[i])
j = i+1
if (i == i0 and j<=(len(data)-1)):
flag = True
if flag:
counter = 1
while (True and j <= (len(data)-1)):
if (abs(T[i0]-T[j])/T[j])> Tcutoff:
i0 = j
#http://www.hs.uni-hamburg.de/DE/Ins/Per/Czesla/PyA/PyA/pyaslDoc/aslDoc/outlier.html
#iin, iout,p , model = pyasl.polyResOutlier(TCalc, C1Calc, deg=3, stdlim=1, controlPlot=True,fullOutput = True)
iin, iout= pyasl.polyResOutlier(TCalc, C1Calc, deg=polyFit, stdlim=stdlim, controlPlot=False)
TOutliered = zeros(len(iin))
C1Outliered = zeros(len(iin))
for idx in range(len(iin)):
print idx, iin
TOutliered[idx] = TCalc[iin[idx]]
C1Outliered[idx] = C1Calc[iin[idx]]
writeIt2.write("%.5f %.10f\n" % (TOutliered[idx],C1Outliered[idx]))
writeIt1.write("%.5f %.5f %.10f %.10f %d\n" % (np.mean(TCalc),np.std(TCalc),np.mean(C1Calc),np.std(C1Calc),len(C1Calc)))
break
flag = False
C1Calc.append(C1[j])
TCalc.append(T[j])
counter = counter +1
if (j == (len(data)-1)):
if(flagDEB and counter!=0): print "save Tav:", i," ",Tav/(counter*1.0)
writeIt1.write("%.5f %.5f %.10f %.10f %d\n" % (np.mean(TCalc),np.std(TCalc),np.mean(C1Calc),np.std(C1Calc),len(C1Calc)))
j = j+1
if(flagDEB): print "Tav: ", Tav, "j: ", j ,"counter: ", counter
writeIt1.close()
writeIt2.close()
print "Process finished"