Sampling data from a 3D model

In this example we’ll see how to sample a 3D model output at arbitrary points within the model domain. This example uses a model with a single timestep, so the time at which data is interpolated doesn’t need to be specified.

First, load the required modules.

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from psipy.data import sample_data
from psipy.model import PLUTOOutput

Load a set of MAS output files, and get the number density variable from the model run. This example works with PLUTO data too. Uncomment the comment block below (and comment out the MAS block!)

"""
from psipy.model import MASOutput
mas_path = sample_data.mas_sample_data()
model = MASOutput(mas_path)
"""
pluto_path = sample_data.pluto_sample_data()
model = PLUTOOutput(pluto_path)
rho = model["rho"]

  0%|                                              | 0.00/18.3k [00:00<?, ?B/s]
  0%|                                              | 0.00/18.3k [00:00<?, ?B/s]
100%|█████████████████████████████████████| 18.3k/18.3k [00:00<00:00, 46.1MB/s]

  0%|                                                | 0.00/200 [00:00<?, ?B/s]
  0%|                                                | 0.00/200 [00:00<?, ?B/s]
100%|██████████████████████████████████████████| 200/200 [00:00<00:00, 484kB/s]

  0%|                                              | 0.00/16.0M [00:00<?, ?B/s]
  0%|                                      | 52.2k/16.0M [00:00<00:45, 347kB/s]
  2%|▌                                      | 244k/16.0M [00:00<00:17, 883kB/s]
  5%|█▋                                    | 731k/16.0M [00:00<00:07, 1.94MB/s]
 16%|█████▉                               | 2.57M/16.0M [00:00<00:02, 5.94MB/s]
 36%|█████████████▏                       | 5.71M/16.0M [00:00<00:00, 11.2MB/s]
 55%|████████████████████▍                | 8.88M/16.0M [00:00<00:00, 14.5MB/s]
 75%|███████████████████████████▊         | 12.1M/16.0M [00:01<00:00, 16.5MB/s]
 95%|███████████████████████████████████▏ | 15.2M/16.0M [00:01<00:00, 17.9MB/s]
  0%|                                              | 0.00/16.0M [00:00<?, ?B/s]
100%|█████████████████████████████████████| 16.0M/16.0M [00:00<00:00, 31.6GB/s]

  0%|                                              | 0.00/16.0M [00:00<?, ?B/s]
  0%|                                      | 51.2k/16.0M [00:00<00:47, 337kB/s]
  2%|▋                                      | 260k/16.0M [00:00<00:16, 935kB/s]
  7%|██▌                                  | 1.13M/16.0M [00:00<00:04, 3.09MB/s]
 23%|████████▌                            | 3.73M/16.0M [00:00<00:01, 8.55MB/s]
 41%|███████████████                      | 6.50M/16.0M [00:00<00:00, 11.9MB/s]
 57%|█████████████████████▏               | 9.19M/16.0M [00:00<00:00, 13.8MB/s]
 77%|████████████████████████████▌        | 12.4M/16.0M [00:01<00:00, 16.0MB/s]
 97%|███████████████████████████████████▉ | 15.5M/16.0M [00:01<00:00, 17.5MB/s]
  0%|                                              | 0.00/16.0M [00:00<?, ?B/s]
100%|█████████████████████████████████████| 16.0M/16.0M [00:00<00:00, 33.5GB/s]

  0%|                                              | 0.00/16.0M [00:00<?, ?B/s]
  0%|▏                                     | 54.3k/16.0M [00:00<00:45, 354kB/s]
  1%|▌                                      | 234k/16.0M [00:00<00:18, 833kB/s]
  7%|██▍                                  | 1.07M/16.0M [00:00<00:05, 2.92MB/s]
 23%|████████▋                            | 3.74M/16.0M [00:00<00:01, 8.57MB/s]
 41%|███████████████                      | 6.51M/16.0M [00:00<00:00, 12.9MB/s]
 53%|███████████████████▍                 | 8.44M/16.0M [00:00<00:00, 14.2MB/s]
 69%|█████████████████████████▋           | 11.1M/16.0M [00:00<00:00, 16.6MB/s]
 82%|██████████████████████████████▎      | 13.1M/16.0M [00:01<00:00, 17.1MB/s]
 98%|████████████████████████████████████▍| 15.8M/16.0M [00:01<00:00, 18.5MB/s]
  0%|                                              | 0.00/16.0M [00:00<?, ?B/s]
100%|█████████████████████████████████████| 16.0M/16.0M [00:00<00:00, 34.5GB/s]

  0%|                                              | 0.00/16.0M [00:00<?, ?B/s]
  0%|                                      | 52.2k/16.0M [00:00<00:46, 344kB/s]
  2%|▋                                      | 261k/16.0M [00:00<00:16, 959kB/s]
  7%|██▋                                  | 1.15M/16.0M [00:00<00:04, 3.18MB/s]
 23%|████████▋                            | 3.75M/16.0M [00:00<00:01, 8.61MB/s]
 43%|███████████████▉                     | 6.91M/16.0M [00:00<00:00, 12.9MB/s]
 63%|███████████████████████▏             | 10.1M/16.0M [00:00<00:00, 15.5MB/s]
 82%|██████████████████████████████▌      | 13.2M/16.0M [00:01<00:00, 17.1MB/s]
  0%|                                              | 0.00/16.0M [00:00<?, ?B/s]
100%|█████████████████████████████████████| 16.0M/16.0M [00:00<00:00, 31.4GB/s]

Choose a set of 1D points to interpolate the model output at.

Here we keep a constant radius, and a set of longitudes that go all the way from 0 to 360 degrees. Then we choose two different, but close latitude values, and plot the results.

As expected, the values at 0 and 360 degrees are the same, and the profiles are similar, but different, due to the small difference in latitude.

fig, ax = plt.subplots()

npoints = 1000
r = 50 * np.ones(npoints) * u.R_sun
lon = np.linspace(0, 360, npoints) * u.deg

for latitude in [0, 1] * u.deg:
    lat = latitude * np.ones(npoints)
    samples = rho.sample_at_coords(lon, lat, r)

    ax.plot(lon, samples, label="lat = " + str(latitude))

ax.legend()
ax.set_xlim(0, 360)
ax.set_ylim(bottom=0)
ax.set_xlabel("Longitude (deg)")
ax.set_ylabel(r"$\rho$ (cm$^{-3}$)")
ax.set_xticks([0, 90, 180, 270, 360])
ax.set_title(f"r={r[0].to_value(u.R_sun):.02f} R_sun")
plt.show()