#Abundances of elements and comparison with abundances from Asplund 1997        

from numpy import *
import matplotlib
from matplotlib.pylab import *
import asciidata
import ppn

element=['H','Li','C','N','O','F','Na','Mg','Al','Si','S','Ca','Ti','Fe','Ni','Zn','Y','Zr','Ba','La']

star=['R CrB','RY Sgr','XX Cam','SU Tau','UX Amt','UV Cas','UW Cen','V482 Cyg','Y Mus','RT Nor','RZ Nor','FH Sct','GU Sgr','RS Tel']

z=[1,3,6,7,8,9,11,12,13,14,16,20,22,26,28,30,39,40,56,57]

h=[6.9,6.9,4.1,7.4,6.8,6.0,6.5,4.8,6.1,6.5,5.2,5.6,0.0,7.4]

li=[2.8,0.0,0.0,2.6,0.0,0.0,3.5,0.0,0.0,0.0,3.5,0.0,0.0,0.0]

c=[9.2,8.9,9.0,8.8,8.9,9.2,8.6,8.9,8.9,8.9,8.9,8.8,8.8,8.9]

n=[8.4,8.5,8.9,8.5,8.3,8.5,8.3,8.8,8.8,9.1,8.7,8.7,8.7,8.8]

o=[9.0,7.9,8.4,8.4,8.8,7.5,7.7,8.1,7.7,8.4,8.9,7.7,8.2,8.3]

f=[6.9,7.0,5.6,7.0,6.2,6.2,7.1,6.6,0.0,0.0,0.0,7.2,7.2,0.0]

na=[6.1,6.3,6.8,5.6,5.8,6.4,6.0,6.3,6.3,6.3,6.4,6.1,6.0,6.0]

mg=[0.0,0.0,0.0,0.0,0.0,0.0,5.7,0.0,6.3,0.0,0.0,0.0,6.9,0.0]

al=[5.8,6.0,0.0,5.2,0.0,6.0,6.3,6.2,6.1,6.2,6.3,5.9,5.7,5.9]

si=[7.2,7.3,7.1,6.7,6.9,7.4,7.0,7.2,7.3,7.4,7.1,7.1,7.2,7.1]

s=[6.8,7.3,6.8,6.5,6.2,7.0,6.7,6.9,6.9,7.7,6.8,7.0,7.0,6.8]

ca=[5.3,5.3,5.4,5.0,5.5,5.6,5.2,5.4,5.3,5.8,5.4,5.1,5.4,5.3]

ti=[0.0,4.1,4.0,3.7,0.0,4.0,4.1,0.0,4.2,0.0,0.0,0.0,0.0,0.0]

fe=[6.5,6.7,6.8,6.1,6.2,6.9,6.3,6.7,6.5,6.8,6.6,6.3,6.3,6.4]

ni=[5.5,5.9,6.1,5.4,5.8,6.2,5.9,5.8,6.0,6.2,5.9,5.8,5.6,5.7]

zn=[0.0,4.5,0.0,3.6,0.0,4.8,4.3,4.4,4.4,4.7,4.4,4.1,4.4,4.3]

y=[1.5,1.9,2.0,1.3,1.5,2.8,1.5,2.6,2.4,3.1,2.0,2.0,2.4,1.9]

zr=[0.0,1.8,0.0,0.0,0.0,0.0,1.8,2.3,2.6,1.8,0.0,2.3,0.0,0.0]

ba=[1.6,1.4,1.5,0.3,1.0,2.1,1.6,2.6,1.5,2.0,1.5,1.6,1.2,1.5]

la=[0.0,1.0,0.0,0.0,0.0,1.5,1.5,2.2,1.3,1.9,1.1,0.0,0.8,0.0]

E=empty([len(element),len(h)])
X_S=empty([len(element),len(h)])


E[0]=h
E[1]=li
E[2]=c
E[3]=n
E[4]=o
E[5]=f
E[6]=na
E[7]=mg
E[8]=al
E[9]=si
E[10]=s
E[11]=ca
E[12]=ti
E[13]=fe
E[14]=ni
E[15]=zn
E[16]=y
E[17]=zr
E[18]=ba
E[19]=la

E_s=[12.0,1.2,8.6,8.1,8.0,4.6,6.3,7.6,6.5,7.6,7.2,6.4,5.0,7.5,6.3,4.6,2.2,2.6,2.1,1.2]

def x_s(e,i):
    return array(e)-array(s)-E_s[i]+7.2

element=['H','Li','C','N','O','F','Na','Mg','Al','Si','S','Ca','Ti','Fe','Ni','Zn','Y','Zr','Ba','La']

X_S[0]=x_s(h,0)
X_S[1]=x_s(li,1)
X_S[2]=x_s(c,2)
X_S[3]=x_s(n,3)
X_S[4]=x_s(o,4)
X_S[5]=x_s(f,5)
X_S[6]=x_s(na,6)
X_S[7]=x_s(mg,7)
X_S[8]=x_s(al,8)
X_S[9]=x_s(si,9)
X_S[10]=x_s(s,10)
X_S[11]=x_s(ca,11)
X_S[12]=x_s(ti,12)
X_S[13]=x_s(fe,13)
X_S[14]=x_s(ni,14)
X_S[15]=x_s(zn,15)
X_S[16]=x_s(y,16)
X_S[17]=x_s(zr,17)
X_S[18]=x_s(ba,18)
X_S[19]=x_s(la,19)

tmp_E=[]
tmp_z=[]
tmp_x=[]


for a in range(len(h)):
    for b in range(len(element)):
        if E[b][a]!=0.0:
            tmp_x.append(X_S[b][a])
            tmp_z.append(z[b])
    figure(1)
    plot(tmp_z,tmp_x,marker='o',color='m',markersize=5,linestyle='None',mec='m',alpha=0.7,mfc='None')
    tmp_x
    tmp_z
    tmp_x=[]
    tmp_z=[]
    print a

axhline(0,color='black',linestyle='solid')
#title('RCB majority stars and AGB model (12th TP) 49901, 2 Ms')
ylabel('[X/S]')
xlabel('Z')
for i in range(len(z)):
    text(z[i]-0.5,max(X_S[i].tolist())+0.09,element[i])



'''
#Outputting multipage pdf plots

from matplotlib.backends.backend_pdf import PdfPages

pp = PdfPages('abundances1.pdf')

for a in range(len(star)):
    for b in range(len(z)):
        if E[a][b]!=0.0:
            tmp_E.append(E[a][b])
            tmp_z.append(z[b])
    c=a%3.0
    if c==0.0:
        print c
        figure(1)
        scatter(z,E_s,marker='^',color='black')
    plot(tmp_z,tmp_E,marker='o',linestyle='--')
    legend()
    xlim((0,60))
    ylim((0,12.15))
    pp.savefig()
    tmp_E=[]
    tmp_z=[]
    print a

pp.close()

tmp_x=[]
tmp_z=[]
c=0

pp = PdfPages('abundances2.pdf')


tmp_x=[]                                                                           
tmp_z=[]                                                                           
c=0       

for a in range(len(star)):
    for b in range(len(z)):
        if E[a][b]!=0.0:
            tmp_x.append(X_S[a][b])
            tmp_z.append(z[b])
    c=a%3.0
    if c==0.0:
        print c
        figure(1)
    plot(tmp_z,tmp_x,markersize=5,marker='s',linestyle='--',color='c')
    legend()
    #xlim((0,60))
    #ylim((0,12.15))
    axhline(0,color='black',linestyle='solid')
    #pp.savefig()
    tmp_x=[]
    tmp_z=[]
    print a

#pp.close()
'''