293 lines
8.5 KiB
Python
293 lines
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)
|
||
|
|
)
|