from pylab import *
import pyfits
from PyAstronomy import pyasl
import scipy
from scipy import interpolate
from scipy import ndimage
from scipy import signal
import pickle
from matplotlib.backends.backend_pdf import PdfPages
import os
#from pyevolve import G1DList
#from pyevolve import GSimpleGA

from multiprocessing import Pool
import time

def download_models(webpage='http://svo2.cab.inta-csic.es/theory/models/coelho/high/data/',dest='../../data/'):
	
	os.system('mkdir '+dest+'/COELHO2014')
	cwd = os.getcwd()
	os.chdir(dest+'/COELHO2014')
	tf = np.arange(6000,10001,250)
	gf = np.arange(2.5,4.6,0.5)
	#gf = np.array([2.5])
	zf = np.array([-1.,-0.5,0.0,0.2])

	for t in tf:
		for g in gf:
			for z in zf:
				modname = get_modname(t,g,z)
				if z<0:
					sz = 'm'
				else:
					sz = 'p'
				sz = sz+str(float(np.absolute(z))).replace('.','')+'p00/'
				os.system('wget ' + webpage+sz+modname+'.fits')
				os.system('wget ' + webpage+sz+modname+'plc.fits')
	os.chdir(cwd)

	return True

def n_Edlen(l):
    sigma = 1e4 / l
    sigma2 = sigma*sigma
    n = 1 + 1e-8 * (8342.13 + 2406030 / (130-sigma2) + 15997/(38.9-sigma2))
    return n

def n_Morton(l):
    sigma = 1e4 / l
    sigma2 = sigma*sigma
    n = 1 + 6.4328e-5 + 2.94981e-2 / (146.-sigma2) + 2.5540e-4/(41.-sigma2)
    return n

def ToAir(l):
    return (l / n_Edlen(l))

def ToVacuum(l):
    cond = 1
    l_prev = l.copy()
    while(cond):
        l_new = n_Edlen(l_prev) * l
        if (max(np.absolute(l_new - l_prev)) < 1e-10): cond = 0
        l_prev = l_new
    return l_prev

def get_modname(t,g,z):
	st = str(int(t))
	if t<10000:
		st = '0'+st
	sg = '+'+str(np.around(g,1))
	if z < 0:
		sz = 'm'
	else:
		sz = 'p'
	z=float(z)
	sz = sz + str(np.around(np.absolute(z),1))
	sz = sz.replace('.','')
	return 't'+st+'_g'+sg+'_'+sz+'p00_hr'

def get_model(t,g,z,model_path='../../data/COELHO2014/'):
	modname = model_path + get_modname(t,g,z)
	try:
		out = pyfits.getdata(modname+'.fits')
	except:
		out = pyfits.getdata(modname+'plc.fits')
	return out

def get_near(x,vec):

	if x == vec[0]:
		mmin = vec[0]
		mmax = vec[1]
	elif x == vec[-1]:
		mmin = vec[-2]
		mmax = vec[-1]
	else:
		tvec = vec - x
		In  = np.where(tvec < 0)[0]
		mmin = tvec[In].max() + x
		Ix = np.where(tvec >= 0)[0]
		mmax = tvec[Ix].min() + x
	return mmin,mmax

def trilinear_interpolation(t,g,z,model_path='../../data/COELHO2014/'):
	teffs = np.arange(6000,10001,250)
	loggs = np.arange(2.5,4.6,0.5)
	fehs  = np.array([-1.,-0.5,0.0,0.2])
	x0,x1 = get_near(t,teffs)
	y0,y1 = get_near(g,loggs)
	z0,z1 = get_near(z,fehs)
	xd = (t-x0)/(x1-x0)
	yd = (g-y0)/(y1-y0)
	zd = (z-z0)/(z1-z0)

	try:
		hd = pyfits.getheader(model_path+get_modname(x0,y0,z0)+'.fits')
	except:
		hd = pyfits.getheader(model_path+get_modname(x0,y0,z0)+'plc.fits')

	c000 = get_model(x0,y0,z0,model_path)
	c001 = get_model(x0,y0,z1,model_path)
	c010 = get_model(x0,y1,z0,model_path)
	c100 = get_model(x1,y0,z0,model_path)
	c110 = get_model(x1,y1,z0,model_path)
	c101 = get_model(x1,y0,z1,model_path)
	c011 = get_model(x0,y1,z1,model_path)
	c111 = get_model(x1,y1,z1,model_path)

	wav = np.arange(len(c111))*hd['CDELT1'] + hd['CRVAL1']

	c00 = c000*(1-xd) + c100*xd
	c01 = c001*(1-xd) + c101*xd
	c10 = c010*(1-xd) + c110*xd
	c11 = c011*(1-xd) + c111*xd

	c0 = c00*(1-yd) + c10*yd
	c1 = c01*(1-yd) + c11*yd

	c = c0*(1-zd) + c1*zd

	return wav,c 

