{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Пример\n",
    "\n",
    "В качестве упражнения смоделируем полёт пушечного ядра. \n",
    "\n",
    "Пушка совершает выстрел сферическим ядром радиуса $r$ и массы $m$ под углом $\\theta_0$ к горизонту с начальной скоростью $v_0$. Смоделировать полёт снаряда. Считать, что на ядро действует [лобовое сопротивление](https://ru.wikipedia.org/wiki/%D0%9B%D0%BE%D0%B1%D0%BE%D0%B2%D0%BE%D0%B5_%D1%81%D0%BE%D0%BF%D1%80%D0%BE%D1%82%D0%B8%D0%B2%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5), \n",
    "т.е. сила сопротивления направлена против скорости движения $\\vec{v}$, её величина пропорциональна характерной площади $S$, плотности среды $\\rho$ и квадрату скорости $v^2$:\n",
    "\n",
    "$$\n",
    "\\vec{F} = - C_{F}\\frac{\\rho v}{2}S \\vec{v}.\n",
    "$$\n",
    "\n",
    "$C_{F}$ — [безразмерный аэродинамический коэффициент сопротивления](https://ru.wikipedia.org/wiki/%D0%9A%D0%BE%D1%8D%D1%84%D1%84%D0%B8%D1%86%D0%B8%D0%B5%D0%BD%D1%82_%D1%81%D0%BE%D0%BF%D1%80%D0%BE%D1%82%D0%B8%D0%B2%D0%BB%D0%B5%D0%BD%D0%B8%D1%8F_%D1%84%D0%BE%D1%80%D0%BC%D1%8B), который для сферы равен $0.47$. $S$ считают равным площади поперечного сечения в случае сферы ($\\pi r^2$). Обычно стандартной величиной плотности воздуха на уровне моря в соответствии с Международной стандартной атмосферой принимается значение 1.2250 кг/м³ (источник - [Wiki](https://ru.wikipedia.org/wiki/%D0%9F%D0%BB%D0%BE%D1%82%D0%BD%D0%BE%D1%81%D1%82%D1%8C_%D0%B2%D0%BE%D0%B7%D0%B4%D1%83%D1%85%D0%B0)).\n",
    "\n",
    "\n",
    "Совместим начало координат с пушкой, ось $Ox$ направим горизонтально в сторону полёта снаряда, а ось $Oy$ вертикально.\n",
    "\n",
    "```{figure} /_static/lecture_specific/hw_3/cannonball.png\n",
    ":scale: 30%\n",
    "```\n",
    "\n",
    "Тогда уравнения движения примут вид\n",
    "\n",
    "$$\n",
    "m \\dfrac{d \\vec{v}}{dt} = m\\vec{g} + \\vec{F}_{conp}\n",
    "$$\n",
    "\n",
    "или\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "m \\dfrac{dv_x}{dt} = - C_{F}\\dfrac{\\rho v}{2}S v_x, \\\\\n",
    "m \\dfrac{dv_y}{dt} = - mg - C_{F}\\dfrac{\\rho v}{2}S v_y. \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Найдите траекторию полёта снаряда."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.constants import g, degree, pi\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook\n",
    "\n",
    "# params\n",
    "theta0 = 45 * degree\n",
    "v0 = 10 # mps\n",
    "m = 1 # kg\n",
    "r = 0.1 # m\n",
    "v0x = v0 * np.cos(theta0)\n",
    "v0y = v0 * np.sin(theta0)\n",
    "\n",
    "\n",
    "# constants\n",
    "C_F = 0.47\n",
    "S = pi * r**2\n",
    "rho = 1.2250 # kg/m^3\n",
    "C_resist = C_F * S * rho / 2. # F_resist = - C_resist * ||v|| * v\n",
    "\n",
    "\n",
    "def resistance_free_solution(t, v0x=v0x, v0y=v0y):\n",
    "    ...\n",
    "    \n",
    "\n",
    "def plot_solution(x_nores, y_nores, x=None, y=None):\n",
    "    fig, ax = plt.subplots(figsize=(8, 4))\n",
    "    ax.plot(x_nores, y_nores, label=\"resistance free\")\n",
    "    ax.set_xlabel(\"distance (m)\")\n",
    "    ax.set_ylabel(\"height (m)\")\n",
    "    if x is not None:\n",
    "        ax.plot(x, y, label=\"with resistance\")\n",
    "    ax.legend()\n",
    "    plt.show(fig)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "8617202e12f254480e1fae3258716b685f1a56bcbf234a446366b4fd3345ed22"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}