read_models_Mo.py

from nugridpy import nugridse as mp
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from collections import OrderedDict
import presolargrains as psg

## range of isotopic ratios
isotopes_x = ['Mo-92','Mo-96']
isotopes_y = ['Mo-94','Mo-96']
#isotopes_y = ['Mo-95','Mo-96']
#isotopes_y = ['Mo-97','Mo-96']
#isotopes_y = ['Mo-98','Mo-96']
#isotopes_y = ['Mo-100','Mo-96']

## References in data
# 'Alexander ApJ 1999', 'Alexander GCA 1993',
# 'Alexander GCA 1994', 'Alexander GCA 1995', 'Amari ApJ 1992',
# 'Amari ApJ 1999', 'Amari ApJ 2001', 'Amari ApJ 2001 (nova)',
# 'Amari unpublished (2003)', 'Amari unpublished 1992',
# 'Avila ApJ 2012', 'Barzyk MAPS 2007', 'Barzyk MAPS 2008',
# 'Besmehn GCA 2003', 'Besmehn Thesis 2001', 'Bose GCA 2012',
# 'Busemann E&P Ltrs 2009', 'Croat Astron. J. 2010',
# 'Floss MAPS 2012', 'Fujiya ApJ Letters 2011',
# 'Fujiya ApJ Letters 2013', 'Gao LPSC 1997', 'Gyngard LPSC 2006',
# 'Gyngard MetSoc 2006', 'Gyngard MetSoc 2010',
# 'Gyngard unpublished 2009', 'Heck ApJ 2007', 'Heck ApJ 2009',
# 'Henkel MAPS 2007', 'Hoppe ApJ 1994', 'Hoppe ApJ 1997',
# 'Hoppe ApJ 2001', 'Hoppe ApJ 2009', 'Hoppe ApJ 2010',
# 'Hoppe ApJ Letters 2012', 'Hoppe GCA 1996', 'Hoppe LPSC 1996',
# 'Hoppe MAPS 2000', 'Hoppe Science 1996', 'Hoppe unpublished',
# 'Hoppe unpublished ', 'Huss GCA 1997', 'Huss MAPS 2007',
# 'Hynes LPSC 2006', 'Hynes LPSC 2008', 'Ireland ApJ 1991',
# 'Jennings LPSC 2002', 'Lin ApJ 2002', 'Liu ApJ 2014',
# 'Lyon LPSC 2006', 'Marhas ApJ 2008', 'Marhas LPSC 2004',
# 'Marhas LPSC 2005', 'Marhas LPSC 2007', 'Marhas MAPS 2007',
# 'Marhas Met Soc 2005', 'Nguyen GCA 2018', 'Nicolussi GCA 1998',
# 'Nicolussi Phys Rev Ltrs 1998', 'Nicolussi Science 1997',
# 'Nittler ApJ 1996', 'Nittler ApJ 2005', 'Nittler GCA 2003',
# 'Nittler Met Soc 2006', 'Nittler Thesis 1996',
# 'Orthous-Daunay LPSC 2012', 'Pellin LPSC 2000', 'Pellin LPSC 2006',
# 'Reference', 'Savina GCA 2003', 'Savina Science 2004',
# 'Stadermann GCA 2006', 'Stroud Met Soc 2004', 'Trappitsch GCA 2018'
# 'Virag GCA 1992', 'Wopenka LPSC 2010', 'Xu Metsoc 2012', 
# 'Yada MAPS 2008','Zhao ApJ 2013', 'Zinner GCA 2007', 'Zinner LPSC 2010',
# 'Zinner MAPS 2003'
reference = 'Barzyk MAPS 2007' # insert here for dataset

## class of presolar grain
# '', 'AB', 'C', 'HM only', 'M', 'MS', 'N', 'Type', 'U', 'U/C', 'X',
# 'X ', 'Y', 'Z', 'uncl.', 'unclass' 
#clas_1 = 'unclass'
clas_4 = 'M'