def normalize_model(w,f):
	ow = w.copy()
	of = f.copy()
	#plot(w,f)
	while True:
		#medflts = scipy.signal.medfilt(f,1001)
		coef = np.polyfit(w,f,6)
		fited = np.polyval(coef,w)
		res = f - fited
		I = np.where(res > -np.sqrt(np.var(res)))[0]
		w,f = w[I],f[I]
		if len(w) < 0.3* len(ow):
			break
	#plot(ow,np.polyval(coef,ow))
	#show()
	return coef

def spec_ccf(sw,sf,mw,mf,vi,vf,dv):
	mf = mf -1
	mf = -mf
	#plot(mw,mf)
	tck = interpolate.splrep(mw,mf,k=1)
	v = vi
	retccf = []
	vels = []
	while v<=vf:
		swt = sw * (1 + v/299792.458)
		mft = interpolate.splev(swt,tck)
		#if v == 0:
		#	plot(swt,mft)
		#	plot(swt,sft)
		#	show()
		mft -= np.mean(mft)
		sft = sf - np.mean(sf)
		#sft = sf.copy()
		#print np.sum(mft**2),np.sum(sft**2)
		retccf.append(np.sum(mft*sft)/np.sqrt(np.sum(mft**2)*np.sum(sft**2)))
		vels.append(v)
		v+=dv

	return np.array(vels),np.array(retccf)

def ccf_fft(swt,sft,mwt,mft):
	mf = mft -1
	mf = -mf

	#plot(mw,mf)
	tck = interpolate.splrep(np.log(mwt),mf,k=1)

	sw = np.log(swt)
	tck2 = interpolate.splrep(sw,sft,k=1)

	nsw = np.linspace(sw[0], sw[-1], 5000)
	sf  =  interpolate.splev(nsw,tck2)

	mf  =  interpolate.splev(nsw,tck)

	sf -= np.mean(sf)
	mf -= np.mean(mf)
	plot(nsw,sf)
	plot(nsw,mf)
	show()
	retccf = np.fft.ifft(np.conj(np.fft.fft(sf))*np.fft.fft(mf))
	retccf = np.hstack((retccf[2500:],retccf[:2500]))
	retvels = np.arange(len(retccf)) - 0.5*len(retccf)
	retvels *= (nsw[1]-nsw[0])
	retvels = 299792.458*(np.exp(retvels)-1.)
	return retvels, retccf

def ccf_simple(sw,sf,mw,mf,rv):
	mf = mf -1
	mf = -mf
	#plot(mw,mf)
	tck = interpolate.splrep(mw,mf,k=1)
	swt = sw * (1 + rv/299792.458)
	mft = interpolate.splev(swt,tck)
	mft -= np.mean(mft)
	sft = sf - np.mean(sf)
	return np.sum(mft*sft)/np.sqrt(np.sum(mft**2)*np.sum(sft**2))

