notebook.community



In [2]:

    
import numpy as np
import pyPLUTO as pp
from astropy.io import ascii
import os
import sys
from ipywidgets import interactive, widgets,fixed
from IPython.display import Audio, display
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.gridspec as gridspec
import matplotlib.ticker as ticker
from matplotlib.animation import FuncAnimation,FFMpegWriter
from matplotlib import rc,rcParams
from scipy.integrate import quad
rc('text', usetex=True)
rcParams['figure.figsize'] = (8, 6.5)
rcParams['ytick.labelsize'],rcParams['xtick.labelsize'] = 17.,17.
rcParams['axes.labelsize']=19.
rcParams['legend.fontsize']=17.
rcParams['text.latex.preamble'] = ['\\usepackage{siunitx}']
import seaborn
seaborn.despine()
seaborn.set_style('white', {'axes.linewidth': 0.5, 'axes.edgecolor':'black'})
seaborn.despine(left=True)
%load_ext autoreload









    



The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload






    





<matplotlib.figure.Figure at 0x7fe9302bb7d0>



In [3]:

    
%autoreload 1



In [4]:

    
%aimport f



In [10]:

    
NCG=np.load('../Data/noCoolingG.npz')
f.quadruple(NCG,np.log10(NCG['RHO']),rows=2,nn=10,tlim=29,Save_Figure='NoCoolGRquad')
f.pprofile(NCG,'RHO',steps=4,itlim=29,tdk='Myrs',Save_Figure='NoCoolGRHOprofile',sc2='log',xlim=[0,20])
f.pprofile(NCG,'Temp',steps=4,itlim=29,tdk='Myrs',yl='Temperature (K)',Save_Figure='NoCoolGTempprofile',sc2='log',xlim=[0,20])
f.pprofile(NCG,'PRS',steps=4,itlim=29,tdk='Myrs',yl='Pressure $(\si{dyne.cm^{-2}})$',y0=f.PRS0,Save_Figure='NoCoolGPRSprofile',sc2='log',xlim=[0,16])

Data from H2Cooling with Gravity



In [5]:

    
#HCG=np.load('../Data/H2CoolingG512.npz')
HCG=np.load('../Data/TabulatedG.npz')
f.quadruple(HCG,np.log10(HCG['RHO']),rows=2,nn=0,tlim=26,Save_Figure='H2CoolGRquad')
f.pprofile(HCG,'RHO',steps=4,itlim=26,tdk='Myrs',Save_Figure='H2CoolGRHOprofile',sc2='log',xlim=[0,20],yprop=256)
f.pprofile(HCG,'Temp',steps=4,itlim=26,tdk='Myrs',yl='Temperature (K)',Save_Figure='H2CoolGTempprofile',sc2='log',xlim=[0,20],yprop=256)
f.pprofile(HCG,'PRS',steps=4,itlim=26,tdk='Myrs',yl='Pressure $(\si{dyne.cm^{-2}})$',y0=f.PRS0,Save_Figure='H2CoolGPRSprofile',yprop=256,sc2='log',xlim=[0,15])



In [8]:

    
(HCG['T'][26]-HCG['T'][25])/1e6









    Out[8]:





0.095213802999999833



In [11]:

    
NCG['RHO'][128,128,:60].argmax()









    Out[11]:





23



In [19]:

    
plt.figure(figsize=(12,6))
plt.plot(NCG['T'][:38]/1e6,NCG['PRS'][128,128,:][:38]*f.Temp0/NCG['RHO'][128,128,:][:38],label='No Cooling Simulation Data')
plt.plot(HCG['T']/1e6,HCG['PRS'][256,256,:]*f.Temp0/HCG['RHO'][256,256,:],label='Cooling Simulation Data')
plt.yscale('log')
plt.ylabel('Temperature (K)')
plt.xlabel('$\si{Myrs}$')
plt.legend()
datafolder='../Document/DataImages/'
plt.savefig(datafolder+'GRcenterTimeTemp.png',bbox_inches='tight',format='png', dpi=100)



In [20]:

    
plt.figure(figsize=(12,6))
plt.plot(NCG['T'][:38]/1e6,NCG['RHO'][128,128,:][:38]*f.RHO0,label='No Cooling Simulation Data')
plt.plot(HCG['T']/1e6,HCG['RHO'][256,256,:]*f.RHO0,label='Cooling Simulation Data')
plt.yscale('log')
plt.ylabel('$ \\rho\, (\si{g.cm^{-3}})$')
plt.xlabel('$\si{Myrs}$')
plt.legend()
datafolder='../Document/DataImages/'
plt.savefig(datafolder+'GRcenterTimeRHO.png',bbox_inches='tight',format='png', dpi=100)



