{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Swept Langmuir Analysis: Floating Potential\n",
    "\n",
    "This notebook covers the use of the [**find_floating_potential()**](../../../api/plasmapy.analysis.swept_langmuir.floating_potential.find_floating_potential.rst#find-floating-potential) function and how it is used to determine the floating potential from a swept Langmuir trace.\n",
    "\n",
    "The floating potential, $V_f$, is defined as the probe bias voltage at which there is no net collected current, $I=0$.  This occurs because the floating potential slows the collected electrons and accelerates the collected ions to a point where the electron- and ion-currents balance each other out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pprint\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "from plasmapy.analysis import swept_langmuir as sla\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = [10.5, 0.56 * 10.5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents:\n",
    "\n",
    "1. [How find_floating_potential() works](#How-find_floating_potential()-works)\n",
    "    1. [Notes about usage](#Notes-about-usage)\n",
    "    1. [Knobs to turn](#Knobs-to-turn)\n",
    "1. [Calculate the Floating Potential](#Calculate-the-Floating-Potential)\n",
    "    1. [Interpreting results](#Interpreting-results)\n",
    "    1. [Plotting results](#Plotting-results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How `find_floating_potential()` works\n",
    "\n",
    "1. The passed current array is scanned for points that equal zero and point-pairs that straddle where the current, $I$, equals zero.  This forms an collection of \"crossing-points.\"\n",
    "1. The crossing-points are then grouped into \"crossing-islands\" based on the `threshold` keyword.\n",
    "    - A new island is formed when a successive crossing-point is more (index) steps away from the previous crossing-point than defined by `threshold`.  For example, if `threshold=4` then an new island is formed if a crossing-point candidate is more than 4 steps away from the previous candidate.\n",
    "    - If multiple crossing-islands are identified, then the function will compare the total span of all crossing-islands to `min_points`.  If the span is greater than `min_points`, then the function is incapable of identifying $V_f$ and will return `numpy.nan` values; otherwise, the span will form one larger crossing-island.\n",
    "1. To calculate the floating potential...\n",
    "    - If the number of points that make up the crossing-island is less than `min_points`, then each side of the \"crossing-island\" is equally padded with the nearest neighbor points until `min_points` is satisfied.\n",
    "    - If `fit_type=\"linear\"`, then a `scipy.stats.linregress` fit is applied to the points that make up the crossing-island.\n",
    "    - If `fit_type=\"exponential\"`, then a `scipy.optimize.curve_fit` fit is applied to the points that make up the crossing-island."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes about usage\n",
    "\n",
    "- The function provides no signal processing.  If needed, the user must smooth, sort, crop, or process the arrays before passing them to the function.\n",
    "- The function requires the voltage array to be monotonically increasing or decreasing.\n",
    "- If the total range spanned by all crossing-islands is less than or equal to `min_points`, then `threshold` is ignored and all crossing-islands are grouped into one island."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Knobs to turn\n",
    "\n",
    "- `fit_type`\n",
    "\n",
    "    There are two types of curves that can be fitted to the identified crossing point data: `\"linear\"` and `\"exponential\"`.  The former will fit a line to the data, whereas, the latter will fit an exponential curve with an offset.  The default curve is `\"exponential\"` since swept Langmuir data is not typically linear as it passes through $I=0$.\n",
    "\n",
    "- `min_points`\n",
    "\n",
    "    This variable specifies the minimum number of points that will be used in the curve fitting.  As mentioned above, the crossing-islands are identified and then padded until `min_points` is satisfied.\n",
    "    \n",
    "    - `min_pints = None` (Default) then the larger of 5 and `factor * array_size` is taken, where `factor = 0.1` for `\"linear\"` and `0.2` for `\"exponential\"`.\n",
    "    - `min_points = numpy.inf` then the entire passed array is fitted.\n",
    "    - `min_points >= 1` then this is the minimum number of points used.\n",
    "    - `0 < min_points < 1` then then the minimum number of points is taken as `min_points * array_size`.\n",
    "\n",
    "\n",
    "- `threshold`\n",
    "\n",
    "    The max allowed index distance between crossing-points before a new crossing-island is formed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate the Floating Potential\n",
    "\n",
    "Below we'll compute the floating potential using the default fitting behavior (`fit_type=\"exponential\"`) and a linear fit (`fit_type=\"linear\"`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "filename = \"Beckers2017_noisy.npy\"\n",
    "filepath = (Path.cwd() / \"..\" / \"..\" / \"langmuir_samples\" / filename).resolve()\n",
    "voltage, current = np.load(filepath)\n",
    "\n",
    "# voltage array needs to be monotonically increasing/decreasing\n",
    "isort = np.argsort(voltage)\n",
    "voltage = voltage[isort]\n",
    "current = current[isort]\n",
    "\n",
    "# get default fit results (exponential fit)\n",
    "results = sla.find_floating_potential(voltage, current, min_points=0.3)\n",
    "\n",
    "# get linear fit results\n",
    "results_lin = sla.find_floating_potential(voltage, current, fit_type=\"linear\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting results\n",
    "\n",
    "The `find_floating_potential()` function returns a six element named tuple, where...\n",
    "\n",
    "- `results[0]` = `results.vf` = the determined floating potential (same units as the pass `voltage` array)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(results[0], results.vf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `results[1]` = `results.vf_err` = the associated uncertainty in the $V_F$ calculation (same units as `results.vf`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(results[1], results.vf_err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `results[2]` = `results.rsq` = the coefficient of determination (r-squared) value of the fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(results[2], results.rsq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `results[3]` = `results.func` = the resulting fitted function\n",
    "\n",
    "    - `results.func` is a callable representation of the fitted function `I = results.func(V)`.\n",
    "    - `resulst.func` is an instance of a sub-class of `AbstractFitFunction`. ([FitFuction classes](../../../api_static/plasmapy.analysis.fit_functions.rst))\n",
    "    - Since `results.func` is a class instance, there are many other attributes available.  For example,\n",
    "        - `results.func.params` is a named tuple of the fitted parameters\n",
    "        - `results.func.param_errors` is a named tuple of the fitted parameter errors\n",
    "        - `results.func.root_solve()` finds the roots of the fitted function. This is how $V_f$ is calculated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    results[3],\n",
    "    results.func,\n",
    "    results.func.params,\n",
    "    results.func.params.a,\n",
    "    results.func.param_errors,\n",
    "    results.func.param_errors.a,\n",
    "    results.func(results.vf),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `results[4]` = `results.islands` = a list of slice objects representing all the indentified crossing-islands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    results[4],\n",
    "    results.islands,\n",
    "    voltage[results.islands[0]],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `results[5]` = `results.indices` = a slice object representing the indices used in the fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    results[5],\n",
    "    results.indices,\n",
    "    voltage[results.indices],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx-thumbnail": {
     "tooltip": "Floating Potential"
    }
   },
   "outputs": [],
   "source": [
    "figwidth, figheight = plt.rcParams[\"figure.figsize\"]\n",
    "figheight = 2.0 * figheight\n",
    "fig, axs = plt.subplots(3, 1, figsize=[figwidth, figheight])\n",
    "\n",
    "# plot original data\n",
    "axs[0].set_xlabel(\"Bias Voltage (V)\", fontsize=12)\n",
    "axs[0].set_ylabel(\"Current (A)\", fontsize=12)\n",
    "\n",
    "axs[0].plot(voltage, current, zorder=10, label=\"Sweep Data\")\n",
    "axs[0].axhline(0.0, color=\"r\", linestyle=\"--\", label=\"I = 0\")\n",
    "axs[0].legend(fontsize=12)\n",
    "\n",
    "# zoom on fit\n",
    "for ii, label, fit in zip([1, 2], [\"Exponential\", \"Linear\"], [results, results_lin]):\n",
    "    # calc island points\n",
    "    isl_pts = np.array([], dtype=np.int64)\n",
    "    for isl in fit.islands:\n",
    "        isl_pts = np.concatenate((isl_pts, np.r_[isl]))\n",
    "\n",
    "    # calc xrange for plot\n",
    "    xlim = [voltage[fit.indices].min(), voltage[fit.indices].max()]\n",
    "    vpad = 0.25 * (xlim[1] - xlim[0])\n",
    "    xlim = [xlim[0] - vpad, xlim[1] + vpad]\n",
    "\n",
    "    # calc data points for fit curve\n",
    "    mask1 = np.where(voltage >= xlim[0], True, False)\n",
    "    mask2 = np.where(voltage <= xlim[1], True, False)\n",
    "    mask = np.logical_and(mask1, mask2)\n",
    "    vfit = np.linspace(xlim[0], xlim[1], 201, endpoint=True)\n",
    "    ifit, ifit_err = fit.func(vfit, reterr=True)\n",
    "\n",
    "    axs[ii].set_xlabel(\"Bias Voltage (V)\", fontsize=12)\n",
    "    axs[ii].set_ylabel(\"Current (A)\", fontsize=12)\n",
    "    axs[ii].set_xlim(xlim)\n",
    "\n",
    "    axs[ii].plot(\n",
    "        voltage[mask], current[mask], marker=\"o\", zorder=10, label=\"Sweep Data\",\n",
    "    )\n",
    "    axs[ii].scatter(\n",
    "        voltage[fit.indices],\n",
    "        current[fit.indices],\n",
    "        linewidth=2,\n",
    "        s=6 ** 2,\n",
    "        facecolors=\"deepskyblue\",\n",
    "        edgecolors=\"deepskyblue\",\n",
    "        zorder=11,\n",
    "        label=\"Points for Fit\",\n",
    "    )\n",
    "    axs[ii].scatter(\n",
    "        voltage[isl_pts],\n",
    "        current[isl_pts],\n",
    "        linewidth=2,\n",
    "        s=8 ** 2,\n",
    "        facecolors=\"deepskyblue\",\n",
    "        edgecolors=\"black\",\n",
    "        zorder=12,\n",
    "        label=\"Island Points\",\n",
    "    )\n",
    "    axs[ii].autoscale(False)\n",
    "    axs[ii].plot(vfit, ifit, color=\"orange\", zorder=13, label=label + \" Fit\")\n",
    "    axs[ii].fill_between(\n",
    "        vfit,\n",
    "        ifit + ifit_err,\n",
    "        ifit - ifit_err,\n",
    "        color=\"orange\",\n",
    "        alpha=0.12,\n",
    "        zorder=0,\n",
    "        label=\"Fit Error\",\n",
    "    )\n",
    "    axs[ii].axhline(0.0, color=\"r\", linestyle=\"--\")\n",
    "    axs[ii].fill_between(\n",
    "        [fit.vf - fit.vf_err, fit.vf + fit.vf_err],\n",
    "        axs[1].get_ylim()[0],\n",
    "        axs[1].get_ylim()[1],\n",
    "        color=\"grey\",\n",
    "        alpha=0.1,\n",
    "    )\n",
    "    axs[ii].axvline(fit.vf, color=\"grey\")\n",
    "    axs[ii].legend(fontsize=12)\n",
    "\n",
    "    # add text\n",
    "    rsq = fit.rsq\n",
    "    txt = f\"$V_f = {fit.vf:.2f} \\\\pm {fit.vf_err:.2f}$ V\\n\"\n",
    "    txt += f\"$r^2 = {rsq:.3f}$\"\n",
    "    txt_loc = [fit.vf, axs[ii].get_ylim()[1]]\n",
    "    txt_loc = axs[ii].transData.transform(txt_loc)\n",
    "    txt_loc = axs[ii].transAxes.inverted().transform(txt_loc)\n",
    "    txt_loc[0] -= 0.02\n",
    "    txt_loc[1] -= 0.26\n",
    "    axs[ii].text(\n",
    "        txt_loc[0],\n",
    "        txt_loc[1],\n",
    "        txt,\n",
    "        fontsize=\"large\",\n",
    "        transform=axs[ii].transAxes,\n",
    "        ha=\"right\",\n",
    "    )"
   ]
  }
 ],
 "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
}