def clean_strong_lines(mw,sc,mode=1):
	if mode==1:
		#""""
		I = np.where((mw>6520)&(mw<6600))[0]
		sc[I] = 1.
		I = np.where((mw>5888)&(mw<5897))[0]
		sc[I] = 1.
		I = np.where((mw>4310)&(mw<4360))[0]
		sc[I] = 1.
		I = np.where((mw>4840)&(mw<4880))[0]
		sc[I] = 1.
		I = np.where((mw>4070)&(mw<4130))[0]
		sc[I] = 1.
		I = np.where((mw>3875)&(mw<3900))[0]
		sc[I] = 1.
		I = np.where((mw>3920)&(mw<3945))[0]
		sc[I] = 1.
		I = np.where((mw>3955)&(mw<3980))[0]
		sc[I] = 1.
		I = np.where(mw<3850)[0]
		sc[I] = 1.
		#"""
	if mode==2:
		#""""
		I = np.where((mw>6550)&(mw<6570))[0]
		sc[I] = 1.
		I = np.where((mw>5888)&(mw<5897))[0]
		sc[I] = 1.
		I = np.where((mw>4320)&(mw<4350))[0]
		sc[I] = 1.
		I = np.where((mw>4850)&(mw<4870))[0]
		sc[I] = 1.
		I = np.where((mw>4090)&(mw<4110))[0]
		sc[I] = 1.
		I = np.where((mw>3875)&(mw<3900))[0]
		sc[I] = 1.
		I = np.where((mw>3920)&(mw<3945))[0]
		sc[I] = 1.
		I = np.where((mw>3955)&(mw<3980))[0]
		sc[I] = 1.
		I = np.where(mw<3850)[0]
		sc[I] = 1.
		#"""
	return sc

def RVforFR(wavs,flxs,teff=6700,logg=4.0,feh=-1.0,vsini=100.,model_path='../../data/COELHO2014/',vmin=-1000.,vmax=1000.,vstep=10.):
	def fitfunc(p,x):
		ret = p[3] + p[0] * np.exp(-.5*((x-p[1])/p[2])**2)
		return ret
	errfunc = lambda p,x,y: np.ravel( (fitfunc(p,x)-y) )

	#sc = get_model(teff,logg,feh)
	#hd = pyfits.getheader(model_path+get_modname(7000,4.5,0.0)+'.fits')
	#wav = np.arange(len(sc))*hd['CDELT1'] + hd['CRVAL1']
	teff = float(teff)

	try:
		sc = get_model(teff,logg,feh)
		hd = pyfits.getheader(model_path+get_modname(7000,4.5,0.0)+'.fits')
		mw = np.arange(len(sc))*hd['CDELT1'] + hd['CRVAL1']
	except:
		mw,sc = trilinear_interpolation(teff,logg,feh,model_path)

	for order in range(len(flxs)):
		flxs[order] = clean_strong_lines(wavs[order],flxs[order])

	sc = clean_strong_lines(mw,sc)

	II = np.where(sc != 1)[0]
	JJ = np.where(sc == 1)[0]

	coef = normalize_model(mw[II],sc[II])
	sc  /= np.polyval(coef,mw)
	sc[JJ] = 1. 
	mw = ToVacuum(mw)
	weis1 = []
	ccftot = []
	for i in range(wavs.shape[0]):
		#plot(wavs[i],flxs[i])
		scf = flxs[i]
		scw = wavs[i]

		J = np.where(scf!=0)[0]
		scw,scf = scw[J],scf[J]
		I = np.where((mw>scw[0]-100) & (mw<scw[-1]+100))
		tmf = pyasl.fastRotBroad(mw[I], sc[I], 0.5, vsini)
		#plot(mw[I],tmf)
		J = np.where(scf!=1)[0]
		if len(J)>100:
			ccv,ccf = spec_ccf(scw,scf,mw[I],tmf,vmin,vmax,vstep)
			#plot(ccv,ccf)
			#show()
			#ccf = np.array(ccf)
			wei1 = len(np.where(scf!=1)[0])**2
			weis1.append(wei1)
			if len(ccftot)==0:
				ccftot = ccf.copy()*wei1
			else:
				ccftot = np.vstack((ccftot,ccf.copy()*wei1))
	#show()
	weis1 = np.array(weis1)
	ccftot = np.sum(ccftot,axis=0)/ np.sum(weis1)

	p0 = [ccftot.min(),ccv[np.argmin(ccftot)],vsini,ccftot[0]]
	p1, success = scipy.optimize.leastsq(errfunc,p0, args=(ccv,ccftot))

	return p1,ccv,ccftot,fitfunc(p1,ccv)
	