In [21]:

    
mid=256
ind= HCG['RHO'][mid,mid,:].argmax()
print 'Maximum Density: {:.2e} g/cm3 ({:.2e} cu) at time {:.2f} Myrs (index: {})'.format(HCG['RHO'][mid,mid:,ind].max()*f.RHO0,HCG['RHO'][mid,mid,:].max(),HCG['T'][ind]/1e6,ind)

xs=HCG['X'][mid:]#/10.
ys=HCG['RHO'][mid,mid:,ind]#*f.RHO0
ps=HCG['PRS'][mid,mid:,ind]
cs=np.sqrt((5./3)*ps*f.PRS0/(ys*f.RHO0))
#cs=np.sqrt(ps*f.PRS0/(ys*f.RHO0))
xspc=xs*10.
xau=xs*206265.
yscgs=ys*f.RHO0
ybe=cs**2/(2.*np.pi*6.67e-8*(xspc*3.086e18)**2)
ax=plt.subplot(111)
ax.plot(xspc,yscgs,label='Cooling Simulation Data')
ax.plot(xspc[xspc<2.],ybe[xspc<2.],label='Bonnor-Ebert Sphere')

ax.plot(NCG['X'][128:]*10.,NCG['RHO'][128,128:,23]*f.RHO0,'--',alpha=1.,label='No Cooling Simulation Data')

ax.set_yscale('log')
ax.set_ylabel('$ \\rho\, (\si{g.cm^{-3}})$')
ax.set_xlabel('$\si{pc}$')
#ax.vlines([0.11267],1e-21,1.3e-17,linewidth=0.5,linestyles='--',label='Critical Radius')
ax.set_xlim(0,5.2)
plt.legend()
datafolder='../Document/DataImages/'
plt.savefig(datafolder+'H2CoolGRHOprofile-BE.png',bbox_inches='tight',format='png', dpi=100)









    



Maximum Density: 1.29e-17 g/cm3 (7.70e+06 cu) at time 2.47 Myrs (index: 26)



In [15]:

    
kmtopc=3.24078e-14
1.3*kmtopc/np.sqrt(4.*np.pi*6.67e-8*1e-19)









    Out[15]:





0.14552079321216815



In [91]:

    
plt.plot(xspc,cs*1e-5)
plt.ylabel('$C_s \, (\si{km.s^{-1}})$')
plt.xlabel('$\si{pc}$')
plt.xlim(0,8)
#plt.yscale('log')









    Out[91]:





(0, 8)



In [92]:

    
ax1=plt.subplot(111)
ax1.plot(xs,ys)
ax1.set_yscale('log')
ax1.set_ylabel('$ \\rho\, (\si{cu})$')
ax1.set_xlabel('$x \, (\si{cu}) $')
ax1.set_xlim(0,0.8)









    Out[92]:





(0, 0.8)



In [18]:

    
regions= [[0.0,0.02],[0.04,0.2],[0.21,0.36],[0.01,0.35]]
regions=[[0.015,0.29]]
for region in regions:
    xmin,xmax=region[0],region[1]
    whr=np.logical_and(xs>xmin,xs<xmax)
    xsl=np.log(xs[whr])
    ysl=np.log(ys[whr])
    from scipy.optimize import curve_fit
    def ff(x,a,b,): return a*x+b 
    p,dp2=curve_fit(ff,xsl,ysl,[-2.1,100.])
    dp=np.sqrt(np.diag(dp2))
    print 'NonLinear Fit: ln(p) = ({:.1f}+-{:.2f})*ln(x)+({:.1f}+-{:.1f})'.format(p[0],dp[0],p[1],dp[1])
    xx=np.linspace(xmin,xmax,100)   
    plt.plot(xx,np.exp(p[1])*xx**(p[0]),'--',label='{:.1f}'.format(p[0]))
plt.plot(xs,ys)
plt.yscale('log')
plt.ylabel('$ \\rho\, (\si{cu})$')
plt.xlabel('$x\, (\si{cu})$')
plt.legend()
plt.xlim(None,0.45)









    



NonLinear Fit: ln(p) = (-2.3+-0.12)*ln(x)+(2.5+-0.3)






    Out[18]:





(-0.095703125, 0.45)

$$ \rho_{BE} (r) = \frac{C_s^2}{4 \pi G r^2} =\frac{L_0^2 \rho _0 t_0}{t_0^2 L_0^2} $$



