#!/usr/bin/python
# -*- Coding: Latin1 -*-
# Written my Jerome Kieffer (kieffer@crans.org) 2001-2004 
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

import os, os.path, GracePlot
from numarray import *

def ReadXY(file):
	"""Read a XY file and retuns the  table. """
	try:
		f=open(file,"r")
	except:
		print "le fichier %s n'existe visiblement pas"%file
		return [],[]
	X=[]
	Y=[]
	for ligne in f.readlines():
		try:
			mots=ligne.split()
			e=float(mots[0])
			i=float(mots[1])
			X.append(e)
			Y.append(i)
		except:
			print ligne
	return X,Y

def ReadFluo(file):
	"""Read a Fluorescence XY file and returns  L, I table. """
	if file[-3:].upper()!=".PS":
		file+=".PS"
	try:
		f=open(file,"r")
	except:
		print "le fichier %s n'existe visiblement pas"%file
		return [],[]
	L=[]
	I=[]
	start=False
	for ligne in f.readlines():
		if ligne.find("#DATA")==0:start=True
		if start:
			try:
				mots=ligne.split()
				L.append(float(mots[0]))
				I.append(float(mots[1]))
			except: 
				print ligne
	return L,I


	
def ReadVoltaMaster(file):
	"""Read a VoltaMaster file and retuns the E and I table. """
	if file[-4:].upper()!=".CRV":
		file+=".CRV"
	try:
		f=open(file,"r")
	except:
		print "le fichier %s n'existe visiblement pas"%file
		return [],[]
	E=[]
	I=[]
	for ligne in f.readlines():
		try:
			mots=ligne.split("	")
			e=float(mots[0])
			i=float(mots[1])
			E.append(e)
			I.append(i)
		except: 
			print ligne
	return E,I
	
def ReadSPEC(file):
	"""open a specter and return 2 arrays : lambda and the number of hits"""
	lo=[]
	cp=[]
	for ligne in open(file+".PRN","r").readlines(): 
#		print ligne
		x,y=ligne.split(",")
		x=float(x.strip())
		#on se limite au domaine du visible ...
		if x>290  and x<840 :
			lo.append(x)
			cp.append(float(y.strip()))
	return lo,cp

def LoadAllSpectra(dir="."):
	"""load all the spectra abs retruns an array of all of them."""
	spectra=[]
	list=[]
	for file in os.listdir(dir):
		head,tail=os.path.splitext(file)
		if tail==".PRN":
			if head[-1] in ["0","1","2","3","4","5","6","7","8","9"]: 
				base,num=SplitNum(head,"")
				list.append(int(num))
			else: continue
		list.sort()
#		list=range(65,100)
	for i in list:
		print "reading Spectre : %s #%s"%(base,i)#,type(openSPEC(base,i))
		spectra.append(ReadSPEC(base+str(i)))		
	return spectra		


def AllGrace(data):
	"""Display a graph with all the data in a Grace window"""
	p=GracePlot.GracePlot()
	datagrace=[]
	for i in data:
		datagrace.append(GracePlot.Data(x=i[0],y=i[1]))
	p.plot(datagrace)

def SaveAllGnuPlot(file,data):
	"""save all the data in the file in a gnuplot format"""
	if file[-4:]!=".dat": file+=".dat"
	mgr=open(file,"w")
	idx=1
	for i in data:
		print "Writing GnuPLot #",idx
        	mgr.write("\n")
        	for j in range(len(i[1])):
                	mgr.write("%s   %s      %s\n"%(i[0][j],idx,i[1][j]))
		idx+=1
	mgr.close()

def AllGnuPlot(fich,data):
	"""save all the data in the file in a gnuplot format then lauch the
	gnuplot program to display the 3D graph"""
	if fich[-4:]!=".dat": fich+=".dat"
	gpl=fich[:-4]+".plt"
	SaveAllGnuPlot(fich,data)
	print gpl
	f=open(gpl,"w")
	f.write("""set pm3d \n set xlabel "Longueur d'onde (nm)" \n set ylabel "Courbe" \n set zlabel "Variation d'Absorbance"\n splot "%s" with pm3d \n"""%fich)
	f.close()


def SplitNum(text,tail):
	"""split the text "abc123" in "abc"+123"""
	if text[-1] in ["0","1","2","3","4","5","6","7","8","9"]:
		tail=text[-1]+tail
		text=text[:-1]
#		print text,tail
		text,tail=SplitNum(text,tail)	
	return text,tail

def F2(a,b,c,d,e,f,g,h,i,j,k):
    """This function is the base for the Savitsky Golay 2nd order smoothing """	
#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 F4(a,b,c,d,e,f,g,h,i,j,k):
    """This function is the base for the Savitsky Golay 4th order smoothing """	
#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 lisse2(x):
    """Does the 2nd order savitsky Golay smoothing"""
    tx = ([ x[0] ] * 5)  + x + ([ x[-1] ] * 5 )
    rx = []
    for i in range(5,len(tx)-5):
		rx.append(apply(F2,tuple(tx[i-5:i+6])))
    return rx

def lisse4(x):
    """Does the 2nd order savitsky Golay smoothing"""
    tx = ([ x[0] ] * 5)  + x + ([ x[-1] ] * 5 )
    rx = []
    for i in range(5,len(tx)-5):
		rx.append(apply(F4,tuple(tx[i-5:i+6])))
    return rx
    
    