#run1=mp.se('/run/media/s1785107/Seagate Backup Plus Drive/eddie3/NuclearPhysics/NuPPN/frames/mppnp/RUN_set1upd_m3z3m2_hCBM_TDUx4p3/H5_surf','surf.h5')
#run2=mp.se('/run/media/s1785107/Seagate Backup Plus Drive/eddie3/NuclearPhysics/NuPPN_Ne22_newer_update1/frames/mppnp/RUN_set1upd_m3z3m2_hCBM_TDUx4p3_newerNe22a_stdstd/H5_surf','surf.h5')
#run3=mp.se('/run/media/s1785107/Seagate Backup Plus Drive/eddie3/NuclearPhysics/NuPPN_Ne22_newer_update2/frames/mppnp/RUN_set1upd_m3z3m2_hCBM_TDUx4p3_newerNe22a_lowstd/H5_surf','surf.h5')
#run4=mp.se('/run/media/s1785107/Seagate Backup Plus Drive/eddie3/NuclearPhysics/NuPPN_Ne22_newer_update3/frames/mppnp/RUN_set1upd_m3z3m2_hCBM_TDUx4p3_newerNe22a_stdhi/H5_surf','surf.h5')
#run5=mp.se('/run/media/s1785107/Seagate Backup Plus Drive/eddie3/np/NuPPN_ll12/frames/mppnp/RUN_set1upd_m3z3m2_hCBM_TDUx4p3_95Zrll12/H5_surf','surf.h5')
#run6=mp.se('/run/media/s1785107/Seagate Backup Plus Drive/eddie3/np/NuPPN_ll23/frames/mppnp/RUN_set1upd_m3z3m2_hCBM_TDUx4p3_95Zrll23/H5_surf','surf.h5')

run5=mp.se('/media/ashley/Seagate Backup Plus Drive/eddie3/NuclearPhysics/NuPPN/frames/mppnp/RUN_set1upd_m3z2m2/H5_surf','surf.h5')
run6=mp.se('/home/ashley/Work/Yield_uncertainties/Pavel/H5_surf','surf.h5')

#run_mass=[run1,run2,run3,run4,run5,run6,run7,run8,run9,run10,run11,run12,run13]
run_mass= [run5,run6]#run1,run2,run3,run4]#,

metallicity_label=[#'m3z2m2','m3z2m2, $0.8x$ $^{96}Mo(n,\gamma)$ kad1','m3z2m2, kad1','m3z2m2, $0.6x$ $^{96}Mo(n,\gamma)$ kad1']
#,'m3z3m2, 0.8x $^{96}$Mo(n,$\gamma$) kad1',#'m3z3m2, $0.8x$ $^{96}Mo(n,\gamma)$ kad1',
#'m2z1m2','m2z2m2','m2z3m2','m3z1m2','m3z2m2','m3z3m2','m2z3p5m2','m3z3p5m2']
#'m2z3p5m2','m2z3p5m2\_hCBM\_TDUx7']
#'m3z2m2','m3z2m2, 0.75x $^{96}$Mo(n,$\gamma$)','m3z2m2, kad1','m3z2m2-hCBM','m3z2m2-hCBM, 0.75x $^{96}$Mo(n,$\gamma$)','RI18']
#'m3z3m2-hCBM-rotmix.st','m3z3m2-hCBM-rotmix.stx2','m3z3m2-hCBM-rotmix.stx1.5','m3z2m2, 0.67x $^{96}$Mo(n,$\gamma$)','m3z3m2-hCBM, 0.67x $^{96}$Mo(n,$\gamma$)','m3z3m2-hCBM-12','m3z2m2, 0.8x $^{96}$Mo(n,$\gamma$) kad1, 1.15x $^{92}$Mo(n,$\gamma$)']
#'m3z3m2-hCBM','m3z3m2-hCBM this work', 'm3z3m2-hCBM this work lostd','m3z3m2-hCBM this work stdhi','m3z3m2 longland histd','m3z3m2 longland stdhi'] ## legend labels
'master','modular2']

