PSP_SPAN-I_FOV_diagnostic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SPAN-I Field-of-View Diagnostic\n",
    "Greetings! Below is a guide to qualitatively assessing the field-of-view (FOV) of the PSP SPAN-I instrument, and subsequent interpretation on data reliability. We note that the SPAN-I instrument has many caveats and the data remains preliminary due to ongoing calibration work done by the instrument team. \n",
    "\n",
    "For optimal scientific interpretation of noteworthy event, the reader is strongly encouraged to contact a member of the instrument team:<br>\n",
    "Roberto Livi rlivi@berkeley.edu <br>\n",
    "Ali Rahmati rahmati@berkeley.edu <br>\n",
    "Michael McManus mdmcmanus@berkeley.edu <br>\n",
    "Davin Larson davin@berkeley.edu <br>\n",
    "\n",
    "This tutorial was brought to you by: <br>\n",
    "Jaye Verniero jaye.l.verniero@nasa.gov <br>\n",
    "Kristoff Paulson kristoff.paulson@cfa.harvard.edu <br>\n",
    "who may serve as points of contact for this tutorial. <br>\n",
    "\n",
    "The paper describing the SPAN-I instrument in more detail can be found here: <br>https://doi.org/10.3847/1538-4357/ac93f5\n",
    "\n",
    "Unlike SPC, the SPAN-I sensor is located behind the TPS heat shield and does not have a direct view of the sun line. That being said, when PSP approaches closer to the Sun, the speed of the spacecraft increases and so does SPAN-I's aberration angle. Therefore, there are times when SPAN-I measures the full proton VDF but we have to be a careful. This walkthrough is intended to teach the data user how to investigate the reliability of SPAN-I's Field-of-View (FOV). The three main culprits to check are (1) Phi angle coverage, (2) Density comparison to Quasithermal Noise (QTN), and (3) magnetic variability. There are other anamolies caused by the theta coverage at high energies. The deflection coverage is reduced at energies above 5 keV (which so far, only effects E13).\n",
    "\n",
    "Check here for updated data release notes.\n",
    "http://sweap.cfa.harvard.edu/Data.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#install packages\n",
    "!pip install wget\n",
    "%pip install --update pyspedas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyspedas\n",
    "from pytplot import tplot, store_data, get_data ,tlimit,xlim,ylim,tplot_options,options,split_vec,cdf_to_tplot,divide,tplot_names,get_timespan\n",
    "from pyspedas import time_string, time_double\n",
    "import wget\n",
    "from warnings import simplefilter \n",
    "simplefilter(action='ignore', category=DeprecationWarning)\n",
    "import cdflib\n",
    "\n",
    "#import math functions \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors\n",
    "import dateutil.parser\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we download data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify time range in the form ['yyyy-mm-dd/hh:mm:ss','yyyy-mm-dd/hh:mm:ss']\n",
    "trange=['2020-01-23','2020-02-03']\n",
    "\n",
    "#specify data type to plot\n",
    "datatype='spi_sf00_l3_mom' #protons\n",
    "spi_vars = pyspedas.psp.spi(trange=trange, datatype=datatype, level='l3', time_clip=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.) Phi Angle Coverage\n",
    "The first (and perhaps most dominant) effect on SPAN-I data reliabiliy is how much the plasma is cut off by the heat shield. This is mostly in the phi direction of the instrument. To see this, we plot the energy flux in the three planes of the instrument: Energy, Theta (anode), and Phi (look direction). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prefix='psp_spi_'\n",
    "tplot([prefix+'SUN_DIST',prefix+'EFLUX_VS_ENERGY',prefix+'EFLUX_VS_THETA',prefix+'EFLUX_VS_PHI'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the fourth panel, we see how much energy flux is captured in each $\\phi$ angle bin. To determine \"good\" coverage by eye, one should observe that there exists a peak in the range of angles, meaning the peak (or core) of the distribution is covered in the phi measurement plane. We also want to make sure this peak is located at the 160 degree mark or lower. Note that the 180 bin will always appear to have lower flux, which may be a deceptive peak.\n",
    "\n",
    "\n",
    "Anytime, the collapsed VDFs (E, theta, phi), approach the edges of the coverage, we have partial loss of VDF information. Theta and Energy coverage are also an important factor to consider, especially in later encounters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define variables\n",
    "eflux_phi_data=get_data(prefix+'EFLUX_VS_PHI')\n",
    "times_unix=eflux_phi_data.times\n",
    "eflux = eflux_phi_data.y\n",
    "phi = eflux_phi_data.v\n",
    "\n",
    "#determine phi angle with max eflux\n",
    "max_phi_ind = np.argmax(eflux, axis=1)\n",
    "max_phi = phi[0, max_phi_ind]\n",
    "print(max_phi)\n",
    "\n",
    "#determine average phi\n",
    "avg_phi = np.average(eflux,axis=1)\n",
    "print(avg_phi)\n",
    "print(avg_phi.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = [dateutil.parser.parse(t) for t in time_string(times_unix)]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(12, 5))\n",
    "p = ax.pcolormesh(times[::100], phi[0, :], eflux[::100].T, norm=matplotlib.colors.LogNorm())\n",
    "plt.colorbar(p, ax=ax, label=f'$(cm^2 \\\\ s \\\\ sr \\\\ eV)^{-1}$')\n",
    "ax.plot(times, max_phi, 'k')\n",
    "ax.set(ylim=(140, 170), xlabel='Time', ylabel=f\"$\\\\phi$ [deg]\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define fov array\n",
    "tlen = times_unix.shape[0]\n",
    "phi_fov=np.ones(tlen)\n",
    "\n",
    "#set threshold 163.125 degrees\n",
    "phi_thresh = phi[0,1]\n",
    "\n",
    "#set fov to 0, meaning reliable, if less than 163.125\n",
    "for x in range(0,tlen):\n",
    "    if max_phi[x] < phi_thresh:\n",
    "        phi_fov[x]=0\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have array of zeros and ones indicating if the phi coverage is sufficient. We can mark the insufficient times with red xs below.\n",
    "\n",
    "In reality, this determination is not quantized or binary. This is a 0th order, HIGHLY rudimentary diagnostic. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fig, ax = plt.subplots(figsize=(12, 5))\n",
    "\n",
    "#print(times.shape)\n",
    "start_tind = 66500\n",
    "stop_tind = 67000\n",
    "p = ax.pcolormesh(times[start_tind:stop_tind], phi[0,:], eflux[start_tind:stop_tind].T, norm=matplotlib.colors.LogNorm())\n",
    "plt.colorbar(p, ax=ax, label=f'$(cm^2 \\\\ s \\\\ sr \\\\ eV)^{-1}$')\n",
    "ax.plot(times[start_tind:stop_tind], max_phi[start_tind:stop_tind], 'k')\n",
    "\n",
    "ax.plot(times[start_tind:stop_tind],phi_fov[start_tind:stop_tind]*160,'r+')    \n",
    "        \n",
    "ax.set(ylim=(140, 170), xlabel='Time', ylabel=f\"$\\\\phi$ [deg]\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quasithermal Noise Comparison \n",
    "\n",
    "The Quasithermal Noise (QTN) data from the FIELDS instrument is presently the most reliable measurement of density. See Moncuquet et al. (DOI: 10.3847/1538-4365/ab5a84) for more details.\n",
    "\n",
    " Since SPAN-I measures partial moments of the VDF, density is highly sensitive to inaccuracies. In addition to the $\\phi$ angular coverage assessment described above, the data user is encouraged to compare the SPAN-I measured density to the density derived from QTN. Note that this measurement is also prone to error, and there is currently no agreed upon value of density. Ongoing calibration efforts by the instrument team continue to converge closer to an answer. Below, we show how to compare the QTN and SPAN-I measured density.\n",
    "\n",
    "Note there is not a defined threshold for offsets, as the answer depends on the physical problem at hand. The user is encouraged to reach out to the SPAN-I instrument team for recommendations if the SPAN-I density in the time period of interest deviates strongly from the QTN.\n",
    "\n",
    "Future additions to this tutorial will show how to also compare with density measured by SPC.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load data with pySPEDAS\n",
    "#note this load routine does not currently work\n",
    "datatype = 'sqtn_rfs_v1v2'\n",
    "vars_qtn = pyspedas.psp.fields(trange=trange, datatype=datatype, level='l3', time_clip=True)\n",
    "print(vars_qtn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can instead use cdflib and a handy routine 'cdf_to_tplot', which converts variable in cdf files to tplot variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "year=2020\n",
    "month=1\n",
    "day=29\n",
    "\n",
    "user_datetime = datetime(year,month,day)\n",
    "\n",
    "def yyyymmdd(dt) : return f\"{dt.year:04d}{dt.month:02d}{dt.day:02d}\"\n",
    "\n",
    "#Import from file directory\n",
    "QTNfile_directoryRemote = f'https://spdf.gsfc.nasa.gov/pub/data/psp/fields/l3/sqtn_rfs_v1v2/{user_datetime.year:04d}/'\n",
    "QTNfile_filename = f'psp_fld_l3_sqtn_rfs_v1v2_{yyyymmdd(user_datetime)}_v1.0.cdf'\n",
    "\n",
    "import os.path\n",
    "\n",
    "#check if file is already downloaded. If so, skip download. If not, download in local directory.\n",
    "if os.path.isfile(QTNfile_filename):\n",
    "    print(f\"File already exists in local directory - [{QTNfile_filename}]\")\n",
    "    QTNfile = QTNfile_filename\n",
    "else:\n",
    "    print(\"File doesn't exist. Downloading ...\")\n",
    "    QTNfile = wget.download(QTNfile_directoryRemote + QTNfile_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "QTN_dat = cdflib.CDF(QTNfile)\n",
    "#print variable names in cdf file\n",
    "print(dat._get_varnames())\n",
    "#check variable formats in cdf file\n",
    "print(QTN_dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#convert variables in cdf files to tplot variables for quick plotting\n",
    "cdf_to_tplot(QTNfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify shorter time range in format 'yyyy-mm-dd hh:mm:ss'\n",
    "trange1=['2020-01-29 00:00:00','2020-01-30 00:00:00']\n",
    "#set limits \n",
    "xlim(trange1[0],trange1[1])\n",
    "help,options\n",
    "options(prefix+'DENS','ylog',0)\n",
    "ylim(prefix+'DENS',0,1500)\n",
    "ylim('electron_density',0,1500)\n",
    "tplot(['electron_density',prefix+'DENS'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#overlay both variables on the same axis by defining a new tplot variable\n",
    "store_data(\"dens_compare\",data=['electron_density',prefix+'DENS'])\n",
    "\n",
    "#set legend\n",
    "options('dens_compare', 'color', ['blue', 'red'])\n",
    "options('dens_compare', 'legend_names', ['QTN', 'SPAN-I'])\n",
    "tplot('dens_compare')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 'divide' function divides two tplot variables and is equivalent 'div_data' in IDL. It is useful since it automatically interpolates the two variables if they have different time cadences.\n",
    "\n",
    "We can see what happens when we divide QTN by SPAN-I density and assess the deviation from 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "divide('electron_density',prefix+'DENS',new_tvar='dens_div')\n",
    "ylim('dens_div',0,3)\n",
    "options('dens_div','ytitle','QTN/SPAN-I Dens')\n",
    "tplot('dens_div')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As evident in the above two plots, the SPAN-I density deviates from QTN roughly in the period between 09:00 - 11:00. Thus, if a data user sees this level of discrepency in a time period of interest, extra support should be sought from the instrument team."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Magnetic Field Variability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A third known culprit of questionable data reliability is due to high magnetic field variability. The SPAN-I instrument samples the plasma at a much lower cadence than the FIELDS magnetometer samples the magnetic field. It is therefore important to assess how much the instanteous magnetic field changes with respect to the averaged magnetic field over the lower SPAN-I time cadence. If the magnetic field is changing at a rapid enough timescale, measurements could be skewed and the user is again encouraged to reach out to the instrument team. As before, the SPAN-I team cannot recommend a threshold value for this variability, as it depends on the physical problem. Below, we show how to calculate the angular fluctuations in the magnetic field during the SPAN-I distribution integration time. As described by Eq. 4 in Woodham et al. 2021 (DOI: 10.1051/0004-6361/202039415), this is quantified by:\n",
    "\n",
    "$\\psi_B = \\frac{1}{N}\\Sigma_{i=1}^N arccos(\\mathbf{\\hat{B}}_i \\cdot \\mathbf{\\hat{B}}_{SPAN})$ \n",
    "\n",
    " where $N$ is the number of instantaneous $\\mathbf{\\hat{B}}$ field measurments during one SPAN-I integration time, $\\mathbf{\\hat{B}}_{SPAN}$ is the average magnetic field direction, and $\\mathbf{\\hat{B}}_i$ is the instantaneous magentic field unit vector.\n",
    " \n",
    " Qualitative assessment of this parameter is demonstrated so the user can gain more intuition on the nuances of interpreting partial plasma moments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#specify time range in the form ['yyyy-mm-dd/hh:mm:ss','yyyy-mm-dd/hh:mm:ss']\n",
    "#time_range=['2021-04-29/05:00:00','2021-04-29/12:00:00']\n",
    "\n",
    "#load data\n",
    "#note we are downloading lower time cadence, 4 samples per cycle data\n",
    "#the full cadence is not necessary and would take time to download\n",
    "fields_vars=pyspedas.psp.fields(datatype='mag_sc_4_sa_per_cyc',trange=trange,time_clip=True)\n",
    "print(trange)\n",
    "#plot data vs time\n",
    "tplot('psp_fld_l2_mag_SC_4_Sa_per_Cyc')\n",
    "get_timespan('psp_fld_l2_mag_SC_4_Sa_per_Cyc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "span_b=get_data(prefix+'MAGF_SC')\n",
    "tspan = span_b.times\n",
    "\n",
    "hires_b=get_data('psp_fld_l2_mag_SC_4_Sa_per_Cyc')\n",
    "thr = hires_b.times\n",
    "bhr = hires_b.y\n",
    "\n",
    "import bisect\n",
    "from pyspedas import time_datetime\n",
    "ts_dt=time_datetime(tspan)\n",
    "thr_dt=time_datetime(thr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = np.size(tspan)\n",
    "\n",
    "alpha_mean=np.zeros(N)\n",
    "for i in range(0,N-2):\n",
    "\n",
    "    tr = [ts_dt[i],ts_dt[i+1]]\n",
    "    #find time indices of B in one span accumulation time\n",
    "    tind_start=bisect.bisect_left(thr_dt,tr[0])\n",
    "    tind_stop=bisect.bisect_left(thr_dt,tr[1])\n",
    "\n",
    "    B0 = bhr[tind_start,:]\n",
    "    B_mag = np.sqrt(B0[0]**2 + B0[1]**2 + B0[2]**2)\n",
    "    B0_hat = B0/B_mag\n",
    "\n",
    "    bt_len = tind_stop-tind_start #length new B time array\n",
    "    cos_alpha = np.zeros(bt_len-1)\n",
    "    \n",
    "    \n",
    "    for j in range(1,bt_len-1):\n",
    "        Bi = bhr[j+tind_start,:]\n",
    "        Bi_mag = np.sqrt(Bi[0]**2 + Bi[1]**2 + Bi[2]**2)\n",
    "        cos_alpha[j-1]=np.sum(B0_hat * Bi)/Bi_mag\n",
    "    print(i)\n",
    "    alpha = np.arccos(cos_alpha)*180/np.pi\n",
    "    alpha_mean[i]=np.mean(alpha)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(alpha_mean.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "fig, ax = plt.subplots(figsize=(12, 5))\n",
    "\n",
    "times_tspan = [dateutil.parser.parse(t) for t in time_string(tspan)]\n",
    "\n",
    "ax.plot(times_tspan, alpha_mean, 'k')\n",
    "ax.set( xlabel='Time', ylabel='mean $\\psi$');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.3 ('pyspedas')",
   "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.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "5877fce38de90c1fb6dffa61cadcfee3c47ad77488053f147fda735568ebc1ad"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}