#!/usr/bin/python """ceci est un programme de voltametrie cyclique, il prend en parametre d'execution le nom d'un fichier .dat contenant les caracteristiques de la volta cyclique. de maniere optionnel il peut integerer spécifiquement des cycles $2 à $3""" import os,string,sys,time,readline,dircache,array,random from math import sqrt #,integre """ Returns coefficients to the regression line "y=ax+b" from x[] and y[]. Basically, it solves Sxx a + Sx b = Sxy Sx a + N b = Sy where Sxy = \sum_i x_i y_i, Sx = \sum_i x_i, and Sy = \sum_i y_i. The solution is a = (Sxy N - Sy Sx)/det b = (Sxx Sy - Sx Sxy)/det where det = Sxx N - Sx^2. In addition, Var|a| = s^2 |Sxx Sx|^-1 = s^2 | N -Sx| / det |b| |Sx N | |-Sx Sxx| s^2 = {\sum_i (y_i - \hat{y_i})^2 \over N-2} = {\sum_i (y_i - ax_i - b)^2 \over N-2} = residual / (N-2) R^2 = 1 - {\sum_i (y_i - \hat{y_i})^2 \over \sum_i (y_i - \mean{y})^2} = 1 - residual/meanerror It also prints to few other data, N, a, b, R^2, s^2, which are useful in assessing the confidence of estimation. """ def linreg(X, Y): if len(X) != len(Y): raise ValueError, 'unequal length' N = len(X) Sx = Sy = Sxx = Syy = Sxy = 0.0 for x, y in map(None, X, Y): Sx = Sx + x Sy = Sy + y Sxx = Sxx + x*x Syy = Syy + y*y Sxy = Sxy + x*y det = Sxx * N - Sx * Sx a, b = (Sxy * N - Sy * Sx)/det, (Sxx * Sy - Sx * Sxy)/det meanerror = residual = 0.0 for x, y in map(None, X, Y): meanerror = meanerror + (y - Sy/N)**2 residual = residual + (y - a * x - b)**2 RR = 1 - residual/meanerror ss = residual / (N-2) Var_a, Var_b = ss * N / det, ss * Sxx / det # print "y=ax+b" # print "N= %d" % N # print "a= %g \\pm t_{%d;\\alpha/2} %g" % (a, N-2, sqrt(Var_a)) # print "b= %g \\pm t_{%d;\\alpha/2} %g" % (b, N-2, sqrt(Var_b)) # print "R^2= %g" % RR # print "s^2= %g" % ss # print "ecart type " +str(sqrt(meanerror)/N) return a, b, ss ############################################################# def F(a,b,c,d,e,f,g,h,i,j,k): #SG4, 11 points #coefs=[0.042,-0.105,-0.023,0.14,0.28,0.333,0.28,0.14,-0.023,-0.105,0.042] #SG2, 11 points coefs=[-0.084, 0.021, 0.103, 0.161, 0.196, 0.207, 0.196, 0.161,0.103, 0.021, -0.084] return (a*coefs[0] + b*coefs[1] + c*coefs[2] + d*coefs[3] + e*coefs[4] + f*coefs[5] + g*coefs[6] + h*coefs[7] + i*coefs[8] +j*coefs[9] + k*coefs[10]) def lisse(x): tx = ([ x[0] ] * 5) + x + ([ x[-1] ] * 5 ) rx = [] for i in range(5,len(tx)-5): rx.append(apply(F,tuple(tx[i-5:i+6]))) return rx ############################################################# def Deriv(exp,I): """ derive la fonction I sur 5 points """ J=[0,0] res=(1.0*abs(exp[6]-exp[5])+abs(exp[6]-exp[7]))/exp[11] print "calcul de la dérivée" # print len(I) for i in range(len(I)-4): J.append((I[i]-8*I[i+1]+8*I[i+3]-I[i+4])/12/res) J.append(0) J.append(0) print "Valeur Maximale de la dérivée %g"%(max(J)) return J def LisFich(file): """lit le fichier 'file' et retourne les potentiels E, les intensités I et la vitesse de ballayage""" TE=[] TI=[] E=[] I=[] nc=1 exp=["user","name",0,"date","Localisation",-200,500,-200,50,1,1,1,0] exp[4]=file[:-4] Fichier=open(file,"r") for ligne in Fichier.readlines(): if ligne[0]=="#": #on cherche la vitesse, mot clé : mV/s if string.find(ligne,"Potentiel")==1 or string.find(ligne,"#E/I")==0: if len(E)>0: TE.append(E) TI.append(I) I=[] E=[] continue print ligne[1:-1] if string.find(ligne,"mV/s")==1: exp[8]=int(string.split(ligne,":")[1]) elif string.find(ligne,"Cycle")==1: exp[9]=int(string.split(ligne,":")[1]) elif string.find(ligne,"nPt")==1: exp[11]=int(string.split(ligne,":")[1]) elif string.find(ligne,"Eini")==1: exp[5]=int(string.split(ligne,":")[1]) elif string.find(ligne,"Emax")==1: exp[6]=int(string.split(ligne,":")[1]) elif string.find(ligne,"Efin")==1: exp[7]=int(string.split(ligne,":")[1]) elif string.find(ligne,"user")==1: exp[0]=string.split(ligne,":")[1][1:-1] elif string.find(ligne,"Liss")==1: exp[10]=int(string.split(ligne,":")[1]) elif string.find(ligne,"exp")==1: exp[2]=int(string.split(ligne,":")[1]) elif string.find(ligne,"name")==1: exp[1]=string.split(ligne,":",1)[1][:-1] elif string.find(ligne,"date")==1: exp[3]=string.split(ligne,":",1)[1][:-1] else: mots=string.split(ligne) E.append(float(mots[0])) I.append(float(mots[1])) TE.append(E) TI.append(I) # print exp # for I in TI: print len(I) # for E in TE: print len(E) # print TI return TE,TI,exp def sortieD(entree,E,J): print "sortie de "+entree[:-4]+"-Deriv.dat" volt=open(entree+"-Deriv.dat","w") volt.write("#Potentiel E (mV) / dif2rentiell du Courant I (mA/V)\n") for i in range(len(E)): volt.write(str(E[i])+" "+str(J[i])+"\n") volt.close() def Peche(J,param,nok): K=[] for j in range(len(J)): if J[j]>param: k=J.index(max(J[j-20:j+20])) #k : position du pic principal d'un massif if K==[]:K.append(k) elif k>K[-1]: K.append(k) L=[] D=[0] for k in K: j=min(J[k-20:k+20]) if j<-0.5*J[k]: l=J.index(j) D.append(abs(k-l)), L.append((k+l)/2) # if len(D)<1 : # print "!!!! Pas de pêche détectée : essayez de baisser la limite !!!" # sys.exit(0) delta=max(D) K=[] for l in L: if abs(l-nok)<20:continue if K==[]:K.append(l) elif l>K[-1]+delta: K.append(l) return K,delta def odd(i): if i % 2 == 0 : p=1 else : p=-1 return(p) def TippEx(E,I,exp,limit): J=Deriv(exp,I) # sortieD(exp[4],E,J) nok=E.index(exp[6]) P,delta=Peche(J,limit,nok) print P,delta for p in P: # if abs(E[p]-exp[6])<5: # print "on ne fait pas de correction à moine de 5mV du bord" # continue # print "marche pas" a,b,s1=linreg(range(-3*delta,-delta),I[p-3*delta:p-delta]) i1=b-a*delta a,b,s2=linreg(range(delta,3*delta),I[p+delta:p+3*delta]) i2=a*delta+b i=0.5*(i1+i2) a=0.5*(i2-i1)/delta s=0.5*(sqrt(s1)+sqrt(s2)) print s for q in range(-delta,+delta): I[p+q]=i+a*q+s*odd(q) return I ###################################################### ##### Module de sortie ###################################################### def sortie(exp,E,I): print "sortie" volt=open(exp[4]+".dat","w") volt.write("#user: "+exp[0]+"\n") volt.write("#name: "+exp[1]+"\n") volt.write("#exp : "+str(exp[2])+"\n") volt.write("#date :"+exp[3]+"\n") volt.write("#Eini: "+str(exp[5])+"\n") volt.write("#Emax: "+str(exp[6])+"\n") volt.write("#Efin: "+str(exp[7])+"\n") volt.write("#mV/s: "+str(exp[8])+"\n") volt.write('#NbCy: '+str(exp[9])+"\n") volt.write('#Liss: '+str(exp[10])+"\n") volt.write("#nPt : "+str(len(E))+"\n") volt.write("#Premier Cycle \n") volt.write("#Potentiel E (mV) / Courant I (µA)\n") for i in range(len(E)): volt.write(str(E[i])+" "+str(I[i])+"\n") volt.close() mgr=open(exp[4]+".agr","w") tmpl=open("/usr/local/share/volta/template.grace","r") for ligne in tmpl.readlines(): mgr.write(ligne) tmpl.close() xmin=min(E) ymin=min(I) ymax=max(I) xmax=max(E) mgr.write('@ title "'+string.upper(exp[0][0])+exp[0][1:]+', '+exp[1]+', n° '+str(exp[2])+'"\n') mgr.write('@ subtitle "'+exp[3]+'"\n') mgr.write('@ s0 legend "'+str(exp[8])+'mV/s"\n') mgr.write('@ world xmin '+str(xmin)+'\n') mgr.write('@ world xmax '+str(xmax)+'\n') mgr.write('@ world ymin '+str(ymin)+'\n') mgr.write('@ world ymax '+str(ymax)+'\n') mgr.write('@target G0.S0 \n@type xy \n') for i in range(len(E)): mgr.write(str(E[i])+" "+str(I[i])+"\n") mgr.write('&') mgr.close() def sortie2(exp,E,I,cy): print "sortie n°"+str(cy+1) volt=open(exp[4]+".dat","a") volt.write("#Cycle : "+str(cy+1)+"\n") volt.write("#Potentiel E (mV) / Courant I (µA)\n") for i in range(len(E)): volt.write(str(E[i])+" "+str(I[i])+"\n") volt.close() mgr=open(exp[4]+".agr","a") mgr.write('@target G0.S'+str(cy)+' \n@type xy \n') for i in range(len(E)): mgr.write(str(E[i])+" "+str(I[i])+"\n") mgr.write('&') mgr.close() ######################################## ##Début de programme principale ######## ######################################## if __name__=="__main__": """programme principale : on lit un nom de fichier, on appelle """ if len(sys.argv)<2: for i in dircache.listdir('.'): print i entree=raw_input("Entrer le nom du fichier contenant la volta .dat : ") else: entree=sys.argv[1] TE,TI,exp=LisFich(entree) print "Vitesse de balayage (en mV/s) : "+str(exp[8]) print "Nombre de cycles : %g"%(exp[9],len(TE),len(TI))