sparcity_sindex=2000   ## sparcity to adopt reading s-process index data
sparcity_isoratio=2000 ## sparcity to adopt reading isotopic ratio data
markers=['tab:red','tab:blue','tab:green','tab:gray','tab:red','tab:blue','tab:green','tab:gray','tab:blue','tab:cyan','tab:purple','xkcd:magenta','tab:red','tab:green','tab:gray','tab:brown','tab:green','tab:blue','tab:red','k', 'tab:gray','tab:cyan','tab:olive','tab:brown', 'y', 'tab:pink', 'tab:orange']## markers to use while plotting s-process indices
symbols=['s','o','h','^','*','D','<','p','>','d','p','s','o','h','^','<','D','^','s','o'] ## linestyles to use while plotting s-process indices
line=['solid','solid','solid','solid','dashed','dashed','dashed','dashed','solid','solid','solid','solid','solid','solid','dashed','dotted']

## IF ISO_RATIO IS TRUE THIS SETS THE ISOTOPIC RATIOS YOU WANT TO PLOT
a = [''.join( isotopes_x[0].split('-')[::-1] ),''.join( isotopes_x[1].split('-')[::-1] )]
b = '/'.join(a)
c = str('d()')
iso_ratio_1 = c[:2] + b + c[2:]
d = [''.join( isotopes_y[0].split('-')[::-1] ),''.join( isotopes_y[1].split('-')[::-1] )]
e = '/'.join(d)
iso_ratio_2 = c[:2] + e + c[2:]

## load data
x,y,xerr,yerr,iso_ratio_111,iso_ratio_222,labels=psg.presolargrains(isotopes_x,isotopes_y,reference,clas_4) 
    
##############################
df1=pd.read_csv('isotopi_m2p0z1m2_000_20180116_152314.txt', sep='\s+',index_col='Isotope',usecols=[0,3,5,6,7,8,9,10,11,12,13],engine='python').T
df2=pd.read_csv('isotopi_m2p0z2m2_000_20180116_152257.txt', sep='\s+',index_col='Isotope',usecols=[0,3,5,6,7,8,9,10,11,12,13],engine='python').T
df3=pd.read_csv('isotopi_m3p0z1m2_000_20180116_152226.txt', sep='\s+',index_col='Isotope',usecols=[0,3,5,6,7,8,9,10,11,12,13,14,15],engine='python').T
df4=pd.read_csv('isotopi_m3p0z2m2_000_20180116_152153.txt', sep='\s+',index_col='Isotope',usecols=[0,3,5,6,7,8,9,10,11,12,13,14,15,16,17,18],engine='python').T

df6=pd.read_csv('isotopi_m2p0zsun_000_20190703_144412.txt', sep='\s+',index_col='Isotope',usecols=[0,3,5,6,7,8,9,10,11,12,13],engine='python').T
df7=pd.read_csv('isotopi_m2p0zsun_010_20190703_144412.txt', sep='\s+',index_col='Isotope',usecols=[0,3,5,6,7,8,9,10,11,12],engine='python').T
df8=pd.read_csv('isotopi_m2p0zsun_030_20190703_144412.txt', sep='\s+',index_col='Isotope',usecols=[0,3,5,6,7,8,9,10,11,12,13],engine='python').T
df9=pd.read_csv('isotopi_m2p0zsun_060_20190703_144412.txt', sep='\s+',index_col='Isotope',usecols=[0,3,5,6,7,8,9,10,11,12],engine='python').T

df5=pd.read_csv('isotopi_m2p0zsun_000_20180720_16736.txt', sep='\s+',index_col='Isotope',usecols=[0,3],engine='python').T

inimo92=df5['Mo92'].values
inimo94=df5['Mo94'].values
inimo95=df5['Mo95'].values
inimo96=df5['Mo96'].values
inimo97=df5['Mo97'].values
inimo98=df5['Mo98'].values
inimo100=df5['Mo00'].values

list1=[df1,df3,df4,df6,df7,df8,df9]
#metallicity=['FRUITY m2z1m2','FRUITY m2z2m2','FRUITY m3z1m2','FRUITY m3z2m2']
metallicity=['FRUITY m2z1m2','FRUITY m3z1m2','FRUITY m3z2m2','FRUITY m2zsun 000','FRUITY m2zsun 010','FRUITY m2zsun 030','FRUITY m2zsun 060']

mo92=[]
for i in list1:
    df_plot = i['Mo92'].values
    mo92.append(df_plot)
