In [1]:
import numpy
import tqdm
import warnings
warnings.simplefilter('ignore',DeprecationWarning)
import wendy
%pylab inline
import matplotlib.animation as animation
from matplotlib import cm
from IPython.display import HTML
_SAVE_GIFS= False
rcParams.update({'axes.labelsize': 17.,
'font.size': 12.,
'legend.fontsize': 17.,
'xtick.labelsize':15.,
'ytick.labelsize':15.,
'text.usetex': _SAVE_GIFS,
'figure.figsize': [5,5],
'xtick.major.size' : 4,
'ytick.major.size' : 4,
'xtick.minor.size' : 2,
'ytick.minor.size' : 2,
'legend.numpoints':1})
import copy
numpy.random.seed(2)
Let's investigate how resilient an exponential disk ($\rho(x) \propto e^{-|x|/H}$) is against non-adiabatic changes to the potential. We will first setup an exponential disk in approximate equilibrium, let it evolve for a while to let it fully equilibrate, and then inject energy into it in both an adiabatic and non-adiabatic manner. First we setup the disk and evolve it:
In [2]:
# Initial disk
N= 30000
totmass= 1. # Sigma above
sigma= 1.
zh= sigma**2./totmass # twopiG = 1. in our units
tdyn= zh/sigma
x= -numpy.log(1.-numpy.random.uniform(size=N))*zh
x[numpy.random.uniform(size=N) < 0.5]*= -1.
v= numpy.random.normal(size=N)*sigma*numpy.sqrt(1.-0.75*numpy.exp(-numpy.fabs(x)/zh))
v-= numpy.mean(v) # stabilize
m= totmass*numpy.ones_like(x)/N
In [3]:
g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)
In [4]:
nt= 3000
xt= numpy.empty((N,nt+1))
vt= numpy.empty((N,nt+1))
Et= numpy.empty((nt+1))
xt[:,0]= x
vt[:,0]= v
Et[0]= wendy.energy(x,v,m)
for ii in tqdm.trange(nt):
tx,tv= next(g)
xt[:,ii+1]= tx
vt[:,ii+1]= tv
Et[ii+1]= wendy.energy(tx,tv,m)
x_start= xt[:,-1]
v_start= vt[:,-1]
The following movie shows that the exponential disk density does not significantly evolve:
In [5]:
figsize(6,4)
fig, ax= subplots()
ii= 0
a= ax.hist(xt[:,ii],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16./N*numpy.ones(N))
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,.65)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\rho(x)$')
ax.set_yscale('log')
ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',
horizontalalignment='right',verticalalignment='top',size=18.)
subsamp= 5
def animate(ii):
ax.clear()
norm= 1./N
a= ax.hist(xt[:,ii*subsamp],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16.*norm*numpy.ones(N))
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,.65)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\rho(x)$')
#ax.set_yscale('log')
ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),
(0.95,0.95),xycoords='axes fraction',
horizontalalignment='right',verticalalignment='top',size=18.)
return a[2]
anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,
frames=nt//subsamp,interval=40,blit=True,repeat=True)
# The following is necessary to just get the movie, and not an additional initial frame
plt.close()
out= HTML(anim.to_html5_video())
plt.close()
out
Out[5]:
Now we inject energy into about 40% of the orbits in an adiabatic manner, by changing the velocities over many dynamical times:
In [6]:
nt= 3000
x= copy.deepcopy(x_start)
v= copy.deepcopy(v_start)
xt= numpy.empty((N,nt+1))
vt= numpy.empty((N,nt+1))
xt[:,0]= x
vt[:,0]= v
launch_v= 2.*sigma
frac_lost= 0.4
time_lost= 100
nlost= int(round(frac_lost*N/4./time_lost))
g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)
for ii in tqdm.trange(nt):
if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \
or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):
idx= numpy.random.permutation(N)[:nlost]#numpy.argpartition(numpy.fabs(x),nlost)[:nlost]
if numpy.random.uniform() < 0.5:
v[idx[:nlost//2]]= launch_v
v[idx[nlost//2:]]= -launch_v
else:
v[idx[:nlost//2]]= -launch_v
v[idx[nlost//2:]]= launch_v
g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)
x,v= next(g)
if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \
or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):
xt[:,ii+1]= x
vt[:,ii+1]= v
vt[idx,ii+1]= numpy.nan
else:
xt[:,ii+1]= x
vt[:,ii+1]= v
The density remains approximately exponential throughout the evolution (energy is injected between $t = 15$ and $t = 35$; note that we remove the particles with inflated energies from the histogram to focus on the response of the system):
In [7]:
figsize(6,4)
fig, ax= subplots()
ii= 0
a= ax.hist(xt[:,ii],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16./N*numpy.ones(N))
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,.65)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\rho(x)$')
ax.set_yscale('log')
ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',
horizontalalignment='right',verticalalignment='top',size=18.)
ax.annotate(r'$\mathrm{adiabatic}$',(0.05,0.95),xycoords='axes fraction',
horizontalalignment='left',verticalalignment='top',size=18.)
subsamp= 5
def animate(ii):
ax.clear()
if ii*subsamp < 300:
norm= 1./N
elif ii*subsamp < 700.:
norm= 1./(N-(ii*subsamp-300.)*nlost)
else:
norm= 1./(N-400.*nlost)
idx= True-numpy.any(numpy.isnan(vt[:,:ii*subsamp]),axis=1)
a= ax.hist(xt[idx,ii*subsamp]-numpy.mean(xt[idx,ii*subsamp]),
bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],
weights=101./16.*norm*numpy.ones(N)[idx])
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,.65)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\rho(x)$')
#ax.set_yscale('log')
ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),
(0.95,0.95),xycoords='axes fraction',
horizontalalignment='right',verticalalignment='top',size=18.)
ax.annotate(r'$\mathrm{adiabatic}$',(0.05,0.95),xycoords='axes fraction',
horizontalalignment='left',verticalalignment='top',size=18.)
return a[2]
anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,
frames=nt//subsamp,interval=40,blit=True,repeat=True)
if _SAVE_GIFS:
anim.save('expdiskadiab_density.gif',writer='imagemagick',dpi=80)
# The following is necessary to just get the movie, and not an additional initial frame
plt.close()
out= HTML(anim.to_html5_video())
plt.close()
out
Out[7]:
Now we inject the same amount of energy non-adiabatically, by doing it in brief episodes rather than over multiple dynamical times:
In [8]:
nt= 3000
x= copy.deepcopy(x_start)
v= copy.deepcopy(v_start)
xt= numpy.empty((N,nt+1))
vt= numpy.empty((N,nt+1))
xt[:,0]= x
vt[:,0]= v
launch_v= 2.*sigma
frac_lost= 0.4
time_lost= 10
nlost= int(round(frac_lost*N/4./time_lost))
g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)
for ii in tqdm.trange(nt):
if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \
or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):
idx= numpy.random.permutation(N)[:nlost]#numpy.argpartition(numpy.fabs(x),nlost)[:nlost]
if numpy.random.uniform() < 0.5:
v[idx[:nlost//2]]= launch_v
v[idx[nlost//2:]]= -launch_v
else:
v[idx[:nlost//2]]= -launch_v
v[idx[nlost//2:]]= launch_v
g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)
x,v= next(g)
if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \
or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):
xt[:,ii+1]= x
vt[:,ii+1]= v
vt[idx,ii+1]= numpy.nan
else:
xt[:,ii+1]= x
vt[:,ii+1]= v
The density now clearly evolves non-adiabatically and has large fluctuations (energy is injected around $t = 15, 20, 25$, and $30$), but in the end the density settles back down to the original exponential distribution. An exponential disk is thus remarkably resilient to energy injection!
In [9]:
figsize(6,4)
fig, ax= subplots()
ii= 0
a= ax.hist(xt[:,ii],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16./N*numpy.ones(N))
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,.65)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\rho(x)$')
ax.set_yscale('log')
ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',
horizontalalignment='right',verticalalignment='top',size=18.)
ax.annotate(r'$\mathrm{non\ adiabatic}$',(0.05,0.95),xycoords='axes fraction',
horizontalalignment='left',verticalalignment='top',size=18.)
subsamp= 5
def animate(ii):
ax.clear()
if ii*subsamp < 300:
norm= 1./N
elif ii*subsamp < 300+time_lost:
norm= 1./(N-(ii*subsamp-300.)*nlost)
elif ii*subsamp < 400.:
norm= 1./(N-time_lost*nlost)
elif ii*subsamp < 400+time_lost:
norm= 1./(N-2.*time_lost*nlost-(ii*subsamp-400.)*nlost)
elif ii*subsamp < 500.:
norm= 1./(N-2.*time_lost*nlost)
elif ii*subsamp < 500+time_lost:
norm= 1./(N-2.*time_lost*nlost-(ii*subsamp-500.)*nlost)
elif ii*subsamp < 600:
norm= 1./(N-3.*time_lost*nlost)
elif ii*subsamp < 600+time_lost:
norm= 1./(N-3.*time_lost*nlost-(ii*subsamp-600.)*nlost)
else:
norm= 1./(N-4.*time_lost*nlost)
idx= True-numpy.any(numpy.isnan(vt[:,:ii*subsamp]),axis=1)
a= ax.hist(xt[idx,ii*subsamp]-numpy.mean(xt[idx,ii*subsamp]),
bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],
weights=101./16.*norm*numpy.ones(N)[idx])
xs= numpy.linspace(-8.,8.,101)
ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)
ax.set_xlim(-8.,8.)
ax.set_ylim(10.**-3.,.65)
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$\rho(x)$')
#ax.set_yscale('log')
ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),
(0.95,0.95),xycoords='axes fraction',
horizontalalignment='right',verticalalignment='top',size=18.)
ax.annotate(r'$\mathrm{non\ adiabatic}$',(0.05,0.95),xycoords='axes fraction',
horizontalalignment='left',verticalalignment='top',size=18.)
return a[2]
anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,
frames=nt//subsamp,interval=40,blit=True,repeat=True)
if _SAVE_GIFS:
anim.save('expdisknonadiab_density.gif',writer='imagemagick',dpi=80)
# The following is necessary to just get the movie, and not an additional initial frame
plt.close()
out= HTML(anim.to_html5_video())
plt.close()
out
Out[9]:
In [ ]: