<jupyter_start><jupyter_text>Data from H2Cooling with Gravity<jupyter_code>#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])
(HCG['T'][26]-HCG['T'][25])/1e6
NCG['RHO'][128,128,:60].argmax()
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)
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)
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)
kmtopc=3.24078e-14
1.3*kmtopc/np.sqrt(4.*np.pi*6.67e-8*1e-19)
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')
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)
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)<jupyter_output>NonLinear Fit: ln(p) = (-2.3+-0.12)*ln(x)+(2.5+-0.3)<jupyter_text>$$\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}$$<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Simulate Profile<jupyter_code>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
A=10.
B=0.002
a=2.3
A/B,B**(1/a)
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()
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)
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)
(P0*T0/rho(r))[0]
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()
plt.plot(r,mass(r),label=-2.3)
plt.ylabel('Mass (cu)')
plt.xlabel('Radius $(\si{pc})$')
plt.legend()
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)<jupyter_output>Total Mass: 3.660e+03 M☉
Total Mass: 1.479808e+02 cu<jupyter_text>$$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)$$<jupyter_code>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)))
integral(10.,1.)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2893170426065163, 'fraction_numerical': 0.08380325814536341, 'mean_word_length': 7.144132653061225, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '438085', 'license': 'MIT', 'n_tokens_mistral': 3521, 'n_tokens_neox': 3221, 'n_words': 306, 'path': 'Notebooks/Density Profile.ipynb', 'repo_name': 'Mixpap/MasterThesis', 'size': '6384'}