numerical experiments

{
 "metadata": {
  "name": "",
  "signature": "sha256:e95905e5a9d37c0820b2b369abe5f4ce1bfa5f18118b9fd71b3439a3e0b89915"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from bash_command.bash_command import bash_command as bc\n",
      "import galaxy\n",
      "from scipy import integrate\n",
      "import time"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gal=galaxy.NukerGalaxy('NGC3605', init={'r1':1.E-2*galaxy.pc, 'r2':1.E4*galaxy.pc})"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import integrate\n",
      "\n",
      "start=time.time()\n",
      "r=2.\n",
      "4.*np.pi*galaxy.G*integrate.quad(lambda r1:gal.rho_stars(r1)*r1, r, gal.rmax_star)[0]\n",
      "print time.time()-start"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "23.8156661987\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/local/anaconda/lib/python2.7/site-packages/scipy/integrate/quadpack.py:321: IntegrationWarning: The occurrence of roundoff error is detected, which prevents \n",
        "  the requested tolerance from being achieved.  The error may be \n",
        "  underestimated.\n",
        "  warnings.warn(msg, IntegrationWarning)\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import integrate\n",
      "\n",
      "start=time.time()\n",
      "r=2.\n",
      "print 4.*np.pi*galaxy.G*integrate.quad(lambda r1:gal.rho_stars(r1)*r1, r, 2.E3)[0]\n",
      "print time.time()-start"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3.78489084003e+33\n",
        "5.24468994141\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import integrate\n",
      "\n",
      "start=time.time()\n",
      "r=2.\n",
      "print 4.*np.pi*galaxy.G*integrate.quad(lambda r1:gal.rho_stars(r1)*r1, r, 2.E3)[0]\n",
      "print time.time()-start"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3.78489084003e+33\n",
        "5.82942199707\n"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import integrate\n",
      "\n",
      "start=time.time()\n",
      "r=2.\n",
      "print 4.*np.pi*galaxy.G*integrate.quad(lambda r1:gal.rho_stars(r1)*r1, r, 2.E3)[0]\n",
      "print time.time()-start"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3.78489084003e+33\n",
        "5.38537812233\n"
       ]
      }
     ],
     "prompt_number": 103
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import integrate\n",
      "\n",
      "start=time.time()\n",
      "r=2.\n",
      "print integrate.romberg(lambda r1:gal.rho_stars(r1)*r1, r, 2.E3)\n",
      "print time.time()-start"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.51381253567e+39\n",
        "19.8025808334\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/local/anaconda/lib/python2.7/site-packages/scipy/integrate/quadrature.py:687: AccuracyWarning: divmax (10) exceeded. Latest difference = 6.205486e+36\n",
        "  AccuracyWarning)\n"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import integrate\n",
      "\n",
      "v_rho_stars=np.vectorize(gal.rho_stars)\n",
      "\n",
      "start=time.time()\n",
      "r=2.\n",
      "print 4.*np.pi*galaxy.G*integrate.quadrature(lambda r1:v_rho_stars(r1)*r1, r, 2.E3)[0]\n",
      "print time.time()-start"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3.78402134925e+33\n",
        "23.9541640282\n"
       ]
      }
     ],
     "prompt_number": 102
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import integrate\n",
      "\n",
      "v_rho_stars=np.vectorize(gal.rho_stars)\n",
      "\n",
      "start=time.time()\n",
      "r=2.\n",
      "print 4.*np.pi*galaxy.G*integrate.fixed_quad(lambda r1:v_rho_stars(r1)*r1, r, 2.E3, n=5)[0]\n",
      "print time.time()-start"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3.58294376606e+33\n",
        "0.148263931274\n"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# def integrand(params,r1):\n",
      "params=gal.params\n",
      "pc=galaxy.pc\n",
      "r=2.\n",
      "\n",
      "beta=params['beta']\n",
      "alpha=params['alpha']\n",
      "Gamma=params['gamma']\n",
      "rb=params['rb']*pc\n",
      "gamma=Gamma\n",
      "\n",
      "integrand=lambda r1:2.**(2.+(beta-gamma)/alpha)*\\\n",
      "    (Gamma*rb**alpha+beta*r1**alpha)*\\\n",
      "    r1**(-1-gamma)*((rb**alpha+r1**alpha)**(-1.-(beta-gamma)/alpha))*\\\n",
      "    (r1**2-(r*pc)**2)**0.5\n",
      "\n",
      "print galaxy.G*params['Ib']*params['Uv']*galaxy.M_sun*rb**params['beta']*integrate.quad(integrand, r*pc, 2.E3*pc)[0]/pc"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3.64707993363e+33\n"
       ]
      }
     ],
     "prompt_number": 101
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pickle\n",
      "import galaxy\n",
      "from scipy import integrate \n",
      "\n",
      "gal=pickle.load(open('/Users/aleksey/Second_Year_Project/hydro/NGC3605_3/vw_200.0/grid.p', 'rb'))\n",
      "r=2.\n",
      "\n",
      "integrand=lambda r1:galaxy.nuker_prime(r1, **gal.params)*(r1**2-r**2)**0.5\n",
      "#-4.*galaxy.G*gal.params['Uv']*integrate.quad(integrand, r, 2.E3)\n",
      "print (-4.*galaxy.G*galaxy.M_sun*gal.params['Uv']*integrate.quad(integrand, r, )[0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "4.8755920088e+33\n"
       ]
      }
     ],
     "prompt_number": 120
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import galaxy\n",
      "\n",
      "gal_dict=galaxy.nuker_params()\n",
      "params=gal_dict['A2052']\n",
      "params['gamma']=-0.5\n",
      "\n",
      "params['Uv']*galaxy.inverse_abel(galaxy.nuker_prime,0.5*params['rb'], **params)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 17,
       "text": [
        "-8.8333134044360833"
       ]
      }
     ],
     "prompt_number": 17
    }
   ],
   "metadata": {}
  }
 ]
}