{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    " <center><h2>Associated Legendre Polynomials </h2></center>\n",
    " \n",
    "    The associated Legendre polynomials  are solutions of the equation:\n",
    "    $$ \\frac{d}{dx}\\Bigl[(1-x^2)\\frac{dP_\\ell^m(x)}{dx}\\Bigr]+\\Bigl[\\ell(\\ell+1) -\\frac{m^2}{1-x^2}\\Bigr] P_\\ell^m(x),  \\ \\ \\ -\\ell \\leq m\\leq \\ell. $$\n",
    "         \n",
    "        We use rk4 to solve this equation for $\\ell = 4, m=2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n",
    "    by RH Landau, MJ Paez, and CC Bordeianu (deceased)\n",
    "    Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n",
    "    C Bordeianu, Univ Bucharest, 2020. \n",
    "    Please respect copyright & acknowledge our work.\"\"\"\n",
    "\n",
    "# Plm.py: Associated Legendre Polynomials via Integration\n",
    "\n",
    "import numpy as np, matplotlib.pylab as plt  \n",
    "#from rk4Algor import rk4Algor if in same directory\n",
    "\n",
    "CosTheta = np.zeros((1999),float)\n",
    "Plm = np.zeros((1999),float)  \n",
    "y = [0]*(2);   dCos = 0.001   \n",
    "el = 4;  m = 2   # m intger  m<=el,   m = 1,2,3,...\n",
    "if el == 0 or el == 2: y[0] = 1\n",
    "if (el>2 and (el)%2  ==  0):      \n",
    "    if m  ==  0: y[0] = -1\n",
    "    elif( m>0 ): y[0] = 1\n",
    "    elif m<0 and abs(m)%2 == 0: y[0] = 1\n",
    "    elif m<0 and abs(m)%2 == 1: y[0] = -1\n",
    "if (el>2 and el%2 == 1) :\n",
    "    if m == 0: y[0] = 1\n",
    "    elif m>0:  y[0] = -1\n",
    "    elif m<0:  y[0] = 1\n",
    "y[1] = 1\n",
    "\n",
    "def f(Cos, y):                           # RHS of equation\n",
    "    rhs = np.zeros(2)                # Declare array dimension\n",
    "    rhs[0] = y[1]\n",
    "    rhs[1] = 2*Cos*y[1]/(1-Cos**2)-(el*(el+1)\n",
    "    \t    -m**2/(1-Cos**2))*y[0]/(1-Cos**2)\n",
    "    return rhs\n",
    "\n",
    "f(0,y)     # Call function for xi = 0 with init conds.\n",
    "i = -1\n",
    "def rk4Algor(t, h, N, y, f):\n",
    "    k1=np.zeros(N); k2=np.zeros(N); k3=np.zeros(N); k4=np.zeros(N);\n",
    "    k1 = h*f(t,y)                             \n",
    "    k2 = h*f(t+h/2.,y+k1/2.)\n",
    "    k3= h*f(t+h/2.,y+k2/2.)\n",
    "    k4= h*f(t+h,y+k3)\n",
    "    y=y+(k1+2*(k2+k3)+k4)/6.\n",
    "    return y    \n",
    "for Cos in np.arange(-0.999999, 1-dCos, dCos):\n",
    "    i = i+1\n",
    "    CosTheta[i] = Cos\n",
    "    y = rk4Algor(Cos, dCos, 2, y, f)  # call runge kutt\n",
    "    Plm[i] =  y[0]   #\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(CosTheta,Plm)\n",
    "plt.grid()\n",
    "plt.title('Unormalized $\\mathbf{P_l^m(Cos)}$')\n",
    "plt.xlabel('cos(theta)')\n",
    "plt.ylabel('$\\mathbf{P_l^m(Cos)}$')\n",
    "plt.show()"
   ]
  }
 ],
 "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": 2
}