def calc_bss2(vels,xc,coef, bot_i=0.15, bot_f=0.4, top_i=0.6, top_f=0.9, dt=0.01):
	try:
		
		I1 = np.where((vels>coef[1]-3*coef[2]) & (vels<coef[1]) )[0]
		I2 = np.where((vels<coef[1]+3*coef[2]) & (vels>coef[1]) )[0]
		I3 = np.where(vels<coef[1]-4*coef[2])[0]
		I4 = np.where(vels>coef[1]+4*coef[2])[0]
		I = np.hstack((I3,I4))
		base = np.median(xc[I])

		xc = base - xc
		xc /= xc.max()


		v1,x1 = vels[I1],xc[I1]
		v2,x2 = vels[I2],xc[I2]
		#plot(v1,x1)
		#plot(v2,x2)
		#show()
		dp = top_f
		vect = []
		while dp >= top_i:
			lb = np.where(x1>dp)[0][0]
			m = (v1[lb] - v1[lb-1])/(x1[lb]-x1[lb-1])
			n = v1[lb] - m*x1[lb]
			bs1 = m*dp+n

			lb = np.where(x2>dp)[0][-1]
			m = (v2[lb] - v2[lb+1])/(x2[lb]-x2[lb+1])
			n = v2[lb] - m*x2[lb]
			bs2 = m*dp+n
			vect.append(0.5*(bs2+bs1))
			dp-=dt
		vect = np.array(vect)

		dp = bot_f
		vecb = []
		while dp >= bot_i:

			lb = np.where(x1>dp)[0][0]
			m = (v1[lb] - v1[lb-1])/(x1[lb]-x1[lb-1])
			n = v1[lb] - m*x1[lb]
			bs1 = m*dp+n

			lb = np.where(x2>dp)[0][-1]
			m = (v2[lb] - v2[lb+1])/(x2[lb]-x2[lb+1])
			n = v2[lb] - m*x2[lb]
			bs2 = m*dp+n
			vecb.append(0.5*(bs2+bs1))
			dp-=dt
		vecb = np.array(vecb)

		return np.median(vecb) - np.median(vect) 
	except:
		return -999.0
