{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "from galpy import potential"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rotation curves\n",
    "## Spherical potentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "kp= potential.KeplerPotential(normalize=1.)\n",
    "ip= potential.IsochronePotential(normalize=1.)\n",
    "lp= potential.LogarithmicHaloPotential(normalize=1.)\n",
    "pp= potential.PowerSphericalPotential(normalize=1.,alpha=0.1)\n",
    "np= potential.NFWPotential(normalize=1.,a=1.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kp.plotRotcurve()\n",
    "ylim(0.,2.)\n",
    "#ip.plotRotcurve(overplot=True)\n",
    "lp.plotRotcurve(overplot=True)\n",
    "#pp.plotRotcurve(overplot=True)\n",
    "np.plotRotcurve(overplot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flattened potentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mp= potential.MiyamotoNagaiPotential(normalize=1.,a=.3,b=0.03)\n",
    "flp= potential.LogarithmicHaloPotential(normalize=1.,q=0.8,core=0.5)\n",
    "dp= potential.DoubleExponentialDiskPotential(normalize=1.,hr=.3,hz=0.03)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mp.plotRotcurve(Rrange=[0.,10.])\n",
    "flp.plotRotcurve(overplot=True,Rrange=[0.,10.])\n",
    "dp.plotRotcurve(overplot=True,Rrange=[0.,10.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Densities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spherical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kp.plotDensity(rmin=0.,rmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lp.plotDensity(log=True,rmin=0.1,rmax=2.,zmin=0.1,zmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.plotDensity(log=True,rmin=0.1,rmax=2.,zmin=0.1,zmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pp.plotDensity(log=True,rmin=0.1,rmax=2.,zmin=0.1,zmax=2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flattened"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mp.plotDensity(log=True,rmin=0.1,rmax=4.,zmin=-2.,zmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dp.plotDensity(log=True,rmin=0.1,rmax=4.,zmin=-2.,zmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "flp.plotDensity(log=True,rmin=0.1,rmax=4.,zmin=0.1,zmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mp.plot(rmin=0.1,rmax=4.,zmin=-2.,zmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dp.plot(rmin=0.1,rmax=4.,zmin=-2.,zmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "flp.plot(rmin=0.1,rmax=4.,zmin=-2.,zmax=2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "flp= potential.LogarithmicHaloPotential(normalize=1.,q=0.5,core=0.)\n",
    "flp.plotDensity(log=True,rmin=0.1,rmax=4.,zmin=0.1,zmax=2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The potential of the Milky Way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from galpy.potential import MWPotential2014, plotPotentials, plotDensities, plotRotcurve, plotEscapecurve, \\\n",
    "    evaluatezforces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plotRotcurve(MWPotential2014)\n",
    "plotRotcurve(MWPotential2014[0],overplot=True)\n",
    "plotRotcurve(MWPotential2014[1],overplot=True)\n",
    "plotRotcurve(MWPotential2014[2],overplot=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plotDensities(MWPotential2014,log=True,rmin=0.1,rmax=2.,zmin=-0.5,zmax=.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plotPotentials(MWPotential2014,rmin=0.1,rmax=2.,zmin=-0.5,zmax=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "zs= numpy.linspace(0.,0.25,1001)\n",
    "plot(zs,[-evaluatezforces(1.,z,MWPotential2014) for z in zs])\n",
    "#plot(zs,[-evaluatezforces(1.,z,MWPotential2014[0]) for z in zs])\n",
    "#plot(zs,[-evaluatezforces(1.,z,MWPotential2014[1]) for z in zs])\n",
    "#plot(zs,[-evaluatezforces(1.,z,MWPotential2014[2]) for z in zs])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Motion in the meridional plane"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's investigate orbits in the flattened logarithmic potential:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from galpy.potential import LogarithmicHaloPotential\n",
    "from galpy.orbit import Orbit\n",
    "lp= LogarithmicHaloPotential(normalize=1.,q=0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rs= numpy.linspace(0.08,0.5,201)\n",
    "Lz= 0.2\n",
    "plot(rs,lp(rs,0.)+Lz**2./2./rs**2.,'b-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has a minimum at 0.2 (=Lz) because the rotation curve is flat. We can further also plot the effective potential as a function of (R,z):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lp.plot(rmin=0.12,rmax=1.,nrs=101,effective=True,Lz=0.2,xrange=[0.,1.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Orbits with the same E and L, but very different orbits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "o1= Orbit([0.35,0.,0.2/0.35,0.,numpy.sqrt(2.*(-0.8-lp(0.35,0.)-(0.2/0.35)**2./2.)),0.])\n",
    "ts= numpy.linspace(0.,50.,1001)\n",
    "o1.integrate(ts,lp)\n",
    "print \"First orbit: E, L = %f,%f\" % (o1.E(),o1.L()[:,2])\n",
    "o2= Orbit([0.35,0.35,0.2/0.35,0.,numpy.sqrt(2.*(-0.8-lp(0.35,0.)-(0.2/0.35)**2./2.-0.35**2./2.)),0.])\n",
    "o2.integrate(ts,lp)\n",
    "print \"Second orbit: E, L = %f,%f\" % (o2.E(),o2.L()[:,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "o1.plot(xrange=[0.,0.5],yrange=[-0.2,0.2])\n",
    "o2.plot(xrange=[0.,0.5],yrange=[-0.2,0.2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can calculate their maximmum height above the z=0 plane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print \"First orbit: zmax = %f\" % (o1.zmax())\n",
    "print \"Second orbit: zmax = %f\" % (o2.zmax())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Surface of section (Poincare section)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def surface_section(Rs,zs,vRs):\n",
    "    #Find points where the orbit crosses z from negative to positive\n",
    "    shiftzs= numpy.roll(zs,-1)\n",
    "    indx= (zs[:-1] < 0.)*(shiftzs[:-1] > 0.)\n",
    "    return (Rs[:-1][indx],vRs[:-1][indx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#First orbit\n",
    "Rs= o1.R(ts); zs= o1.z(ts); vRs= o1.vR(ts)\n",
    "sectRs, sectvRs= surface_section(Rs,zs,vRs)\n",
    "plot(sectRs,sectvRs,'bo')\n",
    "#Second orbit\n",
    "Rs= o2.R(ts); zs= o2.z(ts); vRs= o2.vR(ts)\n",
    "sectRs, sectvRs= surface_section(Rs,zs,vRs)\n",
    "plot(sectRs,sectvRs,'go')\n",
    "xlabel(r'$R$')\n",
    "ylabel(r'$v_R$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ts= numpy.linspace(0.,500.,100001)\n",
    "o1.integrate(ts,lp)\n",
    "o2.integrate(ts,lp)\n",
    "#First orbit\n",
    "Rs= o1.R(ts); zs= o1.z(ts); vRs= o1.vR(ts)\n",
    "sectRs, sectvRs= surface_section(Rs,zs,vRs)\n",
    "plot(sectRs,sectvRs,'bo',mec='none')\n",
    "#Second orbit\n",
    "Rs= o2.R(ts); zs= o2.z(ts); vRs= o2.vR(ts)\n",
    "sectRs, sectvRs= surface_section(Rs,zs,vRs)\n",
    "plot(sectRs,sectvRs,'go',mec='none')\n",
    "xlabel(r'$R$')\n",
    "ylabel(r'$v_R$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": 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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}