update

docs/notebooks/formulary/helio/collisional_analysis.ipynb

@@ -33,6 +33,7 @@
  "import random\n",
  "\n",
  "from astropy import units as u\n",
+ "from hapiclient import hapi\n",
  "\n",
  "from plasmapy.formulary.collisions.helio import collisional_analysis as coal\n",
  "from plasmapy.particles import Particle\n",
@@ -144,20 +145,21 @@
  " params[\"n_1\"][0] = [random.randint(n_p_min, n_p_max)] * u.cm**-3\n",
  " params[\"n_2\"][0] = [random.randint(n_a_min, n_a_max)] * u.cm**-3\n",
  " params[\"v_1\"][0] = [random.randint(v_p_min, v_p_max)] * u.km / u.s\n",
- " params[\"T_1\"][0] = [(1 + random.random()) * T_p] * u.K\n",
- " params[\"T_2\"][0] = [(1 + random.random()) * T_a] * u.K\n",
+ " params[\"T_1\"][0] = [random.random() * T_p] * u.K\n",
+ " params[\"T_2\"][0] = [random.random() * T_a] * u.K\n",
  "\n",
  " return params\n",
  "\n",
  "while params[\"T_2\"][0] / params[\"T_1\"][0] > 15:\n",
  " rand_asgn(params)\n",
  "\n",
+ "\n",
  "# Calculate theta for the given parameters\n",
  "for param in (\"n_1\", \"n_2\", \"v_1\", \"T_1\", \"T_2\"): # Vary certain parameters so we can see their individual affect\n",
  " theta[param] = {}\n",
- " for mm_val in (0.1, 1.5):\n",
+ " for mm_val in (0.8, 10):\n",
  " theta[param][mm_val] = []\n",
- " params[str(param)] = params[str(param)] * mm_val\n",
+ " params[str(param)][0] = params[str(param)][0] * mm_val\n",
  " for radius in radii:\n",
  " params[\"r_n\"] = [radius] * u.au\n",
  " theta[param][mm_val].append(float(coal.temp_ratio(**params)[0]))\n",
@@ -188,7 +190,7 @@
  " fontsize=fs,\n",
  ")\n",
  "\n",
- "plt.ylim(0, 10)\n",
+ "#plt.ylim(0, 15)\n",
  "plt.xlim(0.1, 1)\n",
  "\n",
  "plt.tick_params(\n",
@@ -206,7 +208,7 @@
  "plt.figtext(0.1, -0.08, txt, ha=\"left\", fontsize=fs)\n",
  "plt.legend(loc=\"upper right\", fontsize=fs)\n",
  "plt.grid()\n",
