From 9433a3f7efff19d6cffa03ee5fadf3748d7f35e3 Mon Sep 17 00:00:00 2001 From: jameswilburlewis Date: Wed, 10 Jan 2024 17:58:46 -0800 Subject: [PATCH] First cut at a FOTE (first order taylor expansion) magnetic null finder for MMS, Cluster, or other missions that can operate in tetrahedral formations. --- pyspedas/analysis/find_magnetic_nulls.py | 296 ++++++++++++++++++ .../analysis/tests/test_magnetic_nulls.py | 20 ++ 2 files changed, 316 insertions(+) create mode 100644 pyspedas/analysis/find_magnetic_nulls.py create mode 100644 pyspedas/analysis/tests/test_magnetic_nulls.py diff --git a/pyspedas/analysis/find_magnetic_nulls.py b/pyspedas/analysis/find_magnetic_nulls.py new file mode 100644 index 00000000..75ad424a --- /dev/null +++ b/pyspedas/analysis/find_magnetic_nulls.py @@ -0,0 +1,296 @@ +import numpy as np +from .lingradest import lingradest +from pyspedas import tinterpol +from pytplot import get_data, store_data, time_double, time_string, time_clip, deflag, tnames, options, tplot_options + +def find_magnetic_nulls_fote(positions=None, fields=None, scale_factor=1.0): + """ + Find magnetic null points, given a set of four-point magnetic field observations (e.g. MMS or CLUSTER) using the + First Order Taylor Expansion (FOTE) method. + + This method uses tetrahedral interpolation to estimate the magnetic field (B0) and field gradients at the tetrahedron + barycenter. From the field gradients we can construct a Jacobian matrix (J). The field in that region can be + expressed, to a first order approximation, as + + B = B0 + JV + + where V is the position vector relative to the barycenter. At a magnetic null V_null, all components of B are zero, + so we have + + [0] = B0 + JV_null + + so that + + V_null = J_inv[-B0] + + As long as the four field measurements all differ, a null point will always be found. For it to be credible, + it should be in the neighborhood where the linear gradient approximation is expected to be valid. + + The topology of the field around the null can be inferred from the eigenvalues and eigenvectors of the + estimated Jacobian. + + Parameters + ----------- + + positions: list of str + A 4-element list of tplot variable names representing the probe positions + + fields: list of str + A 4-element list of tplot variable names representing the magnetic field measurements + + scale_factor: float + The scale factor passed lingradest routine to scale some of the distances + Default: 1.0 + + Returns + ------- + A list of tplot variables describing the nulls found + + + Example: + >>> import pyspedas + >>> from pytplot + >>> # load data and ephemeris + >>> pyspedas.find_magnetic_nulls_fote(positions,fields) + """ + + # Input data needs to be sanitized, by removing any 'nans', and finding a time range common + # to all the deflagged variables + + positions_df = [] + fields_df = [] + for i in range(4): + positions_df.append(positions[i]+'_df') + fields_df.append(fields[i]+'_df') + deflag(positions[i], method='remove_nan', new_tvar=positions_df[i]) + deflag(fields[i], method='remove_nan', new_tvar=fields_df[i]) + + print(tnames('*_df')) + d=get_data(fields_df[1]) + print(d.y[-1]) + # Interpolate all deflagged positions and field measurements to the first probe's field measurement times + + # We need to find a time range that's contained in all the position and field variables, to + # prevent interpolation from introducing NaNs and causing exceptions to be thrown in the calculations + # below. + + tmin = 0.0 + tmax = time_double('3000-01-01') + for i in range(4): + d = get_data(positions_df[i]) + dmin = d.times[0] + dmax = d.times[-1] + tmin = np.max([tmin, dmin]) + tmax = np.min([tmax, dmax]) + + d = get_data(fields_df[i]) + dmin = d.times[0] + dmax = d.times[-1] + tmin = np.max([tmin, dmin]) + tmax = np.min([tmax, dmax]) + + # Shave another millisecond off each end just to be sure (because floating point equality testing is weird) + tmin = tmin + .001 + tmax = tmax - .001 + + # We'll interpolate everything to the first set of times, so we'll time_clip that one to the desired time range + time_clip(fields_df[0], tmin, tmax, new_names='f0_tc') + + # Now we can interpolate safely + tinterpol(positions_df,'f0_tc',newname=['pos1_i','pos2_i','pos3_i','pos4_i']) + tinterpol(fields_df, 'f0_tc', newname=['b1_i','b2_i','b3_i','b4_i']) + d1 = get_data('pos1_i') + d2 = get_data('pos2_i') + d3 = get_data('pos3_i') + d4 = get_data('pos4_i') + + # Position variables for MMS include a total distance as the last element, so only take the first three + r1=d1.y[:,0:3] + r2=d2.y[:,0:3] + r3=d3.y[:,0:3] + r4=d4.y[:,0:3] + + b1 = get_data('b1_i') + b2 = get_data('b2_i') + b3 = get_data('b3_i') + b4 = get_data('b4_i') + + # The MMS field variables also include a fourth component with the total field + bx1 = b1.y[:,0] + by1 = b1.y[:,1] + bz1 = b1.y[:,2] + + bx2 = b3.y[:,0] + by2 = b3.y[:,1] + bz2 = b3.y[:,2] + + bx3 = b3.y[:,0] + by3 = b3.y[:,1] + bz3 = b3.y[:,2] + + bx4 = b4.y[:,0] + by4 = b4.y[:,1] + bz4 = b4.y[:,2] + + datapoint_count = len(bx1) + times = b1.times + + # Get the data arrays from the interpolated tplot variables + + # Call the lingradest routine to get the barycenter locations, field at barycenter, and field gradients at barycenter + + # Output of lingradest is a dictionary of string keys and numpy array values + + lingrad_output = lingradest(bx1, bx2, bx3, bx4, + by1, by2, by3, by4, + bz1, bz2, bz3, bz4, + r1,r2,r3,r4, scale_factor=scale_factor) + + # Output arrays + out_pos_null = np.zeros((datapoint_count,3)) + out_pos_bary = np.zeros((datapoint_count,3)) + out_null_bary_dist = np.zeros((datapoint_count)) + out_null_p1_dist = np.zeros((datapoint_count)) + out_null_p2_dist = np.zeros((datapoint_count)) + out_null_p3_dist = np.zeros((datapoint_count)) + out_null_p4_dist = np.zeros((datapoint_count)) + out_null_eta = np.zeros((datapoint_count)) + out_null_xi = np.zeros((datapoint_count)) + + for i in range(datapoint_count): + Rbary = lingrad_output['Rbary'][i] + out_pos_bary[i,:] = Rbary[:] + + # Get field vector at barycenter for i-th point + bxbc = lingrad_output['Bxbc'][i] + bybc = lingrad_output['Bybc'][i] + bzbc = lingrad_output['Bzbc'][i] + b0 = np.array((bxbc, bybc, bzbc)) # Field at barycenter + b0_neg = np.array((-bxbc,-bybc,-bzbc)) # Negative of field at barycenter + + # Form Jacobian matrix from field gradients at barycenter + LGBx = lingrad_output['LGBx'][i] + LGBy = lingrad_output['LGBy'][i] + LGBz = lingrad_output['LGBz'][i] + J = np.stack((LGBx, LGBy, LGBz)) # Jacobian matrix (or its transpose?) + J_inv = np.linalg.inv(J) + + # Solve for null position with respect to barycenter (maybe check for singular matrix first?) + + r_null = np.linalg.solve(J,b0_neg) + r_null_alt = np.matmul(J_inv,b0_neg)# Check that estimated field at null is near 0 + F_null = b0 + np.matmul(J,r_null) + F_null_alt = b0 + np.matmul(J,r_null_alt) + r_null_bary_dist = np.linalg.norm(r_null) + out_null_bary_dist[i] = r_null_bary_dist + + # Translate r_null to origin of probe coordinate system + pos_null = r_null + Rbary + out_pos_null[i,:] = pos_null[:] + + + # Get distances from null to each probe in the tetrahedron + p1_null = r1[i] - pos_null + p1_null_dist = np.linalg.norm(p1_null) + p2_null = r2[i] - pos_null + p2_null_dist = np.linalg.norm(p2_null) + p3_null = r3[i] + pos_null + p3_null_dist = np.linalg.norm(p3_null) + p4_null = r4[i] - pos_null + p4_null_dist = np.linalg.norm(p4_null) + + out_null_p1_dist[i] = p1_null_dist + out_null_p2_dist[i] = p2_null_dist + out_null_p3_dist[i] = p3_null_dist + out_null_p4_dist[i] = p4_null_dist + + # Estimate field at each probe using the estimated linear gradient + # Might need to be negated for correct sense of distances? + # These should all be quite close to the measured fields. + + dR1 = -1.0 * lingrad_output['dR1'][i] + dR2 = -1.0 * lingrad_output['dR2'][i] + dR3 = -1.0 * lingrad_output['dR3'][i] + dR4 = -1.0 * lingrad_output['dR4'][i] + + F1_est = b0 + np.matmul(J,dR1) + F2_est = b0 + np.matmul(J,dR2) + F3_est = b0 + np.matmul(J,dR3) + F4_est = b0 + np.matmul(J,dR4) + + F1_obs = np.array([bx1[i],by1[i],bz1[i]]) + F2_obs = np.array([bx2[i],by2[i],bz2[i]]) + F3_obs = np.array([bx3[i],by3[i],bz3[i]]) + F4_obs = np.array([bx4[i],by4[i],bz4[i]]) + + # Look for opposite signs in each field component, necessary if null is within tetrahedron + bxmax = np.max( (F1_obs[0],F2_obs[0],F3_obs[0],F4_obs[0])) + bxmin = np.min( (F1_obs[0],F2_obs[0],F3_obs[0],F4_obs[0])) + bymax = np.max( (F1_obs[1],F2_obs[1],F3_obs[1],F4_obs[1])) + bymin = np.min( (F1_obs[1],F2_obs[1],F3_obs[1],F4_obs[1])) + bzmax = np.max( (F1_obs[2],F2_obs[2],F3_obs[2],F4_obs[2])) + bzmin = np.min( (F1_obs[2],F2_obs[2],F3_obs[2],F4_obs[2])) + + tstring = time_string(times[i]) + if (bxmax*bxmin < 0.0) and (bymax*bymin < 0.0) and (bzmax*bzmin < 1.0): + print("Null may be in tetrahedron at time ",tstring ) + print("F1 obs:", F1_obs) + print("F2 obs:", F2_obs) + print("F3 obs:", F3_obs) + print("F4 obs:", F4_obs) + + + + + # Get eigenvalues of J to characterize field topology near the null + try: + eigenvalues, eigevectors = np.linalg.eig(J) + except: + print('np.linalg.eig failed') + print('i=',i) + print('J',J) + return + lambdas = eigenvalues + l1_norm = np.linalg.norm(lambdas[0]) + l2_norm = np.linalg.norm(lambdas[1]) + l3_norm = np.linalg.norm(lambdas[2]) + eig_sum_norm = np.linalg.norm(lambdas[0] + lambdas[1] + lambdas[2]) + eig_max_norm = np.max((l1_norm,l2_norm,l3_norm)) + + + # div B should be close to 0 (exactly 0 in theory), so the difference from 0 is a measure of how + # credible any null we've found might be. We normalize it, dividing by the magnitude of the curl, + # to get the statistic eta. eta < 0.40 for a credible null. + + + divB = lingrad_output['LD'][i] + LCx = lingrad_output['LCxB'][i] + LCy = lingrad_output['LCyB'][i] + LCz = lingrad_output['LCzB'][i] + curlB = np.array([LCx,LCy,LCz]) + curlB_norm = np.linalg.norm(curlB) + # eta - | del dot B| / |del x B| + eta = np.abs(divB)/curlB_norm + out_null_eta[i] = eta + + # Similarly, the sum of the eigenvectors of the Jacobian gives another figure of merit. Here, + # we normalize by dividing the norm of the sum by the norm of the largest eigenvalue to yield xi. + # xi < 0.40 if we trust the topology derived from the eigenvalues. Larger values may mean that the type of + # null detected (radial vs. spiral, etc) may be an artifact. + + xi = eig_sum_norm/eig_max_norm + out_null_xi[i] = xi + + store_data('magnull_pos',data={'x':times,'y':out_pos_null}) + store_data('magnull_pos_bary',data={'x':times,'y':out_pos_bary}) + store_data('magnull_null_bary_dist',data={'x':times,'y':out_null_bary_dist}) + options('magnull_null_bary_dist','yrange',[0.0, 5000.0]) + print('Minimum null distance:',np.min(out_null_bary_dist)) + store_data('magnull_null_p1_dist',data={'x':times,'y':out_null_p1_dist}) + store_data('magnull_null_p2_dist',data={'x':times,'y':out_null_p2_dist}) + store_data('magnull_null_p3_dist',data={'x':times,'y':out_null_p3_dist}) + store_data('magnull_null_p4_dist',data={'x':times,'y':out_null_p4_dist}) + store_data('magnull_eta', data={'x':times, 'y':out_null_eta}) + store_data('magnull_xi', data={'x':times, 'y':out_null_xi}) + + diff --git a/pyspedas/analysis/tests/test_magnetic_nulls.py b/pyspedas/analysis/tests/test_magnetic_nulls.py new file mode 100644 index 00000000..2c3edbbf --- /dev/null +++ b/pyspedas/analysis/tests/test_magnetic_nulls.py @@ -0,0 +1,20 @@ + +import unittest +import pyspedas +from pytplot import data_exists, tplot + +class MagNullTestCases(unittest.TestCase): + + + def test_find_magnetic_nulls_fote_mms(self): + from pytplot import tplot + data = pyspedas.mms.fgm(probe=[1, 2, 3, 4], trange=['2015-09-19/07:40', '2015-09-19/07:45'], data_rate='srvy', time_clip=True, varformat='*_gse_*', get_fgm_ephemeris=True) + fields = ['mms'+prb+'_fgm_b_gse_srvy_l2' for prb in ['1', '2', '3', '4']] + positions = ['mms'+prb+'_fgm_r_gse_srvy_l2' for prb in ['1', '2', '3', '4']] + tplot(fields) + curl = pyspedas.find_magnetic_nulls_fote(fields=fields, positions=positions) + tplot('magnull_null_bary_dist') + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file