{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Current surfaces\n\nThis example demonstrates the resulting current surfaces\nfrom the curl-free and divergence-free SEC poles.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nfrom matplotlib.colors import LogNorm\nimport numpy as np\nfrom pysecs import SECS\n\n# Radius of Earth\nR_E = 6378e3\n\n# Pole of the current system at the North Pole\nsec_loc = np.array([90., 0., R_E + 100e3])\n\n# Set up the system with a single SEC\nsystem_df = SECS(sec_df_loc=sec_loc)\nsystem_cf = SECS(sec_cf_loc=sec_loc)\n\n# Fit unit currents since we aren't fitting to any data\nsystem_df.fit_unit_currents()\nsystem_cf.fit_unit_currents()\n\n# Scale the system corresponding to Figure 2 of\n# Amm and Viljannen (1999) doi:10.1186/BF03352247\n# 10 kA\ntotal_current = 10000\nsystem_df.sec_amps *= total_current\nsystem_cf.sec_amps *= total_current\n\n# Create prediction surface\ndelta = 1\n\nlat_min, lat_max = -90, 90\nlon_min, lon_max = -180, 180\nlats = np.arange(lat_min, lat_max, delta)\nlons = np.arange(lon_min, lon_max, delta)\n\nlats2 = np.arange(lat_min-delta/2, lat_max+delta/2, delta)\nlons2 = np.arange(lon_min-delta/2, lon_max+delta/2, delta)\n\nnlat = len(lats)\nnlon = len(lons)\n\nxx, yy = np.meshgrid(lons, lats)\nxx2, yy2 = np.meshgrid(lons2, lats2)\n\npoints = np.zeros((nlat*nlon, 3))\npoints[:, 0] = yy.ravel()\npoints[:, 1] = xx.ravel()\npoints[:, 2] = R_E\n\nshell_points = points.copy()\nshell_points[:, 2] = sec_loc[2]\n\n# Predict the currents\nJ_df = system_df.predict(shell_points, J=True)\nJ_cf = system_cf.predict(shell_points, J=True)\n\nJ_df_x = J_df[:, 0].reshape(nlat, nlon)\nJ_df_y = J_df[:, 1].reshape(nlat, nlon)\nJ_df_phi = np.linalg.norm(J_df[:, :2], axis=1).reshape(nlat, nlon)\nJ_df_r = J_df[:, 2].reshape(nlat, nlon)\n\nJ_cf_x = J_cf[:, 0].reshape(nlat, nlon)\nJ_cf_y = J_cf[:, 1].reshape(nlat, nlon)\nJ_cf_theta = np.linalg.norm(J_cf[:, :2], axis=1).reshape(nlat, nlon)\nJ_cf_r = J_cf[:, 2].reshape(nlat, nlon)\n\nfig, (ax_cf, ax_df) = plt.subplots(figsize=(8, 3), ncols=2,\n                                   constrained_layout=True)\nfig.suptitle('Amm and Viljanen (1999) Figure 1')\n\n# Limit quiver to only a few points\nnstep = 30\n\nnorm = LogNorm(vmin=1e-5, vmax=1e-2)\ncax = ax_cf.pcolormesh(xx2, yy2, J_cf_theta, norm=norm)\nax_cf.quiver(xx[nstep//2::nstep, nstep//2::nstep*2],\n             yy[nstep//2::nstep, nstep//2::nstep*2],\n             J_cf_y[nstep//2::nstep, nstep//2::nstep*2],\n             J_cf_x[nstep//2::nstep, nstep//2::nstep*2],\n             color='w', zorder=9)\nax_cf.set_ylabel('Latitude (deg)')\nax_cf.set_xlabel('Longitude (deg)')\n\ncax = ax_df.pcolormesh(xx2, yy2, J_df_phi, norm=norm)\nax_df.quiver(xx[nstep//2::nstep, nstep//2::nstep*2],\n             yy[nstep//2::nstep, nstep//2::nstep*2],\n             J_df_y[nstep//2::nstep, nstep//2::nstep*2],\n             J_df_x[nstep//2::nstep, nstep//2::nstep*2],\n             color='w', zorder=9)\nax_df.set_ylabel('Latitude (deg)')\nax_df.set_xlabel('Longitude (deg)')\n\nax_cf.set_title('Curl Free')\nax_df.set_title('Divergence Free')\n\nplt.show()"
      ]
    }
  ],
  "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.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}