import numpy as np
import matplotlib.pyplot as plt
import mesa as ms

runmin = 1
runmax = 2

qqq = runmin
plot=1 ### 0=NO; 1=YES
 
for qqq in range(runmin,runmax):

	print '   '
	print qqq
	print '   '

	logmin = 66
	logmax = 70

	iii = logmin
##	st=ms.star_log('./falk2007_doublef_st_'+str(qqq)+'/LOGS/')
	st=ms.star_log('./LOGS/')

	for iii in range(logmin,logmax):
		#	f2+str(iii)=open('/LOGS/log'+str(iii)+'.data','r')
		apro=open('./LOGS/log'+str(iii)+'.data','r')

		all_file=apro.readlines()

		apro.close()

		first_string=all_file[5].split()

		column_name=[]

		column_number=[]

		a=1
		for a in range(1,89):
			column_name.append(first_string[a-1])
			column_number.append(a-1)

		find_column={}
		for aa,bb in zip(column_name,column_number):
			find_column[aa]=bb 

		Mass=[]
		C13=[]
		C12=[]
		He4=[]
		O16=[]
		N14=[]
		H1=[]
		Dmlt=[]
		Dmix=[]

		C13eff=[]
		C13mc=[]
		F20=[]
		C13limit=[]
		C13um=[]

		n=0
		jjj=0
		kkk=0
		C13tot=0
		N14tot=0
		Mstart=0
		Mfinal=0
		string=0

		for string in range(len(all_file)-6):
	#if all_filef2+str(ttt)[stringf2+str(ttt)+6][0] != 'zone':
			Mass.append(float(all_file[string+6].split()[find_column['mass']]))
			C13.append(float(all_file[string+6].split()[find_column['c13']]))
			N14.append(float(all_file[string+6].split()[find_column['n14']]))

			if C13[n] > 0.001:
				if C13[n] > N14[n]:
					jjj=jjj+1
					if jjj== 1.:
						Mstart=(float(all_file[string+6].split()[find_column['mass']]))
			#				C13ef2+str(ttt)=C13ef2+str(ttt)+(C13f2+str(ttt)[n]-N14f2+str(ttt)[n])
			if jjj > 0:
				kkk=kkk+1
				if jjj == kkk: 
					C13tot=C13tot+((C13[n])*((float(all_file[string+5].split()[find_column['mass']]))-(float(all_file[string+6].split()[find_column['mass']]))))
					N14tot=N14tot+((N14[n])*((float(all_file[string+5].split()[find_column['mass']]))-(float(all_file[string+6].split()[find_column['mass']]))))
				if kkk == (jjj+1):
					Mfinal=(float(all_file[string+6].split()[find_column['mass']]))
			C12.append(float(all_file[string+6].split()[find_column['c12']]))
			O16.append(float(all_file[string+6].split()[find_column['o16']]))
			He4.append(float(all_file[string+6].split()[find_column['he4']]))
			H1.append(float(all_file[string+6].split()[find_column['h1']]))
			Dmlt.append(float(all_file[string+6].split()[find_column['log_mlt_D_mix']]))
			Dmix.append(float(all_file[string+6].split()[find_column['log_D_mix']]))
			surface_c12=st.get('surface_c12')
			surface_o16=st.get('surface_o16')
			COratio=(surface_c12*4.)/(surface_o16*3.)
			n=n+1

			if C13tot == 0.0:
				C13tot = 0.000000009
				N14tot = 0.0
				Mstart=0.0000000002
				Mfinal=0.0000000001
				
		print 'Massa C13 efficent log'+str(iii)+': ', (C13tot-N14tot)
		print 'Massa C13 pocket log'+str(iii)+' profile: ', (Mstart - Mfinal)

################ !!! F O R   P L O T T I N G !!! #####################################################################
		if plot == 1:

			plt.figure()
			plt.semilogy(Mass, C13, 'r-', label='C13', lw=1.8)
			plt.semilogy(Mass, N14, 'g-', label='N14', lw=1.8)
			plt.semilogy(Mass, H1, 'b-', label='H1', lw=1.8)
			plt.semilogy(Mass, O16, 'r:', label='O16', lw=1.8)
			plt.semilogy(Mass, C12, 'c-', label='C12', lw=1.8)
			plt.semilogy(Mass, He4, 'm-', label='He4', lw=1.8)
			plt.axis([0.575, 0.581, 0.00001, 1])
			mod_number = 57000 + (26*iii)
			plt.legend(loc=3,prop={'size':30})
			plt.ylabel('Mass Fraction', fontsize=40)
			plt.xlabel('Mass Coordinate (M$_\odot$)', fontsize=40)
			#plt.title('Model '+str(mod_number)+'; D20=10$^{10}$; f2=0.32', fontsize=32)
			size=30
			plt.xticks(size=size)
			plt.yticks(size=size)
			plt.twinx()
			plt.plot(Mass, Dmlt, 'k-', label='Dmlt', lw=1.8)
			plt.plot(Mass, Dmix, 'k:', label='Dmix', lw=1.8)
			plt.legend()
			size=24
			yticks = np.arange(4,17,1)
			plt.yticks(yticks)
			plt.xticks(size=size)
			plt.yticks(size=size)
			plt.ylabel('Log(Dmix) [$cm^{2}$/s]', fontsize=32)
			plt.axis([0.6105, 0.6115, 4., 17.])
			#plt.axis([0.575, 0.581, 0.00001, 1])
			#plt.savefig('./C13poc/C13form_mod'+str(iii)+'.png')


	print 'COratio growth: ', ((COratio[len(COratio)-1]/(COratio[0]) - 1.)*100), '%'