{
 "metadata": {
  "name": "",
  "signature": "sha256:f4cb8887c283a3e99aa23821c93f7e1841c2cb2840af296a0ba5d71a4956eed9"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Formula from [Kimura and Ohta 1969](http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=1212239&tool=pmcentrez&rendertype=abstract)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import romb,trapz, quad"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "integral = quad"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 68
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def integrand1(s,N,x):\n",
      "    if x == 1:\n",
      "        return 2*N*s * (1 - exp(-2*N*s))\n",
      "    return (1 - exp(-2*N*s*x) - exp(-(1-x)) - exp(-2*N*s))/(x*(1-x))\n",
      "\n",
      "def integrand2(s,N,x):\n",
      "    if x == 0:\n",
      "        return 0\n",
      "    return (exp(2*N*s*x) - 1) * (1 - exp(-2*N*s*x))/(x*(1-x))\n",
      "def T_kimura(s,N,x):\n",
      "    J1 = 2/(s*(1 - exp(-2*N*s))) * integral(lambda t: integrand1(s,N,t), x, 1)[0]\n",
      "    u = (1-exp(-2*N*s*x))/(1-exp(-2*N*s))\n",
      "    J2 = 2/(s*(1-exp(-2*N*s))) * integral(lambda t: integrand2(s,N,t), 0, x)[0]\n",
      "    return J1 + ((1-u)/u) * J2\n",
      "T_kimura(0.01, 1e5, 1e-5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}