lisa
404049b068
This code allows you to visualize the position, trajectory of all spacecraft, and Key Science Regions of PMO in 3D. One .npy file is loaded per spacecraft. To plot the trajectory over a time interval, the variables start and end must correspond to the desired time interval. To plot a discrete point in time, set start=None and end=None, and only target_time receives the desired timestamp.
293 lines
No EOL
8.5 KiB
Python
293 lines
No EOL
8.5 KiB
Python
#!/usr/bin/env python3
|
|
# -*- coding: utf-8 -*-
|
|
"""
|
|
Created on Mon Jun 1 09:00:23 2026
|
|
|
|
@author: vieweg
|
|
"""
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from datetime import datetime
|
|
|
|
|
|
# load data
|
|
|
|
def load_data(file):
|
|
all_data = np.load(file, allow_pickle=True).item()
|
|
dates = np.array([datetime.fromisoformat(str(d)) for d in all_data['dates']])
|
|
position = all_data['position'].astype(np.float64)
|
|
|
|
x_loaded = position[:, 0]
|
|
y_loaded = position[:, 1]
|
|
z_loaded = position[:, 2]
|
|
|
|
return dates, x_loaded, y_loaded, z_loaded
|
|
|
|
|
|
# cylinder
|
|
|
|
def plot_cylinder(ax, x0, y0, z0, radius=0.2, height=0.05, color='blue'):
|
|
theta = np.linspace(0, 2*np.pi, 20)
|
|
z = np.linspace(-height/2, height/2, 2)
|
|
theta, z = np.meshgrid(theta, z)
|
|
|
|
x = x0 + radius * np.cos(theta)
|
|
y = y0 + radius * np.sin(theta)
|
|
z = z0 + z
|
|
|
|
ax.plot_surface(x, y, z, color=color)
|
|
|
|
# Deckel + Boden
|
|
theta = np.linspace(0, 2*np.pi, 20)
|
|
r = np.linspace(0, radius, 10)
|
|
theta, r = np.meshgrid(theta, r)
|
|
|
|
x_top = x0 + r * np.cos(theta)
|
|
y_top = y0 + r * np.sin(theta)
|
|
z_top = np.full_like(x_top, z0 + height/2)
|
|
|
|
x_bottom = x0 + r * np.cos(theta)
|
|
y_bottom = y0 + r * np.sin(theta)
|
|
z_bottom = np.full_like(x_bottom, z0 - height/2)
|
|
|
|
ax.plot_surface(x_top, y_top, z_top, color=color)
|
|
ax.plot_surface(x_bottom, y_bottom, z_bottom, color=color)
|
|
return
|
|
|
|
def plot_3d_multi(datasets, start=None, end=None, target_time=None):
|
|
fig = plt.figure()
|
|
ax = fig.add_subplot(projection='3d')
|
|
# Make data
|
|
u = np.linspace(0, 2 * np.pi, 100)
|
|
v = np.linspace(0, np.pi, 100)
|
|
xe = 1 * np.outer(np.cos(u), np.sin(v))
|
|
ye = 1 * np.outer(np.sin(u), np.sin(v))
|
|
ze = 1 * np.outer(np.ones(np.size(u)), np.cos(v))
|
|
|
|
# Plot the surface
|
|
ax.plot_surface(xe, ye, ze, color = 'b', label = 'Earth')
|
|
|
|
#Plot Sun at (10, 0, 0)
|
|
r = 1.5
|
|
u = np.linspace(0, 2 * np.pi, 50)
|
|
v = np.linspace(0, np.pi, 50)
|
|
|
|
x0, y0, z0 = 20, 0 , 0
|
|
|
|
ax.quiver(x0, y0, z0, 1.0,0.0,0.0, color = 'orange', length = 5, normalize = True, linewidth = 2, zorder = 20, label = 'Direction to Sun')
|
|
|
|
|
|
|
|
#plot PMO key science region - Fore Shock
|
|
r1 = 13.5
|
|
r2 = 17
|
|
theta_max = np.deg2rad(45)
|
|
|
|
# Winkel definieren
|
|
theta = np.linspace(0, theta_max, 40)
|
|
phi = np.linspace(0, 2*np.pi, 60)
|
|
theta, phi = np.meshgrid(theta, phi)
|
|
|
|
# outer surface
|
|
x_outer = r2 * np.cos(theta)
|
|
y_outer = r2 * np.sin(theta) * np.cos(phi)
|
|
z_outer = r2 * np.sin(theta) * np.sin(phi)
|
|
|
|
# inner surface
|
|
x_inner = r1 * np.cos(theta)
|
|
y_inner = r1 * np.sin(theta) * np.cos(phi)
|
|
z_inner = r1 * np.sin(theta) * np.sin(phi)
|
|
|
|
# Plot Shell
|
|
ax.plot_surface(x_outer, y_outer, z_outer, color='yellow', alpha=0.3)
|
|
ax.plot_surface(x_inner, y_inner, z_inner, color='yellow', alpha=0.3)
|
|
|
|
r = np.linspace(r1, r2, 30)
|
|
phi = np.linspace(0, 2*np.pi, 60)
|
|
r, phi = np.meshgrid(r, phi)
|
|
|
|
x_cap = r * np.cos(theta_max)
|
|
y_cap = r * np.sin(theta_max) * np.cos(phi)
|
|
z_cap = r * np.sin(theta_max) * np.sin(phi)
|
|
|
|
ax.plot_surface(x_cap, y_cap, z_cap, color='yellow', alpha=0.3, label = r'Fore Shock ($13.5 < R_0 < 17$)')
|
|
|
|
#plot PMO key science region - Shock
|
|
r1 = 12.5
|
|
r2 = 13.5
|
|
theta_max = np.deg2rad(45)
|
|
|
|
# define angle
|
|
theta = np.linspace(0, theta_max, 40)
|
|
phi = np.linspace(0, 2*np.pi, 60)
|
|
theta, phi = np.meshgrid(theta, phi)
|
|
|
|
# outer surface
|
|
x_outer = r2 * np.cos(theta)
|
|
y_outer = r2 * np.sin(theta) * np.cos(phi)
|
|
z_outer = r2 * np.sin(theta) * np.sin(phi)
|
|
|
|
# inner surface
|
|
x_inner = r1 * np.cos(theta)
|
|
y_inner = r1 * np.sin(theta) * np.cos(phi)
|
|
z_inner = r1 * np.sin(theta) * np.sin(phi)
|
|
|
|
# plot shell
|
|
ax.plot_surface(x_outer, y_outer, z_outer, color='orange', alpha=0.3)
|
|
ax.plot_surface(x_inner, y_inner, z_inner, color='orange', alpha=0.3)
|
|
|
|
r = np.linspace(r1, r2, 30)
|
|
phi = np.linspace(0, 2*np.pi, 60)
|
|
r, phi = np.meshgrid(r, phi)
|
|
|
|
x_cap = r * np.cos(theta_max)
|
|
y_cap = r * np.sin(theta_max) * np.cos(phi)
|
|
z_cap = r * np.sin(theta_max) * np.sin(phi)
|
|
|
|
ax.plot_surface(x_cap, y_cap, z_cap, color='orange', alpha=0.3, label = r'Shock ($12.5 < R_0 < 13.5$)')
|
|
|
|
#plot PMO key science region - Magnetosheath
|
|
r1 = 10.5
|
|
r2 = 12.5
|
|
theta_max = np.deg2rad(45)
|
|
|
|
# define angle
|
|
theta = np.linspace(0, theta_max, 40)
|
|
phi = np.linspace(0, 2*np.pi, 60)
|
|
theta, phi = np.meshgrid(theta, phi)
|
|
|
|
# outer surface
|
|
x_outer = r2 * np.cos(theta)
|
|
y_outer = r2 * np.sin(theta) * np.cos(phi)
|
|
z_outer = r2 * np.sin(theta) * np.sin(phi)
|
|
|
|
# inner surface
|
|
x_inner = r1 * np.cos(theta)
|
|
y_inner = r1 * np.sin(theta) * np.cos(phi)
|
|
z_inner = r1 * np.sin(theta) * np.sin(phi)
|
|
|
|
# Plot shell
|
|
ax.plot_surface(x_outer, y_outer, z_outer, color='red', alpha=0.3)
|
|
ax.plot_surface(x_inner, y_inner, z_inner, color='red', alpha=0.3)
|
|
|
|
r = np.linspace(r1, r2, 30)
|
|
phi = np.linspace(0, 2*np.pi, 60)
|
|
r, phi = np.meshgrid(r, phi)
|
|
|
|
x_cap = r * np.cos(theta_max)
|
|
y_cap = r * np.sin(theta_max) * np.cos(phi)
|
|
z_cap = r * np.sin(theta_max) * np.sin(phi)
|
|
|
|
ax.plot_surface(x_cap, y_cap, z_cap, color='red', alpha=0.3, label = r'Magnetosheath ($10.5 < R_0 < 12.5$)')
|
|
|
|
#plot PMO key science region - Magnetopause
|
|
r1 = 10.5
|
|
r2 = 9.5
|
|
theta_max = np.deg2rad(45)
|
|
|
|
# define angle
|
|
theta = np.linspace(0, theta_max, 40)
|
|
phi = np.linspace(0, 2*np.pi, 60)
|
|
theta, phi = np.meshgrid(theta, phi)
|
|
|
|
# outer surface
|
|
x_outer = r2 * np.cos(theta)
|
|
y_outer = r2 * np.sin(theta) * np.cos(phi)
|
|
z_outer = r2 * np.sin(theta) * np.sin(phi)
|
|
|
|
# inner surface
|
|
x_inner = r1 * np.cos(theta)
|
|
y_inner = r1 * np.sin(theta) * np.cos(phi)
|
|
z_inner = r1 * np.sin(theta) * np.sin(phi)
|
|
|
|
# Plot shell
|
|
ax.plot_surface(x_outer, y_outer, z_outer, color='pink', alpha=0.3)
|
|
ax.plot_surface(x_inner, y_inner, z_inner, color='pink', alpha=0.3)
|
|
|
|
r = np.linspace(r1, r2, 30)
|
|
phi = np.linspace(0, 2*np.pi, 60)
|
|
r, phi = np.meshgrid(r, phi)
|
|
|
|
x_cap = r * np.cos(theta_max)
|
|
y_cap = r * np.sin(theta_max) * np.cos(phi)
|
|
z_cap = r * np.sin(theta_max) * np.sin(phi)
|
|
|
|
ax.plot_surface(x_cap, y_cap, z_cap, color='pink', alpha=0.3, label = r'Magnetopause ($9.5 < R_0 < 10.5$)')
|
|
|
|
#plot PMO key science regions - Transition Region (TR)
|
|
ax.bar3d(x=-12, y = -8, z=-7, dx = 7, dy = 16, dz = 14, color = 'green', alpha = 0.4, label = r'Transition Region' )
|
|
|
|
|
|
#plot PMO key science regions - Magnetotail current sheets (CS)
|
|
ax.bar3d(x=-17, y = -11, z=-10, dx = 5, dy = 22, dz = 20, color = 'purple', alpha = 0.4, label = r'Magnetotail Current Sheet' )
|
|
|
|
|
|
|
|
|
|
|
|
colors = ['red', 'green', 'blue', 'purple']
|
|
labels = ['C1', 'C2', 'C3', 'C4']
|
|
for i in range(len(labels)):
|
|
ax.plot([],[],[], color = colors[i], label = labels[i])
|
|
|
|
|
|
for i, data in enumerate(datasets):
|
|
dates, x_loaded, y_loaded, z_loaded = data
|
|
|
|
if start is not None and end is not None:
|
|
mask = (dates >= start) & (dates <= end)
|
|
|
|
ax.plot(
|
|
x_loaded[mask],
|
|
y_loaded[mask],
|
|
z_loaded[mask],
|
|
color=colors[i],
|
|
label=labels[i],
|
|
linewidth=1
|
|
)
|
|
|
|
if target_time is not None:
|
|
idx = min(range(len(dates)), key = lambda i: abs(dates[i] - target_time))
|
|
|
|
plot_cylinder(
|
|
ax,
|
|
x_loaded[idx],
|
|
y_loaded[idx],
|
|
z_loaded[idx],
|
|
radius=0.8,
|
|
height=0.3,
|
|
color=colors[i]
|
|
)
|
|
ax.set_xlabel('X')
|
|
ax.set_ylabel('Y')
|
|
ax.set_zlabel('Z')
|
|
ax.legend(ncols = 2, fontsize = 18)
|
|
ax.set_title(f'Cluster Spacecraft at {target_time} - GSE', fontsize = 26)
|
|
|
|
ax.set_box_aspect([1,1,1])
|
|
|
|
if start is not None:
|
|
ax.legend()
|
|
ax.auto_scale_xyz([-20,20], [-20,20], [-20,20])
|
|
|
|
filename = '/home/asterix/vieweg/Desktop/All_SC_ORBIT_20020906.pdf'
|
|
plt.savefig(filename)
|
|
plt.show()
|
|
|
|
#enter file name for all SC
|
|
data1 = load_data('C1_merged_2002.npy')
|
|
data2 = load_data('C2_merged_2002.npy')
|
|
data3 = load_data('C3_merged_2002.npy')
|
|
data4 = load_data('C4_merged_2002.npy')
|
|
|
|
datasets = [data1, data2, data3, data4]
|
|
|
|
|
|
|
|
#Enter time interval or specific point in time
|
|
plot_3d_multi(
|
|
datasets,
|
|
start=None,
|
|
end=None,
|
|
target_time=datetime(2002, 9, 6, 16, 0)
|
|
) |