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

runmin = 1
runmax = 2

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

	print '   '
	print qqq
	print '   '

	logmin = 50
	logmax = 99
	iii = logmin
##	st=ms.star_log('./falk2007_doublef_st_'+str(qqq)+'/LOGS/')
##	st=ms.mesa_profile('./LOGS/',iii)
	C12=[]
	C12ck=[]
	He4=[]
	He4ck=[]
	O16=[]
	O16ck=[]
	H1=[]
	N14=[]
	zone=[]

#	for iii in range(logmin,logmax):
		#	f2+str(iii)=open('/LOGS/log'+str(iii)+'.data','r')
	jonesmod=str(raw_input("what log number?"))
	apro=open('./LOGS/log'+str(jonesmod)+'.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,98):
		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 

	jjj=0
	lin=1
	n=0

	for lin in range(len(all_file)-6):
	#if all_filef2+str(ttt)[stringf2+str(ttt)+6][0] != 'zone':
		H1.append(float(all_file[lin+6].split()[find_column['h1']]))
		He4ck.append(float(all_file[lin+6].split()[find_column['he4']]))
		C12ck.append(float(all_file[lin+6].split()[find_column['c12']]))
		O16ck.append(float(all_file[lin+6].split()[find_column['o16']]))
		N14.append(float(all_file[lin+6].split()[find_column['n14']]))

		if N14[lin] < 0.0000000000001:
			if He4ck[lin] > 0.10:
				if C12ck[lin] > 0.10:
					#if  He4ck[lin] > (1.1*(O16ck[lin])):
					if (O16ck[lin]+C12ck[lin]+He4ck[lin]) > 0.90:
						jjj=jjj+1
						if jjj== 1.:
							He4.append(float(all_file[lin+6].split()[find_column['he4']]))
							zone.append(float(all_file[lin+6].split()[find_column['zone']]))
							C12.append(float(all_file[lin+6].split()[find_column['c12']]))
							zone.append(float(all_file[lin+6].split()[find_column['zone']]))
							O16.append(float(all_file[lin+6].split()[find_column['o16']]))
							zone.append(float(all_file[lin+6].split()[find_column['zone']]))
							print lin
							print lin, 'He4: ', He4[0] 
							print lin, 'C12: ', C12[0] 
							print lin, 'O16: ', O16[0] 

#	n=n+1


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

	plt.figure()
	tps=[1,2,3,4,5,6,7,8,9,10,11,12,13]
	plt.semilogy(tps, O16, 'r-.', label='O16', lw=1.8)
	plt.semilogy(tps, C12, 'g--', label='C12', lw=1.8)
	plt.semilogy(tps, He4, 'b-', label='He4', lw=1.8)
##			plt.axis([0.6329, 0.634, 0.00001, 1])
	plt.legend(loc=2)
	plt.ylabel('Mass Fraction', fontsize=32)
	plt.xlabel('TDU event number', fontsize=32)
	plt.title('Initial mass: 3Msun; Initial Z: 0.02', fontsize=32)
	size=24
	plt.xticks(size=size)
	plt.yticks(size=size)