#For plotting the HRDs

import nugrid_set as set1

path='/apod2/NuGrid/data/set1ext'

#mesaset=set1.mesa_set('/nfs/rpod3/critter/Results/DATA/Set1_extension/set1.4a/see_wind',multi_dir=['M1.000Z0.0010','M1.650Z0.0010','M2.000Z0.0010','M3.000Z0.0010','M4.000Z0.0010','M5.000Z0.0010','M6.000Z0.0010','M7.000Z0.0010'])#,extra_label=20*[''])
symbs=[]
end_model=20*[-1]
symbs=['o','s','D','+','<','*','p','x']
linestyle=['--','-','-.',':','-','-.','--',':']
#mesaset.set_plot_kipp(startfirstTP=False,xtime=False)


#mesaset=set1.mesa_set('/nfs/rpod3/critter/Results/DATA/Set1_extension/set1.5a/see_wind',multi_dir=['M1.000Z0.0001','M2.000Z0.0001','M5.000Z0.0001','M7.000Z0.0001'])#,extra_label=20*[''])

#mesaset=set1.mesa_set('/nfs/rpod3/critter/Results/DATA/Set1_extension/set1.5a/see_wind',multi_dir=['M3.000Z0.0001'])#,extra_label=20*[''])

#mesaset=set1.mesa_set('/nfs/rpod3/critter/Results/DATA/Set1_extension/set1.3a/see_wind',multi_dir=['M6.000Z0.0060','M7.000Z0.0060'])#,extra_label=20*[''])



#ixaxis='log_time_left'
#ixaxis='time'
#mesaset.set_plot_kipp_cont(fig=32,mix_burn_zones=20,modstart=0,modstop=-1,ixaxis='log_time_left',xres=10000,yres=10000,xlims=[0,0],ylims=[0,0],plot_radius=False,engenPlus=True,engenMinus=False,boundaries=True,savepath='')

#t0_model=4849 for M3



from imp import *
ms=load_source('mesa','/rpod2/critter/PPN/forum.astro.keele.ac.uk_update/utils/pylib/mesa.py')

#######M3
historydata=ms.history_data(path+'/set1.5a/see_wind/M3.000Z0.0001/LOGS')

'''
################version 1
#'model_number'
#'age'
historydata.kip_cont(ifig=1112, modstart=0, modstop=-1, \
                 xlims=[0.,0.], ylims=[0.,0.],outfile='', \
                 xres=10000, yres=20000, ixaxis='log_time_left', \
                 mix_zones=20, burn_zones=20, plot_radius=False, \
                 engenPlus=False, engenMinus=False, \
                 landscape_plot=False, rad_lines=False, profiles=[], \
                 showfig=True, outlines=True, boundaries=True, \
                 c12_boundary=False, rasterise=False, yscale='1.')
ax=plt.gca()
plt.legend( (ax.lines[0],ax.lines[1],ax.lines[2]),('M$_{\star}$','H-free core','He-free core'))
plt.legend(loc=1)
plt.xlabel('log(t$_{final}$ - t)[yr]')
plt.xlim(8.435083281747918,0)
plt.ylabel('$M/M_{\odot}$')
#plt.xlim(8.435083281747918,0)
#plt.xlim(0,100000)
#plt.xlabel('$t-t_0 [yr]$')
plt.xlabel('$log(t_{final}-t)\ [yr]$',fontsize=20,labelpad=1)
plt.legend().set_visible(False)
plt.tick_params(axis='both', which='major', labelsize=16)
plt.savefig('M3_set1_5a_kippcont_preAGB.png')

historydata.kippenhahn(xax='logtimerev',t0_model=0,num_frame=1,plot_star_mass=False)
plt.xlim(8.435083281747918,0)
plt.xlim(8.435083281747918,2.696)
plt.ylabel('$M/M_{\odot}$')
plt.ylabel('$M [M_{\odot]$')

plt.xlabel('$log(t_{final}-t)\ [yr]$',fontsize=20,labelpad=1)
plt.legend().set_visible(False)
plt.ylabel('$M$ $[M_{\odot}]$',fontsize=20,)
plt.legend(loc=5,fontsize=12)
plt.legend(loc=6,fontsize=12)
plt.legend(loc=6,fontsize=14)
plt.savefig('M3_set1_5a_kippcont_preAGB.png')


################### version 2 of figure

#'age'
historydata.kip_cont(ifig=1122, modstart=0, modstop=-1, \
		 t0_model = 4849, \
                 xlims=[0.,0.], ylims=[0.,0.],outfile='', \
                 xres=10000, yres=20000, ixaxis='age', \
                 mix_zones=20, burn_zones=20, plot_radius=False, \
                 engenPlus=False, engenMinus=False, \
                 landscape_plot=False, rad_lines=False, profiles=[], \
                 showfig=True, outlines=True, boundaries=True, \
                 c12_boundary=False, rasterise=False, yscale='1.')

plt.xlim(0,90000)
plt.ylabel('$M-0.836\ [10^{-3}M_{\odot}]$',fontsize=20,labelpad=1)
ax=plt.gca()
plt.tick_params(axis='both', which='major', labelsize=16)
ax.set_xticks([2e4,4e4,6e4,8e4])
plt.ylim(0.836,0.855)
ax.set_xticklabels([2,4,6,8])
ax.set_yticks([0.836,0.84,0.844,0.848,0.852])
ax.set_yticklabels([0,4,8,12,16])
plt.xlabel('$t-t_0\ [10^4yr]$',fontsize=20,labelpad=0)
plt.ylabel('$M-0.836\ [10^{-3}M_{\odot}]$',fontsize=20,labelpad=0)
plt.savefig('M3_set1_5a_kippcont_AGB.pdf')




'''

#important part 3M AGB phase

historydata.kippenhahn(xax='time',t0_model=4849,num_frame=1,plot_star_mass=False)
plt.ylabel('$M-0.836\ [10^{-3}M_{\odot}]$',fontsize=20,labelpad=1)
ax=plt.gca()
plt.tick_params(axis='both', which='major', labelsize=16)
ax.set_xticks([2e4,4e4,6e4,8e4])
#ax.set_xticks([2,4,6,8]) 
plt.ylim(0.836,0.855)
ax.set_xticklabels([2,4,6,8])
ax.set_yticks([0.836,0.84,0.844,0.848,0.852])
ax.set_yticklabels([0,4,8,12,16])
plt.xlabel('$t-t_0\ [10^4yr]$',fontsize=20,labelpad=0)
plt.ylabel('$M-0.836\ [10^{-3}M_{\odot}]$',fontsize=20,labelpad=0)
plt.legend().set_visible(False)
plt.xlim(0,90000)
plt.savefig('M3_set1_5a_kippcont_AGB.png',dpi=300) #changed from pdf