In [1]:
%pylab nbagg
from ppmpy import ppm
import os
import matplotlib.patches as patches
from nugridpy.utils import colourblind as cb
import numpy as np
from matplotlib.pyplot import *
import nugridpy.mesa as ms
import nugridpy.astronomy as ast
import nugridpy.utils as u
from matplotlib import gridspec
In [2]:
# network:
ZH = [0,1,1,2,2,3,5,6,6,21,22,26]
AH = [1,1,2,3,4,7,8,12,13,45,46,56]
Z = { 'n' : 7,
'o' : 8,
'f' : 9,
'ne' : 10,
'na' : 11,
'mg' : 12,
'al' : 13,
'si' : 14,
'p' : 15,
's' : 16,
'cl' : 17,
'ar' : 18,
'k' : 19,
'ca' : 20,
'sc' : 21,
'ti' : 22,
'fe' : 26 }
lines=open('./text_data/network.txt','r').readlines()
for line in lines:
zn = line.split()[0]
a1 = int(line.split()[1])
a2 = int(line.split()[2])
z = Z[zn]
for i in range(a1,a2+1):
ZH.append(z)
AH.append(i)
ZH=np.array(ZH)
AH=np.array(AH)
figure(78)
plot(AH-ZH,ZH,'o',color=cb(5))
# stable nuclei
lines = open( './text_data/stable.dat', 'r' ).readlines()[1:]
z = [ int( l.split()[2] ) for l in lines ]
n = [ int( l.split()[3] ) for l in lines ]
ax = gca()
ax.set_aspect('equal')
for i in range(len(z)):
ax.add_patch(
patches.Rectangle(
(n[i] - 0.5, z[i] - 0.5), # (x,y)
1., # width
1., # height
facecolor = 'None',
lw = 0.5,
zorder = 3,
edgecolor = 'k'
)
)
magic = [1,8,16,20,28,40,50,82]
for j in magic:
axvline( j, ls = '-', lw = 0.25, c = cb(7), zorder = 1)
axhline( j, ls = '-', lw = 0.25, c = cb(7), zorder = 1 )
xticks( magic )
yticks( magic )
xlim(-1,30)
ylim(-1,22)
xlabel('N')
ylabel('Z')
Out[2]:
In [3]:
hd = ms.history_data('/data/ppm_rpod2/Stellar_models/O-shell-M25/M25Z0.02/LOGS')
In [4]:
hd.kip_cont(modstart = 16850, modstop = 36680,
xres = 1000, yres = 1000, ylims = [0.,3.], c12_boundary=True,
outlines=False, ixaxis = 'log_time_left', plot_radius=False)
In [5]:
mesa_model = ms.mesa_profile('/data/ppm_rpod2/Stellar_models/O-shell-M25/M25Z0.02/LOGS_N2b', 29350)
D1 = ppm.yprofile('/data/ppm_rpod2/YProfiles/O-shell-M25/D1')
In [6]:
ppm.energy_comparison(D1,mesa_model, xthing = 'm',ifig = 2, silent = True,
range_conv1 = [1.1593, 1.8553],range_conv2 = [2.0824, 2.50],
xlim = [0.5,2.5] , ylim = [8,13], radbase = 4.1297,
dlayerbot = 0.5, totallum = 20.153)
In [7]:
p = ms.mesa_profile('/data/ppm_rpod2/Stellar_models/O-shell-M25/M25Z0.02/LOGS_N2b', 29350)
m = p.get('mass')
mu = p.get('mu')
T = 10.**p.get('logT')
Pppm = np.load( 'text_data/P_ppmsetup.npy' ) * 1.e19 # convert to barye
Mrppm = np.load( 'text_data/Mr_ppmsetup.npy' ) * 5.025e-07 # convert to Msun
rhoppm = np.load( 'text_data/rho_ppmsetup.npy' ) * 1.e3 # convert to g/cc
f = figure(figsize=(5, 8))
gs = gridspec.GridSpec(2, 1, height_ratios=[1.6, 1])
ax1 = subplot(gs[0])
ax2 = subplot(gs[1], sharex=ax1)
# Fine-tune figure; make subplots close to each other and hide x ticks for
# all but bottom plot.
f.subplots_adjust(hspace=0.1)
setp([a.get_xticklabels() for a in f.axes[:-1]], visible=False)
axes(ax1)
# plot simulation boundaries
rad1 = 3.0 # Mm
rad2 = 9.8 # Mm
rad = 10. ** p.get('logR') * ast.rsun_cm / 1.e8
#print rad
idx1 = np.abs(rad - rad1).argmin()
idx2 = np.abs(rad - rad2).argmin()
m1 = m[idx1]
m2 = m[idx2]
#print m1, m2
cb = u.colourblind
# plot conv zone boundaries
rlconv = 1.1593
ruconv = 1.8553
fill_between( [ rlconv, ruconv ], 0, 100, color='#fffdd8' ) #'f3f3df' )
rlconv = 2.0824
ruconv = 2.50
fill_between( [ rlconv, ruconv ], 0, 100, color='#b39eb5' )
plot(m,p.get('logP'),
label='$P_\mathrm{tot}$',
color=cb(8),
linestyle='-')
prad=ast.Prad(T,mu) # of course, the mu is not really used
plot(m,np.log10(prad),
label='$P_\mathrm{rad}$',
color=cb(0),
linestyle='--')
rho= 10.**p.get('logRho')
pideal = ast.Pgas(rho,T,mu)
plot(m,np.log10(pideal),
label='$P_\mathrm{gas}$',
color=cb(4),
linestyle='-.')
# use the difference:
pdeg = 10.** p.get('logP') - prad - pideal
plot(m,np.log10(pdeg),
label='$P_\mathrm{deg}$',
color=cb(5),
linestyle=':')
plot( Mrppm[1:], np.log10(Pppm[1:]), linestyle='-', \
marker='.', markevery=0.1, label='$P_\mathrm{PPM}$', \
color=cb(4), mfc='w', mec=cb(4), mew=0.4 )
xlim(0.5,2.5)
#xlim(0.8,2.3)
ylim(20.5,24.5)
ylabel('$\log_{10}(P\,/\,\mathrm{Ba})$')
#legend(loc=(0.75,0.6))
legend(loc='lower left')
ax=gca()
ax.yaxis.set_ticks([21,22,23,24,25])
draw()
ylim(21,24.5)
axes(ax2)
plot(m,p.get('logRho'),
#plot(m,10.**p.get('logRho'),
label='$\\rho_\mathrm{MESA}$',
color=cb(8),
linestyle='-')
plot( Mrppm[1:], np.log10(rhoppm[1:]), linestyle='-', \
#plot( Mrppm[1:], rhoppm[1:], linestyle='-', \
marker='.', markevery=0.1,
color=cb(4), mfc='w', mec=cb(4), mew=0.4,
label='$\\rho_\mathrm{PPM}$' )
rlconv = 1.1593
ruconv = 1.8553
fill_between( [ rlconv, ruconv ], 0, 100, color='#fffdd8' ) #'#f3f3df' )
rlconv = 2.0824
ruconv = 2.50
fill_between( [ rlconv, ruconv ], 0, 100, color='#b39eb5' )
#legend(loc=(0.75,0.6))
legend(loc='lower left')
ylim(4.5,7)
yticks(np.linspace(5,7,3))
ylabel('$\log_{10}(\\rho\,/\,\mathrm{g\,cm}^{-3})$')
xlabel('$\mathrm{Mass\,/\,M}_\odot$')
Out[7]:
In [8]:
ppm_path = ppm.set_YProf_path('/data/ppm_rpod2/YProfiles/O-shell-M25/')
mesa_path = '/data/ppm_rpod2/Stellar_models/O-shell-M25/M25Z0.02/LOGS_N2b'
ppm.plot_N2('D1', 0, 132,[3.5,9.,-0.5, 8.],[7.7,8.7,-0.15,1.3], mesa_path, mesa_model_num = 29350)
In [ ]: