pysics

{
 "metadata": {
  "name": "",
  "signature": "sha256:98fcea68ba105712b72f7e888ccceecb26f7e7b0f814f9f6a58c21f28de2453b"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Introduction\n",
      "------------\n",
      "\n",
      "When I started studying physics, I noticed that much could be learned with programming, and sadly many teachers and students don't make use of this wonderful tool. Pysics is a little project I've been doing to show other students and enthusiasts how much we can learn from computational physics using Python. Because I am a beginner, there are no complicated codes, but they are not very optimized, so bear with me.\n",
      "\n",
      "Notice 1: this notebook does not require more packages than the ones you probably already have installed with IPython. So, you are good to go! \n",
      "\n",
      "Notice 2: you can see the line number of the codes by selecting the code cell, pressing Esc and then the key L on the keyboard.\n",
      "\n",
      "Notice 3: run the following code cell in order to see matplotlib plots on the notebook."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Electric charges\n",
      "----------------\n",
      "\n",
      "The following code plots the electric potential of two electric charges in a grayscale. The method used is quite simple: first I define a series of variables (lines 5 to to 15), and then define a function that returns the value of the electric potential created by an electric charge (lines 17 and 18). On lines 20 to 23, the code creates the numpy arrays that will hold the results.\n",
      "\n",
      "The electric potential is defined as\n",
      "$$ \\phi = \\frac{q}{4 \\pi \\epsilon_0 r} $$\n",
      "where $q$ is the electric charge and $r$ is the distance from the charge. $4\\pi \\epsilon_0$ is a constant.\n",
      "\n",
      "From the theory, we know that the resulting electric potential from various charges, in a specific point of space, is just the sum of the potential produced by each charge at that point (see line 35). Maybe this is not very clear on the code, but I defined the origin at the exact middle of the 2-dimensional grid."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from math import sin,log,sqrt,pi,exp\n",
      "import numpy as np\n",
      "from pylab import imshow,show,gray\n",
      "\n",
      "q1 = 1.0            # Electric charge 1, in C\n",
      "q2 = 1.0            # Electric charge 2, in C\n",
      "eps0 = 8.854E-12    # Vacuum permitivity, in F/m**2\n",
      "dx = 0.001          # sample distance in x, in m\n",
      "dy = 0.001          # sample distance in y, in m\n",
      "Dx = 1.0            # grid size in x, in m\n",
      "Dy = 1.0            # grid size in y, in m\n",
      "rx1 = -0.05         # position of q1 in the axis x, in m\n",
      "ry1 = 0.0           # position of q1 in the axis y, in m\n",
      "rx2 = 0.05          # position of q2 in the axis x, in m\n",
      "ry2 = 0.0           # position of q2 in the axis y, in m\n",
      "\n",
      "def phi(q,r):\n",
      "    return q/(4*pi*eps0*r)\n",
      "    \n",
      "Nx = int(Dx/dx)     # number of samples in x\n",
      "Ny = int(Dy/dy)     # number of samples in y\n",
      "EP = np.zeros([Nx,Ny])\n",
      "logEP = np.zeros([Nx,Ny])       # using logarithmic scale makes the plot more eye-appealing... sometimes\n",
      "\n",
      "ry = Dy/2           # starting position in the axis y: bottom         \n",
      "for k in range(Ny):\n",
      "    rx = -Dx/2      # starting position in the axis x: left\n",
      "    for i in range(Nx):\n",
      "        r1 = sqrt((rx1-rx)**2+(ry1-ry)**2)   # Distance from charge 1\n",
      "        r2 = sqrt((rx2-rx)**2+(ry2-ry)**2)   # Distance from charge 2\n",
      "        if abs(r1) < 1E-10 or abs(r2) < 1E-10:\n",
      "            EP[k][i] = None\n",
      "            logEP[k][i] = None\n",
      "        else:\n",
      "            EP[k][i] = phi(q1,r1)+phi(q2,r2)\n",
      "            if abs(EP[k][i]) > 1E-10 and EP[k][i]>0:        # This crap here is to deal with the log of\n",
      "                logEP[k][i] = log(EP[k][i])                 # negative values of electric potential.\n",
      "            elif abs(EP[k][i]) > 1E-10 and EP[k][i]<0:      # If you're using different signal \n",
      "                logEP[k][i] = -log(abs(EP[k][i]))           # charges, it won't look cool anyway.\n",
      "            else:                                           # In this case, you might want to use\n",
      "                logEP[k][i] = None                          # the linear scale instead of logarithmic.\n",
      "        rx += dx\n",
      "    ry -= dy\n",
      "  \n",
      "imshow(logEP)\n",
      "gray()\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Coupled pendula\n",
      "---------------\n",
      "\n",
      "Here's another cool one. Actually, I wrote this code after I got very frustrated with an experiment with springs in the lab. We couldn't get the coupled pendula to oscillate in one of its harmonics, because required a lot of finesse and practice at it. So I decided to just simulate it with a Python program.\n",
      "\n",
      "From Newton's classical mechanics, we can specify exactly how a system will evolve if we have the positions and initial velocities of all the particles that constitute the system, given that we know how the particles interact. And this is what we do from lines 7 to 19: get the initial positions and velocities, and then the masses and spring specifications (which are related to how the pendula interact).\n",
      "\n",
      "This system can be described by a set of differential equations, as follows:\n",
      "\n",
      "$\\ddot{\\theta}_1 = -\\frac{g}{l}\\theta_1 + \\frac{k_1}{m_1}(\\theta_2-\\theta_1)$\n",
      "\n",
      "$\\ddot{\\theta}_2 = -\\frac{g}{l}\\theta_2 + \\frac{k_1}{m_2}(\\theta_2-\\theta_1) + \\frac{k_2}{m_2}(\\theta_3-\\theta_2)$\n",
      "\n",
      "$\\ddot{\\theta}_3 = -\\frac{g}{l}\\theta_3 + \\frac{k_2}{m_3}(\\theta_2-\\theta_3)$\n",
      "\n",
      "where $\\theta_n$ is the position of mass n, $k_1$ is the spring constant of the spring that connects masses $m_1$ and $m_2$, and $k_2$ is for the spring connecting masses $m_2$ and $m_3$, $g$ is the acceleration of gravity and $l$ is the length of the springs when relaxed.\n",
      "\n",
      "We solve the differential equations that rule this system in just a few lines (23 to 28). SciPy's differential equation solver is very weird to implement, so I'll leave that up for you to figure it out. After the calculation, we just plot it!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from math import pi\n",
      "import matplotlib.pyplot as plt\n",
      "from scipy.integrate import odeint\n",
      "from numpy import arange, array\n",
      "\n",
      "# Initial data\n",
      "pos1 = 5.0*pi/180   # Angular initial position of the pendulum 1, in radians\n",
      "vel1 = 0.0          # Angular initial velocity of the pendulum 1, in rad/s\n",
      "pos2 = -5.0*pi/180  # Angular initial position of the pendulum 2, in radians\n",
      "vel2 = 0.0          # Angular initial velocity of the pendulum 2, in rad/s\n",
      "pos3 = 5.0*pi/180   # Angular initial position of the pendulum 3, in radians\n",
      "vel3 = 0.0          # Angular initial velocity of the pendulum 3, in rad/s\n",
      "m1 = 0.1353         # Mass of the pendulum 1\n",
      "m2 = 0.1352         # Mass of the pendulum 2\n",
      "m3 = 0.1349         # Mass of the pendulum 3\n",
      "l = 41*0.01         # Lenght of the strings (in this simulation, they must be equal)\n",
      "g = 9.785           # Gravity acceleration\n",
      "k1 = 2.90           # Spring 1 constant\n",
      "k2 = 2.39           # Spring 2 constant\n",
      "deltat = 10.0       # Time interval for the simulation\n",
      "\n",
      "# Solving the differential equations system\n",
      "def deriv(y, t):\n",
      "    return array([y[1], k1*(y[2]-y[0])/m1-(g/l)*y[0], y[3], k2*(y[4]-y[2])/m2-k1*(y[2]-y[0])/m2-(g/l)*y[2], y[5], -k2*(y[4]-y[2])/m3-(g/l)*y[4]])\n",
      "\n",
      "time = arange(0.0, deltat, 0.01)\n",
      "yinit = array([pos1, vel1, pos2, vel2, pos3, vel3])\n",
      "y = odeint(deriv, yinit, time)\n",
      "\n",
      "# Plotting the graph\n",
      "ax = plt.axes(ylim = (-0.20, 0.20))\n",
      "plt.plot(time, y[:,0], label='Mass 1')\n",
      "plt.plot(time, y[:,2], label='Mass 2')\n",
      "plt.plot(time, y[:,4], label='Mass 3')\n",
      "plt.xlabel('t (s)')\n",
      "plt.ylabel('Amplitude')\n",
      "plt.title('Graph title')\n",
      "plt.grid(True)\n",
      "plt.legend()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There is a version of this code that is capable of simulating three coupled pendula and produce an animation, but sadly IPython notebook does not support matplotlib's animations, as far as I know. You can download it here: http://goo.gl/MBjSfm"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Bessel functions and shaving some for-loops\n",
      "-----------------\n",
      "\n",
      "Python is an amazing language to code, but it has its drawbacks. One of them, and probably the most worrisome, is that its for loops are very slow to compute. But there are ways to optimize it. for instance by using NumPy arrays and SciPy's functions.\n",
      "\n",
      "In the following example, we compute Bessel functions of the first kind using this definition:\n",
      "$$J_m = \\frac{1}{\\pi} \\int_0^\\pi \\cos{(m\\theta - x\\sin{\\theta})}d\\theta$$\n",
      "with $x$ varying from 0 to 20 and $m = 0$ (J0), $m = 1$ (J1) and $m = 2$ (J2). We use SciPy's Simpson integration function to compute the integrals.\n",
      "\n",
      "Keep in mind that this code is a bit slow to run, so be patient. In my computer, it took about 20 seconds."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy.integrate as si\n",
      "import matplotlib.pyplot as plt\n",
      "from timeit import default_timer as timer\n",
      "\n",
      "def func(x,m,theta):\n",
      "    return np.cos(m*theta-x*np.sin(theta))\n",
      "\n",
      "x0 = 0.0        # start\n",
      "x1 = 20.0       # end\n",
      "N = 1000        # number of points\n",
      "x = np.linspace(x0,x1,N)\n",
      "theta = np.linspace(0.0,np.pi,N)\n",
      "\n",
      "start = timer()   # Starting the timer to measure computation time\n",
      "J0 = np.zeros(N,float)\n",
      "for i in range(N):\n",
      "    y0 = np.zeros(N,float)\n",
      "    for k in range(N):\n",
      "        y0[k] = func(x[i],0,theta[k])\n",
      "    J0[i] = si.simps(y=y0,x=theta)    \n",
      "\n",
      "J1 = np.zeros(N,float)\n",
      "for i in range(N):\n",
      "    y1 = np.zeros(N,float)\n",
      "    for k in range(N):\n",
      "        y1[k] = func(x[i],1,theta[k])\n",
      "    J1[i] = si.simps(y=y1,x=theta)\n",
      "    \n",
      "J2 = np.zeros(N,float)\n",
      "for i in range(N):\n",
      "    y2 = np.zeros(N,float)\n",
      "    for k in range(N):\n",
      "        y2[k] = func(x[i],2,theta[k])\n",
      "    J2[i] = si.simps(y=y2,x=theta)\n",
      "    \n",
      "time = timer()-start\n",
      "print \"Calculation time = %f seconds\" % time\n",
      "    \n",
      "plt.plot(x,J0,label=r'$J_0$')\n",
      "plt.plot(x,J1,label=r'$J_1$')\n",
      "plt.plot(x,J2,label=r'$J_2$')\n",
      "plt.xlabel('x')\n",
      "plt.ylabel('J')\n",
      "plt.legend()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It is possible to shorten the previous code by using a bit of numpy array hacking. Actually, we can compute lines 16 to 35 from the previous code in just 3 lines! We can save us the for loops by using list comprehension (see this link: http://goo.gl/5qxntm). We don't get results faster though :("
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy.integrate as si\n",
      "import matplotlib.pyplot as plt\n",
      "from timeit import default_timer as timer\n",
      "\n",
      "def func(x,m,theta):\n",
      "    return np.cos(m*theta-x*np.sin(theta))\n",
      "\n",
      "x0 = 0.0        # start\n",
      "x1 = 20.0       # end\n",
      "N = 1000        # number of points\n",
      "x = np.linspace(x0,x1,N)\n",
      "theta = np.linspace(0.0,np.pi,N)\n",
      "\n",
      "start = timer()\n",
      "J0 = np.array([si.simps(y=np.array([func(xk,0,tk) for tk in theta]),x=theta) for xk in x])\n",
      "J1 = np.array([si.simps(y=np.array([func(xk,1,tk) for tk in theta]),x=theta) for xk in x]) \n",
      "J2 = np.array([si.simps(y=np.array([func(xk,2,tk) for tk in theta]),x=theta) for xk in x]) \n",
      "time = timer()-start\n",
      "print \"Calculation time = %f seconds\" % time\n",
      "    \n",
      "plt.plot(x,J0,label=r'$J_0$')\n",
      "plt.plot(x,J1,label=r'$J_1$')\n",
      "plt.plot(x,J2,label=r'$J_2$')\n",
      "plt.xlabel('x')\n",
      "plt.ylabel('J')\n",
      "plt.legend()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "But if you really want to optimize things and get fast results, there is nothing better than using SciPy's library of special functions. Just import scipy.special (line 3) and use the function jv(v,x). And with this, we compute the Bessel functions in a fraction of a second. Now we are talking!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy.integrate as si\n",
      "import scipy.special as sp\n",
      "import matplotlib.pyplot as plt\n",
      "from timeit import default_timer as timer\n",
      "\n",
      "def func(x,m,theta):\n",
      "    return np.cos(m*theta-x*np.sin(theta))\n",
      "\n",
      "x0 = 0.0        # start\n",
      "x1 = 20.0       # end\n",
      "N = 1000        # number of points\n",
      "x = np.linspace(x0,x1,N)\n",
      "theta = np.linspace(0.0,np.pi,N)\n",
      "\n",
      "start = timer()\n",
      "J0 = np.array([sp.jv(0,xk) for xk in x])\n",
      "J1 = np.array([sp.jv(1,xk) for xk in x])\n",
      "J2 = np.array([sp.jv(2,xk) for xk in x]) \n",
      "time = timer()-start\n",
      "print \"Calculation time = %f seconds\" % time\n",
      "    \n",
      "plt.plot(x,J0,label=r'$J_0$')\n",
      "plt.plot(x,J1,label=r'$J_1$')\n",
      "plt.plot(x,J2,label=r'$J_2$')\n",
      "plt.xlabel('x')\n",
      "plt.ylabel('J')\n",
      "plt.legend()\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Diffraction limit of a telescope\n",
      "----------------\n",
      "\n",
      "I love astronomy, so I couldn't just pass slipping an exercise about it on this notebook. You see, the amount of detail we see in astronomical images are limited by the diffraction of light in the telescope. When a point source of light is seen through a telescope, it is not exactly a dot, but a very bright smudge with concentric rings. It turns out that this diffraction pattern follows the form of a Bessel function of the first kind and of order $m = 1$. Well, we saw on the previous exercises how to compute Bessel functions using various methods. \n",
      "\n",
      "In this example, what was $x$ in the previous exercise will be the distance of a point to the center of the image (which I will call $r$). The intensity of light in function of $r$ is given by the following equation:\n",
      "$$I(r) = \\left( \\frac{J_1(kr)}{kr} \\right)^2$$\n",
      "where $k = 2 \\pi/\\lambda$ and $\\lambda$ is the wavelength of the light being observed. \n",
      "\n",
      "Computing this have some caveats, though. First of all, the wavenlengths $\\lambda$ of visible light are very small (hundreds of nanometers), and so are the sizes $r$ of the diffraction patters (of the order of a few micrometers), and as we all know, computers don't particularly like infinitesimally small numbers, because they tend to treat these numbers as zeroes. The solution to that is to assign an arbitary unit to the dimension of sizes and use more palpable numbers, as long as we are consistent with reality.\n",
      "\n",
      "The other caveat is another limitation of computers: they don't know how to deal with divisions by zero. Well, luckily we are not machines, and we can do some tricks. For instance, when $r \\rightarrow 0$, the computer will complain and yell at you for asking a division by zero. However, mathematics says that\n",
      "$$ lim_{x \\rightarrow 0} \\frac{J_1(x)}{x} = \\frac{1}{2} $$\n",
      "so then we have a way to circumvent the computer's lack of ability to deal with mathematical limits.\n",
      "\n",
      "Okay, now here is the deal. I noticed from my tinkering that the best way to do away with the small numbers is to shove a gigantic multiplicative constant $G$ in the intensity function, like so:\n",
      "$$I'(r) = \\left( \\frac{J_1(Gkr)}{Gkr} \\right)^2$$\n",
      "Notice that the limit when $kr \\rightarrow 0$ still holds. If we assign a small $G$, the computer will complain because it will still think we are dividing by zero. On the other hand, if we use an exceedingly large $G$, the diffraction pattern will be too small. I found that a sweetspot for $G$ is $5 \\times 10^{13}$, when using SI units for $r$ and $k$.\n",
      "\n",
      "Well, without further ado, let's do this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy.special as sp\n",
      "from pylab import imshow,show,gray,hot\n",
      "import matplotlib.pyplot as plt\n",
      "from timeit import default_timer as timer\n",
      "\n",
      "def func(x,m,theta):\n",
      "    return np.cos(m*theta-x*np.sin(theta))\n",
      "\n",
      "wl = 500E-9     # wavelength, in meters\n",
      "size = 2.0      # image size, in micrometers\n",
      "N = 500         # plot resolution\n",
      "G = 5E12        # Gigantic constant\n",
      "\n",
      "start = timer()\n",
      "dists = np.array([[np.sqrt(float((N/2-i)**2+(N/2-j)**2)) for i in range(N)] for j in range(N)])\n",
      "dists = 1E-6*np.sqrt(size)*dists/np.max(dists)     # distances matrix\n",
      "\n",
      "image = np.zeros([N,N],float)\n",
      "for i in range(N):\n",
      "    for j in range(N):\n",
      "        if abs(i-N/2) < 2 and abs(j-N/2)<2:     # This deals with the limit r -> 0\n",
      "            image[i,j] = 0.5**2\n",
      "        else:\n",
      "            image[i,j] = (sp.jv(1,G*2.*np.pi*wl*dists[i,j])/(G*2.*np.pi*wl*dists[i,j]))**2\n",
      "time = timer()-start\n",
      "print \"Calculation time = %f seconds\" % time\n",
      "\n",
      "imshow(image)\n",
      "gray()\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Well, now that's awkward. Where the heck is the diffraction pattern? \n",
      "\n",
      "Here is the thing: it's very subtle. It turns out the brightness of the astronomical objects completely wash out the diffraction pattern. In order to see it, we need to mask out the central brightness. A solution to that is to assign a maximum value for the density plot we are using, and it's very simple. Just use the argument vmax on pylab's imshow. I found that a value of v = 0.01 is a pretty good mask. It's like we have a really tall peak at the center of the image, and simply cut it at a particular height."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(image,vmax=0.01)\n",
      "gray()\n",
      "show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that is a beautiful diffraction pattern (crank up the brightness and constrast of your screen to make the fringes pop-out). If you have a telescope, you'll notice that when you slightly de-focus a point source, it shows a very similar pattern."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}