"""
def lnlike(theta, W, F, Ferr):
	mw,sc = trilinear_interpolation(int(theta[0]),theta[1],theta[2])
	sct  = clean_strong_lines(mw,sc.copy())
	#plot(mw,sc)
	#show()
	coef = normalize_model(mw,sct)
	sc /= np.polyval(coef,mw)

	#print gfd
	mw = ToVacuum(mw)
	mw *= 1 + theta[3]/299792.458

	totD,totM,totE = np.array([]),np.array([]),np.array([])
	for i in range(W.shape[0]):
		scf = F[i]
		scw = W[i]
		scfe = Ferr[i]

		J = np.where(scf!=0)[0]
		scw,scf,scfe = scw[J],scf[J],scfe[J]

		I = np.where((mw>scw[0]-10) & (mw<scw[-1]+10))
		tmf = pyasl.fastRotBroad(mw[I], sc[I], 0.5, theta[4])
		tck = interpolate.splrep(mw[I],tmf,k=1)
		tmf = interpolate.splev(scw,tck)

		tmf = clean_strong_lines(scw,tmf.copy())
		I = np.where(tmf!=1)[0]
		#plot(scw,tmf)
		#plot(scw[I],tmf[I])
		#plot(scw[I],scf[I])

		#show()
		#print gfd

		tmf = tmf[I]
		scf = scf[I]
		scfe = scfe[I]

		tmf /= np.sum(tmf)
		tsf = scf/np.sum(scf)
		tse = scfe*(np.sum(scf)**2)

		totD = np.hstack((totD,tsf))
		totM = np.hstack((totM,tmf))
		totE = np.hstack((totE,tse))
		#plot(scw[I],tsf)
		#plot(scw[I],tmf)
		#plot(scw[I],tsf + 1./np.sqrt(tse))

	#show()
	#print fds
	#print theta
	#show()
	#print gvfd
	#ret = -np.log(2*np.pi) + np.log(np.sum(np.exp(-0.5*((y-model)/yerr)**2)/yerr))
	#ret = -0.5*(np.sum(inv_sigma2*(F-model)**2 - np.log(inv_sigma2)))
	ret = -0.5*(np.sum(totE*(totD-totM)**2 - np.log(totE)))

	#for i in range(len(F)):
	#	errorbar(Y,F[i],yerr=Ferr[i],fmt='b')
	#for j in model:
	#	plot(Y,j,'r')
	#show()
	#print theta, ret
	if np.isnan(ret):
		return -np.inf
	else:
		return ret


def lnprior(theta):
	if 6000 < theta[0] < 9000 and 3.0 < theta[1] < 4.5 and -1 < theta[2] < 0.2 and -500 < theta[3] < 500 and  1. < theta[4] < 500.:
	    return 0.0
	return -np.inf

def lnprob(theta, W,F,Ferr):
	lp = lnprior(theta)
	if not np.isfinite(lp):
		return -np.inf
	return lp + lnlike(theta,W,F,Ferr)
"""
def multiccf(pars):
	teff,logg,feh,vsini=pars[0],pars[1],pars[2],pars[3]
	vmin=-500
	vmax=500.
	vstep=20.
	sc = get_model(teff,logg,feh)
	hd = pyfits.getheader(model_path+get_modname(7000,4.5,0.0)+'.fits')
	wav = np.arange(len(sc))*hd['CDELT1'] + hd['CRVAL1']
	try:
		sc = get_model(teff,logg,feh)
		hd = pyfits.getheader(model_path+get_modname(7000,4.5,0.0)+'.fits')
		mw = np.arange(len(sc))*hd['CDELT1'] + hd['CRVAL1']
	except:
		mw,sc = trilinear_interpolation(teff,logg,feh,model_path)

	sc = clean_strong_lines(mw,sc)

	II = np.where(sc != 1)[0]
	JJ = np.where(sc == 1)[0]

	coef = normalize_model(mw[II],sc[II])
	sc  /= np.polyval(coef,mw)
	sc[JJ] = 1. 
	mw = ToVacuum(mw)
	weis1 = []
	ccftot = []
	for i in range(wavs.shape[0]):
		scf = flxs[i].copy()
		scw = wavs[i].copy()

		J = np.where(scf!=0)[0]
		scw,scf = scw[J],scf[J]
		I = np.where((mw>scw[0]-100) & (mw<scw[-1]+100))
		tmf = pyasl.fastRotBroad(mw[I], sc[I], 0.5, vsini)
		#plot(mw[I],tmf)
		J = np.where(scf!=1)[0]
		if len(J)>100:
			ccv,ccf = spec_ccf(scw,scf,mw[I],tmf,vmin,vmax,vstep)
			#ccv,ccf = ccf_fft(scw,scf,mw[I],tmf)
			#plot(ccv,ccf)
			#show()

			wei1 = len(np.where(scf!=1)[0])**2
			weis1.append(wei1)
			if len(ccftot)==0:
				ccftot = ccf.copy()*wei1
			else:
				ccftot = np.vstack((ccftot,ccf.copy()*wei1))
	weis1 = np.array(weis1)
	ccftot = np.sum(ccftot,axis=0)/ np.sum(weis1)


	#print gfds
	#ccftot = np.mean(ccftot,axis=0)
	#print pars, ccftot.min()
	return ccftot.min()


