import mppnp as mp
import utils
import mesa as ms



logdir='/rpod2/athira/mesa/case11/LOGS'
p=ms.mesa_profile(logdir,100,num_type='model')
cycles_1=p.log_ind.keys()
cycles_1.sort()
m=mp.se('/rpod2/athira/mppnp/case11/H5_out')
#m=mp.se('/rpod2/athira/mppnp/RCB_0.7/H5_out')

#index=[2000]
matplotlib.rc('xtick',labelsize=15,color='k')
matplotlib.rc('ytick',labelsize=12)
matplotlib.rc('ytick.major',pad=15)
matplotlib.rc('axes',labelsize=15,labelcolor='b')

#lines=utils.symbol_list('lines2')
#specs0_l=['He-4']
#specs1a_l=['C-12','C-14','O-15','O-17','Ne-20']
#specs1b_l=['N-14','F-17','Ne-21','Ne-22','Mg-24','Mg-25']
#specs1a_l=['C-13','He-4']
#specs1c_l=['O-15','F-19']
#specs1a_l=['H-1','C-12','C-14','N-13','N-14','N-15','O-15','F-18','Ne-19','F-19']
#specs1b_l=['He-4','C-13','Ne-22','Mg-25']

specs2_l=['$Na$','$Mg$','$Al$','$Si$','$S$','$Ca$','$Ti$']
specs3_l=['$Fe$','$Ni$','$Zn$','$Y$','$Zr$','$Ba$','$La$']
#lines=['k-.','ro-','ks-','r--','k+:','rd-','k.-','rh-','k>-','r:']
lines=['b-.','ro-','bs-','k--','r+:','kd-','b.-','rh-','b>-','r:']
lines_n=['c-']
lines0=['k-.','r-.','b-.']
lines_n=['k-','r-','b-']
lines1a=['k-.','ko-','r-.','ro-.','b-.','bo-.']
lines1b=['ks-','k--','k+:','rs-','r--','r+:','bs-','b--','b+:']
lines1c=['k.-','kd-.','r.-','rd-.','b.-','bd-.']
index_1=arange(1,3000,50)

i=0;j=0;k=0
d=0
close(0)
close(10)
close(20)
close(30)
close(40)

