{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Scientific Computing\n",
    "#### Lecture 06: Matplotlib Review and Animations\n",
    "J.R. Gladden, Spring 2016, Univ. of Mississippi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "As of about a year ago, the following is the preferred method of importing pyplot and numpy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib wx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x = np.linspace(0,20,100)\n",
    "y = np.sin(x)*np.cos(x)\n",
    "plt.plot(x,y,'g-',lw=2,label='This curve')\n",
    "plt.xlabel('X')\n",
    "plt.ylabel('Y')\n",
    "plt.legend()\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "** Reading and writing files in numpy: ** Use 'loadtxt' and 'savetxt' to read and write (respectively) numpy arrays from/to files on the hard drive.  The active directory must contain the files, or you can use a full path to point to the location of the file.  The following is an excercise to illustrate these techniques: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "##First just generate some data and record to a file for practice\n",
    "x=np.linspace(0,10,50)\n",
    "y1=(x-2.)**2*np.sin(x)\n",
    "y2=np.exp( (x-2.) ) * np.sin(x)\n",
    "\n",
    "#These will be ROW vectors, want to record as columns - use transpose\n",
    "data = np.array([x, y1, y2])\n",
    "data=data.transpose()\n",
    "\n",
    "#Save the data to a file with some formatting\n",
    "np.savetxt('mydata.dat',data, fmt='%2.4f', delimiter='\\t\\t')\n",
    "\n",
    "## Now read in and plot\n",
    "\n",
    "plt.close('all')\n",
    "newx,newy1,newy2 = np.loadtxt('mydata.dat',unpack=True)\n",
    "\n",
    "plt.plot(newx,newy1,'g-s',label='Y1 Data')\n",
    "plt.plot(newx,newy2,'k--',lw=2,label='Y2 Data')\n",
    "plt.legend(loc=2)\n",
    "plt.xlabel('X values')\n",
    "plt.ylabel('Y values')\n",
    "plt.savefig('myplot.png')\n",
    "plt.savefig('myplot.svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "** Exercise 2: Flow Profile:** In this bit of code, we compute and plot the radial profile of the flow velocity in a circular pipe for fluids with a range of viscosities (indicated by $\\mu$). Note how we only have to change the values of $\\mu$ in one location in the program, and everything e;se automatically updates!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.close('all')\n",
    "\n",
    "def velocity(r, dP = 0.3, n=1.0,R = 1.0, mu = 0.05):\n",
    "    return (dP/(2*mu))**(1.0/n)*(n/(n+1.0))*( (R**(1+1.0/n) - r**(1+1.0/n)))\n",
    "    \n",
    "dP = 0.3    #driving pressure difference\n",
    "n=1.0      #parameter for fluid type\n",
    "R = 1.0     # Radius of tube\n",
    "\n",
    "mus = [0.1,0.5,1.0]\n",
    "r = np.linspace(0,R,30)\n",
    "\n",
    "for mu in mus:\n",
    "    labelString = 'Viscosity = %2.2f' % mu\n",
    "    v = velocity(r,dP=dP,n=n,R=R,mu=mu)\n",
    "    plt.plot(r,v,'-', label = labelString)\n",
    "\n",
    "plt.xlabel(' Radius (m)')\n",
    "plt.ylabel('Velocity (m/s)')\n",
    "plt.title('Viscous Fluid Flow in a Pipe') \n",
    "plt.legend()\n",
    "\n",
    "# Normalized version in a new figure window\n",
    "\n",
    "plt.figure(2)\n",
    "for mu in mus:\n",
    "    labelString = 'Viscosity = %2.2f' % mu\n",
    "    norm_v = velocity(r,dP=dP,n=n,R=R,mu=mu) / velocity(0.0,dP=dP,n=n,R=R,mu=mu)\n",
    "    plt.plot(r,norm_v,'-', label = labelString)\n",
    "plt.xlabel(' Radius (m)')\n",
    "plt.ylabel('Normalized Velocity')\n",
    "plt.title('Viscous Fluid Flow in a Pipe') \n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Here is a way to do this so each plot is updated in the same loop (as opposed to writing two loops):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "plt.close('all')\n",
    "\n",
    "#Set up Figure 1\n",
    "fig1 = plt.figure(1)\n",
    "plt.xlabel(' Radius (m)')\n",
    "plt.ylabel('Velocity (m/s)')\n",
    "plt.title('Viscous Fluid Flow in a Pipe') \n",
    "plt.legend()\n",
    "\n",
    "#Set up Figure 2\n",
    "fig2 = plt.figure(2)\n",
    "plt.xlabel(' Radius (m)')\n",
    "plt.ylabel('Normalized Velocity')\n",
    "plt.title('Viscous Fluid Flow in a Pipe') \n",
    "plt.legend()\n",
    "\n",
    "for mu in mus:\n",
    "    labelString = 'Viscosity = %2.2f' % mu\n",
    "    v = velocity(r,dP=dP,n=n,R=R,mu=mu)\n",
    "    plt.figure(1)\n",
    "    plt.plot(r,v,'-', label = labelString)\n",
    "    #Now plot normalized curve in other window\n",
    "    norm_v = velocity(r,dP=dP,n=n,R=R,mu=mu) / velocity(0.0,dP=dP,n=n,R=R,mu=mu)\n",
    "    plt.figure(2)\n",
    "    plt.plot(r,norm_v,'-', label = labelString)\n",
    "\n",
    "# Now call the legend() command for each plot\n",
    "plt.figure(1)\n",
    "plt.legend()\n",
    "plt.figure(2)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "#### Animations\n",
    "Computer animations are simply still frames rapidly replaced to give the illusion of motion. Matplotlib has a special library to handle animations.  The most common type of animation is a FuncAnimation inwhich a updating function is called between each frame to change the curve. Saving the animation is a bit tricky.  It has save() function which can write the sequence of images to a movie file, but an external library (not Python) MUST be installed!  A common free on for MPEGs is ffmpeg and is available for all platforms.  An example of this is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "import numpy as np\n",
    "\n",
    "plt.close('all')\n",
    "fig=plt.figure()\n",
    "\n",
    "#Set up the initial plot.  Note 'plot()' returns a line object.\n",
    "#  ALSO NOTE the ',' after the name of the line object - forces it to be a tuple\n",
    "x=np.linspace(0,2*np.pi,100)\n",
    "line,=plt.plot(x,np.sin(x),lw=2)\n",
    "plt.grid()\n",
    "\n",
    "\n",
    "def update(amp):\n",
    "    # Set new values for the y values of the function represented by 'line' and return\n",
    "    # the updated 'line'.  Note the trailing comma...\n",
    "    line.set_ydata(amp*np.sin(x))\n",
    "    return line,\n",
    "\n",
    "def init():\n",
    "    # plot an array of \"not a numbers\" (np.nan) and set the ydata\n",
    "    # This clears the orifinal plot, else it will stay there\n",
    "    line.set_ydata(np.zeros(len(x))+np.nan)\n",
    "    return line,\n",
    "\n",
    "#Create an array of values to change for each frame\n",
    "amps=np.cos(np.linspace(0,2*np.pi,100))\n",
    "\n",
    "#Run the animation\n",
    "ani = animation.FuncAnimation(fig,              #the figure window to animate\n",
    "                              update,           #function to update line\n",
    "                              amps,             #array of parameter to change\n",
    "                              interval=10,       #delay in ms\n",
    "                              init_func=init,   #function to initialize\n",
    "                              blit=True)        #employ blit refreshing\n",
    "\n",
    "##Saves the animation to a MPEG file. Seems a bit flaky, crashes sometimes.\n",
    "## Note this REQUIRES 'mencoder' or 'ffmpeg' to be installed. Open source mpeg encoders.\n",
    "ani.save('animated_sine.mp4', writer='ffmpeg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "A more 'brute force' method is to write a bunch of graphic files to the disk and stitch them together using other software.  The following does that.  Note the file names are updated for each frame in some logical manner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "plt.close('all')\n",
    "basename='frame_'\n",
    "output=True\n",
    "\n",
    "plt.ion()\n",
    "x=np.linspace(-2*np.pi,2*np.pi,200)\n",
    "y=np.sin(x)\n",
    "\n",
    "#Note the ',' after line - plot returns a \"line type\" object, this forces to be a tuple\n",
    "line,=plt.plot(x,y,'-',lw=2,animated=True)\n",
    "#show()\n",
    "\n",
    "amps=np.cos(np.linspace(0,10,30))\n",
    "\n",
    "for i in range(len(amps)):\n",
    "\tif output: \n",
    "\t\tcurrentFile = basename+'%03i.png'%i\n",
    "\t\tplt.savefig(currentFile)\n",
    "\tline.set_ydata(amps[i]*y)\n",
    "\t#t.sleep(0.1)\n",
    "\tplt.draw()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now you can either pull all the images into another program, or even call another command line program from within python to bind the stills together into an animation file.  An animated GIF in this case, using ImageMagick's very handy 'convert' command line program.  NOTE: ImageMagick is standard on Linux and used to be on Mac, but I noticed my new mac does not have it installed.  Hit google to find out how to install it on your machine if you like.  More on 'subprocess' and 'glob' next week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Conversion was not successful!\n"
     ]
    }
   ],
   "source": [
    "from subprocess import call\n",
    "import os, glob\n",
    "\n",
    "# 'call' returns '0' if successful and '1' of there was an error \n",
    "result = call(['convert','-loop', '0', 'frame*.png','standingwave.gif'])\t\n",
    "\n",
    "# 'call' will return '0' if successful and '1' if there was a problem.  We can use that to\n",
    "# check before removing all the input files.\n",
    "if result ==0:\n",
    "    print \"Conversion was succesfull.  Cleaning up input files...\"\n",
    "    framefiles = glob.glob('frame_*.png')\n",
    "    for file in framefiles:\n",
    "        os.remove(file)\n",
    "else: print \"Conversion was not successful!\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}