def sortie0(exp,param):
	"""this prints the headers of the file"""
	file=os.path.normpath("/%s/%s/%s"%(exp["Home"],exp["ExpName"],exp["ExpNumber"]))
	dir=os.path.dirname(file)
	if not os.path.isdir(dir): makedir(dir)
	volt=open(file+".dat","w")
	exp["Date"]=time.strftime("%a, %d %b %Y %H:%M", time.localtime())

	for i in exp:
		volt.write("#%s: %s \n"%(i,exp[i]))
	for i in param:
		volt.write("%"+i)
	volt.close()
	sortie0GRACE(file,exp)
	return file    

def InitGrace(file):
    """prints the headers of the XmGrace file"""
    mgr=open(file+".agr","w")
    tmpl=open(os.getenv("HOME")+"/spectroelectrochem/template.agr","r")
    for ligne in tmpl.readlines(): mgr.write(ligne)
    tmpl.close()
    mgr.write('&')
    mgr.close()
    
def WriteGrace(file,data,curve):

	mgr=open(file+".agr","a")
	print "Writing Curve #%s"%curve
	mgr.write('@target G0.S%s \n@type xy \n'%(curve))
#	mgr.write('@    s%s legend  "Qpic=%0.3fmC;Qsynth=%0.3fmC"\n'%(cy-1,Qp,Qs))
	for i in range(len(data[1])):
		mgr.write("%s	%s \n"%(data[0][i],data[1][i]))
	mgr.write('& \n')
	mgr.write('@    world xmin %i \n'%min(data[0]))
	mgr.write('@    world xmax %i \n'%max(data[0]))
	mgr.write('@    world ymin %i \n'%min(data[1]))
	mgr.write('@    world ymax %i \n'%max(data[1]))
	mgr.close()

def WriteGnuPlot(file,data,curve):
	"""just write a file to be read again with GnuPLot"""
	print "Writing GnuPLot #",curve
	if file[-4:]!=".dat":
		file+=".dat"
	mgr=open(file,"a")
	mgr.write("\n")
	for i in range(len(data[1])):
		mgr.write("%s	%s	%s\n"%(data[0][i],curve,data[1][i]))
	mgr.close()


def WriteSPEC(file,data):
	"""just write a file to be read again  using openSPEC"""
	mgr=open(file+".PRN","w")
	print "Writing Spec "
	for i in range(len(data[1])):
		mgr.write("%s ,	%s \n"%(data[0][i],data[1][i]))
	mgr.close()



def MeanCurves(curves):
	"""calculate the average on all the curves"""
	length=len(curves)
#	if finish<=0:finish=length+finish
#	print finish,start,len(curves)
	lo=curves[0][0]
	do=array([0*len(lo)])
	for i in curves:
		do=do+array(i[1])
	return lo,do/len(curves)
	
def Absorbance(data,blanc):
	"""calculate the absorbance known as log10(I0/I)"""
	return data[0],lisse2(list(log10(array(blanc[1])/array(data[1]))))
	
	
def ReducePoint(data,N=10):
	"""refuce the number of points : one every N nanometere"""
	lo,I=data
	NewLo=[]
	NewAbs=[]
	lowest=max(int(min(lo))+N/2,280) #glass cut-off at 280nm
	highest=int(max(lo))-N/2
	for i in range(lowest,highest,N):
		match=[]
		for idx in range(len(lo)):
				if lo[idx]>i-N/2 and lo[idx]<i+N/2:match.append(idx)
		idx=sum(match)/len(match)
		NewLo.append(idx)
		NewAbs.append(sum(I[min(match):max(match)])/(1.0*len(match)))	
#		print idx
	return NewLo,NewAbs

def SubsBlack(data,noir):
    """substract the black noise to the data, returns 1 if the signal is zero or less"""
    d=array(data[1])-array(noir[1])
    d1=[]
    for i in d:
	if i>0:
	    d1.append(i)
	else:
	    d1.append(1)
    return data[0],d1
    
	
if __name__=='__main__':
	"""main program : just realize a summary of all spectra in the directory"""
	a=LoadAllSpectra(".")
#	for i in a: print len(i[0]),len(i[1])
#	InitGrace("summary")
#	for i in range(len(a)):
#		WriteGrace("summary",a[i],i)		



#	InitGrace("moyenne")
#	b=MeanCurves(a,0,0)
#	WriteGrace("moyenne",b,1)
#	WriteGrace("moyenne",[b[0],lisse2(list(b[1]))],2)
#	WriteGrace("moyenne",[b[0],lisse4(list(b[1]))],4)

#	WriteSPEC("Blanc",[b[0],lisse4(list(b[1]))])


	blanc=ReducePoint(ReadSPEC("Blanc"))
	print type(blanc)	
#	InitGrace("Absorbance")
	open("Absorbance.dat","w")
	idx=65
	for i in a:
		b=Absorbance(i,blanc)
		print len(i),len(b)
				
		WriteGnuPlot("Absorbance",b,idx)
		idx+=1		
#		WriteGrace("Absorbance",b,i)			
#	os.popen("xmgrace Absorbance.agr &")