{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook analyzes the simulations of the vortex expansion. The data should be generated by running commands from the file ``vortex_expansion.py``.\n", "$\\newcommand{d}{\\mathrm{d}}$" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pylab inline\n", "import os.path\n", "import vortex_expansion\n", "import h5py\n", "PREFIX = '_vortex_expansion'\n", "DATA = os.path.join(PREFIX, 'data')\n", "files = !ls $DATA # This is an IPython feature (calling the shell command ls)\n", "files = [_f for _f in files if _f.endswith('.hd5')]\n", "print files\n", "\n", "# Get a representative p object for conversions etc.\n", "p = vortex_expansion._MIT.get_realistic(minimize=False, fast_load=False, dim=2, basis='DVR')[1]" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "#%load https://bitbucket.org/mforbes/notes_ipython/raw/tip/notebook_utils/animation.py" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "%pylab inline\n", "from tempfile import NamedTemporaryFile\n", "import IPython.display\n", "import matplotlib.animation\n", "\n", "VIDEO_TAG = \"\"\"\"\"\"\n", "\n", "class Animate(object):\n", " fps = 20\n", " codec = 'libx264'\n", " extra_args = ['-pix_fmt', 'yuv420p']\n", " keep = False\n", " blit = False\n", "\n", " def __init__(self, keep=False, show=0):\n", " \"\"\"\n", " Parameters\n", " ----------\n", " keep : bool\n", " If `True`, then keep the generated movie file.\n", " show : int\n", " If > 0, then display every `save`th frame in notebook. (Slow!)\n", " \"\"\"\n", " # call the animator. blit=True means only re-draw the parts that have changed.\n", " self._frame = 0\n", " self.keep = keep\n", " self.show = show \n", "\n", " @property\n", " def fig(self):\n", " if not hasattr(self, '_fig'):\n", " self._fig = plt.gcf()\n", " return self._fig\n", " \n", " def _generate_anim(self):\n", " _frames = self.frames()\n", "\n", " # Call this here to allow frames() to set self._fig\n", " _initial_frame = _frames.next()\n", "\n", " def _init():\n", " return _initial_frame\n", "\n", " def _f():\n", " while True:\n", " yield _frames.next()\n", " \n", " self.anim = matplotlib.animation.FuncAnimation(\n", " fig=self.fig, \n", " func=self._animate,\n", " init_func=_init,\n", " frames=_f,\n", " blit=self.blit)\n", " \n", " def _encode_video(self):\n", " with NamedTemporaryFile(suffix='.mp4', delete=not self.keep) as f:\n", " if self.keep:\n", " print(\"Saving movie to file '%s'\" % (f.name,))\n", " self.anim.save(f.name, \n", " fps=self.fps, \n", " codec=self.codec,\n", " extra_args=self.extra_args)\n", " video = open(f.name, \"rb\").read()\n", " self.anim._encoded_video = video.encode(\"base64\")\n", " \n", " def _animate(self, res):\n", " if 0 < self.show and 0 == self._frame % self.show:\n", " IPython.display.clear_output()\n", " display(self.fig)\n", " self._frame += 1\n", " return res\n", " \n", " def _repr_html_(self):\n", " if not hasattr(self, 'anim'):\n", " self._generate_anim()\n", " if not hasattr(self.anim, '_encoded_video'):\n", " self._encode_video()\n", "\n", " res = VIDEO_TAG.format(self.anim._encoded_video)\n", " plt.close(self._fig)\n", " if 0 < self.show:\n", " IPython.display.clear_output()\n", " return res" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are some tools to convert the density profile (in cylindrical coordinates) into the line-of-sight and integrated averages reported in the experiment. For the integrated density, one simply does $\\int_0^\\infty 2\\pi r \\d{r} \\rho(x, r)$. For the line-of-sight density we effect the change of variables $r = z\\cosh(\\alpha)$ to obtain\n", "\n", "$$\n", " \\int_{-\\infty}^\\infty \\d{y}\\; \\rho(x, \\sqrt{y^2+z^2}) = \n", " \\int_{-\\infty}^\\infty \\d{\\alpha}\\; 2 r \\rho(x, r) = \n", " \\int_{-\\infty}^\\infty \\d{\\alpha}\\; 2z\\cosh\\alpha \\rho(x, z\\cosh\\alpha).\n", "$$\n", "\n", "The problem here is that we have the integrand tabulated at a set of discrete $r$. Since the integrand is an even function of $\\alpha$, we instead manipulate the arrays so that we obtain an even extension of the integrand as a function of $\\alpha$ and then apply Simpsons rule." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def get_integrated_densities(x, r, psi, Nz=200):\n", " \"\"\"Return the integrated 1D and 2D densities.\"\"\"\n", " r = r.ravel()\n", " z = np.linspace(-r.max(), r.max(), Nz)\n", " rho = 2.0*abs(psi)**2\n", " integrand = np.zeros(rho.shape + (Nz,))\n", " rho_2d = np.zeros((rho.shape[0], Nz))\n", " Nr = len(r)\n", " for _nz, _z in enumerate(z):\n", " _nr = np.where(r>=abs(_z))[0][0]\n", " integrand = r[_nr:][None,:]*rho[:,_nr:]\n", " alpha = np.arccosh(r[_nr:]/abs(_z))\n", " integrand = np.hstack([integrand[:,::-1], integrand])\n", " alpha = np.hstack([-alpha[::-1], alpha])\n", " rho_2d[:,_nz] = sp.integrate.simps(integrand, alpha[None,:], axis=1)\n", " r = r[None, :]\n", " rho_1d = np.trapz(2*np.pi*r * rho, r, axis=1)\n", " return ((x.ravel(), rho_1d),\n", " ((x, z.ravel()[None, :]), rho_2d))\n", "\n", "def simulate_imaging(x, y, rho_2d, dx=3.0, signal_to_noise=0.03, seed=1):\n", " \"\"\"Simulate the MIT imaging proceedure with 3 micron resolution\n", " and noise with a 3% standard deviation.\"\"\"\n", "\n", " nx0 = dx // np.diff(x.ravel()).mean()\n", " ny0 = dx // np.diff(y.ravel()).mean()\n", " nx1 = x.size//nx0\n", " ny1 = y.size//ny0\n", " Nx = nx0*nx1\n", " Ny = ny0*ny1\n", " print nx0, nx1, ny0, ny1\n", " rho_2d = rho_2d[:Nx, :Ny].reshape((nx1, nx0, ny1, ny0)).mean(axis=3).mean(axis=1)\n", " x = x.ravel()[:Nx].reshape((nx1, nx0)).mean(axis=1)\n", " y = y.ravel()[:Nx].reshape((ny1, ny0)).mean(axis=1)\n", "\n", " np.random.seed(seed)\n", " rho_2d += np.random.normal(0.0, rho_2d.max()*signal_to_noise, rho_2d.shape)\n", " return x, y, rho_2d" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a little test, we compute the following Gaussian which can be integrated analytically." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import sympy.interactive\n", "from sympy import S\n", "sympy.interactive.init_printing(use_latex=True)\n", "rho = S('exp(-(r**2 + x**2)/2)')\n", "_s = dict(r=S('sqrt(y**2+z**2)'))\n", "rho_2d = sympy.integrate(rho.subs(_s), ('z',-S('oo'), S('oo')))\n", "rho_1d = sympy.integrate(rho_2d, ('y',-S('oo'), S('oo')))\n", "display(rho_1d)\n", "display(rho_2d)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "x, r = p.model.basis.xyz\n", "r *= 5/r.max()\n", "x *= 5/x.max()\n", "psi = sqrt(sympy.lambdify(['x', 'r'], rho, numpy)(x=x, r=r)/2)\n", "(_x, _rho_1d), ((_x, _y), _rho_2d) = get_integrated_densities(x, r, psi)\n", "assert np.allclose(_rho_1d, sympy.lambdify(['x'], rho_1d, numpy)(x.ravel()), rtol=0.001)\n", "assert np.allclose(_rho_2d, sympy.lambdify(['x', 'y'], rho_2d, numpy)(x,_y), atol=0.00001, rtol=0.007)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "lam=6.2;k_max_kF=np.pi\n", "_MIT = vortex_expansion._MIT\n", "B_min=580.0\n", "_args = dict(lam=lam, dim=2, basis='DVR', B_min=B_min, minimize=False, fast_load=False)\n", "psi0, p, filename = _MIT.get_realistic(**_args)\n", "self = vortex_expansion.VortexExpansion()\n", "_N = self.get_N(b=p.model.basis, k_max=k_max_kF*p.kF)\n", "psi0, p, filename = _MIT.get_realistic(_N=_N, **_args)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "figure(figsize=(20,5))\n", "plt.subplot(2, len(files), 1)\n", "for _n, _f in enumerate(files):\n", " with h5py.File(os.path.join(DATA, _f)) as f:\n", " psi = np.asarray(f['psis'][-1])\n", " x = np.asarray(f['x'])/f['micron']\n", " r = np.asarray(f['rs'][-1])/f['micron']\n", " subplot(2, len(files), _n+1)\n", " plt.plot(x, (2.0*abs(psi)**2).sum(axis=1))\n", " (_x, _rho_1d), ((_x, _y), _rho_2d) = get_integrated_densities(x, r, psi)\n", " _x, _y, _rho_2d = simulate_imaging(_x, _y, _rho_2d/250)\n", " plt.plot(_x, np.trapz(_rho_2d, _y, axis=1))\n", " \n", " subplot(2, len(files), len(files) + _n + 1)\n", " c = plt.contourf(_x.ravel(), _y.ravel(), _rho_2d.T, 100, cmap='gist_heat')\n", " plt.setp(c.collections, rasterized=True)\n", " title(_f)\n", "savefig('expansion_image.pdf')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "figure(figsize=(5,2))\n", "_f = 'data_R0.20_x0.00_lam6.2.hd5'\n", "_f = 'data_R0.48_x0.00_lam6.2.hd5'\n", "with h5py.File(os.path.join(DATA, _f)) as f:\n", " psi = np.asarray(f['psis'][-1])\n", " x = np.asarray(f['x'])/f['micron']\n", " r = np.asarray(f['rs'][-1])/f['micron']\n", "\n", "#plt.plot(x, (2.0*abs(psi)**2).sum(axis=1))\n", "(_x, _rho_1d), ((x, y), rho_2d) = get_integrated_densities(x, r, psi)\n", "_x, _y, _rho_2d = simulate_imaging(x, y, rho_2d/250)\n", "plt.plot(_x, np.trapz(_rho_2d, _y, axis=1))\n", "plt.xlabel(r'$x\\quad [\\mu\\mathrm{m}]$')\n", "plt.savefig('expansion_image_1D_noise.pdf')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "figure(figsize=(5,2))\n", "c = plt.contourf(_x.ravel(), _y.ravel(), _rho_2d.T, 100, cmap='gist_heat')\n", "plt.setp(c.collections, rasterized=True)\n", "savefig('expansion_image_2D_noise.pdf')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "figure(figsize=(5,2))\n", "c = plt.contourf(x.ravel(), y.ravel(), rho_2d.T, 100, cmap='gist_heat')\n", "plt.setp(c.collections, rasterized=True)\n", "savefig('expansion_image_2D.pdf')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "figure(figsize=(5,2))\n", "plt.plot(x, np.trapz(rho_2d, y, axis=1))\n", "plt.xlabel(r'$x\\quad [\\mu\\mathrm{m}]$')\n", "plt.savefig('expansion_image_1D.pdf')" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "files\n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "import time\n", "import mmf.utils.mmf_plot as mmfplt\n", "class Expansion(object):\n", " def __init__(self, filename=None):\n", " if filename is None:\n", " filename = files[0]\n", " self.img = None\n", " f = self.f = h5py.File(os.path.join(DATA, filename))\n", " self.ts = np.asarray(f['ts'])/p.ms\n", " self.x = np.asarray(f['x'])/f['micron']\n", " self.r = np.asarray(f['rs'][-1])/f['micron']\n", " self.ts = np.asarray(f['ts'])\n", " self.trs = np.asarray(f['trs'])\n", " self.axis = (self.x.min(), self.x.max(), -self.r.max(), self.r.max())\n", " \n", " def __del__(self):\n", " self.f.close()\n", " self.f = None\n", " \n", " def get_frame(self, n):\n", " tic = time.time()\n", " f = self.f\n", " t = self.ts[n]\n", " if t in self.trs and t == self.ts[-1]:\n", " # These are the duplicate expansion points\n", " return None\n", " r = f['rs'][np.where(t <= f['trs'][:])[0][0]]/f['micron']\n", " Nr = (r.max()/r.min() - 1)//2\n", " _r = (np.arange(Nr) + 0.5)/(Nr - 0.5) * r.max()\n", " rho = 2.0*abs(np.asarray(f['psis'][n]))**2\n", " _rho = sp.interpolate.interp1d(r.ravel(), rho)(_r)\n", " _r = np.hstack((-_r[::-1], _r))\n", " _rho = np.hstack((np.fliplr(_rho), _rho))\n", " print(\"%gs to load and interpolate\" % (time.time() - tic,))\n", " tic = time.time()\n", " if self.img is None:\n", " #self.img = self.ax.imshow(2*abs(psi).T, cmap='gist_heat', aspect=1)\n", " #self.img = self.ax.imshow(_rho.T, cmap='gist_heat', aspect=1)\n", " self.img = mmfplt.imcontourf(self.x, _r, _rho, cmap='gist_heat', aspect=1)\n", " else:\n", " #self.img.set_data(2*abs(psi).T)\n", " #self.img.set_data(_rho.T)\n", " self.img = mmfplt.imcontourf(self.x, _r, _rho, cmap='gist_heat', aspect=1)\n", " #ax.set_title(\"t=%gms\" % (_t,))\n", " plt.axis(self.axis)\n", " #plt.xlabel(r'$x [\\mu\\mathrm{m}]$')\n", " #plt.ylabel(r'$y [\\mu\\mathrm{m}]$')\n", " #plt.title(r\"$t=%6.2fms\" % (t/p.ms,))\n", " plt.xlabel('x [micron]')\n", " plt.ylabel('y [micron]')\n", " plt.title(r\"t=%6.2fms\" % (t/p.ms,))\n", " plt.gca().get_frame().set_color('k')\n", " print(\"%gs to draw\" % (time.time() - tic,))\n", " return self.img\n", " \n", "class ExpansionMovie(Animate):\n", " def frames(self):\n", " filename = files[0]\n", " e = Expansion(filename=filename)\n", " for _n, _t in enumerate(e.f['ts']):\n", " print _n\n", " img = e.get_frame(_n)\n", " if img:\n", " yield img" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "from IPython.display import clear_output\n", "filename = 'data_R0.20_x0.00_lam6.2.hd5'\n", "filename = 'data_R0.48_x0.00_lam6.2.hd5'\n", "e = Expansion(filename)\n", "fig = plt.figure(1, figsize=(20, int(20*e.r.max()/e.x.max())))\n", "for n, t in enumerate(e.ts):\n", " print n, t\n", " _filename = os.path.join(PREFIX, 'movie', filename[:-4], 'img_%04i.png' % (n,))\n", " if not os.path.exists(os.path.dirname(_filename)):\n", " os.makedirs(os.path.dirname(_filename))\n", " if not os.path.exists(_filename):\n", " fig.clf()\n", " img = e.get_frame(n)\n", " if not img:\n", " continue\n", " if n % 20 == 0:\n", " clear_output()\n", " display(fig)\n", " tic = time.time()\n", " print(\"%gs to display\" % (time.time() - tic,))\n", " tic = time.time()\n", " fig.savefig(_filename)\n", " print(\"%gs to save\" % (time.time() - tic,))" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "raw", "metadata": {}, "source": [ "ffmpeg -r 80 -i img_%04d.png -vcodec libx264 -pix_fmt yuv420p expansion_movie.mp4" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%connect_info" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "figure(figsize=(20,5))\n", "plt.subplot(1, len(files), 1)\n", "for _n, _f in enumerate(files):\n", " with h5py.File(os.path.join(DATA, _f)) as f:\n", " psi = np.asarray(f['psis'][-1])\n", " x = np.asarray(f['x'])/f['micron']\n", " r = np.asarray(f['rs'][-1])/f['micron']\n", " subplot(1, len(files), _n+1)\n", " plt.plot(x, (2.0*abs(psi)**2).sum(axis=1))\n", " title(_f)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![caption](files/_images/MIT_S1.pdf)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "_MIT = vortex_expansion._MIT\n", "psi0, p, filename = _MIT.get_realistic(\n", " dim=2, basis='DVR', minimize=False, fast_load=False, B_min=702)\n" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "p.t_B_V_image" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Figures for Paper" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now extract the frames used in our paper and exporting these so they can be archived and used elsewhere." ] }, { "cell_type": "code", "collapsed": false, "input": [ "filename = 'data_R0.20_x0.00_lam6.2.hd5'\n", "filename = 'data_R0.48_x0.00_lam6.2.hd5'\n", "outname = 'figure_data.hd5'\n", "if os.path.exists(outname):\n", " os.remove(outname)\n", "with h5py.File(os.path.join(DATA, filename)) as f, h5py.File(outname) as o:\n", " psi0 = np.asarray(f['psis'][1])\n", " psi = np.asarray(f['psis'][-1])\n", " r0 = np.asarray(f['rs'][0])/p.micron\n", " r = np.asarray(f['rs'][-1])/p.micron\n", " x = np.asarray(f['x'])/p.micron\n", " (_x, rho_1d), ((_x, y), rho_2d) = get_integrated_densities(x, r, psi)\n", " (_x, rho_1d0), ((_x, y0), rho_2d0) = get_integrated_densities(x, r0, psi0)\n", " x_noise, y_noise, rho_2d_noise = simulate_imaging(x, y, rho_2d)\n", " o['psi0'] = psi0\n", " o['psi'] = psi\n", " o['rho_1d'] = rho_1d\n", " o['rho_2d'] = rho_2d\n", " o['rho_1d0'] = rho_1d0\n", " o['rho_2d0'] = rho_2d0\n", " o['x'] = x\n", " o['y'] = y\n", " o['y0'] = y0\n", " o['r'] = r\n", " o['rho_2d_noise'] = rho_2d_noise\n", " o['x_noise'] = x_noise\n", " o['y_noise'] = y_noise " ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }