In [ ]:
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.dpi'] = 200
matplotlib.rcParams['font.family'] = 'stixgeneral'
import matplotlib.pyplot as plt
import numpy as np
In [ ]:
g = 5.0/3
d = 1
b = 0.2
rr = np.logspace(-1, 1, 100)
T = rr**(d)
rho = rr**(-3*b)
P = rho*T
plt.plot(rr, P, label='pressure', alpha=0.7)
plt.plot(rr, rho, c='C1', label='density', alpha=0.7)
rho_ad = rr**(-(3*b-d)/g)
plt.plot(rr, rho_ad, '--', c='C1', label='adiabatic density')
plt.plot(rr, T, '-', c='C2', label='temperature', alpha=0.7)
T_ad = P/rho_ad
plt.plot(rr, T_ad, ':', c='C2', label='adiabatic temperature')
plt.grid(alpha=0.2)
plt.semilogy()
plt.semilogx()
plt.xlabel('R')
plt.legend(ncol=2)
In [ ]:
pp = np.linspace(1,4)
vv = 0.5/pp**(1/g)
vv2 = 4/pp**(1/g)
vv3 = 4/pp**(b/a)
plt.plot(vv, pp)
plt.plot(vv2, pp)
plt.plot(vv3, pp, label='cluster profile')
plt.legend()
plt.xlabel('V')
plt.ylabel('P')
In [ ]:
g = 5.0/3
a = 1
b = 2
X1=1
X2=3
rr_ext = np.linspace(X1, X2, 100)
Fg_ext = 1-rr_ext**(-b+a/g)
plt.plot(rr_ext, Fg_ext, color='C0')
rr = np.linspace(X1, X2, 100)
Fg = 1-rr**(-b+a/g)
plt.fill_between(rr, Fg, hatch="//", facecolor='none', edgecolor='C0')
plt.axhspan(0, Fg[-1], 0, rr[-1]/rr_ext[-1], color='C1', alpha=0.4)
#Fg2 = 1-(rr/X2)**(-b+a/g)
#plt.plot(rr, Fg2)
plt.xlim(0,rr_ext[-1])
plt.xlabel('x')
plt.xticks([0, 1, X2], [0, '1', r'$X=\frac{R}{R_0}$'])
plt.ylabel('$MgR_0$')
plt.yticks([])
plt.ylim(0, 0.9)
plt.axhline(0, ls=':', c='k', alpha=0.2)
plt.savefig('Eth_vs_work.pdf')
#plt.semilogx()
In [ ]:
import scipy.integrate as integrate
x_max = 3
d=0.9
colors = {0.0: 'C0', 0.5:'C1', 1.0: 'C2'}
betas = {0.0: [0.1, 0.2, 0.3, 0.5, 1.0],
0.5: [0.2, 0.3, 0.4, 0.5, 0.7, 1.0],
1.0: [0.4, 0.5, 0.6, 0.8, 1.0]}
fig, axes = plt.subplots(1, 3, sharey=True, figsize=(8,2.8), gridspec_kw={'wspace':0})
for iax, d in enumerate(colors.keys()):
color = colors[d]
first = True
axes[iax].semilogx()
axes[iax].semilogy()
axes[iax].grid(alpha=0.5, ls='--')
axes[iax].set_xlim(1, 300)
axes[iax].set_ylim(1E-1,1E2)
for b in betas[d]:
def Fg(xx):
return (1-xx**(-3*b+(3*b-d)/g))*xx**(d-1)
npoints = 200
ratio = np.ones(npoints)
xx = np.logspace(0, np.log10(300), npoints)
for i, x_max in enumerate(xx):
W, err = integrate.quad(Fg, 1, x_max)
Eth = 1/(3*b-d)/(g-1)*x_max*Fg(x_max)
ratio[i] = (Eth)/W
label = r'$\delta$=%.1f' % (d) if first else None
ls = '--' if b == 0.5 else '-'
line, = axes[iax].plot(xx, ratio, label=label, c=color, ls=ls, alpha=min(1.2-b, 1.0))
# Annotate beta values
itext = min(int(npoints*0.9), np.argmin(np.abs(ratio-0.05)[1:]))
text = axes[iax].annotate(r'$\beta$=%.1f' % b, xy=(xx[itext], ratio[itext]), xytext=(0, -2),
textcoords='offset points',
size=8, color=color,
horizontalalignment='center',
verticalalignment='bottom')
sp1 = axes[iax].transData.transform_point((xx[itext-8], ratio[itext-8]))
sp2 = axes[iax].transData.transform_point((xx[itext-4], ratio[itext-4]))
rise = (sp2[1] - sp1[1])
run = (sp2[0] - sp1[0])
slope_degrees = np.degrees(np.arctan2(rise, run))
text.set_rotation(slope_degrees)
first = False
#xi = 1/(3*b-d)/(g-1)*d -1
#plt.axhline(xi, color=line.get_color(), alpha=1.1-b)
#print('beta=%.1f, ratio = %.2f at X=%.1e, ratio_inf = %.2f' % (b, ratio[-1], xx[-1], xi))
axes[iax].legend(fontsize=9, frameon=False, loc='upper center')
axes[iax].axhline(1, ls=':', c='k', alpha=0.7)
if iax == 0:
axes[iax].set_ylabel(r'$\xi_{\rm{max}} \equiv \frac{E_{th}}{W}$')
if iax == 1:
axes[iax].set_xlabel(r'X $\equiv \frac{R}{R_0}$')
plt.savefig('energy_efficiency.pdf', bbox_inches='tight')
In [ ]:
import scipy.integrate as integrate
x_max = 3
d=0.9
colors = {0.0: 'C0', 0.5:'C1', 1.0: 'C2'}
betas = {0.0: [0.1, 0.2, 0.3, 0.5, 1.0],
0.5: [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 1.0],
1.0: [0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]}
fig, axes = plt.subplots(1, 3, sharey=True, figsize=(8,2.8), gridspec_kw={'wspace':0})
for iax, d in enumerate(colors.keys()):
color = colors[d]
first = True
axes[iax].semilogx()
#axes[iax].semilogy()
axes[iax].grid(alpha=0.5, ls='--')
axes[iax].set_xlim(1, 300)
axes[iax].set_ylim(-5,10)
for b in betas[d]:
def Fg(xx):
return (1-xx**(-3*b+(3*b-d)/g))*xx**(d-1)
npoints = 200
ratio = np.ones(npoints)
xx = np.logspace(0, np.log10(300), npoints)
for i, x_max in enumerate(xx):
W, err = integrate.quad(Fg, 1, x_max)
Eth = 1/(3*b-d)/(g-1)*x_max*Fg(x_max)
ratio[i] = (Eth-W)
label = r'$\delta$=%.1f' % (d) if first else None
ls = '--' if b == 0.5 else '-'
line, = axes[iax].plot(xx, ratio, label=label, c=color, ls=ls, alpha=min(1.2-b, 1.0))
# Annotate beta values
if ratio[-1] > 10:
itext = np.argmin(np.abs(ratio-10))-1
xoffset, yoffset = -12, -23
elif ratio[-1] > -4:
itext = len(ratio)-12
xoffset, yoffset = -5, 0
else:
itext = np.argmin(np.abs(ratio+5))-10
xoffset, yoffset = -4, 0
text = axes[iax].annotate(r'$\beta$=%.1f' % b,
xy=(xx[itext], ratio[itext]),
xytext=(xoffset, yoffset),
textcoords='offset points',
size=8, color=color,
horizontalalignment='center',
verticalalignment='bottom')
sp1 = axes[iax].transData.transform_point((xx[itext-8], ratio[itext-8]))
sp2 = axes[iax].transData.transform_point((xx[itext-4], ratio[itext-4]))
rise = (sp2[1] - sp1[1])
run = (sp2[0] - sp1[0])
slope_degrees = np.degrees(np.arctan2(rise, run))
text.set_rotation(slope_degrees)
first = False
#xi = 1/(3*b-d)/(g-1)*d -1
#plt.axhline(xi, color=line.get_color(), alpha=1.1-b)
#print('beta=%.1f, ratio = %.2f at X=%.1e, ratio_inf = %.2f' % (b, ratio[-1], xx[-1], xi))
axes[iax].legend(fontsize=9, frameon=False, loc=3)
axes[iax].axhline(0, ls='-', lw=1, c='k', alpha=0.7)
if iax == 0:
axes[iax].set_ylabel(r'$\Delta \epsilon_{\rm{max}} \equiv \frac{E_{th}-W}{M g_0 R_0}$')
if iax == 1:
axes[iax].set_xlabel(r'X $\equiv \frac{R}{R_0}$')
plt.savefig('energy_gain.pdf', bbox_inches='tight')
In [ ]:
g = 5.0/3
a = 1
X1=1
X2=10
rr = np.linspace(X1, X2, 100)
Fg = (1-rr**(-a+a/g))/rr
plt.fill_between(rr, Fg, hatch="//", facecolor='none', edgecolor='C0')
plt.axhspan(0, Fg[-1], 0, rr[-1], color='C1', alpha=0.4)
#Fg2 = 1-(rr/X2)**(-b+a/g)
#plt.plot(rr, Fg2)
plt.xlim(0,rr_ext[-1])
plt.xlabel('x')
plt.xticks([0, 1, X2], [0, '1', r'$X=\frac{R}{R_0}=%.0f$' % X2])
plt.ylabel('$MgR_0$')
#plt.yticks([])
#plt.ylim(0, 0.9)
plt.axhline(0, ls=':', c='k', alpha=0.2)
plt.savefig('Eth_vs_work.pdf')
#plt.semilogx()