def get_pars_fr(wavst,flxst,model_patht='../../data/COELHO2014/',npools=4,fixG=1.0):
	for order in range(len(flxst)):
		flxst[order] = clean_strong_lines(wavst[order],flxst[order],mode=1)

	t0 = time.time()

	global wavs,flxs
	global model_path

	wavs,flxs=wavst.copy(),flxst.copy()
	model_path=model_patht


	gt = np.array([6000,7000,8000,9000,10000])
	gg = np.array([2.5,3.0,3.5,4.0,4.5])
	if fixG != -1:
		gg = np.array([fixG])
	gz = np.array([-1,-0.5,0.0,0.2])
	gr = np.array([10.,50.,100.,150.,200.,250.,300.])

	#"""
	tr = np.tile(gr,len(gt)*len(gg)*len(gz))
	tg = np.repeat(np.tile(gg,len(gt)),len(gr)*len(gz))
	tz = np.repeat(np.tile(gz,len(gt)*len(gg)),len(gr))
	tt = np.repeat(gt,len(gg)*len(gr)*len(gz))
	tot = np.vstack((tt,tg,tz,tr)).T

	#for pars in tot:
	#	pars = [8000,4.0,-0.5,40.0]
	#	print pars, multiccf(pars)

	p = Pool(npools)
	vals = np.array((p.map(multiccf, list(tot))))
	p.terminate()
	I = np.argmin(vals)
	best_vals = tot[I]
	bt,bg,bz,br = best_vals[0],best_vals[1],best_vals[2],best_vals[3]
	#"""
	t1 = time.time()
	print bt,bg,bz,br, (t1-t0)/60.,'mins'


	#bt,bg,bz,br = 7000.,4.5, 0.2, 100.0
	gt = np.arange(bt-1000,bt+1001,250)
	I = np.where((gt>=6000) & (gt<=10000))[0]
	gt = gt[I]
	gr = np.arange(br-60.,br+61.,20.)
	I = np.where(gr>=10)[0]
	gr = gr[I]

	tr = np.tile(gr,len(gt)*len(gg)*len(gz))
	tg = np.repeat(np.tile(gg,len(gt)),len(gr)*len(gz))
	tz = np.repeat(np.tile(gz,len(gt)*len(gg)),len(gr))
	tt = np.repeat(gt,len(gg)*len(gr)*len(gz))
	tot = np.vstack((tt,tg,tz,tr)).T

	p = Pool(npools)
	vals = np.array((p.map(multiccf, list(tot))))
	p.terminate()
	I = np.argmin(vals)
	best_vals = tot[I]
	bt,bg,bz,br = best_vals[0],best_vals[1],best_vals[2],best_vals[3]
	t2 = time.time()
	print bt,bg,bz,br, (t2-t1)/60.,'mins'
	#np.savetxt('temp_grid.txt',vals)


	if fixG==-1:
		grid = np.reshape(vals,(len(gt),len(gg),len(gz),len(gr)))
		tckt = interpolate.splrep(gt,np.arange(len(gt)),k=1)
		tckg = interpolate.splrep(gg,np.arange(len(gg)),k=1)
		tckz = interpolate.splrep(gz,np.arange(len(gz)),k=1)
		tckr = interpolate.splrep(gr,np.arange(len(gr)),k=1)

		itckt = interpolate.splrep(np.arange(len(gt)),gt,k=1)
		itckg = interpolate.splrep(np.arange(len(gg)),gg,k=1)
		itckz = interpolate.splrep(np.arange(len(gz)),gz,k=1)
		itckr = interpolate.splrep(np.arange(len(gr)),gr,k=1)

		st = np.arange(gt[0],gt[-1]+1,10.)
		sg = np.arange(gg[0],gg[-1]+0.01,0.1)
		sz = np.arange(gz[0],gz[-1]+0.01,0.1)
		sr = np.arange(gr[0],gr[-1]+1.,5.)

		st = interpolate.splev(st,tckt)
		sg = interpolate.splev(sg,tckg)
		sz = interpolate.splev(sz,tckz)
		sr = interpolate.splev(sr,tckr)

		tr2 = np.tile(sr,len(st)*len(sg)*len(sz))
		tg2 = np.repeat(np.tile(sg,len(st)),len(sr)*len(sz))
		tz2 = np.repeat(np.tile(sz,len(st)*len(sg)),len(sr))
		tt2 = np.repeat(st,len(sg)*len(sr)*len(sz))
		tot2 = np.vstack((tt2,tg2,tz2,tr2))

		zi = ndimage.map_coordinates(grid, tot2, order=3, mode='nearest')
		I  = np.argmin(zi)
		minval = tot2[:,I]

		mint = interpolate.splev(minval[0],itckt)
		ming = interpolate.splev(minval[1],itckg)
		minz = interpolate.splev(minval[2],itckz)
		minr = interpolate.splev(minval[3],itckr)

	else:
		grid = np.reshape(vals,(len(gt),len(gz),len(gr)))
		tckt = interpolate.splrep(gt,np.arange(len(gt)),k=1)
		tckz = interpolate.splrep(gz,np.arange(len(gz)),k=1)
		tckr = interpolate.splrep(gr,np.arange(len(gr)),k=1)

		itckt = interpolate.splrep(np.arange(len(gt)),gt,k=1)
		itckz = interpolate.splrep(np.arange(len(gz)),gz,k=1)
		itckr = interpolate.splrep(np.arange(len(gr)),gr,k=1)

		st = np.arange(gt[0],gt[-1]+1,10.)
		sz = np.arange(gz[0],gz[-1]+0.01,0.1)
		sr = np.arange(gr[0],gr[-1]+1.,5.)

		st = interpolate.splev(st,tckt)
		sz = interpolate.splev(sz,tckz)
		sr = interpolate.splev(sr,tckr)

		tr2 = np.tile(sr,len(st)*len(sz))
		tz2 = np.repeat(np.tile(sz,len(st)),len(sr))
		tt2 = np.repeat(st,len(sr)*len(sz))
		tot2 = np.vstack((tt2,tz2,tr2))

		zi = ndimage.map_coordinates(grid, tot2, order=3, mode='nearest')
		I  = np.argmin(zi)
		minval = tot2[:,I]

		mint = interpolate.splev(minval[0],itckt)
		ming = fixG
		minz = interpolate.splev(minval[1],itckz)
		minr = interpolate.splev(minval[2],itckr)

	#d = {'grid':grid, 'zi':zi, 'tot2':tot2, 'gt':gt, 'gg':gg, 'gz':gz, 'gr':gr}
	#pickle.dump(d,open('temp_dict.pkl'))

	return float(mint),float(ming),float(minz),float(minr)

