{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# D-001 -- Euler-Lagrange verification\\n",
    "\\n",
    "Checks that a symbolic Euler-Lagrange routine, applied to the double-pendulum ",
    "Lagrangian `L = T - V`, reproduces the standard textbook equations of motion. ",
    "This is the notebook behind the pass badge on [physics2u.com/compute](https://physics2u.com/compute/euler-lagrange-equation.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build the Lagrangian and generate the equations of motion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import symbols, diff, simplify, sin, cos, sympify\n",
    "from sympy.physics.mechanics import dynamicsymbols\n",
    "\n",
    "# IMPORTANT: use dynamicsymbols' own time symbol, not a separately-created one.\n",
    "# `symbols('t')` and `dynamicsymbols._t` print identically but are different\n",
    "# objects -- differentiating with the wrong one silently returns zero instead\n",
    "# of raising an error. This is the single most common trap in this derivation.\n",
    "t = dynamicsymbols._t\n",
    "m1, m2, l1, l2, g = symbols('m1 m2 l1 l2 g', positive=True)\n",
    "theta1, theta2 = dynamicsymbols('theta1 theta2')\n",
    "theta1d, theta2d = diff(theta1, t), diff(theta2, t)\n",
    "\n",
    "# positions of the two bobs (angles measured from vertical)\n",
    "x1 = l1*sin(theta1);            y1 = -l1*cos(theta1)\n",
    "x2 = x1 + l2*sin(theta2);       y2 = y1 - l2*cos(theta2)\n",
    "x1d, y1d = diff(x1, t), diff(y1, t)\n",
    "x2d, y2d = diff(x2, t), diff(y2, t)\n",
    "\n",
    "# Lagrangian: L = T - V\n",
    "T = sympify(1)/2*m1*(x1d**2 + y1d**2) + sympify(1)/2*m2*(x2d**2 + y2d**2)\n",
    "V = m1*g*y1 + m2*g*y2\n",
    "L = T - V\n",
    "\n",
    "def euler_lagrange(L, q):\n",
    "    \"\"\"d/dt(dL/dq_dot) - dL/dq, for a single generalized coordinate q.\"\"\"\n",
    "    qd = diff(q, t)\n",
    "    return simplify(diff(diff(L, qd), t) - diff(L, q))\n",
    "\n",
    "eom1 = euler_lagrange(L, theta1)   # generated equation of motion, coordinate 1\n",
    "eom2 = euler_lagrange(L, theta2)   # generated equation of motion, coordinate 2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare against the textbook form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# standard textbook form of the double-pendulum equations of motion\n",
    "theta1dd, theta2dd = diff(theta1, t, 2), diff(theta2, t, 2)\n",
    "\n",
    "textbook_eom1 = ((m1+m2)*l1*theta1dd + m2*l2*theta2dd*cos(theta1-theta2)\n",
    "                  + m2*l2*theta2d**2*sin(theta1-theta2) + (m1+m2)*g*sin(theta1))\n",
    "textbook_eom2 = (l2*theta2dd + l1*theta1dd*cos(theta1-theta2)\n",
    "                  - l1*theta1d**2*sin(theta1-theta2) + g*sin(theta2))\n",
    "\n",
    "# The generated equations come out scaled by l1 and l2*m2 respectively --\n",
    "# that's the natural units of \"d/dt(dL/dq_dot) - dL/dq\" for these coordinates,\n",
    "# not an error. Compare like with like:\n",
    "diff1 = simplify(eom1 - l1*textbook_eom1)\n",
    "diff2 = simplify(eom2 - l2*m2*textbook_eom2)\n",
    "\n",
    "assert diff1 == 0 and diff2 == 0, \"generated equations do not match the textbook form\"\n",
    "print(\"PASS -- both generated equations match the textbook form exactly.\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}