In [94]:

    
xs=HCG['X'][mid:]
ysP=HCG['PRS'][mid,mid:,ind]
xscgs=xs*10.
xau=xscgs*206265.
ysPcgs=ysP*f.PRS0
ax=plt.subplot(111)
ax.plot(xscgs,ysPcgs)
ax.set_yscale('log')
ax.set_ylabel('$ P\, (\si{dyne.cm^{-2}})$')
ax.set_xlabel('$\si{pc}$')
ax.set_xlim(None,10)









    Out[94]:





(-0.95703125, 10)

Simulate Profile



In [95]:

    
Radius = 1.0
Density1 = 10.
P0=1e-8
T0=10891304347826.088
def rho(r,a=2.3,rho1=10.,R=1.,rc=0.002): return np.piecewise(r, [r < R , r >= R], [lambda r: rho1/(rc+r**a), 1.])
def mass(r,sm=False): 
    return 4.*np.pi*r**2 *rho(r)*24.73 if sm else 4.*np.pi*r**2 *rho(r)
#def massmo(r,a,rc): return mass(r,a=-2.3,rho0=1e5,rc=0.005,R=10.) *24.73 #(10pc)^3 * hydrogen_mass /cm^3 = 24.73 Mo



In [96]:

    
A=10.
B=0.002
a=2.3



In [97]:

    
A/B,B**(1/a)









    Out[97]:





(5000.0, 0.06707099987476356)



In [98]:

    
r=np.linspace(-2,3,500)
plt.ylabel('Number Density $(\si{cm^{-3}})$')
plt.xlabel('Radius $(\si{pc})$')
plt.yscale('log')
plt.plot(r*10.,rho(r),label=-2.3)
#plt.plot(r*10.,10./r**(2.3),label=-2.3)
#plt.vlines(0.07,1e2,1e5,linewidth=0.5)
plt.legend()









    



/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in power
  """






    Out[98]:





<matplotlib.legend.Legend at 0x7f5e8cf30d10>



In [113]:

    
r=np.linspace(0,3,500)
plt.ylabel('$ \\rho\, (\si{g.cm^{-3}})$')
plt.xlabel('Radius $(\si{pc})$')
plt.yscale('log')
plt.plot(r*10.,rho(r)*f.RHO0,label='{} Density Profile'.format(-2.3))
plt.legend()
plt.xlim(0,16)
datafolder='../Document/DataImages/'
plt.savefig(datafolder+'SimRHOProfile.png',bbox_inches='tight',format='png', dpi=100)



In [114]:

    
plt.ylabel('Temperature$(\si{K})$')
plt.xlabel('Radius $(\SI{10}{pc})$')
plt.yscale('log')
plt.plot(r*10,P0*T0/rho(r),label=-2.3)
#plt.legend()
plt.xlim(0,16)
datafolder='../Document/DataImages/'
plt.savefig(datafolder+'SimTMPProfile.png',bbox_inches='tight',format='png', dpi=100)



In [116]:

    
(P0*T0/rho(r))[0]









    Out[116]:





21.782608695652176



In [101]:

    
plt.plot(r*10,mass(r,sm=True),label=-2.3)
#plt.plot(r,4.*np.pi*r**2 *rho(r),label=-2.3)
plt.ylabel('Mass $(\si{M_\odot})$')
plt.xlabel('Radius $(\si{pc})$')
plt.legend()









    



/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in power
  """






    Out[101]:





<matplotlib.legend.Legend at 0x7f5e8d4d6250>



In [102]:

    
plt.plot(r,mass(r),label=-2.3)
plt.ylabel('Mass (cu)')
plt.xlabel('Radius $(\si{pc})$')
plt.legend()









    



/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in power
  """






    Out[102]:





<matplotlib.legend.Legend at 0x7f5e8d0c79d0>



In [85]:

    
TotalMass=quad(mass,0,1.,args=True)[0]
print u'Total Mass: {:.3e} M☉'.format(TotalMass)
TotalMass=quad(mass,0,1.)[0]
print u'Total Mass: {:e} cu'.format(TotalMass)









    



Total Mass: 3.660e+03 M☉
Total Mass: 1.479808e+02 cu

$$ 4\pi \rho_0 \int _0 ^R \frac{R_c r^2}{R_c +r^2}dr = 4\pi \rho_0 Rc \Big( R-\sqrt{R_c} \arctan\big( \frac{R}{\sqrt{R_C}}\big) \Big) $$



In [139]:

    
def integral(rho1,R):
    return 17.952*rho1*R**0.7
    #return 4.*np.pi*rho0*Rc*(R-np.sqrt(Rc)*np.arctan(R/np.sqrt(Rc)))



In [140]:

    
integral(10.,1.)









    Out[140]:





179.52