|
| 1 | +from pathlib import Path |
| 2 | + |
| 3 | +import matplotlib.pyplot as plt |
| 4 | +from matplotlib.animation import FuncAnimation |
| 5 | + |
| 6 | +import sdf_xarray as sdfxr |
| 7 | + |
| 8 | +input_dir = Path("datasets/2_1_two_stream_instability") |
| 9 | +ds = sdfxr.open_mfdataset( |
| 10 | + input_dir, data_vars=["dist_fn_x_px_Left", "dist_fn_x_px_Right"] |
| 11 | +) |
| 12 | + |
| 13 | +# Rescale coords to account for kilometers |
| 14 | +ds = ds.epoch.rescale_coords(1e-3, "km", ["X_x_px_Left"]) |
| 15 | + |
| 16 | +# Sum phase-space of species "Left" and "Right" in "x_px" distribution function |
| 17 | +# NOTE: We only use the values from the right distribution function as if we inherit |
| 18 | +# the coords we from the right we end up with 4 coords instead of 2 |
| 19 | +total_phase_space = ds["dist_fn_x_px_Left"] + ds["dist_fn_x_px_Right"].values |
| 20 | +total_phase_space.attrs["long_name"] = "Phase Space Distribution" |
| 21 | +total_phase_space.attrs["units"] = "kg.m/s" |
| 22 | + |
| 23 | +fig, ax = plt.subplots() |
| 24 | + |
| 25 | +# construct first frame |
| 26 | +# Data needs to be transposed so that axes match |
| 27 | +total_phase_space = total_phase_space.T |
| 28 | +col_min = total_phase_space.values.min() |
| 29 | +col_max = total_phase_space.values.max() |
| 30 | +gif_frame = ax.pcolormesh( |
| 31 | + total_phase_space["X_x_px_Left"], |
| 32 | + total_phase_space["Px_x_px_Left"], |
| 33 | + total_phase_space.isel(time=0), |
| 34 | + vmin=col_min, |
| 35 | + vmax=col_max, |
| 36 | +) |
| 37 | +ax.set_xlabel( |
| 38 | + f"{total_phase_space['X_x_px_Left'].attrs['long_name']} [{total_phase_space['X_x_px_Left'].attrs['units']}]" |
| 39 | +) |
| 40 | +ax.set_ylabel( |
| 41 | + f"{total_phase_space.attrs['long_name']} [{total_phase_space.attrs['units']}]" |
| 42 | +) |
| 43 | +ax.set_title(f"t = {ds['time'].values[0]:.3f} [{ds['time'].attrs['units']}]") |
| 44 | +cbar = fig.colorbar(gif_frame, ax=ax) |
| 45 | +cbar.set_label("Summed particle weight per bin") |
| 46 | + |
| 47 | + |
| 48 | +# Make GIF |
| 49 | +def update(i): |
| 50 | + gif_frame.set_array(total_phase_space.isel(time=i)) |
| 51 | + ax.set_title(f"t = {ds['time'].values[i]:.3f} [{ds['time'].attrs['units']}]") |
| 52 | + return (gif_frame,) |
| 53 | + |
| 54 | + |
| 55 | +anim = FuncAnimation(fig, update, frames=ds.sizes["time"]) |
| 56 | +anim.save(input_dir / "phase_space.gif") |
0 commit comments