mo94=[]
for i in list1:
    df_plot = i['Mo94'].values
    mo94.append(df_plot)
mo95=[]
for i in list1:
    df_plot = i['Mo95'].values
    mo95.append(df_plot)
mo96=[]
for i in list1:
    df_plot = i['Mo96'].values
    mo96.append(df_plot)
mo97=[]
for i in list1:
    df_plot = i['Mo97'].values
    mo97.append(df_plot)
mo98=[]
for i in list1:
    df_plot = i['Mo98'].values
    mo98.append(df_plot)
mo100=[]
for i in list1:
    df_plot = i['Mo00'].values
    mo100.append(df_plot)
    
ratio2=(((np.array(mo92)/np.array(mo96))/(inimo92/inimo96))-1)*1000
ratio4=(((np.array(mo94)/np.array(mo96))/(inimo94/inimo96))-1)*1000
ratio5=(((np.array(mo95)/np.array(mo96))/(inimo95/inimo96))-1)*1000
ratio7=(((np.array(mo97)/np.array(mo96))/(inimo97/inimo96))-1)*1000
ratio8=(((np.array(mo98)/np.array(mo96))/(inimo98/inimo96))-1)*1000
ratio0=(((np.array(mo100)/np.array(mo96))/(inimo100/inimo96))-1)*1000

df10 = df1[(df1['C12']/df1['O16']>=1.)]
df20 = df2[(df2['C12']/df2['O16']>=1.)]
df30 = df3[(df3['C12']/df3['O16']>=1.)]
df40 = df4[(df4['C12']/df4['O16']>=1.)]
df60 = df6[(df6['C12']/df6['O16']>=1.)]
df70 = df7[(df7['C12']/df7['O16']>=1.)]
df80 = df8[(df8['C12']/df8['O16']>=1.)]
df90 = df9[(df9['C12']/df9['O16']>=1.)]

#list2=[df10,df20,df30,df40]
list2=[df10,df30,df40,df60,df70,df80,df90]

mo92_co=[]
for i in list2:
    df_plot = i['Mo92'].values
    mo92_co.append(df_plot)
mo94_co=[]
for i in list2:
    df_plot = i['Mo94'].values
    mo94_co.append(df_plot)
mo95_co=[]
for i in list2:
    df_plot = i['Mo95'].values
    mo95_co.append(df_plot)
mo96_co=[]
for i in list2:
    df_plot = i['Mo96'].values
    mo96_co.append(df_plot)
mo97_co=[]
for i in list2:
    df_plot = i['Mo97'].values
    mo97_co.append(df_plot)
mo98_co=[]
for i in list2:
    df_plot = i['Mo98'].values
    mo98_co.append(df_plot)
mo100_co=[]
for i in list2:
    df_plot = i['Mo00'].values
    mo100_co.append(df_plot)
    
ratio2_co=(((np.array(mo92_co)/np.array(mo96_co))/(inimo92/inimo96))-1)*1000
ratio4_co=(((np.array(mo94_co)/np.array(mo96_co))/(inimo94/inimo96))-1)*1000
ratio5_co=(((np.array(mo95_co)/np.array(mo96_co))/(inimo95/inimo96))-1)*1000
ratio7_co=(((np.array(mo97_co)/np.array(mo96_co))/(inimo97/inimo96))-1)*1000
ratio8_co=(((np.array(mo98_co)/np.array(mo96_co))/(inimo98/inimo96))-1)*1000
ratio0_co=(((np.array(mo100_co)/np.array(mo96_co))/(inimo100/inimo96))-1)*1000
###############################
        
