2022-08-30_caustique.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# caustiques (2020) / iridiscence (2022)\n",
    "\n",
    "* This repository *https://github.com/NaturalPatterns/2020_caustiques* presents 2 succesive projects\n",
    "  1. One where I developped a modelling of caustics produced by shallow water and which are detailed in https://laurentperrinet.github.io/sciblog/posts/2020-06-19-caustic-optics.html\n",
    "  2. Another, where I build upon that work to decompose the spectrum of light into its different wavelength, producing *iridescent caustics*.\n",
    " * Rendering of the [full notebook](https://naturalpatterns.github.io/2020_caustiques/) in which I explore the role of different parameters in prouding the images.\n",
    "\n",
    "## caustiques (2020)\n",
    "\n",
    "[![DOI](https://zenodo.org/badge/273226625.svg)](https://zenodo.org/badge/latestdoi/273226625)\n",
    "\n",
    "Caustics ([wikipedia](https://en.wikipedia.org/wiki/Caustic_(optics))) are luminous patterns which are resulting from the superposition of smoothly deviated light rays. It is the heart-shaped pattern in your cup of coffee which is formed as the rays of from the sun are reflected on the cup's surface. It is also the wiggly patterns of light that you will see on the floor of a pool as the sun's light is *refracted* at the surface of the water. Here we will simulate this particular physical phenomenon. Simply because they are mesmerizingly beautiful, but also as it is of interest in visual neuroscience. Indeed, it speaks to how images are formed (more on this later), hence how the brain may understand images.\n",
    "\n",
    "In [this post](https://laurentperrinet.github.io/sciblog/posts/2020-06-19-caustic-optics.html), I have developed a simple formalism to generate such patterns, with the paradoxical result that it is *very* simple to code yet generates patterns with great complexity, such as:\n",
    "\n",
    "<BR>\n",
    "<center>\n",
    "<img src=\"caustique.gif\" width=\"100%\"/>\n",
    "</center>\n",
    "<BR>\n",
    "\n",
    "This is joint work with artist [Etienne Rey](https://laurentperrinet.github.io/authors/etienne-rey/), in which I especially follow the ideas put forward in the series [Turbulence](http://ondesparalleles.org/projets/turbulences/).\n",
    "\n",
    "## iridiscence (2022)\n",
    "\n",
    "Upon further observation, one may discover that caustics exhibit some [iridescence](https://en.wikipedia.org/wiki/Iridescence)$^\\ddagger$, that is, that the light pattern which forms the waggling lines of the caustics may decompose into different colors, forming evanescent rainbows. Here, we will simply use a modulation of the [Snell-Descartes law](https://en.wikipedia.org/wiki/Snell's_law) that we used to compute different angle of refraction. This will be put in relation with the dependance of the refraction index with the wavelength of light and the transformation of a [monochromatic light into RGB](https://github.com/laurentperrinet/lambda2color) that we used in a previous post about [colors of the sky](https://laurentperrinet.github.io/sciblog/posts/2020-07-04-colors-of-the-sky.html). The results are close to subjective observations, with the surprising (to me) observation that colors appear more *between* nodes...\n",
    "  \n",
    "<BR>\n",
    "<center>\n",
    "<video src=\"iridiscence.mp4\" width=\"100%\"/>\n",
    "</center>\n",
    "<BR>\n",
    "\n",
    "Note: $\\ddagger$ I use the term iridescence which is improper in the physical sense as it rather concerns the property of an object to exhibit different colors depending on the angle of view. However, in the perspective of the work with [Etienne Rey](https://laurentperrinet.github.io/authors/etienne-rey/) it resonates with our endeavour to show that percepetion, in particular visual perception, is an active process of the observer within its environment.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# source code\n",
    "\n",
    "## defining the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "figpath = '2022-08-30_caustique'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "PRECISION = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "PRECISION = 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import datetime\n",
    "import numpy as np\n",
    "\n",
    "def init(args=[], figpath=figpath, PRECISION=PRECISION):\n",
    "\n",
    "    import argparse\n",
    "    if figpath is None:\n",
    "        import datetime\n",
    "        date = datetime.datetime.now().date().isoformat()\n",
    "        figpath = f'{date}_caustique'\n",
    "\n",
    "    parser = argparse.ArgumentParser()\n",
    "    parser.add_argument(\"--tag\", type=str, default='caustique', help=\"Tag\")\n",
    "    parser.add_argument(\"--figpath\", type=str, default=figpath, help=\"Folder to store images\")\n",
    "    parser.add_argument(\"--vext\", type=str, default='mp4', help=\"Folder to store images\")\n",
    "    parser.add_argument(\"--nx\", type=int, default=8*2**PRECISION, help=\"number of pixels (vertical)\")\n",
    "    parser.add_argument(\"--ny\", type=int, default=8*2**PRECISION, help=\"number of pixels (horizontal)\")\n",
    "    parser.add_argument(\"--nframe\", type=int, default=5*2**PRECISION, help=\"number of frames\")\n",
    "    parser.add_argument(\"--bin_dens\", type=int, default=1, help=\"relative bin density\")\n",
    "    parser.add_argument(\"--bin_spectrum\", type=int, default=12, help=\"bin spacing in spectrum\")\n",
    "    parser.add_argument(\"--seed\", type=int, default=42, help=\"seed for RNG\")\n",
    "    parser.add_argument(\"--H\", type=float, default=6., help=\"depth of the pool\")\n",
    "    parser.add_argument(\"--variation\", type=float, default=.4, help=\"variation of diffraction index: http://www.philiplaven.com/p20.html 1.40 at 400 nm and 1.37 at 700nm makes a 2% variation\")\n",
    "    parser.add_argument(\"--scale\", type=float, default=150, help=\"sf\")\n",
    "    parser.add_argument(\"--B_sf\", type=float, default=0.25, help=\"bandwidth in sf\")\n",
    "    parser.add_argument(\"--V_Y\", type=float, default=0.3, help=\"horizontal speed\")\n",
    "    parser.add_argument(\"--V_X\", type=float, default=0.3, help=\"vertical speed\")\n",
    "    parser.add_argument(\"--B_V\", type=float, default=2.0, help=\"bandwidth in speed\")\n",
    "    parser.add_argument(\"--theta\", type=float, default=2*np.pi*(2-1.61803), help=\"angle with the horizontal\")\n",
    "    parser.add_argument(\"--B_theta\", type=float, default=np.pi/3, help=\"bandwidth in theta\")\n",
    "    parser.add_argument(\"--min_lum\", type=float, default=.2, help=\"diffusion level for the rendering\")\n",
    "    parser.add_argument(\"--gamma\", type=float, default=2.8, help=\"Gamma exponant to convert luminosity to luminance\")\n",
    "    parser.add_argument(\"--fps\", type=float, default=18, help=\"frames per second\")\n",
    "    parser.add_argument(\"--multispectral\", type=bool, default=True, help=\"Compute caustics on the full spectrogram.\")\n",
    "    parser.add_argument(\"--cache\", type=bool, default=True, help=\"Cache intermediate output.\")\n",
    "    parser.add_argument(\"--verbose\", type=bool, default=False, help=\"Displays more verbose output.\")\n",
    "\n",
    "    opt = parser.parse_args(args=args)\n",
    "\n",
    "    if opt.verbose:\n",
    "        print(opt)\n",
    "    return opt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving our simulations in=2022-08-30_caustique\n"
     ]
    }
   ],
   "source": [
    "print(f'Saving our simulations in={figpath}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = init()\n",
    "opt.figpath = figpath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Namespace(tag='caustique', figpath='2022-08-30_caustique', vext='mp4', nx=1024, ny=1024, nframe=640, bin_dens=1, bin_spectrum=12, seed=42, H=6.0, variation=0.4, scale=150, B_sf=0.25, V_Y=0.3, V_X=0.3, B_V=2.0, theta=2.399988291783386, B_theta=1.0471975511965976, min_lum=0.2, gamma=2.8, fps=18, multispectral=True, cache=True, verbose=False)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## utilities\n",
    "\n",
    "Transfoming a sequence of PNG frames into gif or mp4:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_gif(gifname, fnames, fps, do_delete=True):\n",
    "    import imageio\n",
    "\n",
    "    with imageio.get_writer(gifname, mode='I', fps=fps) as writer:\n",
    "        for fname in fnames:\n",
    "            writer.append_data(imageio.imread(fname))\n",
    "\n",
    "    from pygifsicle import optimize\n",
    "    optimize(str(gifname))\n",
    "    if do_delete: \n",
    "        for fname in fnames: os.remove(fname)\n",
    "    return gifname\n",
    "\n",
    "# https://moviepy.readthedocs.io/en/latest/getting_started/videoclips.html#imagesequenceclip\n",
    "def make_mp4(mp4name, fnames, fps, do_delete=True):\n",
    "    import moviepy.editor as mpy\n",
    "    clip = mpy.ImageSequenceClip(fnames, fps=fps)\n",
    "    clip.write_videofile(mp4name, fps=fps, codec='libx264', verbose=False, logger=None)\n",
    "    if do_delete: \n",
    "        for fname in fnames: os.remove(fname)\n",
    "    return mp4name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Utilities to compute the spectrum of the blue sky and convert them later to RGB values (check out https://laurentperrinet.github.io/sciblog/posts/2020-07-04-colors-of-the-sky.html for details)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lambda2color import Lambda2color, xyz_from_xy\n",
    "\n",
    "# borrowed from https://github.com/gummiks/gummiks.github.io/blob/master/scripts/astro/planck.py\n",
    "def planck(wav, T):\n",
    "    import scipy.constants as const\n",
    "    c = const.c # c = 3.0e+8\n",
    "    h = const.h # h = 6.626e-34\n",
    "    k = const.k # k = 1.38e-23\n",
    "    a = 2.0*h*c**2\n",
    "    b = h*c/(wav*k*T)\n",
    "    intensity = a / ( (wav**5) * (np.exp(b) - 1.0) )\n",
    "    return intensity\n",
    "\n",
    "def scattering(wav, a=0.005, p=1.3, b=0.45):\n",
    "    \"\"\"\n",
    "    b is  proportionate  to  the  column  density  of  aerosols\n",
    "    along  the  path  of  sunlight,  from  outside  the  atmosphere\n",
    "    to  the  point  of  observation\n",
    "    \n",
    "    see https://laurentperrinet.github.io/sciblog/posts/2020-07-04-colors-of-the-sky.html for more details\n",
    "\n",
    "    \"\"\"\n",
    "    # converting wav in µm:\n",
    "    intensity = np.exp(-a/((wav/1e-6)**4)) # Rayleigh extinction by nitrogen\n",
    "    intensity *= (wav/1e-6)**-4\n",
    "    intensity *= np.exp(-b/((wav/1e-6)**p)) # Aerosols\n",
    "    return intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## computing the caustics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import MotionClouds as mc\n",
    "class Caustique:\n",
    "    def __init__(self, opt):\n",
    "        \"\"\"\n",
    "        Image coordinates follow 'ij' indexing, that is,\n",
    "        * their origin at the top left,\n",
    "        * the X axis is vertical and goes \"down\",\n",
    "        * the Y axis is horizontal and goes \"right\".\n",
    "\n",
    "        \"\"\"\n",
    "        self.mc = mc\n",
    "        self.ratio = opt.ny/opt.nx # ratio between height and width (>1 for portrait, <1 for landscape)\n",
    "        X = np.linspace(0, 1, opt.nx, endpoint=False) # vertical\n",
    "        Y = np.linspace(0, self.ratio, opt.ny, endpoint=False) # horizontal\n",
    "        self.xv, self.yv = np.meshgrid(X, Y, indexing='ij')\n",
    "        self.opt = opt\n",
    "        # https://stackoverflow.com/questions/16878315/what-is-the-right-way-to-treat-python-argparse-namespace-as-a-dictionary\n",
    "        self.d = vars(opt)\n",
    "        print(self.opt.figpath)\n",
    "        os.makedirs(self.opt.figpath, exist_ok=True)\n",
    "        self.cachepath = os.path.join('/tmp', self.opt.figpath)\n",
    "        if opt.verbose: print(f'{self.cachepath=}')\n",
    "        os.makedirs(self.cachepath, exist_ok=True)\n",
    "\n",
    "        # a standard white:\n",
    "        illuminant_D65 = xyz_from_xy(0.3127, 0.3291), \n",
    "        illuminant_sun = xyz_from_xy(0.325998, 0.335354)\n",
    "        # color conversion class\n",
    "        self.cs_srgb = Lambda2color(red=xyz_from_xy(0.64, 0.33),\n",
    "                               green=xyz_from_xy(0.30, 0.60),\n",
    "                               blue=xyz_from_xy(0.15, 0.06),\n",
    "                               white=illuminant_sun)\n",
    "    def wave(self):\n",
    "        filename = f'{self.cachepath}/{self.opt.tag}_wave.npy'\n",
    "        if os.path.isfile(filename):\n",
    "            z = np.load(filename)\n",
    "        else:\n",
    "            # A simplistic model of a wave using https://github.com/NeuralEnsemble/MotionClouds\n",
    "            fx, fy, ft = mc.get_grids(self.opt.nx, self.opt.ny, self.opt.nframe)\n",
    "            env = mc.envelope_gabor(fx, fy, ft, V_X=self.opt.V_Y, V_Y=self.opt.V_X, B_V=self.opt.B_V,\n",
    "                                    sf_0=1./self.opt.scale, B_sf=self.opt.B_sf/self.opt.scale,\n",
    "                                    theta=self.opt.theta, B_theta=self.opt.B_theta)\n",
    "            z = mc.rectif(mc.random_cloud(env, seed=self.opt.seed))\n",
    "            if self.opt.cache: np.save(filename, z)\n",
    "        return z\n",
    "\n",
    "    def transform(self, z_, modulation=1.):\n",
    "        xv, yv = self.xv.copy(), self.yv.copy()\n",
    "\n",
    "        dzdx = z_ - np.roll(z_, 1, axis=0)\n",
    "        dzdy = z_ - np.roll(z_, 1, axis=1)\n",
    "        xv = xv + modulation * self.opt.H * dzdx\n",
    "        yv = yv + modulation * self.opt.H * dzdy\n",
    "\n",
    "        xv = np.mod(xv, 1)\n",
    "        yv = np.mod(yv, self.ratio)\n",
    "\n",
    "        return xv, yv\n",
    "\n",
    "    def plot(self, z, do_color=True, videoname=None, dpi=50):\n",
    "        #hist = self.do_raytracing(z)\n",
    "        binsx, binsy = self.opt.nx//self.opt.bin_dens, self.opt.ny//self.opt.bin_dens\n",
    "\n",
    "        if videoname is None:\n",
    "            videoname=f'{self.opt.figpath}/{self.opt.tag}.{self.opt.vext}'\n",
    "\n",
    "        subplotpars = matplotlib.figure.SubplotParams(left=0., right=1., bottom=0., top=1., wspace=0., hspace=0.,)\n",
    "\n",
    "        if self.opt.multispectral:\n",
    "            # multiply by the spectrum of the sky\n",
    "            #wavelengths = self.cs_srgb.cmf[::self.opt.bin_spectrum, 0]*1e-9\n",
    "            wavelengths = self.cs_srgb.cmf[:, 0]*1e-9\n",
    "            N_wavelengths = len(wavelengths)\n",
    "            intensity5800 = planck(wavelengths, 5800.)\n",
    "            scatter = scattering(wavelengths)\n",
    "            spectrum_sky = intensity5800 * scatter\n",
    "\n",
    "            #image_rgb = self.cs_srgb.spec_to_rgb(hist)\n",
    "            image_rgb = np.zeros((self.opt.nx//self.opt.bin_dens,  self.opt.ny//self.opt.bin_dens, 3, self.opt.nframe))\n",
    "            for i_frame in range(self.opt.nframe):\n",
    "                for ii_wavelength, i_wavelength in enumerate(range(self.opt.bin_spectrum//2, N_wavelengths, self.opt.bin_spectrum)):\n",
    "                    modulation = 1. + self.opt.variation/2 - self.opt.variation*i_wavelength/N_wavelengths\n",
    "                    # print(i_wavelength, N_wavelengths, modulation)\n",
    "                    xv, yv = self.transform(z[:, :, i_frame], modulation=modulation)\n",
    "                    hist_, edge_x, edge_y = np.histogram2d(xv.ravel(), yv.ravel(),\n",
    "                                                           bins=[binsx, binsy],\n",
    "                                                           range=[[0, 1], [0, self.ratio]],\n",
    "                                                           density=True)\n",
    "\n",
    "                    spec = np.zeros((N_wavelengths))\n",
    "                    spec[i_wavelength] = spectrum_sky[i_wavelength]\n",
    "                    rgb = self.cs_srgb.spec_to_rgb(spec)\n",
    "\n",
    "                    image_rgb[:, :, :, i_frame] += hist_[:, :, None] * rgb[None, None, :]\n",
    "\n",
    "            # image_rgb -= image_rgb.min()\n",
    "            image_rgb /= image_rgb.max()\n",
    "        else:\n",
    "            hist = np.zeros((binsx, binsy, self.opt.nframe))\n",
    "            for i_frame in range(self.opt.nframe):\n",
    "                xv, yv = self.transform(z[:, :, i_frame])\n",
    "                hist_, edge_x, edge_y = np.histogram2d(xv.ravel(), yv.ravel(),\n",
    "                                                       bins=[binsx, binsy],\n",
    "                                                       range=[[0, 1], [0, self.ratio]],\n",
    "                                                       density=True)\n",
    "            #hist /= hist.max()\n",
    "\n",
    "        fnames = []\n",
    "        for i_frame in range(self.opt.nframe):\n",
    "            fig, ax = plt.subplots(figsize=(self.opt.ny/self.opt.bin_dens/dpi, self.opt.nx/self.opt.bin_dens/dpi), subplotpars=subplotpars)\n",
    "            if self.opt.multispectral:\n",
    "                ax.imshow(image_rgb[:, :, :, i_frame] ** (1/self.opt.gamma), vmin=0, vmax=1)\n",
    "            else:\n",
    "                if do_color:\n",
    "                    bluesky = np.array([0.268375, 0.283377]) # xyz\n",
    "                    sun = np.array([0.325998, 0.335354]) # xyz\n",
    "                    # ax.pcolormesh(edge_y, edge_x, hist[:, :, i_frame], vmin=0, vmax=1, cmap=plt.cm.Blues_r)\n",
    "                    # https://en.wikipedia.org/wiki/CIE_1931_color_space#Mixing_colors_specified_with_the_CIE_xy_chromaticity_diagram\n",
    "                    L1 = 1 - hist[:, :, i_frame]\n",
    "                    L2 = hist[:, :, i_frame]\n",
    "                    image_denom = L1 / bluesky[1] + L2 / sun[1]\n",
    "                    image_x = (L1 * bluesky[0] / bluesky[1] + L2 * sun[0] / sun[1]) / image_denom\n",
    "                    image_y = (L1 + L2) / image_denom \n",
    "                    image_xyz = np.dstack((image_x, image_y, 1 - image_x - image_y))\n",
    "                    image_rgb = self.cs_srgb.xyz_to_rgb(image_xyz)\n",
    "                    image_L = self.opt.min_lum + (1-self.opt.min_lum)* L2 ** .61803\n",
    "                    ax.imshow(image_L[:, :, None]*image_rgb, vmin=0, vmax=1)\n",
    "\n",
    "                else:\n",
    "                    ax.imshow(1-image_L, vmin=0, vmax=1)\n",
    "\n",
    "            fname = f'{self.cachepath}/{self.opt.tag}_frame_{i_frame:04d}.png'\n",
    "            fig.savefig(fname, dpi=dpi)\n",
    "            fnames.append(fname)\n",
    "            plt.close('all')\n",
    "        if self.opt.vext == 'gif':\n",
    "            return make_gif(videoname, fnames, fps=self.opt.fps)\n",
    "        else:\n",
    "            return make_mp4(videoname, fnames, fps=self.opt.fps)\n",
    "\n",
    "    def show(self, videoname, width=1024):\n",
    "        if self.opt.vext == 'gif':\n",
    "            from IPython.display import Image, display\n",
    "            return display(Image(url=videoname, width=width))\n",
    "        else:\n",
    "            #import moviepy.editor as mpy\n",
    "            #return mpy.ipython_display(videoname, width=width)\n",
    "            # https://github.com/NeuralEnsemble/MotionClouds/blob/master/MotionClouds/MotionClouds.py#L858\n",
    "            from IPython.display import HTML, display\n",
    "            opts = ' loop=\"1\" autoplay=\"1\" controls '\n",
    "            html = HTML(f'<video {opts} src=\"{videoname}\" type=\"video/{self.opt.vext}\" width={width}\\>')\n",
    "            html.reload()\n",
    "            return display(html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# a simple caustics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## a wave pattern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_wave.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "c = Caustique(opt)\n",
    "wave_videoname = f'{opt.figpath}/{opt.tag}_wave'\n",
    "if not os.path.isfile(f'{wave_videoname}{c.mc.vext}'):\n",
    "    print(f'Doing {wave_videoname}{c.mc.vext}')\n",
    "    z = c.wave()\n",
    "    c.mc.anim_save(z.swapaxes(0, 1), wave_videoname)\n",
    "c.show(f'{wave_videoname}{c.mc.vext}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## generating the caustics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1024, 1024, 640)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = Caustique(opt)\n",
    "z = c.wave()\n",
    "z.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "videoname = f'{opt.figpath}/{opt.tag}.{opt.vext}'\n",
    "if not os.path.isfile(videoname):\n",
    "    c = Caustique(opt)\n",
    "    z = c.wave()\n",
    "    videoname = c.plot(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "c.show(videoname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# exploring parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_scan = 9\n",
    "base = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## water depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "2022-08-30_caustique\n",
      "H = 3.000\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_3.000.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "H = 3.568\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_3.568.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "H = 4.243\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_4.243.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "H = 5.045\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_5.045.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "H = 6.000\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_6.000.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "H = 7.135\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_7.135.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "H = 8.485\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_8.485.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "H = 10.091\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_10.091.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "H = 12.000\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_H_12.000.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "opt = init()\n",
    "opt.figpath = figpath\n",
    "\n",
    "c = Caustique(opt)\n",
    "# compute just once\n",
    "z = c.wave()\n",
    "\n",
    "for H_ in c.opt.H*np.logspace(-1, 1, N_scan, base=base):\n",
    "    opt = init()\n",
    "    opt.figpath = figpath\n",
    "    c = Caustique(opt)\n",
    "\n",
    "    print(f'H = {H_:.3f}')\n",
    "    c.opt.H = H_\n",
    "    c.opt.tag = f'{opt.tag}_H_{H_:.3f}'\n",
    "    videoname = f'{opt.figpath}/{c.opt.tag}.{opt.vext}'\n",
    "    if not os.path.isfile(videoname):\n",
    "        url=c.plot(z, videoname=videoname)\n",
    "    c.show(videoname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## refraction index variation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "2022-08-30_caustique\n",
      "variation = 0.010\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.010.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "variation = 0.017\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.017.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "variation = 0.028\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.028.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "variation = 0.046\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.046.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "variation = 0.077\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.077.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "variation = 0.129\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.129.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "variation = 0.215\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.215.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "variation = 0.359\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.359.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2022-08-30_caustique\n",
      "variation = 0.599\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<video  loop=\"1\" autoplay=\"1\" controls  src=\"2022-08-30_caustique/caustique_variation_0.599.mp4\" type=\"video/mp4\" width=1024\\>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "opt = init()\n",
    "opt.figpath = figpath\n",
    "\n",
    "c = Caustique(opt)\n",
    "z = c.wave()\n",
    "\n",
    "for variation_ in np.logspace(-2, 0, N_scan, base=10, endpoint=False):\n",
    "    opt = init()\n",
    "    opt.figpath = figpath\n",
    "    c = Caustique(opt)\n",
    "    print(f'variation = {variation_:.3f}')\n",
    "    c.opt.variation = variation_\n",
    "    c.opt.tag = f'{opt.tag}_variation_{variation_:.3f}'\n",
    "    videoname = f'{opt.figpath}/{c.opt.tag}.{opt.vext}'\n",
    "    if not os.path.isfile(videoname):\n",
    "        url=c.plot(z, videoname=videoname)\n",
    "    c.show(videoname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## other variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "======scale======\n",
      "2022-08-30_caustique\n",
      "scale=scale(default)*0.500=7.500E+01\n",
      "Doing  2022-08-30_caustique/caustique_scale_modul_0.500.mp4\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'func' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Input \u001b[0;32mIn [19]\u001b[0m, in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     18\u001b[0m         \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mDoing \u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mmcname\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;132;01m{\u001b[39;00mc\u001b[38;5;241m.\u001b[39mmc\u001b[38;5;241m.\u001b[39mvext\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m     19\u001b[0m         c\u001b[38;5;241m.\u001b[39mmc\u001b[38;5;241m.\u001b[39manim_save(z\u001b[38;5;241m.\u001b[39mswapaxes(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m1\u001b[39m), \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mmcname\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m---> 20\u001b[0m     url\u001b[38;5;241m=\u001b[39m\u001b[43mc\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mplot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mz\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mvideoname\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mvideoname\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     21\u001b[0m c\u001b[38;5;241m.\u001b[39mshow(videoname)\n",
      "Input \u001b[0;32mIn [11]\u001b[0m, in \u001b[0;36mCaustique.plot\u001b[0;34m(self, z, do_color, videoname, dpi)\u001b[0m\n\u001b[1;32m    137\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m make_gif(videoname, fnames, fps\u001b[38;5;241m=\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mopt\u001b[38;5;241m.\u001b[39mfps)\n\u001b[1;32m    138\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 139\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mmake_mp4\u001b[49m\u001b[43m(\u001b[49m\u001b[43mvideoname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfnames\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfps\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopt\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfps\u001b[49m\u001b[43m)\u001b[49m\n",
      "Input \u001b[0;32mIn [9]\u001b[0m, in \u001b[0;36mmake_mp4\u001b[0;34m(mp4name, fnames, fps, do_delete)\u001b[0m\n\u001b[1;32m     16\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmoviepy\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01meditor\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mmpy\u001b[39;00m\n\u001b[1;32m     17\u001b[0m clip \u001b[38;5;241m=\u001b[39m mpy\u001b[38;5;241m.\u001b[39mImageSequenceClip(fnames, fps\u001b[38;5;241m=\u001b[39mfps)\n\u001b[0;32m---> 18\u001b[0m \u001b[43mclip\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwrite_videofile\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmp4name\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfps\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfps\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcodec\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mlibx264\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mverbose\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mFalse\u001b[39;49;00m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlogger\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mNone\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m     19\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m do_delete: \n\u001b[1;32m     20\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m fname \u001b[38;5;129;01min\u001b[39;00m fnames: os\u001b[38;5;241m.\u001b[39mremove(fname)\n",
      "File \u001b[0;32m~/opt/anaconda3/envs/default/lib/python3.9/site-packages/decorator.py:232\u001b[0m, in \u001b[0;36mdecorate.<locals>.fun\u001b[0;34m(*args, **kw)\u001b[0m\n\u001b[1;32m    230\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m kwsyntax:\n\u001b[1;32m    231\u001b[0m     args, kw \u001b[38;5;241m=\u001b[39m fix(args, kw, sig)\n\u001b[0;32m--> 232\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mcaller\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mextras\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkw\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/opt/anaconda3/envs/default/lib/python3.9/site-packages/moviepy/decorators.py:54\u001b[0m, in \u001b[0;36mrequires_duration\u001b[0;34m(f, clip, *a, **k)\u001b[0m\n\u001b[1;32m     52\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mAttribute \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mduration\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m not set\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     53\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m---> 54\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[43mclip\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mk\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/opt/anaconda3/envs/default/lib/python3.9/site-packages/decorator.py:232\u001b[0m, in \u001b[0;36mdecorate.<locals>.fun\u001b[0;34m(*args, **kw)\u001b[0m\n\u001b[1;32m    230\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m kwsyntax:\n\u001b[1;32m    231\u001b[0m     args, kw \u001b[38;5;241m=\u001b[39m fix(args, kw, sig)\n\u001b[0;32m--> 232\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mcaller\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mextras\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43m \u001b[49m\u001b[43margs\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkw\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/opt/anaconda3/envs/default/lib/python3.9/site-packages/moviepy/decorators.py:118\u001b[0m, in \u001b[0;36muse_clip_fps_by_default\u001b[0;34m(f, clip, *a, **k)\u001b[0m\n\u001b[1;32m    111\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m clip\u001b[38;5;241m.\u001b[39mfps\n\u001b[1;32m    112\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mAttributeError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mNo \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mfps\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m (frames per second) attribute specified\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    113\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m for function \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m and the clip has no \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mfps\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m attribute. Either\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    114\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m provide e.g. fps=24 in the arguments of the function, or define\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    115\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m the clip\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124ms fps with `clip.fps=24`\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m%\u001b[39m f\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m)\n\u001b[0;32m--> 118\u001b[0m names \u001b[38;5;241m=\u001b[39m inspect\u001b[38;5;241m.\u001b[39mgetfullargspec(\u001b[43mfunc\u001b[49m)\u001b[38;5;241m.\u001b[39margs[\u001b[38;5;241m1\u001b[39m:]\n\u001b[1;32m    120\u001b[0m new_a \u001b[38;5;241m=\u001b[39m [fun(arg) \u001b[38;5;28;01mif\u001b[39;00m (name\u001b[38;5;241m==\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mfps\u001b[39m\u001b[38;5;124m'\u001b[39m) \u001b[38;5;28;01melse\u001b[39;00m arg\n\u001b[1;32m    121\u001b[0m          \u001b[38;5;28;01mfor\u001b[39;00m (arg, name) \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(a, names)]\n\u001b[1;32m    122\u001b[0m new_kw \u001b[38;5;241m=\u001b[39m {k: fun(v) \u001b[38;5;28;01mif\u001b[39;00m k\u001b[38;5;241m==\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mfps\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;28;01melse\u001b[39;00m v\n\u001b[1;32m    123\u001b[0m          \u001b[38;5;28;01mfor\u001b[39;00m (k,v) \u001b[38;5;129;01min\u001b[39;00m k\u001b[38;5;241m.\u001b[39mitems()}\n",
      "\u001b[0;31mNameError\u001b[0m: name 'func' is not defined"
     ]
    }
   ],
   "source": [
    "for variable in ['scale', 'B_sf', 'B_theta', 'V_X', 'B_V', 'gamma', ]: #  'theta', 'V_Y'\n",
    "    print(f'======{variable}======')\n",
    "    for modul in np.logspace(-1, 1, N_scan, base=base):\n",
    "        opt = init()\n",
    "        opt.figpath = figpath\n",
    "\n",
    "        c = Caustique(opt)\n",
    "        c.d[variable] *= modul\n",
    "        c.opt.tag = f'{opt.tag}_{variable}_modul_{modul:.3f}'\n",
    "        videoname = f'{opt.figpath}/{c.opt.tag}.{opt.vext}'\n",
    "\n",
    "        print(f'{variable}={variable}(default)*{modul:.3f}={c.d[variable]:.3E}')\n",
    "        if not os.path.isfile(videoname):\n",
    "            print('Doing ', videoname)\n",
    "            z = c.wave()\n",
    "            mcname = f'{opt.figpath}/{c.opt.tag}'\n",
    "            if False: #not os.path.isfile(f'{mcname}{c.mc.vext}'): \n",
    "                print('Doing ', f'{mcname}{c.mc.vext}')\n",
    "                c.mc.anim_save(z.swapaxes(0, 1), f'{mcname}')\n",
    "            url=c.plot(z, videoname=videoname)\n",
    "        c.show(videoname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## installation\n",
    "\n",
    "Install [dependencies](https://pip.pypa.io/en/stable/user_guide/#requirements-files), then this notebook:\n",
    "\n",
    "```\n",
    "python3 -m pip install -r requirements.txt\n",
    "```"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "tags": []
   },
   "source": [
    "%pip install --upgrade -r requirements.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  },
  "nteract": {
   "version": "0.28.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}