for ind in index_1:
    neut_den_1=[]
    neut=m.get(ind,'Neutron-1')
    mass=m.get(ind,'mass')
    rho= m.get(ind,'rho')
    for a in range(len(neut)):
        res=(neut[a]*6*10**23*rho[a])
        neut_den_1.append(res)

    na23=m.get(ind,'Na-23');
    mg24=m.get(ind,'Mg-24');mg25=m.get(ind,'Mg-25');mg26=m.get(ind,'Mg-26')
    al27=m.get(ind,'Al-27')
    si28=m.get(ind,'Si-28');si29=m.get(ind,'Si-29');si30=m.get(ind,'Si-30')
    s32=m.get(ind,'S-32');s33=m.get(ind,'S-33');s34=m.get(ind,'S-34');s36=m.get(ind,'S-36')
    ca40=m.get(ind,'Ca-40');ca42=m.get(ind,'Ca-42');ca43=m.get(ind,'Ca-43');ca44=m.get(ind,'Ca-44');ca46=m.get(ind,'Ca-46');ca48=m.get(ind,'Ca-48')
    ti46=m.get(ind,'Ti-46');ti47=m.get(ind,'Ti-47');ti48=m.get(ind,'Ti-48');ti49=m.get(ind,'Ti-49');ti50=m.get(ind,'Ti-50')

    na=na23/23.0
    mg=mg24/24.0+mg25/25.0+mg26/26.0
    al=al27/27.0
    si=si28/28.0+si29/29.0+si30/30.0
    s=s32/32.0+s33/33.0+s34/34.0+s36/36.0
    ca=ca40/40.0+ca42/42.0+ca43/43.0+ca44/44.0+ca46/46.0+ca48/48.0
    ti=ti46/46.0+ti47/47.0+ti48/48.0+ti49/49.0+ti50/50.0

    fe56=m.get(ind,'Fe-56');fe57=m.get(ind,'Fe-57');fe58=m.get(ind,'Fe-58')
    ni58=m.get(ind,'Ni-58');ni60=m.get(ind,'Ni-60');ni61=m.get(ind,'Ni-61');ni62=m.get(ind,'Ni-62');ni64=m.get(ind,'Ni-64')
    zn64=m.get(ind,'Zn-64');zn66=m.get(ind,'Zn-66');zn67=m.get(ind,'Zn-67');zn68=m.get(ind,'Zn-68');zn70=m.get(ind,'Zn-70')
    y89=m.get(ind,'Y-89')
    zr90=m.get(ind,'Zr-90');zr91=m.get(ind,'Zr-91');zr92=m.get(ind,'Zr-92');zr94=m.get(ind,'Zr-94');zr96=m.get(ind,'Zr-96')
    ba130=m.get(ind,'Ba-130');ba132=m.get(ind,'Ba-132');ba134=m.get(ind,'Ba-134');ba135=m.get(ind,'Ba-135');ba136=m.get(ind,'Ba-136');ba137=m.get(ind,'Ba-137');ba138=m.get(ind,'Ba-138')
    la138=m.get(ind,'La-138');la139=m.get(ind,'La-139')
    
    fe=fe56/56.0+fe57/57.0+fe58/58.0
    ni=ni58/58.0+ni60/60.0+ni61/61.0+ni62/62.0+ni64/64.0
    zn=zn64/64.0+zn66/66.0+zn67/67.0+zn68/68.0+zn70/70.0
    y=y89/89.0
    zr=zr90/90.0+zr91/91.0+zr92/92.0+zr94/94.0+zr96/96.0
    ba=ba130/130.0+ba132/132.0+ba134/134.0+ba135/135.0+ba136/136.0+ba137/137.0+ba138/138.0
    la=la138/138.0+la139/139.0

    specs2=[na,mg,al,si,s,ca,ti]
    specs3=[fe,ni,zn,y,zr,ba,la]

    close(ind+20)
    ifig1=figure(ind+20,figsize=(10,7))
    ax1=ifig1.add_subplot(111)
    for a in range (len(specs2)):
        ax1.semilogy(mass,specs2[a],lines[a],linewidth=2,label=specs2_l[a],markersize=6.2,alpha=0.65,markevery=15)       
    ylim(1e-6,1e-4)
    box1 = ax1.get_position()
    ax1.set_position([box1.x0, box1.y0, box1.width * 0.90, box1.height])
    ax1.legend(loc='center', bbox_to_anchor=(1.14, 0.23),prop={'size':11.0})
    ax1.set_ylabel('Number fraction')
    ax1.set_xlabel('$\mathrm{Mr/M_{\odot}}$')
    ax2=ax1.twinx()
    ax2.plot(mass,log10(neut_den_1),'orange',linestyle='-',linewidth=1.2)
    ax2.set_ylabel('Log Neutron Density (/cc)')
    xlim(0.44,0.90)
    ylim(0,max(log10(neut_den_1)+1))
    savefig('mppnp_profile2_'+'0'+str(ind)+'.png')
    #savefig('mppnp_profile2_'+'0'+str(ind)+'.pdf')
    #savefig('mppnp_profile2_'+'0'+str(ind)+'.eps')
    close(ind+20)

    close(ind+30)
    ifig1=figure(ind+30,figsize=(10,7))
    ax1=ifig1.add_subplot(111)
    for a in range (len(specs2)):
        ax1.semilogy(mass,specs3[a],lines[a],linewidth=2,label=specs3_l[a],markersize=6.2,alpha=0.65,markevery=15)       
    #ylim(1e-6,1e-4)
    box1 = ax1.get_position()
    ax1.set_position([box1.x0, box1.y0, box1.width * 0.90, box1.height])
    ax1.legend(loc='center', bbox_to_anchor=(1.14, 0.23),prop={'size':11.0})
    ax1.set_ylabel('Number fraction')
    ax1.set_xlabel('$\mathrm{Mr/M_{\odot}}$')
    ax2=ax1.twinx()
    ax2.plot(mass,log10(neut_den_1),'orange',linestyle='-',linewidth=1.2)
    ax2.set_ylabel('Log Neutron Density (/cc)')
    xlim(0.44,0.90)
    ylim(0,max(log10(neut_den_1)+1))
    savefig('mppnp_profile3_'+'0'+str(ind)+'.png')
    #savefig('mppnp_profile3_'+'0'+str(ind)+'.pdf')
    #savefig('mppnp_profile3_'+'0'+str(ind)+'.eps')
    close(ind+30)