# get initial value
initial_isotopic_ratio_x = []
initial_isotopic_ratio_y = []
# sparcity for cycles I am looking at.
sparsity = sparcity_isoratio
isotopic_ratio_x  = []
isotopic_ratio_y = []
isotopic_ratio_x_tps = []
isotopic_ratio_y_tps = []
isotopic_ratio_x_tps_co = []
isotopic_ratio_y_tps_co = [] 
k = 0
for i in run_mass: 
    initial_isotopic_ratio_x.append(float(i.se.get(min(i.se.cycles),'iso_massf',isotopes_x[0]))/float(i.se.get(min(i.se.cycles),'iso_massf',isotopes_x[1])))
    initial_isotopic_ratio_y.append(float(i.se.get(min(i.se.cycles),'iso_massf',isotopes_y[0]))/float(i.se.get(min(i.se.cycles),'iso_massf',isotopes_y[1])))  
    dum_isotopic_ratio_x = []
    dum_isotopic_ratio_y = []
    dumdum_isotopic_ratio_x = []
    dumdum_isotopic_ratio_y = []
    dum_isotopic_ratio_x_tps = []
    dum_isotopic_ratio_y_tps = []
    dum_isotopic_ratio_x_tps_co = []
    dum_isotopic_ratio_y_tps_co = []
    co_ratio=[]
    for j in i.se.cycles[0::sparsity]:
        print(j)
        dum_isotopic_ratio_x.append(float(i.se.get(j,'iso_massf',isotopes_x[0]))/float(i.se.get(j,'iso_massf',isotopes_x[1])))
        dum_isotopic_ratio_y.append(float(i.se.get(j,'iso_massf',isotopes_y[0]))/float(i.se.get(j,'iso_massf',isotopes_y[1])))
        co_ratio.append((float(i.se.get(j,'iso_massf','C-12')*4.))/(float(i.se.get(j,'iso_massf','O-16'))*3))
        if (len(co_ratio)>1):
            if (co_ratio[len(co_ratio)-1]>(co_ratio[len(co_ratio)-2]+0.02)):
                dum_isotopic_ratio_x_tps.append(float(i.se.get(j,'iso_massf',isotopes_x[0]))/float(i.se.get(j,'iso_massf',isotopes_x[1])))
                dum_isotopic_ratio_y_tps.append(float(i.se.get(j,'iso_massf',isotopes_y[0]))/float(i.se.get(j,'iso_massf',isotopes_y[1])))
                if (co_ratio[len(co_ratio)-1]>1.):
                    dum_isotopic_ratio_x_tps_co.append(float(i.se.get(j,'iso_massf',isotopes_x[0]))/float(i.se.get(j,'iso_massf',isotopes_x[1])))
                    dum_isotopic_ratio_y_tps_co.append(float(i.se.get(j,'iso_massf',isotopes_y[0]))/float(i.se.get(j,'iso_massf',isotopes_y[1])))

        dumdum_isotopic_ratio_x=(np.array(dum_isotopic_ratio_x)/initial_isotopic_ratio_x[-1:]-1.)*1000.
        dumdum_isotopic_ratio_y=(np.array(dum_isotopic_ratio_y)/initial_isotopic_ratio_y[-1:]-1.)*1000.
        dumdum_isotopic_ratio_x_tps=(np.array(dum_isotopic_ratio_x_tps)/initial_isotopic_ratio_x[-1:]-1.)*1000.
        dumdum_isotopic_ratio_y_tps=(np.array(dum_isotopic_ratio_y_tps)/initial_isotopic_ratio_y[-1:]-1.)*1000.
        dumdum_isotopic_ratio_x_tps_co=(np.array(dum_isotopic_ratio_x_tps_co)/initial_isotopic_ratio_x[-1:]-1.)*1000.
        dumdum_isotopic_ratio_y_tps_co=(np.array(dum_isotopic_ratio_y_tps_co)/initial_isotopic_ratio_y[-1:]-1.)*1000.
   
    isotopic_ratio_x.append(dumdum_isotopic_ratio_x)
    isotopic_ratio_y.append(dumdum_isotopic_ratio_y)
    isotopic_ratio_x_tps.append(dumdum_isotopic_ratio_x_tps)
    isotopic_ratio_y_tps.append(dumdum_isotopic_ratio_y_tps)
    isotopic_ratio_x_tps_co.append(dumdum_isotopic_ratio_x_tps_co)
    isotopic_ratio_y_tps_co.append(dumdum_isotopic_ratio_y_tps_co)
    k = k+1

mass_label =[]
for i in run_mass:
    mass_label.append(float(i.se.get('mini')))