def plot_CCF_FR(xc_dict,path='XC.pdf'):

	vels       = xc_dict['vels']
	xc_av      = xc_dict['xc_av']
	XCmodelgau = xc_dict['XCmodelgau']
	#refvel     = xc_dict['refvel']
	p1gau      = xc_dict['p1gau']

	f1 = figure()

	pp = PdfPages(path)
	ax1 = f1.add_subplot(111)
	ax1.plot(vels, xc_av,'b.', label='CCF')
	ax1.plot(vels, XCmodelgau,'r-',label='Gaussian fit')

	xlabel('Velocity (km/s)')
	ylabel('XC')
	ax1.axvline(p1gau[1],linestyle=':',color='r')
	ax1.axhline(0.0,linestyle='-')
	title('Average Cross-Correlation Function + Fit')
	handles, labels = ax1.get_legend_handles_labels()
	ax1.legend(handles[::-1], labels[::-1],prop={'size':6})
	pp.savefig()
	pp.close()
	clf()
	pass

"""
def trans_chromosome(chromosome):
	teff = chromosome[0]*100.+chromosome[1]*10.+chromosome[2]
	m = (10000.- 6000.)/999.
	n = 6000. 
	teff = teff*m + n
	logg = chromosome[3] + chromosome[4]*0.1
	m = (4.5 - 3.0)/9.9
	n = 3. 
	logg = logg*m + n

	feh  = chromosome[5] + chromosome[6]*0.1
	m = (0.2 - -1.)/9.9
	n = -1.
	feh = feh*m + n

	vsini  = chromosome[7]*10. + chromosome[8]
	m = (300. - 10.)/99.
	n = 10.
	vsini = vsini*m + n

	return teff, logg, feh, vsini

global wavs, flxs

def find_pars_GA(wavs,flxs,model_path='../../data/COELHO2014/'):

	def eval_func(chromosome):
		print list(chromosome)
		teff, logg, feh, vsini = trans_chromosome(chromosome)
		print teff, logg, feh, vsini
		pt,vels,ccf,mod = RVforFR(wavs,flxs,teff=teff,logg=logg,feh=feh,vsini=vsini,model_path=model_path)
		score = -ccf.min()
		return score

	genome = G1DList.G1DList(9)
	genome.evaluator.set(eval_func)

	ga = GSimpleGA.GSimpleGA(genome, interactiveMode=True)
	ga.setGenerations(40)
	ga.setMutationRate(0.2)
	ga.setPopulationSize(20)
  	#ga.setCrossoverRate(1.0)
	genome.setParams(rangemin=0, rangemax=9)
	#ga.setMultiProcessing(True)
	ga.evolve(freq_stats=10)
	print ga.bestIndividual()
	print trans_chromosome(ga.bestIndividual())

"""