'''
    ifig1=figure(40,figsize=(10,7))
    ax1=ifig1.add_subplot(111)
    for a in range (len(specs1a)):
        #ax1.plot(mass,specs1a[a],lines1a[i], label=specs1a_l[a],linewidth=2,markersize=6.2,alpha=0.65,markevery=15)
        ax1.semilogy(mass,specs1a[a],lines0[i],linewidth=2,markersize=6.2,alpha=0.65,markevery=15)
        i+=1
    #ylim(1e-5,max(c13))
    xlim(0.45,0.75)
    box1 = ax1.get_position()
    #ax1.set_position([box1.x0, box1.y0, box1.width * 0.83, box1.height])
    #ax1.legend(loc='center', bbox_to_anchor=(1.17, 0.17),prop={'size':10.0})
    ax1.set_ylabel('$^{13}C$ Mass fraction')
    ax1.set_xlabel('$Mr/M_{\odot}$')
    ax2=ax1.twinx()
    ax2.plot(mass,log10(neut_den_1),lines_n[d],linewidth=1.2)
    #ylim(max(log10(neut_den_1))/3.0,max(log10(neut_den_1)+1))
    d+=1
    ylim(0,7)
    xlim(0.45,0.75)   
    #savefig('mppnp_profile1.2a_'+'0'+str(ind)+'.png')
    #savefig('mppnp_profile1.2a_'+'0'+str(ind)+'.pdf')
    #savefig('mppnp_profile1.2a_'+'0'+str(ind)+'.eps')
    #close(10)
    #close(20)
ax2.set_ylabel('Log Neutron Density (/cc)')


    ifig2=figure(20,figsize=(10,7))
    ax2=ifig2.add_subplot(111)
    for a in range (len(specs1b)):
        ax2.semilogy(mass,specs1b[a],lines1b[j], label=specs1b_l[a],linewidth=2,markersize=6.2,alpha=0.65,markevery=5)
        j+=1
    box2 = ax2.get_position()
    ax2.set_position([box2.x0, box2.y0, box2.width * 0.83, box2.height])
    ax2.legend(loc='center', bbox_to_anchor=(1.17, 0.17),prop={'size':10.0})
    ax2.set_ylabel('Mass fraction')
    ax2.set_xlabel('$Mr/M_{\odot}$')
    ax3=ax2.twinx()
    ax3.plot(mass,log10(neut_den_1),'c-',linewidth=1.2)
    ylim(max(log10(neut_den_1))/3.0,max(log10(neut_den_1)+1))
    xlim(0.46,0.58)   
'''
'''
specs0_l=['He-4']
specs1a_l=['C-12','N-14']
specs1b_l=['H-1','C-14','N-15']
specs1c_l=['O-15','F-19']
#specs1a_l=['H-1','C-12','C-14','N-13','N-14','N-15','O-15','F-18','Ne-19','F-19']
#specs1b_l=['He-4','C-13','Ne-22','Mg-25']

#specs1b_l=['N-15','O-16','O-18','F-19','Ne-22']
#specs2_l=['Na-23','Mg-24','Al-27','Si-28','S-32','Ca-40','Ti-48']
#specs3_l=['Fe-56','Ni-58','Zn-64','Y-89','Zr-90','Ba-138','La-139']
#lines=['k-.','ro-','ks-','r--','k+:','rd-','k.-','rh-','k>-','r:']
lines0=['k:','r:','b:']
lines_n=['k-','r-','b-']
lines1a=['k-.','ko-','r-.','ro-.','b-.','bo-.']
lines1b=['ks-','k--','k+:','rs-','r--','r+:','bs-','b--','b+:']
lines1c=['k.-','kd-.','r.-','rd-.','b.-','bd-.']
index_1=[300,500,700]

i=0;j=0;k=0
d=0
close(0)
close(10)
close(20)
close(30)



for ind in index_1:
    neut_den_1=[]
    neut=m.get(ind,'Neutron-1')
    mass=m.get(ind,'mass')
    rho= m.get(ind,'rho')
    for a in range(len(neut)):
        res=(neut[a]*6*10**23*rho[a])
        neut_den_1.append(res)
    h1=m.get(ind,'H-1');he4=m.get(ind,'He-4');c12=m.get(ind,'C-12');c13=m.get(ind,'C-13');c14=m.get(ind,'C-14');n13=m.get(ind,'N-13');n14=m.get(ind,'N-14');o15=m.get(ind,'O-15');
    n15=m.get(ind,'N-15');o16=m.get(ind,'O-16');o18=m.get(ind,'O-18');f18=m.get(ind,'F-18');f19=m.get(ind,'F-19');ne19=m.get(ind,'Ne-19');ne22=m.get(ind,'Ne-22')
    na23=m.get(ind,'Na-23');mg24=m.get(ind,'Mg-24');mg25=m.get(ind,'Mg-25');al27=m.get(ind,'Al-27');si28=m.get(ind,'Si-28');s32=m.get(ind,'S-32');ca40=m.get(ind,'Ca-40');ti48=m.get(ind,'Ti-48')
    fe56=m.get(ind,'Fe-56');ni58=m.get(ind,'Ni-58');zn64=m.get(ind,'Zn-64');y89=m.get(ind,'Y-89');zr94=m.get(ind,'Zr-94');ba138=m.get(ind,'Ba-138');la139=m.get(ind,'La-139')
    specs0=[he4]
    specs1a=[c12,n14]
    specs1b=[h1,c14,n15]
    specs1c=[o15,f19]

    ifig0=figure(0,figsize=(10,7))
    ax1=ifig0.add_subplot(111)
    for a in range (len(specs0)):
        ax1.plot(mass,specs0[a],lines0[d], label=specs0_l[a],linewidth=2,markersize=6.2,alpha=0.65,markevery=15)
    ylim(0.3,1.0)
    if ind==300:
        box1 = ax1.get_position()
        ax1.set_position([box1.x0, box1.y0, box1.width * 0.83, box1.height])
        ax1.legend(loc='center', bbox_to_anchor=(1.17, 0.17),prop={'size':10.0})
        ax1.set_ylabel('Mass fraction')
        ax1.set_xlabel('$Mr/M_{\odot}$')
    #ax2=ax1.twinx()
    #ax2.plot(mass,log10(neut_den_1),lines_n[d],linewidth=1.2)
    #ylim(max(log10(neut_den_1))/3.0,max(log10(neut_den_1)+1))
    d+=1

    ifig1=figure(10,figsize=(10,7))
    ax1=ifig1.add_subplot(111)
    for a in range (len(specs1a)):
        ax1.plot(mass,specs1a[a],lines1a[i], label=specs1a_l[a],linewidth=2,markersize=6.2,alpha=0.65,markevery=15)
        i+=1
    #ylim(1e-8,1.1)
    xlim(0.46,0.58)
    if ind==300:
        box1 = ax1.get_position()
        ax1.set_position([box1.x0, box1.y0, box1.width * 0.83, box1.height])
        ax1.legend(loc='center', bbox_to_anchor=(1.17, 0.17),prop={'size':10.0})
        ax1.set_ylabel('Mass fraction')
        ax1.set_xlabel('$Mr/M_{\odot}$')
    xlim(0.46,0.58)   
    #savefig('mppnp_profile1.2a_'+'0'+str(ind)+'.png')
    #savefig('mppnp_profile1.2a_'+'0'+str(ind)+'.pdf')
    #savefig('mppnp_profile1.2a_'+'0'+str(ind)+'.eps')
    #close(10)
    #close(20)
    ifig2=figure(20,figsize=(10,7))
    ax2=ifig2.add_subplot(111)
    for a in range (len(specs1b)):
        ax2.semilogy(mass,specs1b[a],lines1b[j], label=specs1b_l[a],linewidth=2,markersize=6.2,alpha=0.65,markevery=5)
        j+=1
    if ind==300:
        box2 = ax2.get_position()
        ax2.set_position([box2.x0, box2.y0, box2.width * 0.83, box2.height])
        ax2.legend(loc='center', bbox_to_anchor=(1.17, 0.17),prop={'size':10.0})
        ax2.set_ylabel('Mass fraction')
        ax2.set_xlabel('$Mr/M_{\odot}$')

    ifig3=figure(30,figsize=(10,7))
    ax3=ifig3.add_subplot(111)
    for a in range (len(specs1c)):
        ax3.semilogy(mass,specs1c[a],lines1c[k], label=specs1c_l[a],linewidth=2,markersize=6.2,alpha=0.65,markevery=5)
        k+=1
    if ind==300:
        box3 = ax3.get_position()
        ax3.set_position([box3.x0, box3.y0, box3.width * 0.83, box3.height])
        ax3.legend(loc='center', bbox_to_anchor=(1.17, 0.17),prop={'size':10.0})
        ax3.set_ylabel('Mass fraction')
        ax3.set_xlabel('$Mr/M_{\odot}$')
    #ylim(1e-8,1.1)
    xlim(0.46,0.58)   
    #ax2=ax1.twinx()
    #ax2.plot(mass,log10(neut_den_1),linewidth=1.2,color='c')
    #ylim(max(log10(neut_den_1))/3.0,max(log10(neut_den_1)+1))
    xlim(0.46,0.58)
    #savefig('mppnp_profile1.2b_'+'0'+str(ind)+'.png')
    #savefig('mppnp_profile1.2b_'+'0'+str(ind)+'.pdf')
    #savefig('mppnp_profile1.2b_'+'0'+str(ind)+'.eps')
    #close(ind+15)

'''