array_to_plot_x_tps = isotopic_ratio_x_tps
array_to_plot_y_tps = isotopic_ratio_y_tps
array_to_plot_x_tps_co = isotopic_ratio_x_tps_co
array_to_plot_y_tps_co = isotopic_ratio_y_tps_co

plt.rcParams['text.usetex'] = True
plt.rcParams['xtick.major.size'] = 20
plt.rcParams['xtick.major.width'] = 2
plt.rcParams['ytick.major.size'] = 20
plt.rcParams['ytick.major.width'] = 2
plt.rcParams['xtick.labelsize'] = 30
plt.rcParams['ytick.labelsize'] = 30
plt.tick_params(axis='both',pad=5,direction='in')
plt.axhline(y=0.,linewidth=2, color='k')
plt.axvline(x=0.,linewidth=2, color='k')

######################################
#if isotopes_y == ['Mo-94','Mo-96']:
#    for i in range(len(mo94)):
#        plt.plot(ratio2[i],ratio4[i],c=markers[i],ls=line[i],linewidth=2.)
#        plt.plot(ratio2_co[i],ratio4_co[i],marker=symbols[i],ls=line[i],c=markers[i],markersize=20.,label=str(metallicity[i]))
#if isotopes_y == ['Mo-95','Mo-96']:   
#    for i in range(len(mo95)):
#        plt.plot(ratio2[i],ratio5[i],c=markers[i],ls=line[i],linewidth=2.)
#        plt.plot(ratio2_co[i],ratio5_co[i],marker=symbols[i],ls=line[i],c=markers[i],markersize=20.,label=str(metallicity[i]))
#if isotopes_y == ['Mo-97','Mo-96']:
#    for i in range(len(mo97)):
#        plt.plot(ratio2[i],ratio7[i],c=markers[i],ls=line[i],linewidth=2.)
#        plt.plot(ratio2_co[i],ratio7_co[i],marker=symbols[i],ls=line[i],c=markers[i],markersize=20.,label=str(metallicity[i]))
#if isotopes_y == ['Mo-98','Mo-96']:
#    for i in range(len(mo98)):
#        plt.plot(ratio2[i],ratio8[i],c=markers[i],ls=line[i],linewidth=2.)
#        plt.plot(ratio2_co[i],ratio8_co[i],marker=symbols[i],ls=line[i],c=markers[i],markersize=20.,label=str(metallicity[i]))
#if isotopes_y == ['Mo-100','Mo-96']:
#    for i in range(len(mo100)):
#        plt.plot(ratio2[i],ratio0[i],c=markers[i],ls=line[i],linewidth=2.)
#        plt.plot(ratio2_co[i],ratio0_co[i],marker=symbols[i],ls=line[i],c=markers[i],markersize=20.,label=str(metallicity[i]))
#######################################

plt.errorbar(x.astype(np.float64),y.astype(np.float64),xerr=2*xerr.astype(np.float64),yerr=2*yerr.astype(np.float64),marker=symbols[-1],c=markers[-1],ls='None',markersize=8.,label=reference+', '+clas_4)
for k in range(0,len(array_to_plot_x_tps)):
    plt.plot(array_to_plot_x_tps[k],array_to_plot_y_tps[k],c=markers[k],ls=line[k],markersize=10.,linewidth=2.)
    plt.plot(array_to_plot_x_tps_co[k],array_to_plot_y_tps_co[k],marker=symbols[k],c=markers[k],ls=line[k],markersize=20.,linewidth=2.,label=str(metallicity_label[k]))
plt.xlabel(iso_ratio_111, fontsize=40)
plt.ylabel(iso_ratio_222, fontsize=40)

handles, labels = plt.gca().get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys(),prop={'size':25})
plt.tight_layout()
plt.show()                                   

for label, x, y in zip(labels, x, y):
    plt.annotate(
        label,
        xy=(x, y), xytext=(-20, 20),
        textcoords='offset points', #ha='right', va='bottom',
        bbox=dict(boxstyle='round,pad=0.5', fc='white', alpha=0.5),
        arrowprops=dict(arrowstyle = '->', connectionstyle='arc3,rad=0'))
        
#plt.show()