"""
analysis.py
This module handles the analysis derived quantities like USFLUX, DM, Polar Field of SFT (2D) outputs.
Functions:
- calculate_usflux: Computes the total unsigned magnetic flux on the magnetic field map generated by the SFT model.
- calculate_dm: Calculates the axial dipole moment.
- calculate_polar_field: Computes the polar field.
"""
import numpy as np
[docs]
def calculate_usflx(all_br_data,grid,time_duration):
"""
Calculates the total unsigned magnetic flux of the surface magnetic fields for a given time duration.
Parameters:
all_br_data (np.ndarray): 3D array containing the surface magnetic field data [time,latitude,longitude].
grid (dict): Dictionary containing grid information ('theta', 'phi', and their spacings).
time_duration (int): Number of time steps in the simulation. Set this to a range of days from 0 to the number of days as a touple, e.g. [0, 365].
Returns:
float or np.ndarray: Total unsigned flux.
"""
# Read the grid information
colatitude = grid['colatitude']
longitude = grid['longitude']
delta_theta = grid['dtheta']
delta_phi = grid['dphi']
solar_radius = 6.955 * 10**8 # solar radius in meters
temp_bsinth = np.zeros_like(all_br_data[0,:,:])
if time_duration[0] == time_duration[1]:
for i1 in range(longitude.shape[0]):
temp_bsinth[:,i1] = np.abs(all_br_data[int(time_duration[0]),:,i1]) * np.sin(colatitude)
usflx = np.sum(temp_bsinth) * delta_theta * delta_phi * (solar_radius*1e2)**2
else:
time_start = int(time_duration[0])
time_end = int(time_duration[1])
time_total = time_end-time_start+1
usflx = np.zeros(time_end-time_start+1)
for i in range(time_total):
for i1 in range(longitude.shape[0]):
temp_bsinth[:,i1] = np.abs(all_br_data[i,:,i1]) * np.sin(colatitude)
usf_1d = np.sum(temp_bsinth) * delta_theta * delta_phi * (solar_radius*1e2)**2
usflx[i] = usf_1d
return usflx
[docs]
def calculate_dm(all_br_data,grid,time_duration):
"""
Calculates the axial dipole moment for a given time duration.
Parameters:
all_br_data (np.ndarray): 3D array containing the surface magnetic field data [time,latitude,longitude].
grid (dict): Dictionary containing grid information ('theta', 'phi', and their spacings).
time_duration (int): Number of time steps in the simulation. Set this to a range of days from 0 to the number of days as a touple, e.g. [0, 365].
Returns:
float or np.ndarray: Axial dipole moment.
"""
# Read the grid information
colatitude = grid['colatitude']
longitude = grid['longitude']
if time_duration[0] == time_duration[1]:
dipole_moment = compute_axial_dipole_moment(all_br_data[int(time_duration[0]),:,:], colatitude, longitude, 1.0)
else:
time_start = int(time_duration[0])
time_end = int(time_duration[1])
time_total = time_end-time_start+1
dipole_moment = np.zeros(time_end-time_start+1)
for i in range(time_total):
dipole_moment[i] = compute_axial_dipole_moment(all_br_data[i,:,:], colatitude, longitude, 1.0)
return dipole_moment
def compute_axial_dipole_moment(B_field, colatitude, longitude, solar_radius):
"""
Computes the magnetic axial dipole moment from the magnetic field array.
Parameters:
-----------
B_field : 2D numpy array
Radial component of the magnetic field on the solar surface, with shape (num_theta+1, num_phi+1).
colatitude : 1D numpy array
Array of colatitudes (theta), ranging from 0 to pi.
longitude : 1D numpy array
Array of longitudes (phi), ranging from 0 to 2*pi.
solar_radius : float
Radius of the Sun in solar radius units e.g., set 1 for the surface.
Returns:
--------
dipole_moment : float
The computed axial dipole moment.
"""
# Create a meshgrid for the angles
delta_theta = np.abs(colatitude[1] - colatitude[0])
delta_phi = np.abs(longitude[1] - longitude[0])
# Calculate the weight for the dipole moment (cos(theta) * sin(theta))
weight = np.cos(colatitude[:, np.newaxis]) * np.sin(colatitude[:, np.newaxis])
# Integrate the magnetic field weighted by the dipole term over the solar surface
dipole_moment = solar_radius**2 * np.sum(B_field * weight) * delta_theta * delta_phi
return dipole_moment
[docs]
def calculate_polar_field(all_br_data,grid,time_duration):
"""
Calculates the polar magnetic flux for a given time duration.
Parameters:
all_br_data (np.ndarray): 3D array containing the surface magnetic field data [time,latitude,longitude].
grid (dict): Dictionary containing grid information ('theta', 'phi', and their spacings).
time_duration (int): Number of time steps in the simulation. Set this to a range of days from 0 to the number of days as a touple, e.g. [0, 365].
Returns:
float or np.ndarray: Polar magnetic field for north and south pole.
"""
# Read the grid information
colatitude = grid['colatitude']
deg_pol = 70
if time_duration[0] == time_duration[1]:
polar_field_south = np.mean(all_br_data[int(time_duration[0]),np.argwhere(colatitude >= np.deg2rad(180-deg_pol))[0][0]:,:])
polar_field_north = np.mean(all_br_data[int(time_duration[0]),0:np.argwhere(colatitude >= np.deg2rad(deg_pol))[0][0],:])
else:
time_start = int(time_duration[0])
time_end = int(time_duration[1])
time_total = time_end-time_start+1
polar_field_north = np.zeros(time_total)
polar_field_south = np.zeros(time_total)
for i in range(time_total):
polar_field_south[i] = np.mean(all_br_data[i,np.argwhere(colatitude >= np.deg2rad(180-deg_pol))[0][0]:,:])
polar_field_north[i] = np.mean(all_br_data[i,0:np.argwhere(colatitude >= np.deg2rad(deg_pol))[0][0],:])
return polar_field_north, polar_field_south
[docs]
def calculate_polar_flux(all_br_data,grid,time_duration, pol_cap_extent_deg = 20, R_sun=6.98e10):
"""
Compute the total polar flux in Maxwell (Mx) near the polar caps (20-degree extent).
Parameters:
all_br_data (np.ndarray): 3D array containing the surface magnetic field data [time,latitude,longitude].
grid (dict): Dictionary containing grid information ('theta', 'phi', and their spacings).
time_duration (int): Number of time steps in the simulation. Set this to a range of days from 0 to the number of days as a touple, e.g. [0, 365].
R_sun (float): Solar radius in cm.
Returns:
float or np.ndarray: North and South polar flux in Maxwell.
"""
# Define theta and phi grid
theta = grid['colatitude']
dtheta = grid['dtheta']
dphi = grid['dphi']
# Compute area element dA = R_sun^2 * sin(theta) * dtheta * dphi
dA = (R_sun**2) * np.sin(theta)[:, None] * dtheta * dphi
# Polar cap extent
polar_cap_extent = np.deg2rad(pol_cap_extent_deg)
# Define polar cap indices
north_cap = (theta <= polar_cap_extent)
south_cap = (theta >= np.pi - polar_cap_extent)
# Compute flux in both caps
if time_duration[0] == time_duration[1]:
flux_north = np.sum(all_br_data[int(time_duration[0]),north_cap, :] * dA[north_cap, :]) # North polar flux
flux_south = np.sum(all_br_data[int(time_duration[0]),south_cap, :] * dA[south_cap, :]) # South polar flux
else:
time_start = int(time_duration[0])
time_end = int(time_duration[1])
time_total = time_end-time_start+1
flux_north = np.zeros(time_total)
flux_south = np.zeros(time_total)
for i in range(time_total):
flux_north[i] = np.sum(all_br_data[i,north_cap, :] * dA[north_cap, :])
flux_south[i] = np.sum(all_br_data[i,south_cap, :] * dA[south_cap, :])
return flux_north, flux_south