From 41408c36be11efd9f6886818bde60b654957bba0 Mon Sep 17 00:00:00 2001 From: Leonie Date: Tue, 8 Jul 2025 16:18:34 +0300 Subject: [PATCH 1/5] Tasks 1-6 --- B_mean.py | 129 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 129 insertions(+) create mode 100644 B_mean.py diff --git a/B_mean.py b/B_mean.py new file mode 100644 index 0000000..394d5cf --- /dev/null +++ b/B_mean.py @@ -0,0 +1,129 @@ +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from scipy.linalg import norm + +data = np.load("C:/Users/user/Downloads/vtk_field/npy_field/Bnlfffe_NORH_NLFFFE_170904_055842.npy") + +nx, ny, nz = data.shape[:3] +x = np.arange(nx) +y = np.arange(ny) +z = np.arange(nz) +X, Y, Z = np.meshgrid(x, y, z, indexing='ij') +coordinates = np.stack((X, Y, Z), axis=-1) +squared_coords = coordinates ** 2 + + +def calculate_B_mean(data, coordinates, squared_coords, O, R): + term1 = squared_coords[..., 0] - 2 * coordinates[..., 0] * O[0] + O[0] ** 2 + term2 = squared_coords[..., 1] - 2 * coordinates[..., 1] * O[1] + O[1] ** 2 + term3 = squared_coords[..., 2] - 2 * coordinates[..., 2] * O[2] + O[2] ** 2 + squared_dist = term1 + term2 + term3 + mask = squared_dist <= R ** 2 + B_mean = np.mean(data[mask], axis=0) + indices = np.where(mask.ravel())[0] + return B_mean, indices + + +def create_local_frame(B_mean_): + B_n = B_mean_ / norm(B_mean_) + arbitrary_vec = np.array([1, 0, 0]) if abs(B_n[0]) < 0.9 else np.array([0, 1, 0]) + x_prime = np.cross(arbitrary_vec, B_n) + x_prime /= norm(x_prime) + y_prime = np.cross(B_n, x_prime) + return np.vstack([x_prime, y_prime, B_n]), B_n + + +def transform_field(B_global_, rotation_matrix_): + return np.dot(rotation_matrix_, B_global_) + + +def plot_all_2d_components(B_local_, O, R, title_prefix=""): + fig, axes = plt.subplots(1, 3, figsize=(18, 6)) + z_slice = O[2] + + x_min = max(0, int(O[0] - 2 * R)) + x_max = min(B_local_.shape[0], int(O[0] + 2 * R)) + y_min = max(0, int(O[1] - 2 * R)) + y_max = min(B_local_.shape[1], int(O[1] + 2 * R)) + + region = B_local_[x_min:x_max, y_min:y_max, z_slice] + + x_coords = np.arange(x_min, x_max) - O[0] + y_coords = np.arange(y_min, y_max) - O[1] + X, Y = np.meshgrid(y_coords, x_coords) + + vmax = max(np.max(np.abs(region[..., :3])), 1e-6) + + im1 = axes[0].imshow(region[..., 0], cmap='bwr', origin='lower', + vmin=-vmax, vmax=vmax, + extent=[y_coords[0], y_coords[-1], x_coords[-1], x_coords[0]]) + axes[0].set_title(fr"{title_prefix}$B_r$") + plt.colorbar(im1, ax=axes[0]) + + im2 = axes[1].imshow(region[..., 1], cmap='bwr', origin='lower', + vmin=-vmax, vmax=vmax, + extent=[y_coords[0], y_coords[-1], x_coords[-1], x_coords[0]]) + axes[1].set_title(fr"{title_prefix}$B_\phi$") + plt.colorbar(im2, ax=axes[1]) + + im3 = axes[2].imshow(region[..., 2], cmap='bwr', origin='lower', + vmin=-vmax, vmax=vmax, + extent=[y_coords[0], y_coords[-1], x_coords[-1], x_coords[0]]) + axes[2].set_title(fr"{title_prefix}$B_n$") + plt.colorbar(im3, ax=axes[2]) + + for ax in axes: + circle = plt.Circle((0, 0), R, color='black', + fill=False, linestyle='--', linewidth=2) + ax.add_patch(circle) + ax.scatter(0, 0, c='red', s=100, marker='*') + ax.set_xlim(y_coords[0], y_coords[-1]) + ax.set_ylim(x_coords[-1], x_coords[0]) + ax.grid(True, linestyle=':', alpha=0.5) + + plt.tight_layout() + plt.show() + + +def calculate_curl(B): + gradients = np.stack(np.gradient(B, axis=(0, 1, 2)), axis=-1) + curl_x = gradients[:, :, :, 2, 1] - gradients[:, :, :, 1, 2] + curl_y = gradients[:, :, :, 0, 2] - gradients[:, :, :, 2, 0] + curl_z = gradients[:, :, :, 1, 0] - gradients[:, :, :, 0, 1] + return np.stack([curl_x, curl_y, curl_z], axis=-1) + + +O = (203, 144, 11) +R = 15 + +B_mean, indices = calculate_B_mean(data, coordinates, squared_coords, O, R) +print(f"Среднее магнитное поле в окрестности R={R}:") +print(f"Bx = {B_mean[0]:.2f}, By = {B_mean[1]:.2f}, Bz = {B_mean[2]:.2f}") + +points_in_sphere = coordinates.reshape(-1, 3)[indices] +fig = plt.figure(figsize=(10, 8)) +ax = fig.add_subplot(111, projection='3d') +ax.scatter(points_in_sphere[:, 0], points_in_sphere[:, 1], points_in_sphere[:, 2], + s=5, c='blue', alpha=0.3, label=f'Точки в R={R}') +ax.scatter(*O, s=100, c='red', marker='*', label='Центр O') +ax.quiver(*O, *B_mean, length=20, color='black', + arrow_length_ratio=0.2, linewidth=3, label='B_mean') +ax.set_xlabel('Ось X') +ax.set_ylabel('Ось Y') +ax.set_zlabel('Ось Z') +ax.set_title(f'Окрестность точки O={O} радиусом R={R}') +plt.legend() +plt.tight_layout() +plt.show() + +rotation_matrix, B_n = create_local_frame(B_mean) +print(f"B_mean = {B_mean} (норма = {np.linalg.norm(B_mean):.3f})") +print(f"B_n = {B_n} (норма = {np.linalg.norm(B_n):.3f})") + +B_local = np.apply_along_axis(transform_field, 3, data, rotation_matrix) +plot_all_2d_components(B_local, O, R, title_prefix="Локальные компоненты ") + +curl_global = calculate_curl(data) +curl_local = np.apply_along_axis(transform_field, 3, curl_global, rotation_matrix) +plot_all_2d_components(curl_local, O, R, title_prefix="Curl ") \ No newline at end of file From aa3ae8d3279da72330aa56aa87d41b5c2b8e8318 Mon Sep 17 00:00:00 2001 From: Leonie Date: Thu, 10 Jul 2025 14:06:59 +0300 Subject: [PATCH 2/5] Local coordinate system problem was fixed. Linear interpolation was added. --- B_mean.py | 66 ++++++++++++++++++++++++++++++++++--------------------- 1 file changed, 41 insertions(+), 25 deletions(-) diff --git a/B_mean.py b/B_mean.py index 394d5cf..3e21cf0 100644 --- a/B_mean.py +++ b/B_mean.py @@ -14,6 +14,7 @@ squared_coords = coordinates ** 2 + def calculate_B_mean(data, coordinates, squared_coords, O, R): term1 = squared_coords[..., 0] - 2 * coordinates[..., 0] * O[0] + O[0] ** 2 term2 = squared_coords[..., 1] - 2 * coordinates[..., 1] * O[1] + O[1] ** 2 @@ -35,41 +36,54 @@ def create_local_frame(B_mean_): def transform_field(B_global_, rotation_matrix_): - return np.dot(rotation_matrix_, B_global_) + return np.dot(B_global_, rotation_matrix_.T) -def plot_all_2d_components(B_local_, O, R, title_prefix=""): +def plot_all_2d_components(B_local_, O, R, rotation_matrix, title_prefix=""): fig, axes = plt.subplots(1, 3, figsize=(18, 6)) - z_slice = O[2] - x_min = max(0, int(O[0] - 2 * R)) - x_max = min(B_local_.shape[0], int(O[0] + 2 * R)) - y_min = max(0, int(O[1] - 2 * R)) - y_max = min(B_local_.shape[1], int(O[1] + 2 * R)) + x_prime = rotation_matrix[0] + y_prime = rotation_matrix[1] + normal = rotation_matrix[2] + + grid_size = int(2 * R) + u = np.linspace(-R, R, grid_size) + v = np.linspace(-R, R, grid_size) + U, V = np.meshgrid(u, v) + + X_plane = O[0] + U * x_prime[0] + V * y_prime[0] + Y_plane = O[1] + U * x_prime[1] + V * y_prime[1] + Z_plane = O[2] + U * x_prime[2] + V * y_prime[2] + + from scipy.interpolate import RegularGridInterpolator + points = (np.arange(B_local_.shape[0]), np.arange(B_local_.shape[1]), np.arange(B_local_.shape[2])) - region = B_local_[x_min:x_max, y_min:y_max, z_slice] + Br_interp = RegularGridInterpolator(points, B_local_[..., 0], bounds_error=False, fill_value=0) + Bphi_interp = RegularGridInterpolator(points, B_local_[..., 1], bounds_error=False, fill_value=0) + Bn_interp = RegularGridInterpolator(points, B_local_[..., 2], bounds_error=False, fill_value=0) - x_coords = np.arange(x_min, x_max) - O[0] - y_coords = np.arange(y_min, y_max) - O[1] - X, Y = np.meshgrid(y_coords, x_coords) + sample_points = np.stack((X_plane, Y_plane, Z_plane), axis=-1) + Br_plane = Br_interp(sample_points) + Bphi_plane = Bphi_interp(sample_points) + Bn_plane = Bn_interp(sample_points) - vmax = max(np.max(np.abs(region[..., :3])), 1e-6) + vmax = max(np.max(np.abs(Br_plane)), np.max(np.abs(Bphi_plane)), np.max(np.abs(Bn_plane)), 1e-6) - im1 = axes[0].imshow(region[..., 0], cmap='bwr', origin='lower', + im1 = axes[0].imshow(Br_plane, cmap='bwr', origin='lower', vmin=-vmax, vmax=vmax, - extent=[y_coords[0], y_coords[-1], x_coords[-1], x_coords[0]]) + extent=[u[0], u[-1], v[-1], v[0]]) axes[0].set_title(fr"{title_prefix}$B_r$") plt.colorbar(im1, ax=axes[0]) - im2 = axes[1].imshow(region[..., 1], cmap='bwr', origin='lower', + im2 = axes[1].imshow(Bphi_plane, cmap='bwr', origin='lower', vmin=-vmax, vmax=vmax, - extent=[y_coords[0], y_coords[-1], x_coords[-1], x_coords[0]]) + extent=[u[0], u[-1], v[-1], v[0]]) axes[1].set_title(fr"{title_prefix}$B_\phi$") plt.colorbar(im2, ax=axes[1]) - im3 = axes[2].imshow(region[..., 2], cmap='bwr', origin='lower', + im3 = axes[2].imshow(Bn_plane, cmap='bwr', origin='lower', vmin=-vmax, vmax=vmax, - extent=[y_coords[0], y_coords[-1], x_coords[-1], x_coords[0]]) + extent=[u[0], u[-1], v[-1], v[0]]) axes[2].set_title(fr"{title_prefix}$B_n$") plt.colorbar(im3, ax=axes[2]) @@ -77,10 +91,12 @@ def plot_all_2d_components(B_local_, O, R, title_prefix=""): circle = plt.Circle((0, 0), R, color='black', fill=False, linestyle='--', linewidth=2) ax.add_patch(circle) - ax.scatter(0, 0, c='red', s=100, marker='*') - ax.set_xlim(y_coords[0], y_coords[-1]) - ax.set_ylim(x_coords[-1], x_coords[0]) + ax.scatter(0, 0, c='black', s=100, marker='*') + ax.set_xlim(u[0], u[-1]) + ax.set_ylim(v[-1], v[0]) ax.grid(True, linestyle=':', alpha=0.5) + ax.set_xlabel('u') + ax.set_ylabel('v') plt.tight_layout() plt.show() @@ -94,8 +110,8 @@ def calculate_curl(B): return np.stack([curl_x, curl_y, curl_z], axis=-1) -O = (203, 144, 11) -R = 15 +O = (200, 200, 35) +R = 10 B_mean, indices = calculate_B_mean(data, coordinates, squared_coords, O, R) print(f"Среднее магнитное поле в окрестности R={R}:") @@ -122,8 +138,8 @@ def calculate_curl(B): print(f"B_n = {B_n} (норма = {np.linalg.norm(B_n):.3f})") B_local = np.apply_along_axis(transform_field, 3, data, rotation_matrix) -plot_all_2d_components(B_local, O, R, title_prefix="Локальные компоненты ") +plot_all_2d_components(B_local, O, R, rotation_matrix, title_prefix="Локальные компоненты ") curl_global = calculate_curl(data) curl_local = np.apply_along_axis(transform_field, 3, curl_global, rotation_matrix) -plot_all_2d_components(curl_local, O, R, title_prefix="Curl ") \ No newline at end of file +plot_all_2d_components(curl_local, O, R, rotation_matrix, title_prefix="Curl ") \ No newline at end of file From c0bd94b9125b71dedced4c88bcde30c1424e6e40 Mon Sep 17 00:00:00 2001 From: Leonie Date: Thu, 10 Jul 2025 14:08:35 +0300 Subject: [PATCH 3/5] Local coordinate system problem was fixed. Linear interpolation was added. --- B_mean.py | 1 - 1 file changed, 1 deletion(-) diff --git a/B_mean.py b/B_mean.py index 3e21cf0..1a4a79b 100644 --- a/B_mean.py +++ b/B_mean.py @@ -14,7 +14,6 @@ squared_coords = coordinates ** 2 - def calculate_B_mean(data, coordinates, squared_coords, O, R): term1 = squared_coords[..., 0] - 2 * coordinates[..., 0] * O[0] + O[0] ** 2 term2 = squared_coords[..., 1] - 2 * coordinates[..., 1] * O[1] + O[1] ** 2 From b2d374b6e75972bf9bc67f356a1645b7080cb26f Mon Sep 17 00:00:00 2001 From: Leonie Date: Fri, 11 Jul 2025 02:47:24 +0300 Subject: [PATCH 4/5] The function create_local_frame was altered, stricter conditions were enforced. The problem with cylindrical coordinates was fixed. --- B_mean.py | 45 ++++++++++++++++++++++++++++----------------- 1 file changed, 28 insertions(+), 17 deletions(-) diff --git a/B_mean.py b/B_mean.py index 1a4a79b..250c9bd 100644 --- a/B_mean.py +++ b/B_mean.py @@ -27,10 +27,14 @@ def calculate_B_mean(data, coordinates, squared_coords, O, R): def create_local_frame(B_mean_): B_n = B_mean_ / norm(B_mean_) - arbitrary_vec = np.array([1, 0, 0]) if abs(B_n[0]) < 0.9 else np.array([0, 1, 0]) + basis_vectors = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) + dots = np.abs(basis_vectors @ B_n) + arbitrary_vec = basis_vectors[np.argmin(dots)] x_prime = np.cross(arbitrary_vec, B_n) x_prime /= norm(x_prime) y_prime = np.cross(B_n, x_prime) + y_prime /= norm(y_prime) + return np.vstack([x_prime, y_prime, B_n]), B_n @@ -57,14 +61,21 @@ def plot_all_2d_components(B_local_, O, R, rotation_matrix, title_prefix=""): from scipy.interpolate import RegularGridInterpolator points = (np.arange(B_local_.shape[0]), np.arange(B_local_.shape[1]), np.arange(B_local_.shape[2])) - Br_interp = RegularGridInterpolator(points, B_local_[..., 0], bounds_error=False, fill_value=0) - Bphi_interp = RegularGridInterpolator(points, B_local_[..., 1], bounds_error=False, fill_value=0) - Bn_interp = RegularGridInterpolator(points, B_local_[..., 2], bounds_error=False, fill_value=0) + Bx_interp = RegularGridInterpolator(points, B_local_[..., 0], bounds_error=False, fill_value=0) + By_interp = RegularGridInterpolator(points, B_local_[..., 1], bounds_error=False, fill_value=0) + Bz_interp = RegularGridInterpolator(points, B_local_[..., 2], bounds_error=False, fill_value=0) sample_points = np.stack((X_plane, Y_plane, Z_plane), axis=-1) - Br_plane = Br_interp(sample_points) - Bphi_plane = Bphi_interp(sample_points) - Bn_plane = Bn_interp(sample_points) + Bx_plane = Bx_interp(sample_points) + By_plane = By_interp(sample_points) + Bz_plane = Bz_interp(sample_points) + + r = np.sqrt(U ** 2 + V ** 2) + phi = np.arctan2(V, U) + + Br_plane = Bx_plane * np.cos(phi) + By_plane * np.sin(phi) + Bphi_plane = -Bx_plane * np.sin(phi) + By_plane * np.cos(phi) + Bn_plane = Bz_plane vmax = max(np.max(np.abs(Br_plane)), np.max(np.abs(Bphi_plane)), np.max(np.abs(Bn_plane)), 1e-6) @@ -110,10 +121,10 @@ def calculate_curl(B): O = (200, 200, 35) -R = 10 +R = 8 B_mean, indices = calculate_B_mean(data, coordinates, squared_coords, O, R) -print(f"Среднее магнитное поле в окрестности R={R}:") +print(f"Mean magnetic field in a spherical region of radius R={R}:") print(f"Bx = {B_mean[0]:.2f}, By = {B_mean[1]:.2f}, Bz = {B_mean[2]:.2f}") points_in_sphere = coordinates.reshape(-1, 3)[indices] @@ -121,23 +132,23 @@ def calculate_curl(B): ax = fig.add_subplot(111, projection='3d') ax.scatter(points_in_sphere[:, 0], points_in_sphere[:, 1], points_in_sphere[:, 2], s=5, c='blue', alpha=0.3, label=f'Точки в R={R}') -ax.scatter(*O, s=100, c='red', marker='*', label='Центр O') +ax.scatter(*O, s=100, c='red', marker='*', label='Center O') ax.quiver(*O, *B_mean, length=20, color='black', arrow_length_ratio=0.2, linewidth=3, label='B_mean') -ax.set_xlabel('Ось X') -ax.set_ylabel('Ось Y') -ax.set_zlabel('Ось Z') -ax.set_title(f'Окрестность точки O={O} радиусом R={R}') +ax.set_xlabel('Axis X') +ax.set_ylabel('Axis Y') +ax.set_zlabel('Axis Z') +ax.set_title(f'Vicinity of point O={O} with radius R={R}') plt.legend() plt.tight_layout() plt.show() rotation_matrix, B_n = create_local_frame(B_mean) -print(f"B_mean = {B_mean} (норма = {np.linalg.norm(B_mean):.3f})") -print(f"B_n = {B_n} (норма = {np.linalg.norm(B_n):.3f})") +print(f"B_mean = {B_mean} (magnitude = {np.linalg.norm(B_mean):.3f})") +print(f"B_n = {B_n} (magnitude = {np.linalg.norm(B_n):.3f})") B_local = np.apply_along_axis(transform_field, 3, data, rotation_matrix) -plot_all_2d_components(B_local, O, R, rotation_matrix, title_prefix="Локальные компоненты ") +plot_all_2d_components(B_local, O, R, rotation_matrix, title_prefix="Local components: ") curl_global = calculate_curl(data) curl_local = np.apply_along_axis(transform_field, 3, curl_global, rotation_matrix) From 31be94f26dfab08fed0c1ca617c10d7671745d6e Mon Sep 17 00:00:00 2001 From: Leonie Date: Fri, 18 Jul 2025 06:15:00 +0300 Subject: [PATCH 5/5] Linear extrapolation was added to avoid invalid results for points close to borders of the cube. All warnings were fixed. --- B_mean.py | 194 ++++++++++++++++------------------------------- requirements.txt | 1 + 2 files changed, 68 insertions(+), 127 deletions(-) diff --git a/B_mean.py b/B_mean.py index 250c9bd..df50f31 100644 --- a/B_mean.py +++ b/B_mean.py @@ -1,155 +1,95 @@ import numpy as np import matplotlib.pyplot as plt -from mpl_toolkits.mplot3d import Axes3D from scipy.linalg import norm +from scipy.ndimage import map_coordinates +from typing import Tuple + + +def linear_extrapolate(data_, points): + return map_coordinates(data_, points.T, order=1, mode='nearest') + data = np.load("C:/Users/user/Downloads/vtk_field/npy_field/Bnlfffe_NORH_NLFFFE_170904_055842.npy") nx, ny, nz = data.shape[:3] -x = np.arange(nx) -y = np.arange(ny) -z = np.arange(nz) +x, y, z = np.arange(nx), np.arange(ny), np.arange(nz) X, Y, Z = np.meshgrid(x, y, z, indexing='ij') coordinates = np.stack((X, Y, Z), axis=-1) -squared_coords = coordinates ** 2 - - -def calculate_B_mean(data, coordinates, squared_coords, O, R): - term1 = squared_coords[..., 0] - 2 * coordinates[..., 0] * O[0] + O[0] ** 2 - term2 = squared_coords[..., 1] - 2 * coordinates[..., 1] * O[1] + O[1] ** 2 - term3 = squared_coords[..., 2] - 2 * coordinates[..., 2] * O[2] + O[2] ** 2 - squared_dist = term1 + term2 + term3 - mask = squared_dist <= R ** 2 - B_mean = np.mean(data[mask], axis=0) - indices = np.where(mask.ravel())[0] - return B_mean, indices - - -def create_local_frame(B_mean_): - B_n = B_mean_ / norm(B_mean_) - basis_vectors = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) - dots = np.abs(basis_vectors @ B_n) - arbitrary_vec = basis_vectors[np.argmin(dots)] - x_prime = np.cross(arbitrary_vec, B_n) - x_prime /= norm(x_prime) - y_prime = np.cross(B_n, x_prime) - y_prime /= norm(y_prime) - - return np.vstack([x_prime, y_prime, B_n]), B_n - - -def transform_field(B_global_, rotation_matrix_): - return np.dot(B_global_, rotation_matrix_.T) - - -def plot_all_2d_components(B_local_, O, R, rotation_matrix, title_prefix=""): - fig, axes = plt.subplots(1, 3, figsize=(18, 6)) - - x_prime = rotation_matrix[0] - y_prime = rotation_matrix[1] - normal = rotation_matrix[2] - grid_size = int(2 * R) - u = np.linspace(-R, R, grid_size) - v = np.linspace(-R, R, grid_size) - U, V = np.meshgrid(u, v) - X_plane = O[0] + U * x_prime[0] + V * y_prime[0] - Y_plane = O[1] + U * x_prime[1] + V * y_prime[1] - Z_plane = O[2] + U * x_prime[2] + V * y_prime[2] +def calculate_b_mean(data_, coordinates_, center_, radius_): + squared_dist = np.sum((coordinates_ - np.array(center_)) ** 2, axis=-1) + mask = np.array(squared_dist <= radius_**2) + return np.mean(data_[mask], axis=0), np.where(mask.ravel())[0] - from scipy.interpolate import RegularGridInterpolator - points = (np.arange(B_local_.shape[0]), np.arange(B_local_.shape[1]), np.arange(B_local_.shape[2])) - Bx_interp = RegularGridInterpolator(points, B_local_[..., 0], bounds_error=False, fill_value=0) - By_interp = RegularGridInterpolator(points, B_local_[..., 1], bounds_error=False, fill_value=0) - Bz_interp = RegularGridInterpolator(points, B_local_[..., 2], bounds_error=False, fill_value=0) +# noinspection PyUnreachableCode +def create_local_frame(b_mean_: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + b_n_ = b_mean_ / norm(b_mean_) + basis = np.eye(3) + x_prime = np.cross(basis[np.argmin(np.abs(basis @ b_n_))], b_n_) + x_prime /= norm(x_prime) + y_prime = np.cross(b_n_, x_prime) + y_prime /= norm(y_prime) + return np.vstack([x_prime, y_prime, b_n_]), b_n_ - sample_points = np.stack((X_plane, Y_plane, Z_plane), axis=-1) - Bx_plane = Bx_interp(sample_points) - By_plane = By_interp(sample_points) - Bz_plane = Bz_interp(sample_points) - r = np.sqrt(U ** 2 + V ** 2) - phi = np.arctan2(V, U) +def transform_field(b_, rotation_matrix_): + return np.dot(b_, rotation_matrix_.T) - Br_plane = Bx_plane * np.cos(phi) + By_plane * np.sin(phi) - Bphi_plane = -Bx_plane * np.sin(phi) + By_plane * np.cos(phi) - Bn_plane = Bz_plane - vmax = max(np.max(np.abs(Br_plane)), np.max(np.abs(Bphi_plane)), np.max(np.abs(Bn_plane)), 1e-6) +def plot_all_2d_components(b_local_, center_, radius_, rotation_matrix_, title_prefix=""): + fig, axes = plt.subplots(1, 3, figsize=(18, 6)) + x_prime, y_prime = rotation_matrix_[0], rotation_matrix_[1] + + u = v = np.linspace(-radius_, radius_, int(2 * radius_)) + u_grid, v_grid = np.meshgrid(u, v) + x_plane = center_[0] + u_grid * x_prime[0] + v_grid * y_prime[0] + y_plane = center_[1] + u_grid * x_prime[1] + v_grid * y_prime[1] + z_plane = center_[2] + u_grid * x_prime[2] + v_grid * y_prime[2] + sample_points = np.stack((x_plane, y_plane, z_plane), axis=-1) + + b_x = linear_extrapolate(b_local_[..., 0], sample_points) + b_y = linear_extrapolate(b_local_[..., 1], sample_points) + b_z = linear_extrapolate(b_local_[..., 2], sample_points) + + phi = np.arctan2(v_grid, u_grid) + b_r = b_x * np.cos(phi) + b_y * np.sin(phi) + b_phi = -b_x * np.sin(phi) + b_y * np.cos(phi) + + vmax = max(np.max(np.abs(b_r)), np.max(np.abs(b_phi)), np.max(np.abs(b_z)), 1e-6) + for ax, component in zip(axes, [b_r, b_phi, b_z]): + im = ax.imshow(component, cmap='bwr', vmin=-vmax, vmax=vmax, + extent=[u[0], u[-1], v[-1], v[0]], origin='lower') + plt.colorbar(im, ax=ax) + ax.add_patch(plt.Circle((0, 0), radius, color='k', fill=False, linestyle='--')) + ax.scatter(0, 0, c='k', s=100, marker='*') + ax.grid(True, linestyle=':', alpha=0.5) - im1 = axes[0].imshow(Br_plane, cmap='bwr', origin='lower', - vmin=-vmax, vmax=vmax, - extent=[u[0], u[-1], v[-1], v[0]]) axes[0].set_title(fr"{title_prefix}$B_r$") - plt.colorbar(im1, ax=axes[0]) - - im2 = axes[1].imshow(Bphi_plane, cmap='bwr', origin='lower', - vmin=-vmax, vmax=vmax, - extent=[u[0], u[-1], v[-1], v[0]]) axes[1].set_title(fr"{title_prefix}$B_\phi$") - plt.colorbar(im2, ax=axes[1]) - - im3 = axes[2].imshow(Bn_plane, cmap='bwr', origin='lower', - vmin=-vmax, vmax=vmax, - extent=[u[0], u[-1], v[-1], v[0]]) axes[2].set_title(fr"{title_prefix}$B_n$") - plt.colorbar(im3, ax=axes[2]) - - for ax in axes: - circle = plt.Circle((0, 0), R, color='black', - fill=False, linestyle='--', linewidth=2) - ax.add_patch(circle) - ax.scatter(0, 0, c='black', s=100, marker='*') - ax.set_xlim(u[0], u[-1]) - ax.set_ylim(v[-1], v[0]) - ax.grid(True, linestyle=':', alpha=0.5) - ax.set_xlabel('u') - ax.set_ylabel('v') - plt.tight_layout() plt.show() -def calculate_curl(B): - gradients = np.stack(np.gradient(B, axis=(0, 1, 2)), axis=-1) - curl_x = gradients[:, :, :, 2, 1] - gradients[:, :, :, 1, 2] - curl_y = gradients[:, :, :, 0, 2] - gradients[:, :, :, 2, 0] - curl_z = gradients[:, :, :, 1, 0] - gradients[:, :, :, 0, 1] - return np.stack([curl_x, curl_y, curl_z], axis=-1) - - -O = (200, 200, 35) -R = 8 - -B_mean, indices = calculate_B_mean(data, coordinates, squared_coords, O, R) -print(f"Mean magnetic field in a spherical region of radius R={R}:") -print(f"Bx = {B_mean[0]:.2f}, By = {B_mean[1]:.2f}, Bz = {B_mean[2]:.2f}") - -points_in_sphere = coordinates.reshape(-1, 3)[indices] -fig = plt.figure(figsize=(10, 8)) -ax = fig.add_subplot(111, projection='3d') -ax.scatter(points_in_sphere[:, 0], points_in_sphere[:, 1], points_in_sphere[:, 2], - s=5, c='blue', alpha=0.3, label=f'Точки в R={R}') -ax.scatter(*O, s=100, c='red', marker='*', label='Center O') -ax.quiver(*O, *B_mean, length=20, color='black', - arrow_length_ratio=0.2, linewidth=3, label='B_mean') -ax.set_xlabel('Axis X') -ax.set_ylabel('Axis Y') -ax.set_zlabel('Axis Z') -ax.set_title(f'Vicinity of point O={O} with radius R={R}') -plt.legend() -plt.tight_layout() -plt.show() - -rotation_matrix, B_n = create_local_frame(B_mean) -print(f"B_mean = {B_mean} (magnitude = {np.linalg.norm(B_mean):.3f})") -print(f"B_n = {B_n} (magnitude = {np.linalg.norm(B_n):.3f})") - -B_local = np.apply_along_axis(transform_field, 3, data, rotation_matrix) -plot_all_2d_components(B_local, O, R, rotation_matrix, title_prefix="Local components: ") +def calculate_curl(b_): + gradients = np.stack(np.gradient(b_, axis=(0, 1, 2)), axis=-1) + curl = np.stack([gradients[..., 2, 1] - gradients[..., 1, 2], + gradients[..., 0, 2] - gradients[..., 2, 0], + gradients[..., 1, 0] - gradients[..., 0, 1]], axis=-1) + return curl + + +center = (203, 10, 1) +radius = 15 + +b_mean, indices = calculate_b_mean(data, coordinates, center, radius) +rotation_matrix, B_n = create_local_frame(b_mean) +b_local = np.apply_along_axis(transform_field, 3, data, rotation_matrix) + +plot_all_2d_components(b_local, center, radius, rotation_matrix, "Local components: ") curl_global = calculate_curl(data) curl_local = np.apply_along_axis(transform_field, 3, curl_global, rotation_matrix) -plot_all_2d_components(curl_local, O, R, rotation_matrix, title_prefix="Curl ") \ No newline at end of file +plot_all_2d_components(curl_local, center, radius, rotation_matrix, "Curl ") diff --git a/requirements.txt b/requirements.txt index bb28ce7..26dd036 100644 --- a/requirements.txt +++ b/requirements.txt @@ -16,6 +16,7 @@ idna==3.10 isodate==0.7.2 joblib==1.5.0 lxml==5.4.0 +matplotlib==3.10.3 multidict==6.4.3 numpy==2.2.5 packaging==25.0