{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"source": [
"import numpy\n",
"import tqdm\n",
"import warnings\n",
"warnings.simplefilter('ignore',DeprecationWarning)\n",
"import wendy\n",
"%pylab inline\n",
"import matplotlib.animation as animation\n",
"from matplotlib import cm\n",
"from IPython.display import HTML\n",
"_SAVE_GIFS= False\n",
"rcParams.update({'axes.labelsize': 17.,\n",
" 'font.size': 12.,\n",
" 'legend.fontsize': 17.,\n",
" 'xtick.labelsize':15.,\n",
" 'ytick.labelsize':15.,\n",
" 'text.usetex': _SAVE_GIFS,\n",
" 'figure.figsize': [5,5],\n",
" 'xtick.major.size' : 4,\n",
" 'ytick.major.size' : 4,\n",
" 'xtick.minor.size' : 2,\n",
" 'ytick.minor.size' : 2,\n",
" 'legend.numpoints':1})\n",
"import copy\n",
"numpy.random.seed(2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Adiabatic versus non-adiabatic changes to an exponential disk\n",
"\n",
"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:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Initial disk\n",
"N= 30000\n",
"totmass= 1. # Sigma above\n",
"sigma= 1.\n",
"zh= sigma**2./totmass # twopiG = 1. in our units\n",
"tdyn= zh/sigma\n",
"x= -numpy.log(1.-numpy.random.uniform(size=N))*zh\n",
"x[numpy.random.uniform(size=N) < 0.5]*= -1.\n",
"v= numpy.random.normal(size=N)*sigma*numpy.sqrt(1.-0.75*numpy.exp(-numpy.fabs(x)/zh))\n",
"v-= numpy.mean(v) # stabilize\n",
"m= totmass*numpy.ones_like(x)/N"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 3000/3000 [00:28<00:00, 104.44it/s]\n"
]
}
],
"source": [
"nt= 3000\n",
"xt= numpy.empty((N,nt+1))\n",
"vt= numpy.empty((N,nt+1))\n",
"Et= numpy.empty((nt+1))\n",
"xt[:,0]= x\n",
"vt[:,0]= v\n",
"Et[0]= wendy.energy(x,v,m)\n",
"for ii in tqdm.trange(nt):\n",
" tx,tv= next(g)\n",
" xt[:,ii+1]= tx\n",
" vt[:,ii+1]= tv\n",
" Et[ii+1]= wendy.energy(tx,tv,m)\n",
"x_start= xt[:,-1]\n",
"v_start= vt[:,-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following movie shows that the exponential disk density does not significantly evolve:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"figsize(6,4)\n",
"fig, ax= subplots()\n",
"ii= 0\n",
"a= ax.hist(xt[:,ii],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16./N*numpy.ones(N))\n",
"xs= numpy.linspace(-8.,8.,101)\n",
"ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)\n",
"ax.set_xlim(-8.,8.)\n",
"ax.set_ylim(10.**-3.,.65)\n",
"ax.set_xlabel(r'$x$')\n",
"ax.set_ylabel(r'$\\rho(x)$')\n",
"ax.set_yscale('log')\n",
"ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',\n",
" horizontalalignment='right',verticalalignment='top',size=18.)\n",
"subsamp= 5\n",
"def animate(ii):\n",
" ax.clear()\n",
" norm= 1./N\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))\n",
" xs= numpy.linspace(-8.,8.,101)\n",
" ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)\n",
" ax.set_xlim(-8.,8.)\n",
" ax.set_ylim(10.**-3.,.65)\n",
" ax.set_xlabel(r'$x$')\n",
" ax.set_ylabel(r'$\\rho(x)$')\n",
" #ax.set_yscale('log')\n",
" ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),\n",
" (0.95,0.95),xycoords='axes fraction',\n",
" horizontalalignment='right',verticalalignment='top',size=18.)\n",
" return a[2]\n",
"anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,\n",
" frames=nt//subsamp,interval=40,blit=True,repeat=True)\n",
"# The following is necessary to just get the movie, and not an additional initial frame\n",
"plt.close()\n",
"out= HTML(anim.to_html5_video())\n",
"plt.close()\n",
"out"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we inject energy into about 40% of the orbits in an *adiabatic manner*, by changing the velocities over many dynamical times:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 3000/3000 [00:32<00:00, 91.25it/s] \n"
]
}
],
"source": [
"nt= 3000\n",
"x= copy.deepcopy(x_start)\n",
"v= copy.deepcopy(v_start)\n",
"xt= numpy.empty((N,nt+1))\n",
"vt= numpy.empty((N,nt+1))\n",
"xt[:,0]= x\n",
"vt[:,0]= v\n",
"launch_v= 2.*sigma\n",
"frac_lost= 0.4\n",
"time_lost= 100\n",
"nlost= int(round(frac_lost*N/4./time_lost))\n",
"g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)\n",
"for ii in tqdm.trange(nt):\n",
" if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \\\n",
" or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):\n",
" idx= numpy.random.permutation(N)[:nlost]#numpy.argpartition(numpy.fabs(x),nlost)[:nlost]\n",
" if numpy.random.uniform() < 0.5:\n",
" v[idx[:nlost//2]]= launch_v\n",
" v[idx[nlost//2:]]= -launch_v\n",
" else:\n",
" v[idx[:nlost//2]]= -launch_v\n",
" v[idx[nlost//2:]]= launch_v\n",
" g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)\n",
" x,v= next(g)\n",
" if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \\\n",
" or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):\n",
" xt[:,ii+1]= x\n",
" vt[:,ii+1]= v\n",
" vt[idx,ii+1]= numpy.nan\n",
" else:\n",
" xt[:,ii+1]= x\n",
" vt[:,ii+1]= v"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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):"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"figsize(6,4)\n",
"fig, ax= subplots()\n",
"ii= 0\n",
"a= ax.hist(xt[:,ii],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16./N*numpy.ones(N))\n",
"xs= numpy.linspace(-8.,8.,101)\n",
"ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)\n",
"ax.set_xlim(-8.,8.)\n",
"ax.set_ylim(10.**-3.,.65)\n",
"ax.set_xlabel(r'$x$')\n",
"ax.set_ylabel(r'$\\rho(x)$')\n",
"ax.set_yscale('log')\n",
"ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',\n",
" horizontalalignment='right',verticalalignment='top',size=18.)\n",
"ax.annotate(r'$\\mathrm{adiabatic}$',(0.05,0.95),xycoords='axes fraction',\n",
" horizontalalignment='left',verticalalignment='top',size=18.)\n",
"subsamp= 5\n",
"def animate(ii):\n",
" ax.clear()\n",
" if ii*subsamp < 300:\n",
" norm= 1./N\n",
" elif ii*subsamp < 700.:\n",
" norm= 1./(N-(ii*subsamp-300.)*nlost)\n",
" else:\n",
" norm= 1./(N-400.*nlost)\n",
" idx= True-numpy.any(numpy.isnan(vt[:,:ii*subsamp]),axis=1)\n",
" a= ax.hist(xt[idx,ii*subsamp]-numpy.mean(xt[idx,ii*subsamp]),\n",
" bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],\n",
" weights=101./16.*norm*numpy.ones(N)[idx])\n",
" xs= numpy.linspace(-8.,8.,101)\n",
" ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)\n",
" ax.set_xlim(-8.,8.)\n",
" ax.set_ylim(10.**-3.,.65)\n",
" ax.set_xlabel(r'$x$')\n",
" ax.set_ylabel(r'$\\rho(x)$')\n",
" #ax.set_yscale('log')\n",
" ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),\n",
" (0.95,0.95),xycoords='axes fraction',\n",
" horizontalalignment='right',verticalalignment='top',size=18.)\n",
" ax.annotate(r'$\\mathrm{adiabatic}$',(0.05,0.95),xycoords='axes fraction',\n",
" horizontalalignment='left',verticalalignment='top',size=18.)\n",
" return a[2]\n",
"anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,\n",
" frames=nt//subsamp,interval=40,blit=True,repeat=True)\n",
"if _SAVE_GIFS:\n",
" anim.save('expdiskadiab_density.gif',writer='imagemagick',dpi=80)\n",
"# The following is necessary to just get the movie, and not an additional initial frame\n",
"plt.close()\n",
"out= HTML(anim.to_html5_video())\n",
"plt.close()\n",
"out"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we inject the same amount of energy *non-adiabatically*, by doing it in brief episodes rather than over multiple dynamical times:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"100%|██████████| 3000/3000 [00:19<00:00, 156.08it/s]\n"
]
}
],
"source": [
"nt= 3000\n",
"x= copy.deepcopy(x_start)\n",
"v= copy.deepcopy(v_start)\n",
"xt= numpy.empty((N,nt+1))\n",
"vt= numpy.empty((N,nt+1))\n",
"xt[:,0]= x\n",
"vt[:,0]= v\n",
"launch_v= 2.*sigma\n",
"frac_lost= 0.4\n",
"time_lost= 10\n",
"nlost= int(round(frac_lost*N/4./time_lost))\n",
"g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)\n",
"for ii in tqdm.trange(nt):\n",
" if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \\\n",
" or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):\n",
" idx= numpy.random.permutation(N)[:nlost]#numpy.argpartition(numpy.fabs(x),nlost)[:nlost]\n",
" if numpy.random.uniform() < 0.5:\n",
" v[idx[:nlost//2]]= launch_v\n",
" v[idx[nlost//2:]]= -launch_v\n",
" else:\n",
" v[idx[:nlost//2]]= -launch_v\n",
" v[idx[nlost//2:]]= launch_v\n",
" g= wendy.nbody(x,v,m,0.05*tdyn,approx=True,nleap=1)\n",
" x,v= next(g)\n",
" if (ii >= 300 and ii < 300+time_lost) or (ii >= 400 and ii < 400+time_lost) \\\n",
" or (ii >= 500 and ii < 500+time_lost) or (ii >= 600 and ii < 600+time_lost):\n",
" xt[:,ii+1]= x\n",
" vt[:,ii+1]= v\n",
" vt[idx,ii+1]= numpy.nan\n",
" else:\n",
" xt[:,ii+1]= x\n",
" vt[:,ii+1]= v"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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!"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"figsize(6,4)\n",
"fig, ax= subplots()\n",
"ii= 0\n",
"a= ax.hist(xt[:,ii],bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],weights=101./16./N*numpy.ones(N))\n",
"xs= numpy.linspace(-8.,8.,101)\n",
"ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)\n",
"ax.set_xlim(-8.,8.)\n",
"ax.set_ylim(10.**-3.,.65)\n",
"ax.set_xlabel(r'$x$')\n",
"ax.set_ylabel(r'$\\rho(x)$')\n",
"ax.set_yscale('log')\n",
"ax.annotate(r'$t=0$',(0.95,0.95),xycoords='axes fraction',\n",
" horizontalalignment='right',verticalalignment='top',size=18.)\n",
"ax.annotate(r'$\\mathrm{non\\ adiabatic}$',(0.05,0.95),xycoords='axes fraction',\n",
" horizontalalignment='left',verticalalignment='top',size=18.)\n",
"subsamp= 5\n",
"def animate(ii):\n",
" ax.clear()\n",
" if ii*subsamp < 300:\n",
" norm= 1./N\n",
" elif ii*subsamp < 300+time_lost:\n",
" norm= 1./(N-(ii*subsamp-300.)*nlost)\n",
" elif ii*subsamp < 400.:\n",
" norm= 1./(N-time_lost*nlost)\n",
" elif ii*subsamp < 400+time_lost:\n",
" norm= 1./(N-2.*time_lost*nlost-(ii*subsamp-400.)*nlost)\n",
" elif ii*subsamp < 500.:\n",
" norm= 1./(N-2.*time_lost*nlost)\n",
" elif ii*subsamp < 500+time_lost:\n",
" norm= 1./(N-2.*time_lost*nlost-(ii*subsamp-500.)*nlost)\n",
" elif ii*subsamp < 600:\n",
" norm= 1./(N-3.*time_lost*nlost)\n",
" elif ii*subsamp < 600+time_lost:\n",
" norm= 1./(N-3.*time_lost*nlost-(ii*subsamp-600.)*nlost)\n",
" else:\n",
" norm= 1./(N-4.*time_lost*nlost)\n",
" idx= True-numpy.any(numpy.isnan(vt[:,:ii*subsamp]),axis=1)\n",
" a= ax.hist(xt[idx,ii*subsamp]-numpy.mean(xt[idx,ii*subsamp]),\n",
" bins=101,histtype='step',lw=1.,color='k',range=[-8.,8.],\n",
" weights=101./16.*norm*numpy.ones(N)[idx])\n",
" xs= numpy.linspace(-8.,8.,101)\n",
" ax.plot(xs,totmass/2./zh*numpy.exp(-numpy.fabs(xs)/2./zh)**2.,'b--',lw=2.,zorder=0)\n",
" ax.set_xlim(-8.,8.)\n",
" ax.set_ylim(10.**-3.,.65)\n",
" ax.set_xlabel(r'$x$')\n",
" ax.set_ylabel(r'$\\rho(x)$')\n",
" #ax.set_yscale('log')\n",
" ax.annotate(r'$t=%.0f$' % (ii*subsamp/20.),\n",
" (0.95,0.95),xycoords='axes fraction',\n",
" horizontalalignment='right',verticalalignment='top',size=18.)\n",
" ax.annotate(r'$\\mathrm{non\\ adiabatic}$',(0.05,0.95),xycoords='axes fraction',\n",
" horizontalalignment='left',verticalalignment='top',size=18.)\n",
" return a[2]\n",
"anim = animation.FuncAnimation(fig,animate,#init_func=init_anim_frame,\n",
" frames=nt//subsamp,interval=40,blit=True,repeat=True)\n",
"if _SAVE_GIFS:\n",
" anim.save('expdisknonadiab_density.gif',writer='imagemagick',dpi=80)\n",
"# The following is necessary to just get the movie, and not an additional initial frame\n",
"plt.close()\n",
"out= HTML(anim.to_html5_video())\n",
"plt.close()\n",
"out"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.5.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}