HW2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot\n",
    "from matplotlib.colors import ColorConverter as cc\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Challenge 1.3</h1>\n",
    "<h4>Note: I have used the code without the drag coefficient as the homework said only to use a position dependent gravity without drag</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class particle(object):\n",
    "    \n",
    "    def __init__(self, mass=1., y=0., v=0.):\n",
    "        self.mass = mass\n",
    "        self.y = y\n",
    "        self.v = v\n",
    "        \n",
    "    def euler(self, f, dt):\n",
    "        self.y = self.y + self.v*dt\n",
    "        self.v = self.v + f/self.mass*dt\n",
    "        \n",
    "    def euler_cromer(self, f, dt):\n",
    "        self.v = self.v + f/self.mass*dt\n",
    "        self.y = self.y + self.v*dt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exercise_1_3(y0, v0, vt, dt):\n",
    "    g = 9.8            # g acceleration\n",
    "    mass = 0.01        # mass of the particle\n",
    "    #y0 = 300.          # initial position\n",
    "    #v0 = 0.            # initial velocity\n",
    "    #vt = 30.           # terminal velocity\n",
    "\n",
    "    #dt = 0.5           # time step\n",
    "    \n",
    "    R = 6.37e6\n",
    "    gforce = g*mass    # weight\n",
    "\n",
    "    p = particle(mass, y0, v0)\n",
    "    m = p.mass\n",
    "\n",
    "    y = [y0] # since we do not know the size of the arrays, we define first a python list\n",
    "    v = [v0] # the append method is more efficient for lists than arrays\n",
    "    t = [0.]\n",
    "\n",
    "    while p.y > 0.:\n",
    "        fy = -( g*m )/( 1 + (y[-1]/R) )**2\n",
    "        p.euler(fy, dt)\n",
    "        y.append(p.y)\n",
    "        v.append(p.v)\n",
    "        t.append(t[-1]+dt)\n",
    "\n",
    "    t_data = np.array(t) # we convert the list into a numpy array for plotting\n",
    "    y_data = np.array(y)\n",
    "    v_data = np.array(v)\n",
    "    return t_data, y_data, v_data\n",
    "\n",
    "    #for i in range(0,t_data.size):\n",
    "    #    print i,t_data[i], y_data[i], v_data[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare(y0, v0, vt, dt):\n",
    "    analytical_vf = -np.sqrt( v0**2 + 2*9.8*y0 )\n",
    "    t, y, v = exercise_1_3(y0, v0, vt, dt)\n",
    "    vf = v[-1]\n",
    "    percent_diff = abs ( (analytical_vf - vf)/analytical_vf )*100\n",
    "    return percent_diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.06261987145503226\n",
      "29.288593900491772\n",
      "At y0 =  129941.4876461029 Results are off by  1.0\n"
     ]
    }
   ],
   "source": [
    "y0 = 300.          # initial position\n",
    "v0 = 0.            # initial velocity\n",
    "vt = 30.           # terminal velocity\n",
    "dt = 0.01          # time step, decreased from the example given\n",
    "\n",
    "tol = 1e-4         # When searching for a value within 1%, just look for equivalence within this tolerance\n",
    "\n",
    "# Search for the \n",
    "\n",
    "percent_diff = compare(y0, v0, vt, dt)\n",
    "print(percent_diff)\n",
    "\n",
    "# percent_diff is much lower than 1%, find something drastically higher\n",
    "y0 = 6.37e6       # Using R as the first value\n",
    "percent_diff = compare(y0, v0, vt, dt)\n",
    "print(percent_diff)\n",
    "\n",
    "# 1% difference lies somewhere between these two values\n",
    "upper_y0 = 6.37e6\n",
    "lower_y0 = 300\n",
    "check = True\n",
    "res = ()\n",
    "while check:\n",
    "    y_check = (upper_y0 + lower_y0)/2\n",
    "    percent_diff = compare(y_check, v0, vt, dt)\n",
    "    \n",
    "    if percent_diff > 1:\n",
    "        upper_y0 = y_check\n",
    "    elif percent_diff < 1:\n",
    "        lower_y0 = y_check\n",
    "    \n",
    "    check = not abs(1 - percent_diff) < tol\n",
    "    if not check:\n",
    "        res = (y_check, percent_diff)\n",
    "\n",
    "print(\"At y0 = \", res[0], \"Results are off by \", round(res[1], 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At about y0 = 1.299e5 m, the results differ by 1%. Any higher, and the difference is greater. Constant acceleration is a really good approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Challenge 1.4</h1>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class particle2(object):\n",
    "    \n",
    "    def __init__(self, mass=1., x=0., y=0., vx=0., vy=0.):\n",
    "        self.mass = mass\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        self.vx = vx\n",
    "        self.vy = vy\n",
    "       \n",
    "    def euler(self, fx, fy, dt):\n",
    "        self.vx = self.vx + fx/self.mass*dt\n",
    "        self.vy = self.vy + fy/self.mass*dt\n",
    "        self.x = self.x + self.vx*dt\n",
    "        self.y = self.y + self.vy*dt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Max X (m)')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3iV5f3H8fc3m4SwA4ZlANkbgoBa6hYnSN1K3bTWOmut41dbd1vbqnVV6h5FLWq1St1FQUWIyAYZYc+EyEgCmd/fHzm0lAIJIec8Jzmf13XlCjk5Oc8HLvLJk/u5n/s2d0dERGJHXNABREQkslT8IiIxRsUvIhJjVPwiIjFGxS8iEmMSgg5QE61atfKsrKygY4iI1Ctff/11vrtn7Pl4vSj+rKwscnJygo4hIlKvmNnKvT2uoR4RkRij4hcRiTEqfhGRGKPiFxGJMSp+EZEYo+IXEYkxKn4RkRhTL+bxS/2xo7SCLTtKyd9eSn5hCTvKKigoKmVHaQVllZWUllcSb0Z5pZOUEEdyQhxxZqSnJNA8NYnmaYm0apxMmyYppCTGB/3XEWmQVPxywMorKllVUMyC9dv4dsN2Vm4uZmVBMeu27CBve0mdHadxcgIZ6clktUylU6vGdMpIo3OrNLofkk6rxsl1dhyRWKPil2rlF5YwLXcz36zawpw1W5i7dis7yyoBiI8zMpumcGjLVI7pnkHHFqk0T0uiVeNkMtKTaZQYT/PUJNKS40mIqzrDr3An3oyyykpKyiupqHC27yznu+JSvisuJb+wlI3bdrI59D43v4hpuQXsKKv4d6ZWjZPpmZnOoI7NGd6lJQM6NNNvCCI1ZPVhB67s7GzXkg2RU1hSztQleXy6OJ/pyzezLK8IgJTEOHpmNmFgh+b0zEynZ2YTurZpTHJC+Au3stLZtL2EpZsKWbRhG4s2bGf+um18u2EblQ5JCXEM6tiMYZ1bMrxzSwZ0bBaRXCLRzMy+dvfs/3k8XMVvZs8ApwGb3L3Pbo9fA/wUKAfedfebq3stFX/4fVdUyocLNvLe/A1MXZJPaUUljZMTGNqpBdlZLRjepSV92jYhIT665gNs3VHGjOUFTMvdzLTlm5m/bhvukJwQ9+/fBoZ1bsnAjs1IjLLsIuEWRPGPAAqBF3YVv5kdA9wOnOruJWbW2t03VfdaKv7wKC4t5/35G3h71jqmLMmnvNJp37wRJ/U+hBN7tWHQoc3rXVluLS5j+ooCvly2mWm5m1m4oeoHQYu0JEb2OYTT+mUytFNL4uMs6KgiYRfx4g8dNAt4Z7fifw0Y7+4fHcjrqPjrjrvz9crveHXGaibNXU9RaQXtmjXitP6ZnN6vLb3bNsGs4ZTiluJSvly2mUnzNvDRgo3sKKsgIz2ZU/tmclq/TAZ1bE6cfghIAxUtxT8LeAsYCewEbnL3Gfv42nHAOICOHTsOXrlyr6uLSg2VV1Qyad4Gnp6Sy+w1W0lLiufUfpn8YFB7hmS1iInyKy4t55NFm3hn9nr+9e0mSsorads0hbMGt+fcwzvSrlmjoCOK1KloKf55wCfAdcAQ4FWgs1cTQmf8tVdYUs6Er1bx3BcrWLtlB51apXHZUZ0YM7AdacmxO6mrsKScjxZs5M1v1vLZkjwMOLp7ay44vCNHd8+IumsZIrWxr+KP9Hf+GuCNUNFPN7NKoBWQF+EcDd7OsgpemraSxycvo6ColGGdW3DnGb05tkfrmDi7r07j5ARGD2zH6IHtWF1QzKszVvNqzmqueCGHzKYpnJPdgQuHdaR1ekrQUUXqXKTP+H8MtHX3O8ysG/Ax0FFn/HWntLyS13JW88gnS9i4rYTvdW3FjSd0Y2DH5kFHi3plFZV8vHAjL3+1iilL8kmKj2PMoHZcOaIzXTIaBx1P5IBF/IzfzCYARwOtzGwN8CvgGeCZ0JBPKXBxdaUvNVNZ6bw1ey0PfriEVQXFZB/anIfPG8iwzi2DjlZvJMbHMbJPJiP7ZLI8v4inpuQy8es1vJqzmuN7tuFHIzqTndUi6JgiB003cDUAizdu57Y35pKz8jt6t23CTSd15+huGQ1qdk5Q8gtLeOHLlbzw5Qq2FJcxJKs51x/fjSO6tNS/r0S9QC7u1hUV/97tKK3gT58s4S+f5ZKeksBtp/TkB4Paaww/DIpLy3ltxmr+/GkuG7bt1A8AqRdU/A3M5G838cu35rG6YAdnD27Praf0pEVaUtCxGryS8gpem7Gax/617N8/AG4e2YMhGgKSKKTibyA2bdvJne8s4N056+mSkca9Z/bVOH4Adv0AePRfS9m4rYSTerfh1pN7ktUqLehoIv+m4q/nKiudl6ev4nf/XERJRSXXHHMY477fWQuRBay4tJynpiznz58uo6yikh8Oz+LaY7vSNDUx6GgiKv76bHNhCTdPnMPHizZx1GGtuHt0HzrpzDKqbNq2kz98sJjXvl5Nk5RErjuuK2OHH1rv1jqShkXFX099tjiPn/1tNluLy7jtlB5cfESWLiZGsYXrt3HvuwuZujSfrq0bc9eoPgzvoqE4Cca+il+nI1GqstJ55OMlXPzsdJqnJvLWT4/kkiM7qfSjXM/MJrx4+eE89cNsdpZXcP5fpnHja7PYXFh3O5OJHKzYXawlim0pLuXaV2bx2eI8Rg9oy/1j+tEoSWP59YWZcXyvNhzVtRWPfrKUJz9bxieLNnHryT04e3AHTbeVwOmMP8osXL+NUY99zrRlm7nvzL48eO4AlX49lZIYz00ndWfStd+jW5t0fvH6XC54ahor8ouCjiYxTsUfRd6bt54xj3/BjtIKJowbygVDO2popwHo2iadV8cN4zdj+jJ/3TZOeugzxn+2jPKKyqCjSYxS8UcBd+eJycu46uWZ9MxM551rj2LwobohqCExM847vCMf3fh9RnTL4L5JixjzxBcs3rg96GgSg1T8ASuvqOS2N+fy2/cWcVq/tvz1ymFaCrgBa9MkhfFjB/PI+QNZ890OTntkKuM/W0ZFZfTPrpOGQ8UfoKKSci5/PocJ01dz9TFd+NN5A0hJ1Hh+Q2dmnN6/LR/cMIKjQ2f/543/ktUFxUFHkxih4g/Itp1lXPzMdKYsyeP+MX35+Uk9NJ4fY1o1TubJsYP54zn9WbR+Oyc/PIU3v1kTdCyJASr+AGzavpNzn5zGrNVbePSCQZx/eMegI0lAzIwxg9oz6brv0TMznRtenc11r3zD9p1lQUeTBkzFH2Gbtu3kvPHTWLm5iGcuGcIpfTODjiRRoEOLVF4ZN5wbT+jGO3PWc+qfpjJr9ZagY0kDpeKPoDXfFXPOk1+yYetOnr/scEZ0ywg6kkSR+Djj2uO68uq4YVRUOmc98QVPT11OfVhWReoXFX+ErN2yg3OfnEZBUSkvXj5U67fLPmVntWDStd/j2B6tufudBVz915ka+pE6peKPgA1bd3LhX6axbWcZf71yGIMP1cbnsn9NUxN5cuxgbjulB+/P38gZj37Oog3bgo4lDYSKP8w2F5Zw4VPTyC8s5fnLDqdPu6ZBR5J6wswYN6ILE64cRlFJOWc+9gVvzVobdCxpAFT8YbR1RxkXPzudNd/t4OmLsxnUUWf6cuAO79SCd649ij7tmnDdK7O46x8LKNNyD3IQVPxhsrOsgitfyOHbDdv580WDGartEeUgtE5P4a9XDuPSI7N45vPlXPTUV1rqWWpNxR8GlZXOz/42m+nLC/j92f05pkfroCNJA5AYH8evTu/Ng+f2Z9bqLZzx6OcsWKdxfzlwKv465u7cFdoM/daTezBqQLugI0kDc+bA9vztx8OpqHR+8MQXfLhgY9CRpJ5R8dexF6et5LkvVnDFUZ0YN6Jz0HGkgerXvhlvX3Mk3Q5J50cv5vDStJVBR5J6JGzFb2bPmNkmM5u3l8/dZGZuZq3CdfwgTFmSx53/WMDxPVtz6yk9tfaOhFXr9BQmXDmUY7q35v/+Po8H3l+km72kRsJ5xv8cMHLPB82sA3ACsCqMx464FflFXP3yTA7LaMxD5w0kXtvrSQSkJiXw5NjBXDC0I4/9axk/e202peWa8SP7F7bid/fPgIK9fOpB4GagwZyaFJWUM+7FHOLijKcuzqZxsrYylshJiI/j3tF9uOnEbrzxzVoue26G7vSV/YroGL+ZnQGsdffZNXjuODPLMbOcvLy8CKSrHXfn5tfnsHRTIY+eP4gOLVKDjiQxyMz46bFd+f3Z/ZmWu5lznpzGxm07g44lUSpixW9mqcDtwB01eb67j3f3bHfPzsiI3sXMxn+Wy7tz1nPzyB4c1bVBXbKQeuiswe155pIhrNpcxJmPfc4Sbe0oexHJM/4uQCdgtpmtANoDM83skAhmqFNf5W7md+9/yyl9D+FHmsEjUWJEtwxe/dFwykLTPb/K3Rx0JIkyESt+d5/r7q3dPcvds4A1wCB33xCpDHUpb3sJ10z4hkNbpPLbH/TTDB6JKn3aNeWNq44gIz2ZsU9P590564OOJFEknNM5JwBfAt3NbI2ZXR6uY0VaRaVz3SvfsG1nGY9fNIj0lMSgI4n8jw4tUnn9qiPo174pP50wk6em5AYdSaJE2KafuPv51Xw+K1zHDrc/fPAtXyzbzANn9aPHIU2CjiOyT81Sk3jpiqHc8Oos7nl3Ieu27OT/Tu1JnKYbxzTduXuAPlqwkccnL+P8wztwdnaHoOOIVCslMZ5HLxjEJUdULfB2zSvfaK5/jNOE8wOQX1jCL16fQ8/MJvz6jN5BxxGpsfg441en9yKzaQr3/3MRxSXlPHHRYFIS44OOJgHQGX8NuTu3vD6X7SXlPHTuAJIT9A0j9YuZ8aPvd+G+M/syeXEelzw7ncKS8qBjSQBU/DX06ozVfLRwIzef1J3uh6QHHUek1i4Y2pEHzxnAjBXfcckzKv9YpOKvgZWbi7jrnQUc0aUllx3ZKeg4Igdt9MB2PHL+QL5ZvYVLn51Okco/pqj4q1FZ6fx84hzi44zfn91fsyGkwTilbyZ/Om8gM1dt4dJnZ6j8Y4iKvxoTZ65h+vICbj+lJ22bNQo6jkidOrVfJg+dO4CclQVc9twMiktV/rFAxb8fBUWl3D9pIdmHNuccTd2UBur0/m158NwBzFhRVf47SiuCjiRhpuLfj/smLWT7znLuG9NXQzzSoI0a0I4Hzx3A9OUFXP68yr+hU/Hvw7TczUz8eg1XjuhMtzaaxSMN36gB7fjDOf35MnczV7wwg51lKv+GSsW/FyXlFdz+5lw6tGjEtcd2DTqOSMScObA9vz+rP18s28yVL+So/BsoFf9ejP80l2V5Rdw1qg+NknSjlsSWHwxuzwNn9Wfq0nyVfwOl4t/DivwiHvnXUk7tm8kx3VsHHUckEGcNbs9vf9CPqUvz+dGLX6v8GxgV/27cnV++NY/k+DjuOL1X0HFEAnVOdgd+M6Yvny7O49oJ31BeoYXdGgoV/27enr2OKUvyuemk7rRpkhJ0HJHAnTukI78+vRcfLNjIbW/Oxd2DjiR1QKtzhmzdUcbd7yykX/umXDTs0KDjiESNS47sREFRKX/6ZCnN05K49eSeQUeSg6TiD3nww8UUFJXw3KVDiNecfZH/csMJ3SgoLuXJT3NpmZbEuBFdgo4kB0HFT9UF3ZemreTcIR3p065p0HFEoo6ZcecZfdhSXMZ9kxbRPDVJGxHVYyp+qu7QTUqI44YTNGdfZF/i44w/njOArTvKuOWNuTRLTeKEXm2CjiW1EPMXd79ctpkPFmzkJ0d3oXW6LuiK7E9SQhx/vmgwfdo15eq/zuSr3M1BR5JaiOnir6x07p20gHbNGnHF9zoHHUekXkhLTuDZS4bQoXkjrng+h/nrtgYdSQ5QTBf/P+dtYN7abfzsxG7ae1TkALRIS+LFy4eSnpLAxc9MZ3VBcdCR5ADEbPGXVVTywPuL6N4mnVED2gUdR6TeadusES9cPpSyCueSZ6eztbgs6EhSQzFb/BOmr2LF5mJuObmHpm+K1NJhrRszfuxgVhfsYNyLOZSUa2mH+iBsxW9mz5jZJjObt9tjD5jZIjObY2ZvmlmzcB1/fwpLynn4oyUM69yCo7tnBBFBpMEY2rklD5zdj6+WF/CLiXN0d289EM4z/ueAkXs89iHQx937AYuBW8N4/H168cuVbC4q5Rcje2Cms32RgzVqQDt+flJ3/j5rHX/8cHHQcaQaYSt+d/8MKNjjsQ/cfdemntOA9uE6/r4Ul5bz1JRcRnTLYGDH5pE+vEiD9ZOju3DekA488slSXpuxOug4sh9BjvFfBvxzX580s3FmlmNmOXl5eXV20F1n+9cee1idvaaIVN3de/foPnyvaytufXMuny2uu+9bqVuBFL+Z3Q6UAy/v6znuPt7ds909OyOjbsbht+8s48+fLuP73TLIzmpRJ68pIv+RGB/H4xcOomvrxvzk5ZksXL8t6EiyFzUqfjOLM7OBZnaqmR1rZrW+T9vMLgZOAy70CF8FeuHLlXxXXMaNJ3SL5GFFYkp6SiLPXjqEtOR4Ln12Bhu37Qw6kuxhv8VvZl3MbDywFPgNcD7wE+BDM5tmZpeaWY1/azCzkcAvgDPcPaJ3fOwsq+DZz5czolsG/TsEMplIJGZkNm3EM5cMYdvOMsZp+8aoU11p3wO8BHRx95Pc/SJ3Pys0K+cMoCkwdm9faGYTgC+B7ma2xswuBx4F0qn6wTHLzP5cZ3+Takz8eg35haX8+PtamkEkEnq3bcpD5w5gztqt3KxpnlFlv6tzuvv5+/ncJuChA/zap2sere6UV1Tylym59G/flOGdWwYRQSQmndj7EG46sTsPvP8t3Q9J5+pjNKkiGtRoWWYziwdOBbJ2/xp3/2N4YtWtt2evY+XmYm4bO1jz9kUi7CdHd+HbDdt54P1v6dq6MSf2PiToSDGvpuPz/wAuAVpSNVSz6y3qVVQ6j36ylJ6ZTThRa4eLRJyZ8buz+tG/fVOuf3WWZvpEgZpuxNI+NK5f73y4YCO5+UU8esFAne2LBCQlMZ7xP8zmjEencsXzObz90yNp2Tg56Fgxq6Zn/P80sxPDmiRMnp6aS/vmjRipXy9FAtWmSQrjx2aTX1jCVS/NpLS8MuhIMaumxT8NeNPMdpjZNjPbbmZR//va7NVbmLHiOy49shMJ8TG7EKlI1OjfoRm/O6sf01cUcM+7C4KOE7NqOtTzB2A4MDfSN10djOe/WEF6cgLnZEd8SSAR2YdRA9oxd81Wnpq6nMGHNtd+GAGoafEvAebVp9IH+OVpvRgzqD3pKYlBRxGR3fzi5B7MXrOFW9+YS6/MJnRtUy/mijQYNR3/WA9MNrNbzezGXW/hDFYXmqclcVTXVkHHEJE9JMbH8egFg0hNiueql2dSVFJe/RdJnalp8S8HPgaSqGfTOUUkOrVpksKfzhtIbl4ht74xV3f2RlCNhnrc/c5wBxGR2HPEYa34WejO3uys5vxweFbQkWJCdYu0jTezvvv4XJqZXWZmF4YnmojEgqu+34Vje7Tm7ncWMGv1lqDjxITqhnoeB35pZgvN7G9m9nhoL90pwBdUDfdMDHtKEWmw4uKMP57Tn9bpKVz98ky+KyoNOlKDt9/id/dZ7n4OMAR4DJgCvA1c4e793f1hdy+JQE4RacCapSbxxEWDyNtewvWvzqKyUuP94VSji7vuXujuk919grv/3d2/DXcwEYkt/do3447Te/Hp4jwe/dfSoOM0aLqdVUSixoVDOzJ6QFse/GgxU5fkBx2nwVLxi0jUMDPuG9OXrq0bc+0r37B+646gIzVI1c3q2efnzUz7F4pInUtNSuCJiwZTUlbB1S/PpKxCi7nVterO+HPMbOieD5rZFcDM8EQSkVjXJaMxvz2rHzNXbeH+SYuCjtPgVFf81wLjzewvZtbCzAaa2ZfAScCI8McTkVh1Wr+2XHJEFs98vpxJc9cHHadBqW4651RgELARWEbVVM5fufvZ7r4mAvlEJIbddkpPBnZsxs0T55CbVxh0nAajJhd3zwbOB56garG2c82sRVhTiYgASQlxPHbBIBLjjatemsmO0oqgIzUI1V3c/Qi4EDje3W8DhgKzgBlmNi4C+UQkxrVt1oiHzxvI4k3buf3vWsytLlR3xv+Yu5/u7ssBvMojwJHA98OeTkQEGNEtg+uO68obM9fyyozVQcep96ob439zH49vcHctziYiEXPtsV0Z0S2DX709n/nrtgYdp14L2w1cocXcNpnZvN0ea2FmH5rZktD75uE6vog0LHFxxkPnDqB5aiLXTviG4lJt3lJb4bxz9zlg5B6P3QJ87O5dqdrY5ZYwHl9EGpgWaUk8eM4AcvOLuPsdbdZeWzUqfjNL2ctj+93T0N0/Awr2eHgU8Hzoz88Do2tyfBGRXY44rBU//n4XJkxfzT81v79WanrGP8PMhu36wMx+QNV6/AeqjbuvBwi9b72vJ5rZODPLMbOcvLy8WhxKRBqqG0/oRv/2Tbnljbms26L1fA5UTYv/AuARM3vAzF4GrgSODV8scPfx7p7t7tkZGRnhPJSI1DOJ8XE8fN5Ayisquf7VWVRo/f4DUtP1+OcC9wI/Bo4BflrLO3c3mlkmQOj9plq8hogIWa3SuGtUH6YvL+Bxrd9/QGo6xv80cD3QD7gU+IeZXV2L470NXBz688XAW7V4DRERAMYMaseoAW156OMlfL3yu6Dj1Bs1HeqZBxzj7svd/X1gGFVr+OyTmU0AvgS6m9kaM7sc+A1wgpktAU4IfSwiUitmxt2j+5DZNIUbX5tFUYmmeNaE1Yfbn7Ozsz0nJyfoGCISpb7K3cx5f5nGhUM7cs/ovkHHiRpm9rW7Z+/5eE2Herqa2UQzW2Bmubve6j6miMiBG9q5JVcc1YmXpq3i08WaBVidmg71PEvV6pzlVF3cfQF4MVyhREQO1M9O7E7X1o25eeJsthaXBR0nqtW0+Bu5+8dUDQ2tdPdfE+bpnCIiByIlMZ4/njOAzYWl3PH2vOq/IIbVtPh3hvbfXWJmPzWzM9nPzVciIkHo274p1xzblbdmrePdObqrd19qWvzXA6lUbcU4GBjLf6ZliohEjZ8c04X+7Zvyf3+fy6btO4OOE5VqegPXDHcvdPc17n6pu49x92nhDicicqAS4+P4wzkDKC6t4NbXtXHL3iTs75Nm9vb+Pu/uZ9RtHBGRg3dY68bcPLIHd7+zgNdyVnPukI5BR4oq+y1+YDiwGpgAfAVY2BOJiNSBS4/I4sMFG7jrHws4oksrOrRIDTpS1KhuqOcQ4DagD/AwVXfb5rv7p+7+abjDiYjUVlyc8fuz+2Nm3PS32VRqIbd/q27rxQp3f8/dL6ZqmYalwGQzuyYi6UREDkL75qnccXovvlpewDOfLw86TtSo9uKumSWb2RjgJeBq4E/AG+EOJiJSF84e3J7je7bmd+9/y5KN24OOExX2W/xm9jxVG64MAu509yHufre7r41IOhGRg2Rm3D+mH42TE7jxtdmUVVQGHSlw1Z3xjwW6AdcBX5jZttDbdjPbFv54IiIHLyM9mXtH92Hu2q08prX7qx3jj3P39NBbk93e0t29SaRCiogcrJP7ZnLmwHY88slS5qzZEnScQNX0zl0RkXrv12f0JqNxMjf9bTYl5RVBxwmMil9EYkbTRoncP6YvizcW8tgnsTvko+IXkZhyTI/WjBnYjscnL2PButi8VKniF5GYc8fpvWiWmsTPJ8bmLB8Vv4jEnGapSdwzujfz121j/Gext5mgil9EYtLIPpmc2jeThz9awtJNsXVjl4pfRGLWr8/oTWpyPD+fOIeKGFrLR8UvIjErIz2ZX5/em29WbeHZGFrLR8UvIjFt1IC2HNujNX/8cDFrvisOOk5EqPhFJKaZGXeN6o073PHW/JjYsSuQ4jezG8xsvpnNM7MJZpYSRA4REahavvlnJ3bjk0WbmDR3Q9Bxwi7ixW9m7ajatD3b3fsA8cB5kc4hIrK7S47Iok+7Jvz6H/PZuqMs6DhhFdRQTwLQyMwSgFRgXUA5REQASIiP4zdj+rG5sITfvrco6DhhFfHiD63l/3tgFbAe2OruH0Q6h4jInvq0a8plR3bir1+tYsaKgqDjhE0QQz3NgVFAJ6AtkGZmF+3leePMLMfMcvLy8iIdU0Ri1A0ndKNds0bc+sbcBruCZxBDPccDy909z93LqNrG8Yg9n+Tu4909292zMzIyIh5SRGJTWnIC94zuw9JNhTz5acNcziGI4l8FDDOzVDMz4DhgYQA5RET26pgerTmtXyaPfrKUZXmFQcepc0GM8X8FTARmAnNDGcZHOoeIyP7ccXovUhLjuP3NuQ1ubn8gs3rc/Vfu3sPd+7j7WHcvCSKHiMi+tE5P4ZaTezItt4C/fb0m6Dh1Snfuiojsw3lDOjAkqzn3vruQ/MKGc36q4hcR2Ye4OOP+MX0pLi3nnncWBB2nzqj4RUT247DW6Vx19GH8fdY6Pl3cMKaWq/hFRKrxk6O70DkjjV/+fR47y+r/3H4Vv4hINVIS47lnVB9WFRTz+ORlQcc5aCp+EZEaOOKwVowa0JY/T17G8vyioOMcFBW/iEgN3X5KT5IT4rjjrXn1em6/il9EpIZaN0nhZyd2Y8qS/Hq9br+KX0TkAIwdnkXvtk246535bN9ZP9ftV/GLiByA+DjjntF92LS9hIc+WhJ0nFpR8YuIHKCBHZtz/uEdee6LFSxcvy3oOAdMxS8iUgs3n9SdJikJ/KoebtCu4hcRqYVmqUncPLIH01cU8Pbs+rV7rIpfRKSWzsnuQN92Tblv0kKKSsqDjlNjKn4RkVqKjzPuHNWbjdtKeOSTpUHHqTEVv4jIQRjUsTlnDW7P01Nz681uXSp+EZGD9IuRPUhJiOfOfyyoFxd6VfwiIgcpIz2Z60/oxmeL8/hwwcag41RLxS8iUgd+OPxQurZuzN3vLoj6pZtV/CIidSAxPo47z+jN6oIdPPlpbtBx9kvFLyJSR444rBWn9s3k8clLWV1QHHScfVLxi4jUodtO7UmcGfe+uzDoKPuk4hcRqUPtmjXi6mO68N78DUxZEp179Kr4RUTq2BXf68yhLVP59dvzKS2vDDrO/1Dxi4jUsZTEeO44rRfL8op44adBVzsAAAk9SURBVMsVQcf5H4EUv5k1M7OJZrbIzBaa2fAgcoiIhMuxPVozolsGD3+8hIKi0qDj/JegzvgfBt5z9x5AfyB6r4KIiNSCmfHLU3tSXFrBgx8uDjrOf4l48ZtZE2AE8DSAu5e6+5ZI5xARCbeubdK5aGhHXv5qJd9u2B50nH8L4oy/M5AHPGtm35jZU2aWtueTzGycmeWYWU5eXnReGRcRqc71x3ejcXIC97wbPev4BFH8CcAg4Al3HwgUAbfs+SR3H+/u2e6enZGREemMIiJ1onlaEtcf340pS/L5ZNGmoOMAwRT/GmCNu38V+ngiVT8IREQapLHDD6VzRhr3vrswKqZ3Rrz43X0DsNrMuoceOg5YEOkcIiKRkhgfxy9P7UVufnRM7wxqVs81wMtmNgcYANwXUA4RkYg4untG1EzvDKT43X1WaPy+n7uPdvfvgsghIhIpZsb/Rcn0Tt25KyISId3apHNhFEzvVPGLiETQDVEwvVPFLyISQdEwvVPFLyISYbtP7yyriPz0ThW/iEiEJcbHcfspPcnNL2LC9FURP76KX0QkAMf2aM2wzi146KMlbNtZFtFjq/hFRAJgZtx+Si8Kikp5YvKyiB5bxS8iEpC+7Zty5sB2PD11OWu37IjYcVX8IiIBuumkqtVr/vD+txE7popfRCRA7Zo14vKjOvHGN2uZt3ZrRI6p4hcRCdhVR3ehRVpSxG7qUvGLiASsSUoi1x3XlWm5BRG5qUvFLyISBS4Y2pFOrdK4/5+LKA/zTV0qfhGRKJAYH8ctJ/dg6aZCXs1ZHdZjqfhFRKLEib3acHhWCx78cDGFJeVhO46KX0QkSpgZt53ak/zCUsZ/Gr6bulT8IiJRZECHZpzevy3jp+SyYevOsBxDxS8iEmVuPqk7lZXwhw/Cc1OXil9EJMp0aJHKD4cfyusz14Rlp66EOn9FERE5aFcfcxjfbtwelvX6VfwiIlGoeVoSL14+NCyvraEeEZEYo+IXEYkxKn4RkRgTWPGbWbyZfWNm7wSVQUQkFgV5xn8dsDDA44uIxKRAit/M2gOnAk8FcXwRkVgW1Bn/Q8DNwD4nqJrZODPLMbOcvLy8yCUTEWngIl78ZnYasMndv97f89x9vLtnu3t2RkZGhNKJiDR8Foltvv7rgGb3A2OBciAFaAK84e4X7edr8oCVtTxkKyC/ll8bTsp1YKI1F0RvNuU6MA0x16Hu/j9nzhEv/v86uNnRwE3ufloYj5Hj7tnhev3aUq4DE625IHqzKdeBiaVcmscvIhJjAl2rx90nA5ODzCAiEmti4Yx/fNAB9kG5Dky05oLozaZcByZmcgU6xi8iIpEXC2f8IiKyGxW/iEiMabDFb2bPmNkmM5sXdJbdmVkHM/uXmS00s/lmdl3QmQDMLMXMppvZ7FCuO4POtLtoXNTPzFaY2Vwzm2VmOUHn2cXMmpnZRDNbFPp/NjwKMnUP/TvtettmZtcHnQvAzG4I/Z+fZ2YTzCwl6EwAZnZdKNP8uv63arBj/GY2AigEXnD3PkHn2cXMMoFMd59pZunA18Bod18QcC4D0ty90MwSganAde4+Lchcu5jZjUA20CSc930cCDNbAWS7e1Td9GNmzwNT3P0pM0sCUt19S9C5djGzeGAtMNTda3tjZl1laUfV//Ve7r7DzF4DJrn7cwHn6gO8AhwOlALvAVe5+5K6eP0Ge8bv7p8BBUHn2JO7r3f3maE/b6dqhdJ2waYCr1IY+jAx9BYVZwVa1K/mzKwJMAJ4GsDdS6Op9EOOA5YFXfq7SQAamVkCkAqsCzgPQE9gmrsXu3s58ClwZl29eIMt/vrAzLKAgcBXwSapEhpOmQVsAj5096jIRQ0W9QuIAx+Y2ddmNi7oMCGdgTzg2dDQ2FNmlhZ0qD2cB0wIOgSAu68Ffg+sAtYDW939g2BTATAPGGFmLc0sFTgF6FBXL67iD4iZNQZeB653921B5wFw9wp3HwC0Bw4P/boZqJou6heQI919EHAycHVoeDFoCcAg4Al3HwgUAbcEG+k/QkNPZwB/CzoLgJk1B0YBnYC2QJqZ7XPdsEhx94XAb4EPqRrmmU3V+mZ1QsUfgNAY+uvAy+7+RtB59hQaGpgMjAw4CsCRwBmh8fRXgGPN7KVgI1Vx93Wh95uAN6kajw3aGmDNbr+tTaTqB0G0OBmY6e4bgw4Scjyw3N3z3L0MeAM4IuBMALj70+4+yN1HUDVsXSfj+6Dij7jQRdSngYXu/seg8+xiZhlm1iz050ZUfUMsCjYVuPut7t7e3bOoGiL4ZH8ruUaKmaWFLs4TGko5kapfzwPl7huA1WbWPfTQcUCgEwf2cD5RMswTsgoYZmapoe/N44iSnQHNrHXofUdgDHX47xboWj3hZGYTgKOBVma2BviVuz8dbCqg6gx2LDA3NJ4OcJu7TwowE0Am8HxoxkUc8Jq7R83UySjUBnizqitIAP7q7u8FG+nfrgFeDg2r5AKXBpwHgNBY9QnAj4LOsou7f2VmE4GZVA2lfEP0LN3wupm1BMqAq939u7p64QY7nVNERPZOQz0iIjFGxS8iEmNU/CIiMUbFLyISY1T8IiIxRsUvshszO9PM3Mx6HMRrXGJmj9ZlLpG6pOIX+W/nU7Va43lBBxEJFxW/SEho/aQjgcsJFb+ZHW1mk3db3/7l0B2emNkpocemmtmf9rZXQOiO6NfNbEbo7ciI/qVE9qLB3rkrUgujgffcfbGZFZjZrjVuBgK9qVqu93PgyNDGK08CI9x9eehO8b15GHjQ3aeGbr1/n6old0UCo+IX+Y/zqVoCGqoWhDsfeBeY7u5rAELLbGRRtclPrrsvDz1/ArC3pZmPB3qFfkkAaGJm6aG9GEQCoeIXAUJrohwL9DEzB+KpWm9/ElCy21MrqPq+sf95kb2LA4a7+446jCtyUDTGL1LlLKq26TzU3bPcvQOwHDhqH89fBHQObaYDcO4+nvcB8NNdH5jZgLqJK1J7Kn6RKudTtab+7l4HLtjbk0Nn8D8B3jOzqcBGYOtennotkG1mc8xsAfDjuossUjtanVOklsyscWhzegMeA5a4+4NB5xKpjs74RWrvytDF3vlAU6pm+YhEPZ3xi4jEGJ3xi4jEGBW/iEiMUfGLiMQYFb+ISIxR8YuIxJj/B7+cZNsWzFMYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "g = 9.8            # g acceleration\n",
    "v0 = 30.           # initial velocity\n",
    "\n",
    "dt = 0.001           # time step\n",
    "\n",
    "colors = ['red','orange','yellow','green','magenta','cyan','blue','purple','black']\n",
    "\n",
    "highest_y = []\n",
    "highest_x = []\n",
    "angles = np.arange( 1, 9, 0.01 )\n",
    "\n",
    "for angle in angles:\n",
    "    x = [0]                                  # we need to initialize the arrays for each value of the angle\n",
    "    y = [0]\n",
    "    vx = [math.cos(angle*0.1*math.pi/2.)*v0] \n",
    "    vy = [math.sin(angle*0.1*math.pi/2.)*v0] \n",
    "    t = [0.]\n",
    "    \n",
    "    p = particle2(1., 0., 0., vx[0], vy[0])\n",
    "    m = p.mass\n",
    "    k = 0.1*m                             # k/m = 0.1 so k = 0.1*m\n",
    "    while p.y >= 0.:\n",
    "        curr_vx = p.vx\n",
    "        curr_vy = p.vy\n",
    "        curr_v = np.sqrt( curr_vx**2 + curr_vy**2 )\n",
    "        fx = -k*curr_v*curr_vx\n",
    "        fy = -m*g - k*curr_v*curr_vy\n",
    "        p.euler(fx, fy, dt)\n",
    "        x.append(p.x)\n",
    "        y.append(p.y)\n",
    "        vx.append(p.vx)\n",
    "        vy.append(p.vy)\n",
    "        t.append(t[-1]+dt)\n",
    "    \n",
    "    t_data = np.array(t) # we convert the list into a numpy array for plotting\n",
    "    x_data = np.array(x)\n",
    "    y_data = np.array(y)\n",
    "    vx_data = np.array(vx)\n",
    "    vy_data = np.array(vy)\n",
    "    highest_x.append(max(x))\n",
    "    highest_y.append(max(y))\n",
    "\n",
    "    #my_plot = pyplot.plot(x_data, y_data, color=(colors[angle]), ls='-', lw=3, label = str(angle*0.1))\n",
    "    #pyplot.legend()\n",
    "\n",
    "pyplot.plot(angles, highest_x)\n",
    "pyplot.xlabel(\"Angle\")\n",
    "pyplot.ylabel(\"Max X (m)\")\n",
    "#pyplot.xlabel('position x(m)')\n",
    "#pyplot.ylabel('position y(m)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum range of:  16.373520312554408\n",
      "At an angle of:  0.367 *(Pi/2) radians\n"
     ]
    }
   ],
   "source": [
    "ind = 0\n",
    "for i in range(len(highest_x)):\n",
    "    if highest_x[i] == max(highest_x):\n",
    "        ind = i\n",
    "        break\n",
    "\n",
    "print(\"Maximum range of: \", highest_x[ind])\n",
    "print(\"At an angle of: \", round(angles[ind]*0.1, 4), \"*(Pi/2) radians\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}