- "plt.show()"
+ "plt.show()\n"
  ]
  },
  {
@@ -259,17 +261,15 @@
  "outputs": [],
  "source": [
  "# Get Dst index from CDAWeb HAPI server\n",
- "from hapiclient import hapi\n",
- "\n",
  "server = 'https://cdaweb.gsfc.nasa.gov/hapi'\n",
- "dataset = 'PSP_SWP_SPC_L3I'\n",
- "start = \"2018-10-30T23:59:59\" #2018-10-30T23:59:59Z\n",
- "stop = '2018-11-04T22:39:56' #2023-01-31T22:39:56Z\n",
+ "dataset = 'HELIOS2_40SEC_MAG-PLASMA'\n",
+ "start = \"1976-09-12T01:21:38Z\" #1974-12-12T01:21:38Z\n",
+ "stop = '1976-11-18T01:21:38Z' #1985-09-04T09:52:10Z\n",
  "\n",
  "# Get data\n",
  "time = hapi(server, dataset, 'Time', start, stop)\n",
- "t_p = hapi(server, dataset, 'wp_fit', start, stop)\n",
- "t_a = hapi(server, dataset, 'wa_fit', start, stop)\n",
+ "t_p = hapi(server, dataset, 'Tp', start, stop)\n",
+ "t_a = hapi(server, dataset, 'T_a', start, stop)\n",
  "\n",
  "# Calculate and clean up theta values\n",
  "theta = []\n",
@@ -307,7 +307,8 @@
  "\n",
  "plt.figtext(0.25, -0.04, txt, ha=\"left\", fontsize=fs)\n",
  "plt.grid()\n",
- "plt.show()\n"
+ "plt.show()\n",
+ "\n"
  ]
  },
  {
@@ -327,58 +328,120 @@
  },
  "outputs": [],
  "source": [
+ "r_n = hapi(server, dataset, 'R_Helio', start, stop)\n",
+ "n_p = hapi(server, dataset, 'Np', start, stop)\n",
+ "n_a = hapi(server, dataset, 'N_a', start, stop)\n",
+ "v_p = hapi(server, dataset, 'Vp', start, stop)\n",
+ "\n",
+ "theta_pred = []\n",
+ "L = len(r_n[0])\n",
+ "\n",
+ "for i in range(L):\n",
+ " params = {\n",
+ " \"r_0\": [r_n[0][i][1]] * u.au,\n",
+ " \"r_n\": [1.0] * u.au,\n",
+ " \"n_1\": [n_p[0][i][1]] * u.cm**-3,\n",
+ " \"n_2\": [n_a[0][i][1]] * u.cm**-3,\n",
+ " \"v_1\": [v_p[0][i][1]] * u.km / u.s,\n",
+ " \"T_1\": [abs(t_p[0][i][1])] * u.K,\n",
+ " \"T_2\": [abs(t_a[0][i][1])] * u.K,\n",
+ " \"ions\": [\"p+\", \"He-4++\"],\n",
+ " }\n",
+ "\n",
+ " theta_pred.append(float(coal.temp_ratio(**params)[0]))\n",
+ " print('\\r', f\"{(i / L) * 100:.2f} %\", end=\"\")\n",
+ "\n",
+ "theta_pred = [x for x in theta_pred if str(x) != 'nan']\n",
+ "print(theta_pred)\n",
+ "# Plot the results\n",
+ "fs = 14 # Default font size\n",
+ "figwidth, figheight = plt.rcParams[\"figure.figsize\"]\n",
+ "figheight = 1.6 * figheight\n",
+ "fig = plt.figure(figsize=[figwidth, figheight])\n",
  "\n",
- "n_p = hapi(server, dataset, 'P_DENS', start, stop)\n",
- "n_a = hapi(server, dataset, 'A_DENS', start, stop)\n",
- "v_p_gse = hapi(server, dataset, 'P_VELS', start, stop)\n",
- "\n",
- "print(n_p)\n",
- "\n",
- "v_p = []\n",
- "\n",
- "for i in range(len(v_p_gse[0])):\n",
- " arg = v_p_gse[0][i][1][0]**2 + v_p_gse[0][i][1][1]**2 + v_p_gse[0][i][1][2]**2\n",
- " v_p.append(np.sqrt(arg))\n",
- "\n",
- "params = {\n",
- " \"r_0\": [r_start] * u.au,\n",
- " \"r_n\": [r_end] * u.au,\n",
- " \"n_1\": n_p * u.cm**-3,\n",
- " \"n_2\": n_a * u.cm**-3,\n",
- " \"v_1\": v_p * u.km / u.s,\n",
- " \"T_1\": t_p * u.K,\n",
- " \"T_2\": t_a * u.K,\n",
- " \"ions\": [\"p+\", \"He-4++\"],\n",
- " }\n",
- "\n",
- "print(params)\n",
- "\n",
- "coal.temp_ratio()"
+ "plt.hist(theta_pred, range=(0, 10), density=1)\n",
+ "plt.grid()\n",
+ "plt.show()"
  ]
  },
  {
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "## Comparison\n",
  "\n",
  "here we introduce the wind data\n",
  "\n",
  "pyspedas, spacepy, pystat\n",
  "\n",
- "https://github.com/hapi-server/client-python/blob/master/README.md"
+ "https://github.com/hapi-server/client-python/blob/master/README.md\n",
+ "\n",
+ "\n",
+ "what if instead of the alpha particles being heated prefentially, but that the protons are under heated.\n",
+ "by this i mean that because there are so many of them, some are not heated and this disparity between the heated and underheated ions is what is causing the behaviour."
  ]
  },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "[maruca2013]: https://www.doi.org/10.1103/PhysRevLett.111.241101\n",
+ "\n",
+ "## Comparision\n",
+ "\n",
+ "Now that we have made a prediction, we can compare it with the actual observed values from the Wind spacecraft."
+ ],
+ "metadata": {
+ "collapsed": false
+ }
+ },
  {
  "cell_type": "code",
  "execution_count": null,
+ "outputs": [],
+ "source": [
+ "# Get data from Wind\n",
+ "dataset = 'WI_PM_3DP'\n",
+ "start = \"1996-11-01T01:21:38\" #1994-11-15T12:53:43Z\n",
+ "stop = '1996-11-05T01:21:38' #2023-04-24T23:59:58Z\n",
+ "\n",
+ "# Get data\n",
+ "time = hapi(server, dataset, 'Time', start, stop)\n",
+ "t_p = hapi(server, dataset, 'P_TEMP', start, stop)\n",
+ "t_a = hapi(server, dataset, 'A_TEMP', start, stop)\n",
+ "\n",
+ "theta = []\n",
+ "for i in range(len(t_p[0])):\n",
+ " theta.append(t_a[0][i][1]/t_p[0][i][1])\n",
+ "\n",
+ "theta = [x for x in theta if str(x) != 'nan']\n",
+ "\n",
+ "# Plot the results\n",
+ "fs = 14 # Default font size\n",
+ "figwidth, figheight = plt.rcParams[\"figure.figsize\"]\n",
+ "figheight = 1.6 * figheight\n",
+ "fig = plt.figure(figsize=[figwidth, figheight])\n",
+ "\n",
+ "plt.hist(theta, range=(0,10), density=1)\n",
+ "\n",
+ "plt.show()\n"
+ ],
  "metadata": {
+ "collapsed": false,
  "pycharm": {
  "is_executing": true
  }
- },
+ }
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
  "outputs": [],
- "source": []
+ "source": [],
+ "metadata": {
+ "collapsed": false,
+ "pycharm": {
+ "is_executing": true
+ }
+ }
  }
  ],
  "metadata": {