From a958a36f1d9a9f672b0f16cb9155dfb3cdc6d607 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 23 Jun 2021 16:51:41 -0400 Subject: [PATCH] Experimenting with f2py module. Orbital elements conversion code from Swiftest is now embedded in swiftestio which is used to generate elements when user outputs XV --- examples/symba_mars_disk/param.in | 2 +- examples/symba_mars_disk/test_orbel.ipynb | 1901 +++++++++++++++++++++ python/swiftestio/mars-scatter.ipynb | 1147 ------------- python/swiftestio/orbel.f90 | 160 -- python/swiftestio/swiftestio.py | 70 +- src/modules/module_interfaces.f90 | 2 +- src/orbel/orbel_xv2el.f90 | 2 +- src/python/orbel.f90 | 32 + 8 files changed, 2003 insertions(+), 1313 deletions(-) create mode 100644 examples/symba_mars_disk/test_orbel.ipynb delete mode 100644 python/swiftestio/mars-scatter.ipynb delete mode 100644 python/swiftestio/orbel.f90 create mode 100644 src/python/orbel.f90 diff --git a/examples/symba_mars_disk/param.in b/examples/symba_mars_disk/param.in index b96460880..b4b2ce6b7 100644 --- a/examples/symba_mars_disk/param.in +++ b/examples/symba_mars_disk/param.in @@ -10,7 +10,7 @@ ISTEP_DUMP 1 BIN_OUT bin.dat PARTICLE_FILE particle.dat OUT_TYPE REAL8 -OUT_FORM EL +OUT_FORM XV OUT_STAT REPLACE J2 0.0 J4 0.0 diff --git a/examples/symba_mars_disk/test_orbel.ipynb b/examples/symba_mars_disk/test_orbel.ipynb new file mode 100644 index 000000000..a53e45d0c --- /dev/null +++ b/examples/symba_mars_disk/test_orbel.ipynb @@ -0,0 +1,1901 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 76, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The autoreload extension is already loaded. To reload it, use:\n", + " %reload_ext autoreload\n" + ] + } + ], + "source": [ + "%load_ext autoreload\n", + "%autoreload 1\n", + "%aimport swiftestio" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "metadata": {}, + "outputs": [], + "source": [ + "import swiftestio as swio" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading Swiftest file param.in\n" + ] + } + ], + "source": [ + "config = swio.read_swiftest_config(\"param.in\")" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Reading in time 6.000e+03\n", + "Creating Dataset\n", + "\n", + "Adding particle info to Dataset\n", + "Successfully converted 11 output frames.\n" + ] + } + ], + "source": [ + "ds = swio.swiftest2xr(config)" + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (id: 1521, time: 11)\n",
+       "Coordinates:\n",
+       "  * id           (id) int64 1 2 3 4 5 6 7 ... 1515 1516 1517 1518 1519 1520 1521\n",
+       "  * time         (time) float64 0.0 600.0 1.2e+03 ... 4.8e+03 5.4e+03 6e+03\n",
+       "Data variables:\n",
+       "    px           (time, id) float64 0.0 -2.358e+06 ... 2.857e+06 2.775e+06\n",
+       "    py           (time, id) float64 0.0 8.604e+06 ... -8.703e+06 -8.67e+06\n",
+       "    pz           (time, id) float64 0.0 1.252e+04 ... 5.26e+04 5.588e+04\n",
+       "    vx           (time, id) float64 0.0 -2.134e+03 ... 2.049e+03 2.034e+03\n",
+       "    vy           (time, id) float64 0.0 -591.2 -1.533e+03 ... 645.5 698.3 711.1\n",
+       "    vz           (time, id) float64 0.0 -1.465 -0.1413 ... 4.414 5.343 2.951\n",
+       "    a            (time, id) float64 0.0 9.113e+06 ... 9.176e+06 8.984e+06\n",
+       "    e            (time, id) float64 0.0 0.02115 0.008838 ... 0.01143 0.0297\n",
+       "    inc          (time, id) float64 0.0 0.00155 0.001431 ... 0.006276 0.006327\n",
+       "    capom        (time, id) float64 0.0 6.112 2.23 ... 3.898 3.874 3.696\n",
+       "    omega        (time, id) float64 0.0 2.144 0.9395 4.552 ... 4.107 2.584 3.389\n",
+       "    capm         (time, id) float64 0.0 6.154 0.6746 3.227 ... 3.298 4.877 4.273\n",
+       "    mass         (time, id) float64 4.284e+13 9.907e+04 ... 1.168e+05 1.168e+05\n",
+       "    radius       (time, id) float64 3.39e+06 7.076e+03 ... 5.073e+03 5.073e+03\n",
+       "    rot_x        (time, id) float64 0.0 0.0 0.0 ... -0.0003309 -0.0003309\n",
+       "    rot_y        (time, id) float64 0.0 0.0 0.0 ... -0.0003242 -0.0003242\n",
+       "    rot_z        (time, id) float64 0.0 0.0 0.0 ... -8.195e-05 -8.195e-05\n",
+       "    Ip_x         (time, id) float64 0.4 0.4 0.4 0.4 0.4 ... 0.4 0.4 0.4 0.4 0.4\n",
+       "    Ip_y         (time, id) float64 0.4 0.4 0.4 0.4 0.4 ... 0.4 0.4 0.4 0.4 0.4\n",
+       "    Ip_z         (time, id) float64 0.4 0.4 0.4 0.4 0.4 ... 0.4 0.4 0.4 0.4 0.4\n",
+       "    npl          (time) int32 1501 1501 1501 1501 1501 ... 1518 1518 1518 1517\n",
+       "    ntp          (time) int32 0 0 0 0 0 0 0 0 0 0 0\n",
+       "    time_origin  (id) float64 0.0 0.0 0.0 0.0 ... 2.4e+03 2.4e+03 2.4e+03\n",
+       "    px_origin    (id) float64 0.0 -2.358e+06 ... -3.477e+06 -3.501e+06\n",
+       "    py_origin    (id) float64 0.0 8.604e+06 -6.936e+06 ... -8.553e+06 -8.542e+06\n",
+       "    pz_origin    (id) float64 0.0 1.252e+04 1.515e+04 ... 2.495e+04 3.359e+04\n",
+       "    vx_origin    (id) float64 0.0 -2.134e+03 1.312e+03 ... 1.999e+03 1.975e+03\n",
+       "    vy_origin    (id) float64 0.0 -591.2 -1.533e+03 ... -821.3 -788.3 -781.6\n",
+       "    vz_origin    (id) float64 0.0 -1.465 -0.1413 0.2818 ... 9.078 12.49 11.65\n",
+       "    origin_type  (id) <U32 'Central body' ... 'Supercatastrophic'
" + ], + "text/plain": [ + "\n", + "Dimensions: (id: 1521, time: 11)\n", + "Coordinates:\n", + " * id (id) int64 1 2 3 4 5 6 7 ... 1515 1516 1517 1518 1519 1520 1521\n", + " * time (time) float64 0.0 600.0 1.2e+03 ... 4.8e+03 5.4e+03 6e+03\n", + "Data variables:\n", + " px (time, id) float64 0.0 -2.358e+06 ... 2.857e+06 2.775e+06\n", + " py (time, id) float64 0.0 8.604e+06 ... -8.703e+06 -8.67e+06\n", + " pz (time, id) float64 0.0 1.252e+04 ... 5.26e+04 5.588e+04\n", + " vx (time, id) float64 0.0 -2.134e+03 ... 2.049e+03 2.034e+03\n", + " vy (time, id) float64 0.0 -591.2 -1.533e+03 ... 645.5 698.3 711.1\n", + " vz (time, id) float64 0.0 -1.465 -0.1413 ... 4.414 5.343 2.951\n", + " a (time, id) float64 0.0 9.113e+06 ... 9.176e+06 8.984e+06\n", + " e (time, id) float64 0.0 0.02115 0.008838 ... 0.01143 0.0297\n", + " inc (time, id) float64 0.0 0.00155 0.001431 ... 0.006276 0.006327\n", + " capom (time, id) float64 0.0 6.112 2.23 ... 3.898 3.874 3.696\n", + " omega (time, id) float64 0.0 2.144 0.9395 4.552 ... 4.107 2.584 3.389\n", + " capm (time, id) float64 0.0 6.154 0.6746 3.227 ... 3.298 4.877 4.273\n", + " mass (time, id) float64 4.284e+13 9.907e+04 ... 1.168e+05 1.168e+05\n", + " radius (time, id) float64 3.39e+06 7.076e+03 ... 5.073e+03 5.073e+03\n", + " rot_x (time, id) float64 0.0 0.0 0.0 ... -0.0003309 -0.0003309\n", + " rot_y (time, id) float64 0.0 0.0 0.0 ... -0.0003242 -0.0003242\n", + " rot_z (time, id) float64 0.0 0.0 0.0 ... -8.195e-05 -8.195e-05\n", + " Ip_x (time, id) float64 0.4 0.4 0.4 0.4 0.4 ... 0.4 0.4 0.4 0.4 0.4\n", + " Ip_y (time, id) float64 0.4 0.4 0.4 0.4 0.4 ... 0.4 0.4 0.4 0.4 0.4\n", + " Ip_z (time, id) float64 0.4 0.4 0.4 0.4 0.4 ... 0.4 0.4 0.4 0.4 0.4\n", + " npl (time) int32 1501 1501 1501 1501 1501 ... 1518 1518 1518 1517\n", + " ntp (time) int32 0 0 0 0 0 0 0 0 0 0 0\n", + " time_origin (id) float64 0.0 0.0 0.0 0.0 ... 2.4e+03 2.4e+03 2.4e+03\n", + " px_origin (id) float64 0.0 -2.358e+06 ... -3.477e+06 -3.501e+06\n", + " py_origin (id) float64 0.0 8.604e+06 -6.936e+06 ... -8.553e+06 -8.542e+06\n", + " pz_origin (id) float64 0.0 1.252e+04 1.515e+04 ... 2.495e+04 3.359e+04\n", + " vx_origin (id) float64 0.0 -2.134e+03 1.312e+03 ... 1.999e+03 1.975e+03\n", + " vy_origin (id) float64 0.0 -591.2 -1.533e+03 ... -821.3 -788.3 -781.6\n", + " vz_origin (id) float64 0.0 -1.465 -0.1413 0.2818 ... 9.078 12.49 11.65\n", + " origin_type (id) \n", - "\n", - "\n", - "Show/Hide data repr\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "Show/Hide attributes\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "
xarray.DataArray
'radmarker'
  • 1.0
    array(1.)
    • ntp
      ()
      int32
      0
      array(0, dtype=int32)
" - ], - "text/plain": [ - "\n", - "array(1.)\n", - "Coordinates:\n", - " ntp int32 0" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "disk['radmarker'] = disk['radmarker'] / 1\n", - "disk['radmarker'].max()" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "41" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "disk.sizes['time']" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "data": { - "text/html": [ - "
\n", - "\n", - "\n", - "Show/Hide data repr\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "Show/Hide attributes\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "\n", - "
xarray.Dataset
    • id: 1595
    • time: 41
    • ntp
      ()
      int32
      0
      array(0, dtype=int32)
    • id
      (id)
      float64
      2.0 3.0 4.0 ... 1.009e+04 1.01e+04
      array([2.0000e+00, 3.0000e+00, 4.0000e+00, ..., 1.0093e+04, 1.0094e+04,\n",
-       "       1.0095e+04])
    • time
      (time)
      float64
      0.0 500.0 1e+03 ... 1.95e+04 2e+04
      array([    0.,   500.,  1000.,  1500.,  2000.,  2500.,  3000.,  3500.,  4000.,\n",
-       "        4500.,  5000.,  5500.,  6000.,  6500.,  7000.,  7500.,  8000.,  8500.,\n",
-       "        9000.,  9500., 10000., 10500., 11000., 11500., 12000., 12500., 13000.,\n",
-       "       13500., 14000., 14500., 15000., 15500., 16000., 16500., 17000., 17500.,\n",
-       "       18000., 18500., 19000., 19500., 20000.])
    • npl
      (time)
      int32
      1501 1501 1501 ... 1565 1563 1563
      array([1501, 1501, 1501, 1501, 1501, 1501, 1509, 1509, 1509, 1509, 1509,\n",
-       "       1509, 1509, 1509, 1509, 1520, 1520, 1520, 1519, 1521, 1524, 1524,\n",
-       "       1524, 1524, 1532, 1531, 1530, 1530, 1530, 1538, 1538, 1538, 1546,\n",
-       "       1549, 1548, 1556, 1567, 1567, 1565, 1563, 1563], dtype=int32)
    • a
      (time, id)
      float64
      9.113e+06 1.067e+07 ... 1.119e+07
      array([[ 9112793.22504467, 10670485.73389117, 10699047.60825503, ...,\n",
-       "                      nan,               nan,               nan],\n",
-       "       [ 9112720.74850665, 10670512.89055318, 10698898.55338521, ...,\n",
-       "                      nan,               nan,               nan],\n",
-       "       [ 9112648.67658152, 10670538.29463011, 10698749.72478348, ...,\n",
-       "                      nan,               nan,               nan],\n",
-       "       ...,\n",
-       "       [ 9110990.16558072, 10670636.15293078, 10696488.33257031, ...,\n",
-       "        11085474.15997991, 11124729.12096808, 11178785.37557426],\n",
-       "       [ 9110911.37693149, 10670660.55610576, 10696418.61110883, ...,\n",
-       "        11085477.51062212, 11126343.28086734, 11182760.49522788],\n",
-       "       [ 9110842.49945745, 10670689.10881169, 10696357.85284678, ...,\n",
-       "        11085274.15893465, 11127832.86831051, 11187105.80328711]])
    • e
      (time, id)
      float64
      0.02115 0.008838 ... 0.0227 0.01982
      array([[0.02115088, 0.00883771, 0.00992666, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       [0.02114299, 0.00883922, 0.00994076, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       [0.02113537, 0.00884057, 0.00995482, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       ...,\n",
-       "       [0.02123603, 0.00884585, 0.00999627, ..., 0.02403232, 0.02299171,\n",
-       "        0.02036968],\n",
-       "       [0.0212364 , 0.00884677, 0.0099895 , ..., 0.02398331, 0.02284192,\n",
-       "        0.02009945],\n",
-       "       [0.02123559, 0.00884774, 0.00998352, ..., 0.02394456, 0.02270122,\n",
-       "        0.01981633]])
    • inc
      (time, id)
      float64
      0.00155 0.001431 ... 0.002293
      array([[0.00155003, 0.00143134, 0.00188344, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       [0.00155009, 0.00143134, 0.00188338, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       [0.00155009, 0.00143134, 0.00188351, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       ...,\n",
-       "       [0.00156436, 0.00142959, 0.00189326, ..., 0.00207236, 0.00262489,\n",
-       "        0.00219447],\n",
-       "       [0.00156431, 0.00142955, 0.00189368, ..., 0.00208184, 0.00264587,\n",
-       "        0.00223855],\n",
-       "       [0.00156421, 0.00142947, 0.0018942 , ..., 0.00209115, 0.00266912,\n",
-       "        0.00229255]])
    • capom
      (time, id)
      float64
      6.112 2.23 0.001328 ... 5.949 5.816
      array([[6.11248453e+00, 2.22989742e+00, 1.32835148e-03, ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       [6.11240957e+00, 2.22988740e+00, 1.83660995e-04, ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       [6.11240231e+00, 2.22988034e+00, 6.28223460e+00, ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       ...,\n",
-       "       [6.10471297e+00, 2.22905963e+00, 6.27242920e+00, ...,\n",
-       "        5.87491241e+00, 5.97081975e+00, 5.85889905e+00],\n",
-       "       [6.10470257e+00, 2.22911638e+00, 6.27190993e+00, ...,\n",
-       "        5.86962179e+00, 5.95960204e+00, 5.83691993e+00],\n",
-       "       [6.10467553e+00, 2.22919215e+00, 6.27141008e+00, ...,\n",
-       "        5.86529699e+00, 5.94948284e+00, 5.81565410e+00]])
    • omega
      (time, id)
      float64
      2.144 0.9395 4.552 ... 1.562 1.652
      array([[2.1440744 , 0.93950815, 4.55219084, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       [2.14408855, 0.93975764, 4.5534041 , ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       [2.14399145, 0.93998871, 4.5546904 , ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       ...,\n",
-       "       [2.15107552, 0.93680791, 4.57405186, ..., 1.74643297, 1.54470362,\n",
-       "        1.63401881],\n",
-       "       [2.15150441, 0.93642049, 4.57444671, ..., 1.75497467, 1.55380541,\n",
-       "        1.64414596],\n",
-       "       [2.15189707, 0.93598786, 4.57481251, ..., 1.76276262, 1.56176275,\n",
-       "        1.65167588]])
    • capm
      (time, id)
      float64
      6.154 0.6746 3.227 ... 3.943 3.985
      array([[6.1537018 , 0.6746454 , 3.22702473, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       [6.27272262, 0.76829546, 3.32047233, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       [0.10860316, 0.8619602 , 3.41383948, ...,        nan,        nan,\n",
-       "               nan],\n",
-       "       ...,\n",
-       "       [4.39231767, 4.24597955, 0.48843303, ..., 3.65611071, 3.7630202 ,\n",
-       "        3.78386363],\n",
-       "       [4.51089387, 4.34020076, 0.58210914, ..., 3.74133853, 3.8532184 ,\n",
-       "        3.88328671],\n",
-       "       [4.62952449, 4.43444762, 0.67579585, ..., 3.82636003, 3.94343154,\n",
-       "        3.98450411]])
    • Mass
      (time, id)
      float64
      6.612e-06 2.119e-06 ... 5.007e-06
      array([[6.61213283e-06, 2.11868146e-06, 2.66008238e-05, ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       [6.61213283e-06, 2.11868146e-06, 2.66008238e-05, ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       [6.61213283e-06, 2.11868146e-06, 2.66008238e-05, ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       ...,\n",
-       "       [6.61213283e-06, 2.11868146e-06, 2.66008238e-05, ...,\n",
-       "        5.00653973e-06, 5.00653973e-06, 5.00653973e-06],\n",
-       "       [6.61213283e-06, 2.11868146e-06, 2.66008238e-05, ...,\n",
-       "        5.00653973e-06, 5.00653973e-06, 5.00653973e-06],\n",
-       "       [6.61213283e-06, 2.11868146e-06, 2.66008238e-05, ...,\n",
-       "        5.00653973e-06, 5.00653973e-06, 5.00653973e-06]])
    • Radius
      (time, id)
      float64
      7.076e+03 4.842e+03 ... 4.377e+03
      array([[ 7076.43092   ,  4842.34399   , 11254.653     , ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       [ 7076.43092   ,  4842.34399   , 11254.653     , ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       [ 7076.43092   ,  4842.34399   , 11254.653     , ...,\n",
-       "                   nan,            nan,            nan],\n",
-       "       ...,\n",
-       "       [ 7076.43092   ,  4842.34399   , 11254.653     , ...,\n",
-       "         4376.86324968,  4376.86324968,  4376.86324968],\n",
-       "       [ 7076.43092   ,  4842.34399   , 11254.653     , ...,\n",
-       "         4376.86324968,  4376.86324968,  4376.86324968],\n",
-       "       [ 7076.43092   ,  4842.34399   , 11254.653     , ...,\n",
-       "         4376.86324968,  4376.86324968,  4376.86324968]])
    • radmarker
      (time, id)
      float64
      0.541 0.3702 ... 0.3346 0.3346
      array([[0.54095131, 0.37016857, 0.86035169, ..., 0.        , 0.        ,\n",
-       "        0.        ],\n",
-       "       [0.54095131, 0.37016857, 0.86035169, ..., 0.        , 0.        ,\n",
-       "        0.        ],\n",
-       "       [0.54095131, 0.37016857, 0.86035169, ..., 0.        , 0.        ,\n",
-       "        0.        ],\n",
-       "       ...,\n",
-       "       [0.54095131, 0.37016857, 0.86035169, ..., 0.33458532, 0.33458532,\n",
-       "        0.33458532],\n",
-       "       [0.54095131, 0.37016857, 0.86035169, ..., 0.33458532, 0.33458532,\n",
-       "        0.33458532],\n",
-       "       [0.54095131, 0.37016857, 0.86035169, ..., 0.33458532, 0.33458532,\n",
-       "        0.33458532]])
" - ], - "text/plain": [ - "\n", - "Dimensions: (id: 1595, time: 41)\n", - "Coordinates:\n", - " ntp int32 0\n", - " * id (id) float64 2.0 3.0 4.0 5.0 ... 1.009e+04 1.009e+04 1.01e+04\n", - " * time (time) float64 0.0 500.0 1e+03 1.5e+03 ... 1.9e+04 1.95e+04 2e+04\n", - " npl (time) int32 1501 1501 1501 1501 1501 ... 1567 1565 1563 1563\n", - "Data variables:\n", - " a (time, id) float64 9.113e+06 1.067e+07 ... 1.113e+07 1.119e+07\n", - " e (time, id) float64 0.02115 0.008838 0.009927 ... 0.0227 0.01982\n", - " inc (time, id) float64 0.00155 0.001431 ... 0.002669 0.002293\n", - " capom (time, id) float64 6.112 2.23 0.001328 ... 5.865 5.949 5.816\n", - " omega (time, id) float64 2.144 0.9395 4.552 5.005 ... 1.763 1.562 1.652\n", - " capm (time, id) float64 6.154 0.6746 3.227 ... 3.826 3.943 3.985\n", - " Mass (time, id) float64 6.612e-06 2.119e-06 ... 5.007e-06 5.007e-06\n", - " Radius (time, id) float64 7.076e+03 4.842e+03 ... 4.377e+03 4.377e+03\n", - " radmarker (time, id) float64 0.541 0.3702 0.8604 ... 0.3346 0.3346 0.3346" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "disk" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([ nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, nan,\n", - " nan, nan, 9016553.13090581,\n", - " 9016401.50991297, 9016149.65351738, 9015785.0788279 ,\n", - " 9015304.20639912, 9014713.28750873, 9014031.77210284,\n", - " 9013304.3685138 , 9012603.71454923])" - ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "disk['a'].isel(id=1563).values" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEWCAYAAABWn/G6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU9b3/8dc7CUmAAFkIAcKSsKgsFpCIgGKtWkSvLYi1tZtYW60t3qs/6622996fdrn91S7a2kVriy22datWxdq6FC3WBTAsYgCBsG9CgAQIAUKSz++P+caOMQtLZiYJn+fjMc6Zz/l+z/mccZhPzjnfOUdmhnPOORcvSYlOwDnn3MnFC49zzrm48sLjnHMurrzwOOeciysvPM455+LKC49zzrm48sLjXBsjqVLSoCbmXS3p1Xjn5Fxr8sLjOixJ50nakqB1j5a0SFJVeB59tH3NLMPM1sUyv2iSUiU9LmmDJJN0XoP5knSnpN3hcackRc1vclub6yupp6TXQrxC0huSzm4mzzRJD0jaJ+ldSTfH4O1wceCFx7lWJikVeBr4A5AFzAaeDvG26lXgc8C7jcy7DpgGjAI+BHwM+DIc1bY22ReoBK4BckPfO4FnJKU0keMdwFBgIPAR4OuSphzX1rqE8sLj2rXwV/otkpZJ2ivpUUnpkroCfwP6hkNXlZL6Sroj/HX/qKT9khZLGtXKaZ0HpAA/MbPDZnYPIOD8o9wmkzQkTOdImhP+yl8IDG7lXDGzajP7iZm9CtQ20mQG8GMz22JmW4EfA1eHeefR/LY22dfMDpnZKjOrC31qiRSg7CZSnQF8x8zKzWwl8OuoPFw74oXHdQSfBKYAhUT+qr7azA4AFwPbwqGrDDPbFtpPBf5E5AvuIeApSZ0aW3AoaBVNPH7ZRD4jgGX2/utRLQvxY/UL4BDQh8jewTXNNW4m1wpJtx3H+iGS91tRr9/iX9vS0rY217c+52VEtnEO8Bsz29nIdmUReQ+aXZZrH5rapXWuPbmnvqhIegZo6XzKIjN7PLS/C/gaMB74Z8OGZvah48gnA9jbILYX6HYsC5GUDFwOnB4KaYmk2cC5TfUxs8xjzPVoNNyevUBGOFfT0rY22be+WJnZhySlA5cBTR2OzIjq39h6XDviezyuI4g+L1HFv76kmrK5fiIc5tkC9G3FfCqB7g1i3YH9x7icXCJ/HG6Oim08gbyOV8Pt6Q5UhsLR0rY21/c94bDbw8BtTRz6rIzq39h6XDvihcd1ZE1der1//YSkJKAfsK2xhpKWR50javi4r4nlLwc+FD3yi8ghwOXHmH8ZUBOdLzCguQ7N5Fop6ZvHuP56y4kMDqg3in9tS0vb2lzfxnQCPjCU3MzKge3HuCzXRnnhcR3ZDiBHUo8G8bGSpofRUzcBh4H5jS3AzEZEnSNq+Li+ifX+g8iJ8v8IQ4BvCPGX4L3f4mxoKXkzqwX+DNwhqYuk4UROsDfXp6lcM8zse031C3mmh5epYYBGfTF5ELhZUr6kvkQOTf7uaLa1ub6Sxks6R5Hh3J0l3QrkAQuaSPNB4L8lZUk6Dbg2Kg/XjnjhcR2Wmb0DPAysCyfX6w+nPQ18CigHPg9MN7MjrbjeaiJDiK8CKogMCJgW4hDZg3ntKBd3A5FDh+8S+ZL9bWvl2cAq4CCQDzwfpgeGeb8CngHeBkqAZ0PsaLa1yb5AGpHBE7uBrcAlwL9Fna/7rKToPZrbgbVEDjfOA35oZs+12jvg4kZ+Izh3MpF0BzDEzD6XwBxeAG4MQ4KdO+n4qDbn4szMJic6B+cSyQ+1Oeeciys/1Oaccy6ufI/HOedcXPk5nhb07NnTCgoKEp2Gc861K4sWLdplZrmNzfPC04KCggKKi4sTnYZzzrUrkpq8yoYfanPOORdXXnicc87FlRce55xzceWFxznnXFx54XHOORdXXnicc87FlRce55xzceWFJ0Zq64zv/XUlm/dUJToV55xrU7zwxMjG3Qd4ZOEmpt/7OiVbG96S3jnnTl5eeGJkUG4Gj39lIp2SxKd+9QavrC5LdErOOdcmeOGJoVPyuvHnr55N/+wuXPO7N3l80ZZEp+SccwnnhSfGevdI57HrJ3DWoGxu+dNb/OLlUvxWFM65k5kXnjjont6J3149jmmj+/LD51fx30+VUFNbl+i0nHMuIfzq1HGSmpLEXZ8cTZ/Mztz7j7Xs2HeIn336DDqnJic6Neeciyvf44mjpCRx65TT+PbUEcx9Zyef/vV8dlceTnRazjkXV154EuCqCQXc97mxrNy+j8vvfZ0Nuw4kOiXnnIubmBYeSTdKKpG0XNJNIXZFeF0nqSiqbYGkg5KWhsd9jSxvjqSSqNfZkl6UtCY8Z4W4JN0jqVTSMklnRPWZEdqvkTQjltvfnItG9Oaha8ez9+ARpv7iNV4r3ZWoVJxzLq5iVngkjQSuBcYBo4BLJQ0BSoDpwCuNdFtrZqPD4/oGy5sOVDZofxsw18yGAnPDa4CLgaHhcR1wb1hGNnA7cFbI6/b6YpUIYwdm8fTMc+jdPZ2rHljI715b7yPenHMdXiz3eIYBC8ysysxqgHnAdDNbaWarjmVBkjKAm4HvNpg1FZgdpmcD06LiD1rEfCBTUh/gIuBFM9tjZuXAi8CU49m41jIgpwtPfHUiHzm1F3c8s4JvPvk21TU+4s0513HFsvCUAJMk5UjqAlwC9G+hT6GkJZLmSZoUFf8O8GOg4YXP8sxse5h+F8gL0/nA5qh2W0Ksqfj7SLpOUrGk4rKy2F9xICMthfs/P5aZHxnMwws387lZC3zQgXOuw4pZ4TGzlcCdwAvAc8BSoLaZLtuBAWY2hsjezUOSuksaDQw2sydbWJ8BrXKcyszuN7MiMyvKzc1tjUW2KClJ/OdFp/HTK0fz1uYKpv7iNVZu3xeXdTvnXDzFdHCBmc0ys7Fmdi5QDqxupu1hM9sdphcBa4FTgAlAkaQNwKvAKZL+EbrtCIfQCM87Q3wr79+76hdiTcXbjKmj83nsyxM4UlvH5fe+zvPL3010Ss4516piPaqtV3geQGRAwUPNtM2VlBymBxEZGLDOzO41s75mVgCcA6w2s/NCtzlA/ci0GcDTUfGrwui28cDecEjueWCypKwwqGByiLUpo/pn8swN5zA0rxtf/v0ifjZ3jQ86cM51GLH+Hc8TklYAzwAzzaxC0mWSthDZk3lWUv0X/7nAMklLgceB681sTwvL/z7wUUlrgAvDa4C/AuuAUuDXwFcBwvK+A7wZHt8+inUkRK/u6Tx63XguG5PPj19czZdmF7PnQHWi03LOuRMm/0u6eUVFRVZcXJyw9ZsZs1/fwPf++g5ZXTvxk0+NYcLgnITl45xzR0PSIjMramyeX7mgjZPE1WcX8uTMiXRNS+Ezv5nPXS+s8ouMOufaLS887cSIvj145oZzuPyMftzzUimf/vV8tlYcTHRazjl3zLzwtCNd01L40RWj+OmVo1m5fT+X/PSfPFfio96cc+2LF552aOrofJ79j3MYmNOF6/+wiP9+6m0OHWnuJ1LOOdd2eOFppwbmdOXx6ydy3bmD+MP8TUz7xWss2VSe6LTipq7OKD9QTenOSrb5IUfn2hUf1daCRI9qOxr/WLWTrz++jJ37D/OJsf34+pRT6dUtPdFpnZDDNbW8XrqbxZvK2X2gmj2V1eypqmbPgWrKD1RTXlVNXdRH95whPfnc+IFcOKwXKcn+95RzidbcqDYvPC1oD4UHoPJwDT9/qZRZr64jLSWZfz9/CF84u5DUlPbzJXzoSC3/WFXGcyXbmbtyJ/sP15AkyOqSSnbXVLK6ppLTNTId/di0u4qHF25i295D9O6ezpXj+vPpcQPI696+i69z7ZkXnhPQXgpPvfW7DvDdv6xg7js7GdSzK/9z6XA+clqvRKfVpMrDNbz0zk6eK9nOy++UcfBILZldOjF5eB4Xj+zDxCE5pKW0fHvwmto6Xl5Vxu/nb+SV1WUkJ4mPDsvj8xMGMnFwDpLisDXOuXpeeE5Aeys89V5etZPvPLOCdbsOcP5pvfifS4dT2LNrotMCIns2L6zYwZyl23hlTRnVNXX0zEjjohF5XHJ6H84qzD6hw2Ubdx/goQWbeKx4M+VVRxjUsys3XjiUqaM/cCFy51yMeOE5Ae218ABU19Qx+/UN/HTuGg7X1PKFswuZMbGA/MzOcc+lrs5YuGEPTy7eyl/f3s7+wzX06ZHOlJG9uXhkH8YOzCI5qXX3Sg4dqeWvb2/n1/9cz5od+3n11vPp3cMPvzkXD154TkB7Ljz1du4/xA+fW8WfFm0BYPygbKaP6cfFp/emW3qnmK57XVklTy7Zyp8Xb2VrxUG6piZz8el9mH5GPuMLc0hq5WLTmE27qzjvRy/z5Q8P5tYpp8V8fc45LzwnpCMUnnqb91Tx5JKtPLlkK+t3HSAtJYmPDs9j+hn5TBqaS6dWGA1WV2ds3FPFq6W7+PPiLSzZVEGS4OwhPbn8jH5MHpFHl9SUVtiaY/OVPyzitdJdvPGNC+iaFv/1O3ey8cJzAjpS4alnZizZXMFTS7byzFvbKK86Qk7XVD42qi8fG9WH/lldyOyS2uKIODNjS/lBlm3Zy7KtFby9ZS9vb93L/kM1AJySl8HlZ/Rj6uj8hB/iWrRxD5ff+wbfnjqCqyYUJDQX504GXnhOQEcsPNGqa+r4x6qdPLlkK3NX7qQ66uKj3dJSyArDmLO7dCK7axrZXTuRnJTEiu37eHtLBeVVRwDolCyG9enO6fk9+FC/Hozun8UpeRltZjSZmXHZL1+noqqauV87r9XPJznn3q+5wuPHHE5yqSlJTB7Rm8kjerO36ghvrNvFrsrIjzR3hx9q7jlQTVnlYVbvqGTPgWqqa+s4Ja8bk4f35vR+kUJzau9uRzXsOVEk8aVJhdzw0BLmrtzB5BG9E52ScyctLzzuPT26dGLKyD4ttqurs7gMCmhtU0b0Jj+zM795db0XHucSqP38rN21Ge2x6ACkJCfxhbMLWLh+D8u2VCQ6HedOWjEtPJJulFQiabmkm0LsivC6TlJRVNsCSQclLQ2P+0K8i6RnJb0T+n0/qk+apEcllUpaIKkgat43QnyVpIui4lNCrFTSbbHcftf2fOrM/mSkpfCbf65PdCrOnbRiVngkjQSuBcYBo4BLJQ0BSoDpwCuNdFtrZqPD4/qo+I/M7DRgDHC2pItD/ItAuZkNAe4G7gzrHg5cCYwApgC/lJQsKRn4BXAxMBz4dGjrThLd0jtx5Zn9efbt7X5Va+cSJJZ7PMOABWZWZWY1wDxgupmtNLNVR7uQ0P/lMF0NLAb6hdlTgdlh+nHgAkWGUU0FHjGzw2a2HiglUgDHAaVmti4s65HQ1p1Erj67ADNj9usbEp2KcyelWBaeEmCSpBxJXYBLgP4t9CmUtETSPEmTGs6UlAl8DJgbQvnAZoBQ3PYCOdHxYEuINRVvuJ7rJBVLKi4rK2t5S1270i+rCxef3oeHFm6i8nBNotNx7qQTs8JjZiuJHPp6AXgOWAo0d5vM7cAAMxsD3Aw8JKl7/UxJKcDDwD1mti5WeQOY2f1mVmRmRbm5ubFclUuQL51TyP5DNfypeHPLjZ1zrSqmgwvMbJaZjTWzc4FyYHUzbQ+b2e4wvQhYC5wS1eR+YI2Z/SQqtpWwFxUKUw9gd3Q86BdiTcXdSWbMgCzGDszigdfWU1vnP6J2Lp5iPaqtV3geQGRAwUPNtM0NJ/+RNAgYCqwLr79LpKjc1KDbHGBGmP4E8JJFLsUwB7gyjHorDMtaCLwJDJVUKCmVyACEOa2xra79+dI5hWzec5AXV7yb6FScO6nE+nc8T0haATwDzDSzCkmXSdoCTACelfR8aHsusEzSUiIDBa43sz2S+gH/RWQU2uIw1PpLoc8sIEdSKZHDc7cBmNly4DFgBZHDfDPNrDacB7oBeB5YCTwW2rqT0OQRvemf3dmHVjsXZ36tthZ09Gu1nex++9p6vvXMCp786kTGDMhKdDrOdRjNXavNr1zgTmpXFPWnW3oKs171vR7n4sULjzupZaSl8JlxA/hbybtsKa9KdDrOnRS88LiT3oyJBQD87rUNCc3DuZOFFx530uub2ZmLRuTx1NJt+DlP52LPC49zwDlDctlVeZgNu/1wm3Ox5oXHOWBcYWRE25vr9yQ4E+c6Pi88zgGDczPI7prKwg1eeJyLNS88zhG5NXbRwCze9MLjXMx54XEuGFeYzcbdVezcdyjRqTjXoXnhcS44syAbwA+3ORdjXnicC0b07U6X1GQfYOBcjHnhcS5ISU7ijAFZLNxQnuhUnOvQvPA4F+XMgmzeeXcfew8eSXQqznVYXnici3JmYRZmsHij7/U4FyteeJyLMqZ/Fp2SxQI/z+NczHjhcS5K59RkRub38N/zOBdDsb719Y2SSiQtl3RTiF0RXtdJKopqWyDpYLjD6FJJ90XNGyvpbUmlku6RpBDPlvSipDXhOSvEFdqVSlom6YyoZc0I7ddImoFzDYwryGbZlgoOHalNdCrOdUgxKzySRgLXAuOAUcClkoYAJcB04JVGuq01s9HhcX1U/N6wrKHhMSXEbwPmmtlQYG54DXBxVNvrQn8kZQO3A2eFvG6vL1bO1TuzIJsjtcbSzRWJTsW5DimWezzDgAVmVmVmNcA8YLqZrTSzVUe7EEl9gO5mNt8i16x/EJgWZk8FZofp2Q3iD1rEfCAzLOci4EUz22Nm5cCL/KuIOQdAUYFfMNS5WIpl4SkBJknKkdQFuATo30KfQklLJM2TNCnE8oEtUW22hBhAnpltD9PvAnlRfTY30qep+PtIuk5SsaTisrKyFlJ2HU1ml1ROzevmVzBwLkZiVnjMbCVwJ/AC8BywFGjuoPl2YICZjQFuBh6S1P0Y1mdAq9zFy8zuN7MiMyvKzc1tjUW6dubMwiwWbyynprYu0ak41+HEdHCBmc0ys7Fmdi5QDqxupu1hM9sdphcBa4FTgK1Av6im/UIMYEc4hFZ/SG5niG/l/XtX9X2aijv3PmcWZHOgupaV2/cnOhXnOpxYj2rrFZ4HEBlQ8FAzbXMlJYfpQUQGBqwLh9L2SRofRrNdBTwdus0B6kemzWgQvyqMbhsP7A3LeR6YLCkrDCqYHGLOvc+4Qr9gqHOxEuvf8TwhaQXwDDDTzCokXSZpCzABeFZS/Rf/ucAySUuBx4Hrzaz+X/1Xgd8ApUT2hP4W4t8HPippDXBheA3wV2BdaP/r0J+wvO8Ab4bHt6PW4dx7+vToTL+szj7AwLkYUOTUiGtKUVGRFRcXJzoNlwA3P7qUeavLKP7vCwk/HXPOHSVJi8ysqLF5fuUC55owrjCb3QeqWbfrQKJTca5D8cLjXBPODOd5/HCbc63LC49zTRjUsys9M1J9gIFzrcwLj3NNkETRwGy/YKhzrcwLj3PNOLMwm817DvLu3kOJTsW5DsMLj3PNGFfgv+dxrrV54XGuGcP6dKNrarIPMHCuFXnhca4ZKclJnDEwi4VeeJxrNV54nGvBuIJsVu3YT0VVdaJTca5D8MLjXAvqf89TvKE8wZk41zF44XGuBaP7Z9IpWT6s2rlW4oXHuRakd0rmQ/0yfWSbc63EC49zR+HMgmze3rKXg9XN3cvQOXc0vPA4dxTGFWZRU2cs2ezneZw7UV54nDsKYwdmI8Gb673wOHeivPA4dxR6dO7EqXndWLhhd6JTca7di/Wtr2+UVCJpuaSbQuyK8LpO0gduEiRpgKRKSbdExf5P6FMi6WFJ6SFeKGmBpFJJj0pKDfG08Lo0zC+IWtY3QnyVpItiuf2uYxk/KIdFG8s5XOPneZw7ETErPJJGAtcC44BRwKWShgAlwHTglSa63sW/bm2NpHzgP4AiMxsJJANXhtl3Aneb2RCgHPhiiH8RKA/xu0M7JA0PfUcAU4BfSkpulQ12Hd6EwTkcOlLHW5v3JjoV59q1WO7xDAMWmFmVmdUA84DpZrbSzFY11kHSNGA9sLzBrBSgs6QUoAuwTZF7EZ8PPB7azAamhemp4TVh/gWh/VTgETM7bGbrgVIihdG5Fo0vzEGC19fuSnQqzrVrsSw8JcAkSTmSugCXAP2baiwpA7gV+FZ03My2Aj8CNgHbgb1m9gKQA1SEogawBcgP0/nA5tC/Btgb2r8Xb6RPdC7XSSqWVFxWVnZMG+06rh5dOjGib3feWOvneZw7ETErPGa2ksghrheA54ClQHMHx+8gctisMjooKYvInkoh0BfoKulzsci5npndb2ZFZlaUm5sby1W5dmbi4J4s2VTBoSN+nse54xXTwQVmNsvMxprZuUTOwaxupvlZwA8kbQBuAr4p6QbgQmC9mZWZ2RHgz8BEYDeQGQ6/AfQDtobprYS9qzC/R2j/XryRPs61aMLgHKpr61i00YdVO3e8Yj2qrVd4HkBkQMFDTbU1s0lmVmBmBcBPgO+Z2c+JHGIbL6lLOE9zAbDSzAx4GfhEWMQM4OkwPSe8Jsx/KbSfA1wZRr0VAkOBha22wa7DO7Mgm+Qk+Xke505ASstNTsgTknKAI8BMM6uQdBnwMyAXeFbSUjNrclizmS2Q9DiwGKgBlgD3h9m3Ao9I+m6IzwrxWcDvJZUCewij4MxsuaTHgBVhWTPNzI+ZuKOWkZbCqH49eN3P8zh33BTZEXBNKSoqsuLi4kSn4dqQHz7/DvfNW8dbt08mIy3Wf7s51z5JWmRmH/itJviVC5w7ZhMH96S2zvx22M4dJy88zh2jsQOzSE1O4o11frjNuePhhce5Y5TeKZkxAzJ9gIFzx8kLj3PHYeLgnizfto+9VUcSnYpz7Y4XHueOw4TBOZjB/PV+uM25Y3XUQ3LCFQSGAun1MTNr6kKfznVoo/tnkt4piTfW7uaiEb0TnY5z7cpRFR5JXwJuJPJL/6XAeOANIhfpdO6kk5qSxJkF2X7dNueOw9EearsROBPYaGYfAcYAFTHLyrl2YMLgHFbt2M+uysOJTsW5duVoC88hMzsEkZusmdk7wKmxS8u5tm/i4J4AzPdh1c4dk6MtPFskZQJPAS9KehrYGLu0nGv7RvbtTkZail8+x7ljdFTneMzssjB5h6SXiVzt+bmYZeVcO5CSnMRZhdnM98Lj3DE55uHUZjbPzOaYWXUsEnKuPZkwOId1uw7w7t5DiU7FuXbDf8fj3AmYMDgHgDfW+VUMnDtaXnicOwHDencns0snXi/1w23OHS0vPM6dgKQkMb4wxwcYOHcMvPA4d4ImDslha8VBNu+pSnQqzrULXnicO0ETBkXO8/jVqp07OjEtPJJulFQiabmkm0LsivC6TtIH7k4naYCkSkm3RMUyJT0u6R1JKyVNCPFsSS9KWhOes0Jcku6RVCppmaQzopY1I7RfI2lGLLffnRyG9MqgZ0aaXz7HuaMUs8IjaSRwLTAOGAVcKmkIUAJMB5q6wOhdwN8axH4KPGdmp4VlrQzx24C5ZjYUmBteA1xM5IKmQ4HrgHtDTtnA7cBZIa/b64uVc8dLEhMHR87z+K3knWtZLPd4hgELzKzKzGqAecB0M1tpZqsa6yBpGrAeWB4V6wGcC8wCMLNqM6u/TtxUYHaYng1Mi4o/aBHzgUxJfYCLgBfNbI+ZlQMvAlNab5PdyWrC4Bx27j/Mul0HEp2Kc21eLAtPCTBJUo6kLsAlQP+mGkvKAG4FvtVgViFQBvxW0hJJv5HUNczLM7PtYfpdIC9M5wObo5axJcSaijfM5TpJxZKKy8rKjmJT3clu4uD68zx+uM25lsSs8JjZSuBO4AUil9dZCtQ20+UO4G4zq2wQTwHOAO41szHAAf51SC16fQa0ynEOM7vfzIrMrCg3N7c1Fuk6uAHZXejbI503fICBcy2K6eACM5tlZmPN7FygHFjdTPOzgB9I2gDcBHxT0g1E9kq2mNmC0O5xIoUIYEc4hEZ43hniW3n/3lW/EGsq7twJkcSEwT2Zv24PdXV+nse55sR6VFuv8DyAyICCh5pqa2aTzKzAzAqAnwDfM7Ofm9m7wGZJ9bdhuABYEabnAPUj02YAT0fFrwqj28YDe8MhueeByZKywqCCySHm3AmbODiHPQeqWbVjf6JTca5NO+pbXx+nJyTlAEeAmWZWIeky4GdALvCspKVmdlELy/l34I+SUoF1wBdC/PvAY5K+SOQ2DZ8M8b8SOadUClTVtzezPZK+A7wZ2n3bzPa0xoY6995129buZlif7gnOxrm2Sz78s3lFRUVWXFyc6DRcO3HeD19mUG4GD1x9ZqJTcS6hJC0ysw/8VhP8ygXOtaoLh+Xx6ppd7K06kuhUnGuzvPA414qmjs6nuraOv5Vsb7mxcycpLzzOtaKR+d0ZlNuVp5b6YEnnmuKFx7lWJImpo/JZsH4P2/ceTHQ6zrVJXnica2VTR/fFDJ55a1uiU3GuTfLC41wrK+jZldH9M3lqiRce5xrjhce5GJg6ui8rtu9jjf+Y1LkP8MLjXAxc+qG+JAmeXup7Pc415IXHuRjI7ZbG2UN68vRbW/0ePc414IXHuRiZNjqfzXsOsnhTeaJTca5N8cLjXIxcNLI3aSlJfrjNuQa88DgXIxlpKVw4PI+/LNvOkdq6RKfjXJvhhce5GJo2Op89B6p5dY3fIM65el54nIuhD5+SS4/OnXjaL6Hj3Hu88DgXQ6kpSVxyeh9eWLGDquqaRKfjXJvghce5GJs6ui9V1bW8uGJHolNxrk2I9a2vb5RUImm5pJtC7Irwuk7SB24SJGmApEpJtzSIJ0taIukvUbFCSQsklUp6NNyhFElp4XVpmF8Q1ecbIb5KUkt3PnXuhI0ryKZvj3Qf3eZcELPCI2kkcC0wDhgFXCppCFACTAdeaaLrXcDfGonfCKxsELsTuNvMhgDlwBdD/ItAeYjfHdohaThwJTACmAL8UlLycW2gc0cpKUl8bHRfXlldxp4D1YlOx7mEi+UezzBggZlVmVkNMA+YbmYrzWxVYx0kTQPWA8sbxPsB/wb8Jiom4Hzg8RCaDUwL01PDa8L8C0L7qcAjZnbYzNYDpUQKo3MxNXVUPjV1xrPLfK/HuVgWnhJgkqQcSV2AS4D+TTWWlAHcCnyrkc3+p5IAABeLSURBVNk/Ab4ORP8YIgeoCEUNYAuQH6bzgc0AYf7e0P69eCN9onO5TlKxpOKysrKWttO5Fg3r041T8jL8cJtzxLDwmNlKIoe4XgCeA5YCtc10uYPIYbPK6KCkS4GdZrYoRql+gJndb2ZFZlaUm5sbr9W6DkwSU0fnU7yxnM17qhKdjnMJFdPBBWY2y8zGmtm5RM7BrG6m+VnADyRtAG4CvinpBuBs4OMh/ghwvqQ/ALuBTEkpoX8/oP7HElsJe1dhfo/Q/r14I32ci6mPj+oLwBy/QZw7ycV6VFuv8DyAyICCh5pqa2aTzKzAzAqIHFr7npn93My+YWb9QvxK4CUz+5xFLvn7MvCJsIgZwNNhek54TZj/Umg/B7gyjHorBIYCC1tvi51rWv/sLhQNzOKpJX7Fandyi/XveJ6QtAJ4BphpZhWSLpO0BZgAPCvp+RNY/q3AzZJKiZzDmRXis4CcEL8ZuA3AzJYDjwEriBz+m2lmzR3+c65VTR3dlzU7K1m53W8Q505e8r+8mldUVGTFxcWJTsN1EHsOVDPuf//ONecU8s1LhiU6HediRtIiM/vAbzXBr1zgXFxld01l8og8HlqwyX/T405aXnici7P/c+EpVFXXcO8/ShOdinMJ4YXHuTgbmteNy8b0Y/YbG9m+92Ci03Eu7rzwOJcAN104FAx++vc1iU7FubjzwuNcAvTP7sJnzhrAnxZtYW1ZZcsdnOtAvPA4lyA3nD+EtJQk7nqhud9VO9fxeOFxLkF6ZqTxpXMKefbt7ZRs3ZvodJyLGy88ziXQl84dRGaXTvzg+UYv2O5ch+SFx7kE6p7eiZnnDeGV1WW8sXZ3otNxLi688DiXYJ+fMJDe3dP5wfPv+DXc3EnBC49zCZbeKZkbLxzKkk0VvLhiR6LTcS7mvPA41wZcMbYfg3p25UcvrKK2zvd6XMfmhce5NiAlOYmbJ5/C6h2VPL3UbxHlOjYvPM61EZeM7MOIvt2568XVVNfUtdzBuXbKC49zbURSkvj6lNPYUn6QhxduSnQ6zsWMFx7n2pBzh/bkrMJsfvbSGg4crkl0Os7FRKxvfX2jpBJJyyXdFGJXhNd1kj5wkyBJAyRVSrolvO4v6WVJK0K/G6PaZkt6UdKa8JwV4pJ0j6RSScsknRHVZ0Zov0bSjIbrdy6RpMhez67Kan71yrpEp+NcTMSs8EgaCVwLjANGAZdKGgKUANOBV5roehfwt6jXNcDXzGw4MB6YKWl4mHcbMNfMhgJzw2uAi4Gh4XEdcG/IKRu4HTgr5HV7fbFyrq0YOzCLj4/qy89fWsPrpbsSnY5zrS6WezzDgAVmVmVmNcA8YLqZrTSzRq8PImkasB5YXh8zs+1mtjhM7wdWAvlh9lRgdpieDUyLij9oEfOBTEl9gIuAF81sj5mVAy8CU1pvk51rHd+bfjqDcjO44eElbKvwe/a4jiWWhacEmCQpR1IX4BKgf1ONJWUAtwLfaqZNATAGWBBCeWa2PUy/C+SF6Xxgc1TXLSHWVLzheq6TVCypuKysrKl0nIuZjLQUfvX5sVTX1PGVPyzi0JHaRKfkXKuJWeExs5XAncALwHPAUqC5fz13AHebWaM3JwmF6QngJjPb18j6DGiVX96Z2f1mVmRmRbm5ua2xSOeO2eDcDH50xSje2rKXbz2zvOUOzrUTMR1cYGazzGysmZ0LlAPN3XjkLOAHkjYANwHflHQDgKRORIrOH83sz1F9doRDaITnnSG+lffvXfULsabizrVJU0b25qvnDebhhZt59E0fYu06hliPausVngcQGVDwUFNtzWySmRWYWQHwE+B7ZvZzSQJmASvN7K4G3eYA9SPTZgBPR8WvCqPbxgN7wyG554HJkrLCoILJIeZcm/W1yacyaWhP/ufp5by1uSLR6Th3wmL9O54nJK0AngFmmlmFpMskbQEmAM9KaumL/2zg88D5kpaGxyVh3veBj0paA1wYXgP8FVgHlAK/Br4KYGZ7gO8Ab4bHt0PMuTYrOUn89Mox5Gak8ZU/LGJ35eFEp+TcCZFfhr15RUVFVlxcnOg0nOPtLXu5/L7XObMgi9lfGEdKsv/+28XOpt1VdE5NJrdb2nH1l7TIzD7wW03wKxc4126c3q8H3502ktdKd/OjF5o7XerciTlcU8tX/riIz/5mPnUxuFq6Fx7n2pFPFvXnM2cN4L55a3muZHvLHZw7Dj9+YTXLt+3jlsmnkpSkVl++Fx7n2pnbPzac0f0z+dpjb1G6c3+i03EdzKtrdnH/K+v47FkDmDyid0zW4YXHuXYmLSWZez93Bp1Tk7ny/vks2+Ij3Vzr2HOgmpsfW8qQXhn8978Nb7nDcfLC41w71KdHZx65bgLpnZL51K/mM3el3zLbnRgz4+uPL6Oi6gg/vXI0nVOTY7YuLzzOtVNDemXw569OZEivDK59sJg/LtiY6JRcO/bHBZv4+8odfH3KqYzo2yOm6/LC41w71qtbOo9cN57zTu3Ffz1Zwg+eeycmo5Bcx7Zmx36+++wKJg3tyTVnF8Z8fV54nGvnuqalcP/nx/KZswbwy3+s5ebHlvqts91RO1xTy388spSuqSn8+JOjYjKKraGUmK/BORdzKclJ/O+0keRnduaHz69ix77D3Pf5sfTo3CnRqbk27gfPrWLl9n3MmlFEr27pcVmn7/E410FIYuZHhvCTT42meOMerrjvdbb6vXxcM+atLmPWq+u5asJALhiW13KHVuKFx7kOZtqYfGZfM47tew8x/ZevsWRTeaJTcm3QrsrDfO2xtzglL4NvXjIsruv2wuNcBzRxcE8ev34iKUlJTL/3df7nqRL2HjyS6LRcG2Fm3Pr4MvYdOsI9nx5DeqfYDZ1ujBce5zqoU3t347mbJjFjQgF/XLCRC348j6eXbsUvDHxyMzN++Y+1zH1nJ9+4+DRO69097jl44XGuA+uW3ok7Pj6Cp2eeQ9/MdG58ZClXPbCQ9bsOJDo1lwD7Dh1h5kOL+eHzq/i30/tw9cSChOTht0Vogd8WwXUUtXXGHxds5IfPreJwbR0zzxvC9ecNIi0lvodZXGK8tbmCGx5ezLaKQ3z9olO5dtKgmA6d9tsiOOdIThJXTShg7tc+zEUjenP331dz8U/+yWuluxKdmoshM2PWq+v5xH2vU1cHj315Al/+8OC4/F6nKbG+9fWNkkokLZd0U4hdEV7XSfpANZQ0QFKlpFuiYlMkrZJUKum2qHihpAUh/qik1BBPC69Lw/yCqD7fCPFVki6K5fY71xb16p7Ozz49hgevGUetGZ/9zQJmPLCQf6za6Vc96GAqqqq59sFFfOcvKzjv1F48+x/nMHZgVqLTil3hkTQSuBYYB4wCLpU0BCgBpgOvNNH1LuBvUctJBn4BXAwMBz4tqf6yqXcCd5vZEKAc+GKIfxEoD/G7QztCvyuBEcAU4Jdh+c6ddM49JZfnbzqX/7zoVFZs38fVv32Tj949j9/P30hVdU2i03MnaNHGcv7tnleZt3on//fS4dz/+bFkdklNdFpAbPd4hgELzKzKzGqAecB0M1tpZqsa6yBpGrAeWB4VHgeUmtk6M6sGHgGmShJwPvB4aDcbmBamp4bXhPkXhPZTgUfM7LCZrQdKw/KdOymld0pm5keG8Nqt53P3p0bRNS2F/3mqhPHfm8v/++tKtpRXJTpFd4zq6oz75q3lk796g+Qk8cRXJnLNOYVEvgLbhlheMqcE+F9JOcBB4BKgybP0kjKAW4GPArdEzcoHNke93gKcBeQAFaGo1cfzG/YxsxpJe0P7fGB+g2Xl49xJLjUlicvG9GPa6HwWbyrngdc28JtX1/Prf67johG9+cLZhZxZkNWmvrzc+1XX1PHs29t44NUNvL11L5ec3pvvX/4huqe3vcsmxazwmNlKSXcCLwAHgKVAbTNd7iBy2Kwy0R9uSdcB1wEMGDAgobk4F0+SGDswm7EDs9lWcZAH39jIwws38beSd+mX1ZkLh+Xx0eF5jCvMplOyj01qC3ZVHuahBZv4/fyNlO0/zODcrvzoilFcfkZ+m/1DIaYXCTWzWcAsAEnfI7KH0ZSzgE9I+gGQCdRJOgQsAvpHtesHbAV2A5mSUsJeT32c8Nwf2CIpBegR2m9tYlkN874fuB8iw6mPZZud6yj6ZnbmtotP48YLhvLMW9t4YcW7PLxwE797fQPd0lM479ReXDisF+ed2ssvRpoAK7bt47evrefpt7ZRXVPHh0/J5ZorCpk0pGdCR6wdjZgWHkm9zGynpAFEBhSMb6qtmU2K6ncHUGlmPw+FY6ikQiJF4krgM2Zmkl4GPkHkvM8M4OmwiDnh9Rth/kuh/RzgIUl3AX2BocDCVt1o5zqYzqnJfPLM/nzyzP4crK7l1dJd/H3FDua+s4Nn3tpGSpIYV5jNBcPyGD8om1PzupHie0MxUVNbx9x3dvLb19Yzf90eOndK5pNF/bh6YiFDemUkOr2jFuvbIjwRzvEcAWaaWYWky4CfAbnAs5KWmlmTw5rDOZobgOeBZOABM6sffHAr8Iik7wJLCHtX4fn3kkqBPUSKFWa2XNJjwAqgJuTU3OE/51yUzqnJfHR45HBbXZ2xdEtFpAit3Ml3/rICgC6pyYzql8kZAzM5Y0AWYwZkkd21bYymao827j7AK2t28eqaMl5fu5v9h2rIz+zMNy4+jSvPHECPLu1vb9OvXNACv3KBc0dnS3kVizaWs3hjOYs3VbBi+z5qw++CCnt2ZcyASCEa1qcbQ/O6tcmT3m1BRVU1r6/dzT/X7OLV0jI274nc2iI/szPnDOnJR06LHOJs63uVzV25wAtPC7zwOHd8DlbXsmxLBYs3VbB4UzlLNpWzq7L6vfm9u6czNC+DU/K6MbRXBkPzujE0L+OkKkjlB6pZvWM/q3dWsmbHft7aXMGyrXsxg25pKYwfnMOkoT05Z0hPCnt2bbODBRrjhecEeOFxrnWYGVvKD0a+aHdEvmhX79xP6c5KDh351626e3dPp392Z/pm/uuRn5lOnx6R6e7pKe3qC/hIbR27Kg+ztfwgq3dUsnrHftbs3M+qdyvZVXn4vXYZaSkM69ONs4f0ZNLQnozql9nm92qa44XnBHjhcS626uqiClIoRFvLD7Jt70G2VxyipsFlfDLSUujTI53srqlkdUklq2snMrukktUl8pwdYj06p9I5NZn0lCTSOyWT3imZ5BMc7XWkto6qw7UcqK6hqrqGA2F638EayvYfYse+w+zcf4id+w+zY99hyvYfYveBaqK/ZrukJr+3h3dKXv1zN/r2SG9XBbUlzRWeWA8ucM65ZiUliQE5XRiQ04ULh7//9su1dRbZW6iIFKFtFQcj03sPUn7gCGvLKinfeISKquoPFKjGdEoW6SnJpKcmk94pidSwR2HhP0ZkzyzyDIZRV0ek0Byupbq2rpmlQ5KgZ0Yaed3T6dsjndH9e9CrWzq9uqfRt0dnhvTKID+zc5sf7hxrXnicc21WcpLI655OXvd0aOa33GZG5eEaKqqOUF5VTXlVpBgdOlLLoSN1/3quqQ3TkdfVNaGQCJIkBEiE58jrpCTRNTWZLmkpkefUFLqmvf85Iy2FXt3TyOmadsJ7VScDLzzOuXZPEt3SO9EtvRP9s7skOh3XgvZ75so551y75IXHOedcXHnhcc45F1deeJxzzsWVFx7nnHNx5YXHOedcXHnhcc45F1deeJxzzsWVX6utBZLKgI1NzO4J7IpjOieqveULnnO8eM6x197yhRPLeaCZ5TY2wwvPCZBU3NRF8Nqi9pYveM7x4jnHXnvLF2KXsx9qc845F1deeJxzzsWVF54Tc3+iEzhG7S1f8JzjxXOOvfaWL8QoZz/H45xzLq58j8c551xceeFxzjkXV154okjqL+llSSskLZd0Y4jfIWmrpKXhcUlUn29IKpW0StJFUfEpIVYq6bYY571B0tsht+IQy5b0oqQ14TkrxCXpnpDXMklnRC1nRmi/RtKMGOV6atT7uFTSPkk3tbX3WNIDknZKKomKtdp7Kmls+H9WGvqe8G0rm8j5h5LeCXk9KSkzxAskHYx6v+9rKbemtj8GObfaZ0FSoaQFIf6opNQY5fxoVL4bJC0N8YS/z2r6ey1xn2cz80d4AH2AM8J0N2A1MBy4A7ilkfbDgbeANKAQWAskh8daYBCQGtoMj2HeG4CeDWI/AG4L07cBd4bpS4C/Ebm773hgQYhnA+vCc1aYzorx+50MvAsMbGvvMXAucAZQEov3FFgY2ir0vThGOU8GUsL0nVE5F0S3a7CcRnNravtjkHOrfRaAx4Arw/R9wFdikXOD+T8G/m9beZ9p+nstYZ9n3+OJYmbbzWxxmN4PrATym+kyFXjEzA6b2XqgFBgXHqVmts7MqoFHQtt4mgrMDtOzgWlR8QctYj6QKakPcBHwopntMbNy4EVgSoxzvABYa2ZNXRmiPt+4v8dm9gqwp5FcTvg9DfO6m9l8i/yrfTBqWa2as5m9YGY14eV8oF9zy2ght6a2v1VzbsYxfRbCX93nA4/HK+ewzk8CDze3jHi+z818ryXs8+yFpwmSCoAxwIIQuiHsdj4QteubD2yO6rYlxJqKx4oBL0haJOm6EMszs+1h+l0gL0y3lZwBruT9/0Db8nsMrfee5ofphvFYu4bIX6P1CiUtkTRP0qQQay63prY/Flrjs5ADVEQV3ni8z5OAHWa2JirWZt7nBt9rCfs8e+FphKQM4AngJjPbB9wLDAZGA9uJ7Eq3JeeY2RnAxcBMSedGzwx/hbSpcfPhWPvHgT+FUFt/j9+nLb6nzZH0X0AN8McQ2g4MMLMxwM3AQ5K6H+3yYrz97eqz0MCnef8fU23mfW7key0m6zkaXngakNSJyP+cP5rZnwHMbIeZ1ZpZHfBrIrv2AFuB/lHd+4VYU/GYMLOt4Xkn8GTIb0fYBa7frd/ZlnImUiQXm9mOkHubfo+D1npPt/L+Q14xzV3S1cClwGfDFwzhcNXuML2IyDmSU1rIrantb1Wt+FnYTeQwUUoj29LqwnqmA4/Wx9rK+9zY91oz64n559kLT5RwfHYWsNLM7oqK94lqdhlQP5plDnClpDRJhcBQIifZ3gSGhhE1qUQOKc2JUc5dJXWrnyZyMrkkrK9+1MkM4OmonK8KI1fGA3vD7vbzwGRJWeHQxuQQi5X3/WXYlt/jKK3ynoZ5+ySND5+5q6KW1aokTQG+DnzczKqi4rmSksP0ICLv67oWcmtq+1s751b5LIQi+zLwiVjnHFwIvGNm7x12agvvc1Pfa82sJ/af5+ZGHpxsD+AcIruby4Cl4XEJ8Hvg7RCfA/SJ6vNfRP6KWUXUSI7Qb3WY918xzHkQkVE8bwHL69dF5Pj2XGAN8HcgO8QF/CLk9TZQFLWsa4icsC0FvhDDnLsS+Wu0R1SsTb3HRIriduAIkWPWX2zN9xQoIvKFuhb4OeEqIjHIuZTIcfn6z/N9oe3l4fOyFFgMfKyl3Jra/hjk3GqfhfDvY2F4H/4EpMUi5xD/HXB9g7YJf59p+nstYZ9nv2SOc865uPJDbc455+LKC49zzrm48sLjnHMurrzwOOeciysvPM455+LKC49zbYykTElfDdN9JT3eUh/n2hMfTu1cGxOup/UXMxuZ4FSci4mUlps45+Ls+8BgRe7psgYYZmYjw6VvphH5Ae5Q4EdEbgPweeAwcImZ7ZE0mMgPAHOBKuBaM3sn/pvhXOP8UJtzbc9tRG4XMRr4zwbzRhK5HtiZwP8CVRa5AOUbRC5VAnA/8O9mNha4BfhlXLJ27ij5Ho9z7cvLFrmnyn5Je4FnQvxt4EPhCsQTgT/pXzeBTIt/ms41zQuPc+3L4ajpuqjXdUT+PScRuQfN6Hgn5tzR8kNtzrU9+4ncoviYWeQ+K+slXQGRKxNLGtWayTl3orzwONfGWOT+La9JKgF+eByL+CzwRUn1VyyP923XnWuWD6d2zjkXV77H45xzLq688DjnnIsrLzzOOefiyguPc865uPLC45xzLq688DjnnIsrLzzOOefi6v8DoAFoPOX8WFoAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "#import seaborn as sns\n", - "#sns.lmplot('a', 'e', data=disk, fit_reg=False)\n", - "disk['a'].isel(id=1502).plot.line(x='time')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [ - { - "ename": "SyntaxError", - "evalue": "positional argument follows keyword argument (, line 3)", - "output_type": "error", - "traceback": [ - "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m3\u001b[0m\n\u001b[0;31m disk.isel(time=0).plot.scatter(x='a', y='e', markersize='radmarker', size_norm)\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m positional argument follows keyword argument\n" - ] - } - ], - "source": [ - "fig = plt.figure(1, figsize=(10,6))\n", - "ax = fig.add_subplot(111)\n", - "disk.isel(time=0).plot.scatter(x='a', y='e', markersize='radmarker', size_norm)\n", - "plt.rcParams.update({'font.size': 18})\n", - "ax.set_yscale('log')\n", - "plt.ylim((1e-4,1.0))" - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Reading Swiftest file /Users/daminton/work/Projects/Swiftest/Elliott_Performance/high_high_1500_1/param.in\n" - ] - }, - { - "ename": "NameError", - "evalue": "name 'FortranFile' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mconfig\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mswio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_swiftest_config\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconfig_file_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mconfig\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'BIN_OUT'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mworkingdir\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mconfig\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'BIN_OUT'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0;32mwith\u001b[0m \u001b[0mFortranFile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mconfig\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'BIN_OUT'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 6\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnpl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mplid\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpvec\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mplab\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mntp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtpid\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtvec\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtlab\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mswiftest_stream\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mconfig\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mNameError\u001b[0m: name 'FortranFile' is not defined" - ] - } - ], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "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.8.6" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/python/swiftestio/orbel.f90 b/python/swiftestio/orbel.f90 deleted file mode 100644 index 83e52db1c..000000000 --- a/python/swiftestio/orbel.f90 +++ /dev/null @@ -1,160 +0,0 @@ -! Created by on 6/22/21. - -module orbel - implicit none - -contains - - subroutine xv2el(mu, x, v, a, e, inc, capom, omega, capm) - !! author: David A. Minton - !! - !! Compute osculating orbital elements from relative Cartesian position and velocity - !! All angular measures are returned in radians - !! If inclination < TINY, longitude of the ascending node is arbitrarily set to 0 - !! - !! If eccentricity < sqrt(TINY), argument of pericenter is arbitrarily set to 0 - !! - !! References: Danby, J. M. A. 1988. Fundamentals of Celestial Mechanics, (Willmann-Bell, Inc.), 201 - 206. - !! Fitzpatrick, P. M. 1970. Principles of Celestial Mechanics, (Academic Press), 69 - 73. - !! Roy, A. E. 1982. Orbital Motion, (Adam Hilger, Ltd.), 75 - 95 - !! - !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2el.f90 - !! Adapted from Martin Duncan's Swift routine orbel_xv2el.f - implicit none - ! Arguments - real*8, intent(in) :: mu - real*8, dimension(:), intent(in) :: x, v - real*8, intent(out) :: a, e, inc, capom, omega, capm - !f2py intent(in) mu - !f2py intent(in) x - !f2py intent(in) v - !f2py intent(out) a - !f2py intent(out) e - !f2py intent(out) inc - !f2py intent(out) capom - !f2py intent(out) omega - !f2py intent(out) capm - - ! Internals - - integer :: iorbit_type - real*8 :: r, v2, h2, h, rdotv, energy, fac, u, w, cw, sw, face, cape, tmpf, capf - real*8, dimension(3) :: hvec - integer, parameter :: ELLIPSE = -1 !! Symbolic names for orbit types - ellipse - integer, parameter :: PARABOLA = 0 !! Symbolic names for orbit types - parabola - integer, parameter :: HYPERBOLA = 1 !! Symbolic names for orbit types - hyperbola - integer, parameter :: DP = kind(1.d0) - real*8, parameter :: VSMALL = 4.0E-15_DP - real*8, parameter :: PI = 3.141592653589793238462643383279502884197_DP !! Definition of /(\pi\) - real*8, parameter :: TWOPI = 6.283185307179586476925286766559005768394_DP !! Definition of 2 \pi - - a = 0.0_DP - e = 0.0_DP - inc = 0.0_DP - capom = 0.0_DP - omega = 0.0_DP - capm = 0.0_DP - r = sqrt(dot_product(x(:), x(:))) - v2 = dot_product(v(:), v(:)) - hvec = crossproduct(x(:), v(:)) - h2 = dot_product(hvec(:), hvec(:)) - h = sqrt(h2) - if (h2 == 0.0_DP) return - rdotv = dot_product(x(:), v(:)) - energy = 0.5_DP * v2 - mu / r - fac = hvec(3) / h - if (fac < -1.0_DP) then - inc = PI - else if (fac < 1.0_DP) then - inc = acos(fac) - end if - fac = sqrt(hvec(1)**2 + hvec(2)**2) / h - if (fac**2 < VSMALL) then - u = atan2(x(2), x(1)) - if (hvec(3) < 0.0_DP) u = -u - else - capom = atan2(hvec(1), -hvec(2)) - u = atan2(x(3) / sin(inc), x(1) * cos(capom) + x(2) * sin(capom)) - end if - if (capom < 0.0_DP) capom = capom + TWOPI - if (u < 0.0_DP) u = u + TWOPI - if (abs(energy * r / mu) < sqrt(VSMALL)) then - iorbit_type = PARABOLA - else - a = -0.5_DP * mu / energy - if (a < 0.0_DP) then - fac = -h2 / (mu * a) - if (fac > VSMALL) then - iorbit_type = HYPERBOLA - else - iorbit_type = PARABOLA - end if - else - iorbit_type = ELLIPSE - end if - end if - select case (iorbit_type) - case (ELLIPSE) - fac = 1.0_DP - h2 / (mu * a) - if (fac > VSMALL) then - e = sqrt(fac) - cape = 0.0_DP - face = (a - r) / (a * e) - if (face < -1.0_DP) then - cape = PI - else if (face < 1.0_DP) then - cape = acos(face) - end if - if (rdotv < 0.0_DP) cape = TWOPI - cape - fac = 1.0_DP - e * cos(cape) - cw = (cos(cape) - e) / fac - sw = sqrt(1.0_DP - e**2) * sin(cape) / fac - w = atan2(sw, cw) - if (w < 0.0_DP) w = w + TWOPI - else - cape = u - w = u - end if - capm = cape - e * sin(cape) - case (PARABOLA) - a = 0.5_DP * h2 / mu - e = 1.0_DP - w = 0.0_DP - fac = 2 * a / r - 1.0_DP - if (fac < -1.0_DP) then - w = PI - else if (fac < 1.0_DP) then - w = acos(fac) - end if - if (rdotv < 0.0_DP) w = TWOPI - w - tmpf = tan(0.5_DP * w) - capm = tmpf * (1.0_DP + tmpf * tmpf / 3.0_DP) - case (HYPERBOLA) - e = sqrt(1.0_DP + fac) - tmpf = max((a - r) / (a * e), 1.0_DP) - capf = log(tmpf + sqrt(tmpf**2 - 1.0_DP)) - if (rdotv < 0.0_DP) capf = -capf - fac = e * cosh(capf) - 1.0_DP - cw = (e - cosh(capf)) / fac - sw = sqrt(e * e - 1.0_DP) * sinh(capf) / fac - w = atan2(sw, cw) - if (w < 0.0_DP) w = w + TWOPI - capm = e * sinh(capf) - capf - end select - omega = u - w - if (omega < 0.0_DP) omega = omega + TWOPI - - return - contains - function crossproduct(A, B) result(C) - implicit none - real*8, dimension(:), intent(in) :: A, B - real*8, dimension(3) :: C - C(1) = A(2) * B(3) - A(3) * B(2) - C(2) = A(3) * B(1) - A(1) * B(3) - C(3) = A(1) * B(2) - A(2) * B(1) - return - end function crossproduct - end subroutine xv2el - -end module orbel \ No newline at end of file diff --git a/python/swiftestio/swiftestio.py b/python/swiftestio/swiftestio.py index 97f47841c..473c2bd5d 100644 --- a/python/swiftestio/swiftestio.py +++ b/python/swiftestio/swiftestio.py @@ -6,6 +6,7 @@ import astropy.constants as const import datetime import sys +from swiftestpy import orbel # Constants in SI units AU2M = np.longdouble(const.au.value) @@ -304,6 +305,12 @@ def make_swiftest_labels(config): tlab.append('vx') tlab.append('vy') tlab.append('vz') + tlab.append('a') + tlab.append('e') + tlab.append('inc') + tlab.append('capom') + tlab.append('omega') + tlab.append('capm') elif config['OUT_FORM'] == 'EL': tlab.append('a') tlab.append('e') @@ -394,16 +401,72 @@ def swiftest_stream(f, config): plab, tlab = make_swiftest_labels(config) + if config['OUT_FORM'] == 'XV': + mu = np.empty_like(p1) + mu[0] = Mpl[0] + mu[1:] = Mpl[0] + Mpl[1:] + p7 = [] + p8 = [] + p9 = [] + p10 = [] + p11 = [] + p12 = [] + for i in range(mu.size): + elem = orbel.xv2el(mu[i], p1[i], p2[i], p3[i], p4[i], p5[i], p6[i]) + p7.append(elem[0]) + p8.append(elem[1]) + p9.append(elem[2]) + p10.append(elem[3]) + p11.append(elem[4]) + p12.append(elem[5]) + p7 = np.array(p7) + p8 = np.array(p8) + p9 = np.array(p9) + p10 = np.array(p10) + p11 = np.array(p11) + p12 = np.array(p12) + if ntp[0] > 0: + mu = np.full_like(t1,Mpl[0]) + t7 = [] + t8 = [] + t9 = [] + t10 = [] + t11 = [] + t12 = [] + for i in range(mu.size): + elem = orbel.xv2el(mu[i], t1[i], t2[i], t3[i], t4[i], t5[i], t6[i]) + t7.append(elem[0]) + t8.append(elem[1]) + t9.append(elem[2]) + t10.append(elem[3]) + t11.append(elem[4]) + t12.append(elem[5]) + t7 = np.array(t7) + t8 = np.array(t8) + t9 = np.array(t9) + t10 = np.array(t10) + t11 = np.array(t11) + t12 = np.array(t12) + if npl > 0: if config['ROTATION'] == 'YES': - pvec = np.vstack([p1,p2,p3,p4,p5,p6,Mpl,Rpl,rot_x,rot_y,rot_z,Ip_1,Ip_2,Ip_3]) + if config['OUT_FORM'] == 'EL': + pvec = np.vstack([p1,p2,p3,p4,p5,p6,Mpl,Rpl,rot_x,rot_y,rot_z,Ip_1,Ip_2,Ip_3]) + else: + pvec = np.vstack([p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, Mpl, Rpl, rot_x, rot_y, rot_z, Ip_1, Ip_2, Ip_3]) else: - pvec = np.vstack([p1,p2,p3,p4,p5,p6,Mpl,Rpl]) + if config['OUT_FORM'] == 'EL': + pvec = np.vstack([p1,p2,p3,p4,p5,p6,Mpl,Rpl]) + else: + pvec = np.vstack([p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, Mpl, Rpl]) else: pvec = np.empty((8,0)) plid = np.empty(0) if ntp > 0: - tvec = np.vstack([t1,t2,t3,t4,t5,t6]) + if config['OUT_FORM'] == 'EL': + tvec = np.vstack([t1,t2,t3,t4,t5,t6]) + else: + tvec = np.vstack([t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12]) else: tvec = np.empty((6,0)) tpid = np.empty(0) @@ -771,3 +834,4 @@ def swiftest_xr2_infile(ds, config, framenum=-1): swiftestdat = swiftest2xr(config) + diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index 12ec7d02b..c26078b60 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -673,7 +673,7 @@ END SUBROUTINE orbel_xv2aqt END INTERFACE INTERFACE - SUBROUTINE orbel_xv2el(x, v, mu, a, e, inc, capom, omega, capm) + pure SUBROUTINE orbel_xv2el(x, v, mu, a, e, inc, capom, omega, capm) USE swiftest_globals IMPLICIT NONE REAL(DP), INTENT(IN) :: mu diff --git a/src/orbel/orbel_xv2el.f90 b/src/orbel/orbel_xv2el.f90 index 38e69b36a..284a25e3e 100644 --- a/src/orbel/orbel_xv2el.f90 +++ b/src/orbel/orbel_xv2el.f90 @@ -40,7 +40,7 @@ ! Roy, A. E. 1982. Orbital Motion, (Adam Hilger, Ltd.), 75 - 95. ! !********************************************************************************************************************************** -SUBROUTINE orbel_xv2el(x, v, mu, a, e, inc, capom, omega, capm) +pure SUBROUTINE orbel_xv2el(x, v, mu, a, e, inc, capom, omega, capm) ! Modules USE swiftest diff --git a/src/python/orbel.f90 b/src/python/orbel.f90 new file mode 100644 index 000000000..4dd26a68d --- /dev/null +++ b/src/python/orbel.f90 @@ -0,0 +1,32 @@ +module orbel + use swiftest + private + public :: xv2el +contains + pure elemental subroutine xv2el(mu, px, py, pz, vx, vy, vz, a, e, inc, capom, omega, capm) + use module_interfaces, only : orbel_xv2el + implicit none + ! Arguments + real*8, intent(in) :: mu, px, py, pz, vx, vy, vz + real*8, intent(out) :: a, e, inc, capom, omega, capm + !$f2py intent(in) mu + !$f2py intent(in) px + !$f2py intent(in) py + !$f2py intent(in) pz + !$f2py intent(in) vx + !$f2py intent(in) vy + !$f2py intent(in) vz + !$f2py intent(out) a + !$f2py intent(out) e + !$f2py intent(out) inc + !$f2py intent(out) capom + !$f2py intent(out) omega + !$f2py intent(out) capm + ! Internals + real*8, dimension(3) :: x, v + x = [px, py, pz] + v = [vx, vy, vz] + call orbel_xv2el(x, v, mu, a, e, inc, capom, omega, capm) + return + end subroutine xv2el +end module orbel \ No newline at end of file