read_models_FeNi.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 = ['Fe-58','Fe-56']
#isotopes_y = ['Fe-54','Fe-56']
#isotopes_y = ['Fe-57','Fe-56']
#isotopes_y = ['Ni-64','Ni-58']

isotopes_x = ['Ni-64','Ni-58']
isotopes_y = ['Ni-60','Ni-58']
#isotopes_y = ['Ni-61','Ni-58']
#isotopes_y = ['Ni-62','Ni-58']

## 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', 'Liu ApJ 2015',
# '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 = 'Trappitsch GCA 2018' # insert here for dataset
reference1 = 'Marhas ApJ 2008'

## 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_6 = 'MS'
clas_7 = 'M'

clas_1 = 'Ni58'
clas_2 = 'Ni60'
clas_3 = 'Ni61'
clas_4 = 'Ni62'
clas_5 = 'Ni64'

run1=mp.se('/run/media/s1785107/Seagate Backup Plus Drive/eddie3/NuclearPhysics/NuPPN/frames/mppnp/RUN_set1upd_m3z3m2/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')

#run_mass=[run1,run2,run3,run4,run5,run6,run7,run8,run9,run10,run11,run12,run13]
run_mass=[run1,run2,run3,run4,run5,run6]
metallicity_label=[#'$M1p65.z1m2$','$M1p65.z2m2$']
#'$M2.z2m2$','$M2.z3m2$','$M3.z2m2$','$M3.z3m2$',
'$M3.z2m2$ $diffusion/30$','$M3.z2m2$ $diffusion/100$','$M3.z2m2$ $diffusion/250$','$M3.z2m2$ $diffusion/500$','$M3.z2m2$ $diffusion/1000$','$M3.z2m2$ $diffusion/2000$']
#'$M3.z1m2$']#,'$M3.z1m2$ $KADoNiS$ $1.0$','$M3.z1m2$ $KADoNiS$ $1.0,$ $Ni$','$M3.z1m2$ $KADoNiS$ $1.0,$ $FeNi$']
#'$M3.z2m2$ $KADoNiS$ $0.3$ $Asplund$ $2009$','$M3.z2m2$ $KADoNiS$ $0.3$ $2xNi58$ $iniab$','$M3.z2m2$ $KADoNiS$ $0.3$ $0.5xNi58$ $iniab$','$M3.z2m2$ $KADoNiS$ $0.3$ $Asplund$ $2005$']
#'$M2.z1m2$','$M2.z2m2$','$M3.z1m2$','$M3.z2m2$','$M2.z3m2$','$M3.z3m2$']
#'$M2.z6m3\_hCBM$','$M3.z6m3\_hCBM$','$M2.z1m2\_hCBM$','$M2.z2m2\_hCBM$','$M3.z2m2\_hCBM$']
#'$M1p86.z2m2\_hCBM$','$M1p86.z2m2\_hCBM,$ $Ne22$','$M1p86.z2m2,$ $Ne22,$ $3x$ $Fe56(n,\gamma)$','$M2.z2m2$','$M2.z2m2\_hCBM$']
#'$M2.z1m2$','$M2.z2m2,$ $Ne22$ $\uparrow(\\alpha,\gamma)\downarrow(\\alpha,n)$','$M2.z1m2,$ $Ne22,$ $3x$ $Fe56(n,\gamma)$','$M2.z1m2,$ $Ne22,$ $2x$ $Fe56(n,\gamma)$','$M2.z1m2,$ $Ne22,$ $1.5x$ $Fe56(n,\gamma)$',\
#'$M2.z1m2\_hCBM$','$M2.z1m2\_hCBM,$ $Ne22$ $\uparrow(\\alpha,\gamma)\downarrow(\\alpha,n)$','$M2.z1m2\_hCBM,$ $Ne22,$ $3x$ $Fe56(n,\gamma)$','$M2.z1m2\_hCBM,$ $Ne22,$ $2x$ $Fe56(n,\gamma)$','$M2.z1m2\_hCBM,$ $Ne22,$ $`1.5x$ $Fe56(n,\gamma)$']
#'$M2.z1m2,$ $Ne22,$ $1.5x$ $Fe56(n,\gamma)$','$M2.z1m2,$ $Ne22,$ $2x$ $Fe56(n,\gamma)$','$M2.z1m2,$ $1.15x$ $Fe56(n,\gamma)$','$M2.z1m2,$ $Ne22,$ $3x$ $Fe56(n,\gamma)$']
#'$M3.z2m2\_he07$','$M3.z2m2\_he07,$ $Longland$ $13$','$M3.z2m2\_he07,$ $Longland$ $22$','$M3.z2m2\_he07,$ $Longland$ $31$']
#'$M2.z1m2\_he07$','$M2.z1m2\_he07,$ $Longland$ $13$','$M2.z1m2\_he07,$ $Longland$ $22$','$M2.z1m2\_he07,$ $Longland$ $31$','$M2.z1m2\_he07,$ $new$ $Ne22 + \\alpha$ $13$','$M2.z1m2\_he07,$ $new$ $Ne22 + \\alpha$ $22$','$M2.z1m2\_he07,$ $new$ $Ne22 + \\alpha$ $31$']
#'$M2.z1m2$','$M2.z1m2,$ $Longland$ $22$','$M2.z1m2\_he07,$ $new$ $Ne22 + \\alpha$ $22$','$M2.z1m2,$ $Longland$ $31$']
#'$M1p65.z2m2$','$M2.z2m2$','$M3.z2m2$']
#'$M2.z2m2$','$M2.z2m2,$ $Ritter$ $et$ $al.$ $2017$']
#'$M1p65.z2m2$','$M1p65.z2m2,$ $Ritter$ $et$ $al.$ $2017$','$M3.z1m2\_lowres$','$M3.z1m2$','$M3.z1m2,$ $Ritter$ $et$ $al.$ $2017$']
#'$M2.z1m2\_he07$','$M2.z1m2\_R17$','$M2.z2m2\_he07$','$M2.z2m2\_R17$','$M3.z1m2\_he07$','$M3.z1m2\_R17$','$M3.z2m2\_he07$','$M3.z2m2\_R17$'] ## legend labels
sparcity_sindex=2000   ## sparcity to adopt reading s-process index data
sparcity_isoratio=2000 ## sparcity to adopt reading isotopic ratio data
markers=['tab:green','tab:blue','k', 'tab:cyan', 'tab:olive', 'tab:gray', 'tab:brown', 'y', 'tab:purple', 'tab:pink', 'tab:orange','tab:red']## markers to use while plotting s-process indices
symbols=['o','s','h','^','>','D','<','p','d','^','s','o'] ## linestyles to use while plotting s-process indices
line=['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_6)  

df = pd.read_csv('isotopi_m3p0zsun_000_20180116_151709.txt', sep='\s+')
ini = df['INI'].values
ni58_ini=ini[88]
ni60_ini=ini[90]
ni61_ini=ini[91]
ni62_ini=ini[92]
ni64_ini=ini[94]

df1 = pd.read_csv('isotopi_m2p0z1m2_000_20180116_152314.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12,13])
df2 = pd.read_csv('isotopi_m2p0z1m2_000_20180116_152314.txt', sep='\s+')
df_plot1 = df1[(df2['Isotope']==clas_1)]
df_plot2 = df1[(df2['Isotope']==clas_2)]
df_plot3 = df1[(df2['Isotope']==clas_3)]
df_plot4 = df1[(df2['Isotope']==clas_4)]
df_plot5 = df1[(df2['Isotope']==clas_5)]
ni58 = df_plot1.values
ni60 = df_plot2.values
ni61 = df_plot3.values
ni62 = df_plot4.values
ni64 = df_plot5.values
d_ni6458 = (((ni64/ni58)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058 = (((ni60/ni58)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158 = (((ni61/ni58)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258 = (((ni62/ni58)/(ni62_ini/ni58_ini))-1)*1000

df3 = pd.read_csv('isotopi_m2p0z2m2_000_20180116_152257.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12,13])
df_plot6 = df3[(df2['Isotope']==clas_1)]
df_plot7 = df3[(df2['Isotope']==clas_2)]
df_plot8 = df3[(df2['Isotope']==clas_3)]
df_plot9 = df3[(df2['Isotope']==clas_4)]
df_plot11 = df3[(df2['Isotope']==clas_5)]
ni58_1 = df_plot6.values
ni60_1 = df_plot7.values
ni61_1 = df_plot8.values
ni62_1 = df_plot9.values
ni64_1 = df_plot11.values
d_ni6458_1 = (((ni64_1/ni58_1)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_1 = (((ni60_1/ni58_1)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_1 = (((ni61_1/ni58_1)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_1 = (((ni62_1/ni58_1)/(ni62_ini/ni58_ini))-1)*1000

df4 = pd.read_csv('isotopi_m3p0z1m2_000_20180116_152226.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12,13,14,15])
df5 = pd.read_csv('isotopi_m3p0z1m2_000_20180116_152226.txt', sep='\s+')
df_plot12 = df4[(df5['Isotope']==clas_1)]
df_plot13 = df4[(df5['Isotope']==clas_2)]
df_plot14 = df4[(df5['Isotope']==clas_3)]
df_plot15 = df4[(df5['Isotope']==clas_4)]
df_plot16 = df4[(df5['Isotope']==clas_5)]
ni58_2 = df_plot12.values
ni60_2 = df_plot13.values
ni61_2 = df_plot14.values
ni62_2 = df_plot15.values
ni64_2 = df_plot16.values
d_ni6458_2 = (((ni64_2/ni58_2)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_2 = (((ni60_2/ni58_2)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_2 = (((ni61_2/ni58_2)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_2 = (((ni62_2/ni58_2)/(ni62_ini/ni58_ini))-1)*1000

df6 = pd.read_csv('isotopi_m3p0z2m2_000_20180116_152153.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12,13,14,15])
df_plot17 = df6[(df5['Isotope']==clas_1)]
df_plot18 = df6[(df5['Isotope']==clas_2)]
df_plot19 = df6[(df5['Isotope']==clas_3)]
df_plot20 = df6[(df5['Isotope']==clas_4)]
df_plot21 = df6[(df5['Isotope']==clas_5)]
ni58_3 = df_plot17.values
ni60_3 = df_plot18.values
ni61_3 = df_plot19.values
ni62_3 = df_plot20.values
ni64_3 = df_plot21.values
d_ni6458_3 = (((ni64_3/ni58_3)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_3 = (((ni60_3/ni58_3)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_3 = (((ni61_3/ni58_3)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_3 = (((ni62_3/ni58_3)/(ni62_ini/ni58_ini))-1)*1000

df7 = pd.read_csv('isotopi_m2p0zsun_010_20180720_153729.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12])
df_plot22 = df7[(df2['Isotope']==clas_1)]
df_plot23 = df7[(df2['Isotope']==clas_2)]
df_plot24 = df7[(df2['Isotope']==clas_3)]
df_plot25 = df7[(df2['Isotope']==clas_4)]
df_plot26 = df7[(df2['Isotope']==clas_5)]
ni58_4 = df_plot22.values
ni60_4 = df_plot23.values
ni61_4 = df_plot24.values
ni62_4 = df_plot25.values
ni64_4 = df_plot26.values
d_ni6458_4 = (((ni64_4/ni58_4)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_4 = (((ni60_4/ni58_4)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_4 = (((ni61_4/ni58_4)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_4 = (((ni62_4/ni58_4)/(ni62_ini/ni58_ini))-1)*1000

df8 = pd.read_csv('isotopi_m2p0zsun_000_20180720_16736.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12,13])
df_plot27 = df8[(df2['Isotope']==clas_1)]
df_plot28 = df8[(df2['Isotope']==clas_2)]
df_plot29 = df8[(df2['Isotope']==clas_3)]
df_plot30 = df8[(df2['Isotope']==clas_4)]
df_plot31 = df8[(df2['Isotope']==clas_5)]
ni58_5 = df_plot27.values
ni60_5 = df_plot28.values
ni61_5 = df_plot29.values
ni62_5 = df_plot30.values
ni64_5 = df_plot31.values
d_ni6458_5 = (((ni64_5/ni58_5)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_5 = (((ni60_5/ni58_5)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_5 = (((ni61_5/ni58_5)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_5 = (((ni62_5/ni58_5)/(ni62_ini/ni58_ini))-1)*1000

df9 = pd.read_csv('isotopi_m2p0zsun_T00_20180830_18519.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12,13])
df_plot32 = df9[(df2['Isotope']==clas_1)]
df_plot33 = df9[(df2['Isotope']==clas_2)]
df_plot34 = df9[(df2['Isotope']==clas_3)]
df_plot35 = df9[(df2['Isotope']==clas_4)]
df_plot36 = df9[(df2['Isotope']==clas_5)]
ni58_6 = df_plot32.values
ni60_6 = df_plot33.values
ni61_6 = df_plot34.values
ni62_6 = df_plot35.values
ni64_6 = df_plot36.values
d_ni6458_6 = (((ni64_6/ni58_6)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_6 = (((ni60_6/ni58_6)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_6 = (((ni61_6/ni58_6)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_6 = (((ni62_6/ni58_6)/(ni62_ini/ni58_ini))-1)*1000

df10 = pd.read_csv('isotopi_m2p0z1m2_T00_20180830_18255.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12,13])
df_plot37 = df10[(df2['Isotope']==clas_1)]
df_plot38 = df10[(df2['Isotope']==clas_2)]
df_plot39 = df10[(df2['Isotope']==clas_3)]
df_plot40 = df10[(df2['Isotope']==clas_4)]
df_plot41 = df10[(df2['Isotope']==clas_5)]
ni58_7 = df_plot37.values
ni60_7 = df_plot38.values
ni61_7 = df_plot39.values
ni62_7 = df_plot40.values
ni64_7 = df_plot41.values
d_ni6458_7 = (((ni64_7/ni58_7)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_7 = (((ni60_7/ni58_7)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_7 = (((ni61_7/ni58_7)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_7 = (((ni62_7/ni58_7)/(ni62_ini/ni58_ini))-1)*1000

df11 = pd.read_csv('isotopi_m2p0zsun_030_20180831_123823.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12,13])
df_plot42 = df11[(df2['Isotope']==clas_1)]
df_plot43 = df11[(df2['Isotope']==clas_2)]
df_plot44 = df11[(df2['Isotope']==clas_3)]
df_plot45 = df11[(df2['Isotope']==clas_4)]
df_plot46 = df11[(df2['Isotope']==clas_5)]
ni58_8 = df_plot42.values
ni60_8 = df_plot43.values
ni61_8 = df_plot44.values
ni62_8 = df_plot45.values
ni64_8 = df_plot46.values
d_ni6458_8 = (((ni64_8/ni58_8)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_8 = (((ni60_8/ni58_8)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_8 = (((ni61_8/ni58_8)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_8 = (((ni62_8/ni58_8)/(ni62_ini/ni58_ini))-1)*1000

df12 = pd.read_csv('isotopi_m2p0zsun_060_20180831_123823.txt', sep='\s+',usecols=[3,4,5,6,7,8,9,10,11,12])
df_plot47 = df12[(df2['Isotope']==clas_1)]
df_plot48 = df12[(df2['Isotope']==clas_2)]
df_plot49 = df12[(df2['Isotope']==clas_3)]
df_plot50 = df12[(df2['Isotope']==clas_4)]
df_plot51 = df12[(df2['Isotope']==clas_5)]
ni58_9 = df_plot47.values
ni60_9 = df_plot48.values
ni61_9 = df_plot49.values
ni62_9 = df_plot50.values
ni64_9 = df_plot51.values
d_ni6458_9 = (((ni64_9/ni58_9)/(ni64_ini/ni58_ini))-1)*1000
d_ni6058_9 = (((ni60_9/ni58_9)/(ni60_ini/ni58_ini))-1)*1000
d_ni6158_9 = (((ni61_9/ni58_9)/(ni61_ini/ni58_ini))-1)*1000
d_ni6258_9 = (((ni62_9/ni58_9)/(ni62_ini/ni58_ini))-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['axes.linewidth'] = 5
plt.rcParams['xtick.major.size'] = 20
plt.rcParams['xtick.major.width'] = 2#3
plt.rcParams['ytick.major.size'] = 20
plt.rcParams['ytick.major.width'] = 2#3
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')

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_7)
#plt.errorbar(x1,y1,xerr=2*xerr1,yerr=2*yerr1,marker=symbols[-2],c=markers[-2],ls='None',markersize=8.,label=reference2+', '+clas_7)
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)

#for i in range(len(d_ni6058_2)):
#    if isotopes_y == ['Ni-60','Ni-58']:
#        #plt.plot(d_ni6458[i],d_ni6058[i],'k*-',markersize=10.,label='M2.z1m2 FRUITY')
#        #plt.plot(d_ni6458_1[i],d_ni6058_1[i],'bs-',markersize=10.,label='M2.z2m2 FRUITY')
#        plt.plot(d_ni6458_5[i],d_ni6058_5[i],'gd-',markersize=10.,label='M2.zsun FRUITY')
#        plt.plot(d_ni6458_4[i],d_ni6058_4[i],'r^-',markersize=10.,label='M2.zsun rotation10 FRUITY')
#        #plt.plot(d_ni6458_2[i],d_ni6058_2[i],'ch-',markersize=10.,label='M3.z1m2 FRUITY')
#        #plt.plot(d_ni6458_3[i],d_ni6058_3[i],'mo-',markersize=10.,label='M3.z2m2 FRUITY')
#        #plt.plot(d_ni6458_6[i],d_ni6058_6[i],'yp-',markersize=10.,label='M2.zsun FRUITY extended pocket')
#        #plt.plot(d_ni6458_7[i],d_ni6058_7[i],'c8-',markersize=10.,label='M2.z1m2 FRUITY extended pocket')
#        plt.plot(d_ni6458_8[i],d_ni6058_8[i],'y>-',markersize=10.,label='M2.zsun rotation30 FRUITY')
#        plt.plot(d_ni6458_9[i],d_ni6058_9[i],'bo-',markersize=10.,label='M2.zsun rotation60 FRUITY')
#    elif isotopes_y == ['Ni-61','Ni-58']:
#        #plt.plot(d_ni6458[i],d_ni6158[i],'k*-',markersize=10.,label='M2.z1m2 FRUITY')
#        #plt.plot(d_ni6458_1[i],d_ni6158_1[i],'bs-',markersize=10.,label='M2.z2m2 FRUITY')
#        plt.plot(d_ni6458_5[i],d_ni6158_5[i],'gd-',markersize=10.,label='M2.zsun FRUITY')
#        plt.plot(d_ni6458_4[i],d_ni6158_4[i],'r^-',markersize=10.,label='M2.zsun rotation10 FRUITY')
#        #plt.plot(d_ni6458_2[i],d_ni6158_2[i],'ch-',markersize=10.,label='M3.z1m2 FRUITY')
#        #plt.plot(d_ni6458_3[i],d_ni6158_3[i],'mo-',markersize=10.,label='M3.z2m2 FRUITY')
#        #plt.plot(d_ni6458_6[i],d_ni6158_6[i],'yp-',markersize=10.,label='M2.zsun FRUITY extended pocket')
#        #plt.plot(d_ni6458_7[i],d_ni6158_7[i],'c8-',markersize=10.,label='M2.z1m2 FRUITY extended pocket')
#        plt.plot(d_ni6458_8[i],d_ni6158_8[i],'y>-',markersize=10.,label='M2.zsun rotation30 FRUITY')
#        plt.plot(d_ni6458_9[i],d_ni6158_9[i],'bo-',markersize=10.,label='M2.zsun rotation60 FRUITY')
#    elif isotopes_y == ['Ni-62','Ni-58']:
#        #plt.plot(d_ni6458[i],d_ni6258[i],'k*-',markersize=10.,label='M2.z1m2 FRUITY')
#        #plt.plot(d_ni6458_1[i],d_ni6258_1[i],'bs-',markersize=10.,label='M2.z2m2 FRUITY')
#        plt.plot(d_ni6458_5[i],d_ni6258_5[i],'gd-',markersize=10.,label='M2.zsun FRUITY')
#        plt.plot(d_ni6458_4[i],d_ni6258_4[i],'r^-',markersize=10.,label='M2.zsun rotation10 FRUITY')
#        #plt.plot(d_ni6458_2[i],d_ni6258_2[i],'ch-',markersize=10.,label='M3.z1m2 FRUITY')
#        #plt.plot(d_ni6458_3[i],d_ni6258_3[i],'mo-',markersize=10.,label='M3.z2m2 FRUITY')
#        #plt.plot(d_ni6458_6[i],d_ni6258_6[i],'yp-',markersize=10.,label='M2.zsun FRUITY extended pocket')
#        #plt.plot(d_ni6458_7[i],d_ni6258_7[i],'c8-',markersize=10.,label='M2.z1m2 FRUITY extended pocket')
#        plt.plot(d_ni6458_8[i],d_ni6258_8[i],'y>-',markersize=10.,label='M2.zsun rotation30 FRUITY')
#        plt.plot(d_ni6458_9[i],d_ni6258_9[i],'bo-',markersize=10.,label='M2.zsun rotation60 FRUITY')

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.show()                                   

labels = df_plot10['Grain Label'].values
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()