#!/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) )