flux_by_heading_interactive_3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import altair as alt\n",
    "import matplotlib.transforms\n",
    "import matplotlib.cm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from ipywidgets import interactive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "code_folding": [],
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def gs_from_wind(windspeed):\n",
    "    #Given the magnitude and sign of the wind in the direction \n",
    "    #of the intended velocity (signed wind_par),\n",
    "    #return the adjusted magnitude of the\n",
    "    #fly's intended velocity\n",
    "    if -0.8<=windspeed<=4.:\n",
    "        return 1.6\n",
    "    elif windspeed<-0.8:\n",
    "        return windspeed+2.4\n",
    "    elif windspeed>4.:\n",
    "        return windspeed-2.4\n",
    "    \n",
    "def signed_wind_par_cartesian(wind_par,original_heading):\n",
    "    #version 1: inputs are in cartesian coordinates\n",
    "    #Given the component of wind parallel to the intended velocity,\n",
    "    #(wind_par vector) and the intended velocity vector,\n",
    "    #return the signed magnitude of the parallel wind \n",
    "    #(+ sign = same direction as intended velocity, - sign = opposite)\n",
    "    #(prep function for gs_from_wind)\n",
    "    sign = np.sign(np.dot(wind_par,original_heading))\n",
    "    return sign*np.sqrt(original_heading.dot(original_heading))\n",
    "\n",
    "def sigmoid(x,x_0,L,y_0,k):\n",
    "    return (x-x_0)+(L/2)+y_0 - L/(np.exp(-k*(x-x_0))+1)\n",
    "\n",
    "def speed_sigmoid_func(x):\n",
    "    x_0a = -0.4\n",
    "    x_0b = 3.6\n",
    "    L = 0.8\n",
    "    k = 4.\n",
    "    y_0 = 1.6\n",
    "    output = np.zeros_like(x)\n",
    "    output[(x>=x_0a)&(x<=x_0b)] = 1.6\n",
    "    output[x<x_0a] = sigmoid(x[x<x_0a],x_0a,L,y_0,k)\n",
    "    output[x>x_0b] = sigmoid(x[x>x_0b],x_0b,L,y_0,k)\n",
    "    return output\n",
    "        \n",
    "def signed_wind_par_polar(theta_0,theta,w):\n",
    "    #Similar to above except start with wind and intended velocity angle/magnitude\n",
    "    #theta_0: wind angle\n",
    "    #theta: intended velocity angle\n",
    "    #w: wind magnitude\n",
    "    return w*np.cos(theta_0-theta)\n",
    "def pp_angle_map_array(theta_0,thetas,r_1s,w):\n",
    "    wind_par_mags = signed_wind_par_polar(theta_0,thetas,w)\n",
    "    thetas_adjusted = np.copy(thetas)\n",
    "#     r_1s[wind_par_mags<-0.8] = wind_par_mags[wind_par_mags<-0.8]+2.4\n",
    "#     r_1s[wind_par_mags>4.] = wind_par_mags[wind_par_mags>4.]-2.4\n",
    "    r_1s = speed_sigmoid_func(wind_par_mags)\n",
    "    sign_change_inds = (r_1s<0.)\n",
    "    thetas_adjusted[sign_change_inds] = (thetas[sign_change_inds]+np.pi)%(2*np.pi)\n",
    "    return (np.arctan(w*np.sin(theta_0-thetas_adjusted)/np.abs(r_1s))+thetas_adjusted)%(2*np.pi)\n",
    "\n",
    "# def pp_angle_map_array_0(theta_0,thetas,r_1s,w):\n",
    "    \n",
    "\n",
    "def pp_angle_map_array_1(theta_0,thetas,r_1s,w):\n",
    "#     signed_wind_par_mags = signed_wind_par_polar(theta_0,thetas,w)\n",
    "    signed_wind_par_mags = w*np.cos(theta_0-thetas)\n",
    "    thetas_adjusted = np.copy(thetas)\n",
    "    adjusted_mag = speed_sigmoid_func(signed_wind_par_mags)\n",
    "  #  sign_change_inds = (adjusted_mag<0.)\n",
    "   # thetas_adjusted[sign_change_inds] = (thetas[sign_change_inds]+np.pi)%(2*np.pi)\n",
    "       \n",
    "#     output = (np.arctan(w*np.sin(theta_0-thetas_adjusted)/np.abs(r_1s))+thetas_adjusted)%(2*np.pi)\n",
    "\n",
    "        \n",
    "    intended_heading_unit_vectors = np.vstack((\n",
    "        np.cos(thetas),np.sin(thetas)))\n",
    "#     intended_heading_unit_vectors =\n",
    "    intended_heading_vectors = adjusted_mag*intended_heading_unit_vectors\n",
    "    w_vec = np.array([w*np.cos(theta_0),w*np.sin(theta_0)])\n",
    "    wind_par_vec = (np.dot(\n",
    "        w_vec,intended_heading_unit_vectors))*intended_heading_unit_vectors\n",
    "    w_perp_vec = w_vec[:,None] - wind_par_vec\n",
    "    heading_final_vec = intended_heading_vectors+w_perp_vec\n",
    "    dispersing_speeds = np.sqrt(np.sum(heading_final_vec**2,axis=0))\n",
    "    output = np.arctan2(heading_final_vec[1],heading_final_vec[0])\n",
    "  \n",
    "    return output\n",
    "    \n",
    "\n",
    "def pp_angle_map(theta_0,theta,w):\n",
    "    wind_par_mag = signed_wind_par_polar(theta_0,theta,w)\n",
    "    r_1 = gs_from_wind(wind_par_mag)\n",
    "    if r_1<0:\n",
    "        theta = (theta+np.pi)%(2*np.pi)\n",
    "    return (np.arctan(w*np.sin(theta_0-theta)/np.abs(r_1))+theta)%(2*np.pi)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(w=1.6,theta_d=1.5):#1.*np.pi/4.):\n",
    "    n = 10000 \n",
    "#     n=10\n",
    "    theta_0 = np.arctan(-1.)\n",
    "    theta_0 = 7*np.pi/8.\n",
    "    thetas = np.linspace(360./n,360,n)*np.pi/180\n",
    "    r_1s = 1.6*np.ones_like(thetas)\n",
    "   \n",
    "    \n",
    "    plt.figure()\n",
    "    time = 5*60.\n",
    "    wind_angle = theta_0\n",
    "    wind_mag = w\n",
    "    signed_wind_par_mags = w*np.cos(wind_angle-thetas)\n",
    "#     print(signed_wind_par_mags) \n",
    "    adjusted_mag = speed_sigmoid_func(signed_wind_par_mags)\n",
    "#     print('-------------')\n",
    "#     print(adjusted_mag)\n",
    "    intended_heading_unit_vectors = np.vstack((\n",
    "        np.cos(thetas),np.sin(thetas)))\n",
    "#     intended_heading_unit_vectors =\n",
    "#     print('-------------')\n",
    "#     print(intended_heading_unit_vectors)\n",
    "    intended_heading_vectors = adjusted_mag*intended_heading_unit_vectors\n",
    "    w_vec = np.array([wind_mag*np.cos(wind_angle),wind_mag*np.sin(wind_angle)])\n",
    "#     print('-------------')\n",
    "#     print(w_vec)\n",
    "    wind_par_vec = (np.dot(\n",
    "        w_vec,intended_heading_unit_vectors))*intended_heading_unit_vectors\n",
    "#     print('-------------')\n",
    "#     print(wind_par_vec)\n",
    "    w_perp_vec = w_vec[:,None] - wind_par_vec\n",
    "    heading_final_vec = intended_heading_vectors+w_perp_vec\n",
    "    dispersing_speeds = np.sqrt(np.sum(heading_final_vec**2,axis=0))\n",
    "    track_heading_angles = np.arctan2(heading_final_vec[1],heading_final_vec[0])\n",
    "\n",
    "#     output = pp_angle_map_array_1(theta_0,thetas,r_1s,w)\n",
    "    output = pp_angle_map_array(theta_0,thetas,r_1s,w)\n",
    "    \n",
    "    mag = time*dispersing_speeds\n",
    "    mag = dispersing_speeds\n",
    "    plt.scatter(mag*np.cos(output),mag*np.sin(output),label='output')\n",
    "   \n",
    "    plt.scatter(mag*np.cos(track_heading_angles),mag*np.sin(track_heading_angles),label='track_heading_angles')\n",
    "    plt.scatter(heading_final_vec[0],heading_final_vec[1],label='heading_final_vec')\n",
    "    plt.axes().set_aspect('equal', 'datalim')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8d3ba932a5bf42ea94ea9c3ed50c1a5b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.6, description='w', max=4.0, step=0.2), FloatSlider(value=1.5, descr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interactive_plot = interactive(f, w=(0, 4.0,0.2),theta_d=(0,2*np.pi,np.pi/8),continuous_update=False)\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout.height = '1000px'\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "be103",
   "language": "python",
   "name": "be103"
  },
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}