#For plotting the HRDs

import nugrid_set as set1
reload(set1)

mppnp_dir=['M4.000Z0.0001','M5.000Z0.0001','M6.000Z0.0001','M7.000Z0.0001']
mppnp_dir=['M4.000Z0.0010','M5.000Z0.0010','M6.000Z0.0010','M7.000Z0.0010']
mppnp_dir=['M4.000Z0.0060','M5.000Z0.0060','M6.000Z0.0060','M7.000Z0.0060']

#mppnp_dir=['M4.000Z0.0001','M5.000Z0.0001','M6.000Z0.0001','M7.000Z0.0001','M4.000Z0.0060','M5.000Z0.0060','M6.000Z0.0060','M7.000Z0.0060','M4.000Z0.0010','M5.000Z0.0010','M6.000Z0.0010','M7.000Z0.0010']

mppnp_dir=['M4.000Z0.0010','M5.000Z0.0010','M6.000Z0.0010','M7.000Z0.0010']

mppnp_dir=['M4.000Z0.0001']

mesa_dir=mppnp_dir
path_data='/apod2/NuGrid/data/set1ext'
species_mppnp=['C-12','C-13','N-14','O-16']
species_mesa=['c12','c13','n14','o16']


xaxis='age'
#xaxis='model'

if xaxis=='model':
	t0_model=0
	startfirstTP=False
else:
	t0_model=4000
	startfirstTP=True

for k in range(len(mppnp_dir)):

        if 'Z0.0001' in mppnp_dir[k]:
                path=path_data+'/set1.5a/'
        if 'Z0.0010' in mppnp_dir[k]:
                path=path_data+'/set1.4a/'
        if 'Z0.0060' in mppnp_dir[k]:
                path=path_data+'/set1.3a/'

	fig = mppnp_dir[k]
	plt.figure(fig)
	#for set1.1
	matplotlib.rc('xtick', labelsize=16)
	matplotlib.rc('ytick', labelsize=16)
	plt.minorticks_on()

	mesaset=set1.mesa_set(path+'/see_wind',multi_dir=[mesa_dir[k]])
	symbs=[None,None,None,None,None]
        end_model=20*[-1]
        setse=set1.mppnp_set(path+'/ppd_wind',multi_dir=[mppnp_dir[k]])
        setse.set_plot_surface_abu(fig=fig,species=species_mppnp,xaxis=xaxis,ratio=False,samefigureall=True,samefigure=True,sparsity=20,linestyle=[':','-','-.','--',':']*5,marker=symbs,color=['r','b','g','k','m']*5,markevery=100,t0_model=t0_model,age_years=True)

        symbs=['^','s','D','x','<','*','p','x']
        end_model=20*[-1]
        linestyle=20*[''] #,'-','-.',':','-','-.','--',':']
        #mesaset.set_plot_surface_abu(num_frame=fig,species=species_mesa,ratio=False, xax='model', t0_model=0,log=False,20*label=['standard'],marker=symbs,linestyle=linestyle,markevery=200)

	mesaset.set_plot_surface_abu_prof(num_frame=fig,isos=species_mesa,xax='model',label=20*['Mesa'],marker=symbs,linestyle=linestyle,markevery=1,markersize=6,color=['r','b','g','k','m'],xaxis=xaxis,startfirstTP=startfirstTP)

 

        ax=plt.gca()
        #ax.invert_xaxis()

	if False:
		ax.lines[0].set_label('$^{12}$C, nested')
		ax.lines[1].set_label('$^{14}$N, nested')
		ax.lines[2].set_label('$^{16}$O, nested')
		ax.lines[3].set_label('$^{12}$C, MESA')
		ax.lines[4].set_label('$^{14}$N, MESA')
		ax.lines[5].set_label('$^{16}$O, MESA')
	else:
		#ax.lines[0].set_label('$^{7}$Li, nested')
		ax.lines[0].set_label('$^{12}$C, nested')
		ax.lines[1].set_label('$^{13}$C, nested')
		ax.lines[2].set_label('$^{14}$N, nested')
		ax.lines[3].set_label('$^{16}$O, nested')
		#ax.lines[4].set_label('$^{7}$Li, MESA')
		ax.lines[4].set_label('$^{12}$C, MESA')
		ax.lines[5].set_label('$^{13}$C, MESA')
		ax.lines[6].set_label('$^{14}$N, MESA')
		ax.lines[7].set_label('$^{16}$O, MESA')


        #plt.title(mppnp_dir[k])
        ll=plt.legend()
	plt.legend(loc=4)
	if mppnp_dir[k] == 'M7.000Z0.0001':
		plt.ylim(1e-7,1e-2)

ax.set_ylabel('surface X$_i$')
lines=ll.get_lines()
ax2=ax.twinx()
if False:
	ax2.plot([],[],markersize=8,color='r',linestyle='',marker='^',label='$^{12}$C')
	ax2.plot([],[],markersize=8,color='b',linestyle='',marker='s',label='$^{14}$N')
	ax2.plot([],[],markersize=8,color='g',linestyle='',marker='D',label='$^{16}$O')
else:

	#ax2.plot([],[],markersize=8,color='r',linestyle='',marker='^',label='$^{7}$Li')
	ax2.plot([],[],markersize=8,color='r',linestyle='',marker='s',label='$^{12}$C')
	ax2.plot([],[],markersize=8,color='b',linestyle='',marker='D',label='$^{13}$C')
	ax2.plot([],[],markersize=8,color='g',linestyle='',marker='x',label='$^{14}$N')
	ax2.plot([],[],markersize=8,color='k',linestyle='',marker='<',label='$^{16}$O')


ax2.legend(loc=2,numpoints=1,title='MESA',prop={'size':16})
if False:
	ax.legend((lines[0],lines[1],lines[2]),('$^{12}$C','$^{14}$N','$^{16}$O'),title='nested',loc=4,prop={'size':16})
else:
	ax.legend((lines[0],lines[1],lines[2],lines[3]),('$^{12}$C','$^{13}$C','$^{14}$N','$^{16}$O'),title='nested',loc=4,prop={'size':16})


ax2.set_yscale('log');ax2.set_ylim(1e-6,1e-2)
ax2.yaxis.set_ticklabels([])


plt.xlim(0,1e5)

plt.minorticks_on()

plt.xlim(0,100000)
fig=plt.gcf()

if xaxis=='age':
	ax.set_xlabel('t-t$_0$ [yr]',fontsize=20,labelpad=1)
else:
	ax.set_xlabel('model number',fontsize=20,labelpad=1)


matplotlib.rc('xtick', labelsize=16)
matplotlib.rc('xtick', labelsize=12)
matplotlib.rc('xtick', labelsize=10)
ax.set_ylabel('surface $X_i$',fontsize=20,labelpad=1)
#fig.set_size_inches(10,8,forward=True)
ax.yaxis.set_tick_params(width=2,size=8)
ax.yaxis.set_tick_params(width=2,size=8)
ax.xaxis.set_tick_params(width=2,size=8)
plt.minorticks_on()
ax.axes.tick_params(which='minor',length='5')
fig.savefig('hbb_surface_abundance_M4Z0p0001.pdf')