From 63169edc6a237f03836fc4531dd519505980af0c Mon Sep 17 00:00:00 2001 From: jtgasparik Date: Tue, 10 Feb 2026 17:25:41 -0500 Subject: [PATCH 1/8] created NDM module --- pysp2/util/normalized_derivative_method.py | 30 ++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 pysp2/util/normalized_derivative_method.py diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py new file mode 100644 index 0000000..2c9b297 --- /dev/null +++ b/pysp2/util/normalized_derivative_method.py @@ -0,0 +1,30 @@ +import numpy as np +from scipy.optimize import curve_fit +import xarray as xr + +def central_difference_scheme(S, min_signal=1e-12): + """ + Compute normalized derivative S'(t) / S(t) using the + central difference scheme (Moteki & Kondo, Eq. A.2). + - Fourth-order central difference + Edge cases: + - Fourth-order forward difference at the beginning + - Fourth-order backward difference at the end + + Parameters + ---------- + S: xarray Dataset + The scattering signal dataset. + min_signal: float + Minimum S(t) signal value to avoid division by zero. + + Returns + ------- + dSdt : ndarray + Fourth-order numerical derivative S'(t). + norm_deriv: xarray Dataset + Normalized derivative S'(t) / S(t). + """ + spectra = ds.isel( + + return dSdt, norm_deriv \ No newline at end of file From efac3d797913ddca3361da772151d48288af2194 Mon Sep 17 00:00:00 2001 From: jtgasparik Date: Wed, 11 Feb 2026 17:36:33 -0500 Subject: [PATCH 2/8] not building cleanly yet but some progress --- pysp2/util/__init__.py | 1 + pysp2/util/normalized_derivative_method.py | 66 ++++++++++++++++++++-- tests/test_ndm.py | 14 +++++ 3 files changed, 77 insertions(+), 4 deletions(-) create mode 100644 tests/test_ndm.py diff --git a/pysp2/util/__init__.py b/pysp2/util/__init__.py index 68bba37..0a224ee 100644 --- a/pysp2/util/__init__.py +++ b/pysp2/util/__init__.py @@ -18,3 +18,4 @@ from .particle_properties import calc_diams_masses, process_psds from .deadtime import deadtime from .leo_fit import beam_shape,leo_fit +from .normalized_derivative_method import central_difference diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py index 2c9b297..6f22082 100644 --- a/pysp2/util/normalized_derivative_method.py +++ b/pysp2/util/normalized_derivative_method.py @@ -2,7 +2,7 @@ from scipy.optimize import curve_fit import xarray as xr -def central_difference_scheme(S, min_signal=1e-12): +def central_difference(S, min_signal=1e-12, num_records=None): """ Compute normalized derivative S'(t) / S(t) using the central difference scheme (Moteki & Kondo, Eq. A.2). @@ -17,14 +17,72 @@ def central_difference_scheme(S, min_signal=1e-12): The scattering signal dataset. min_signal: float Minimum S(t) signal value to avoid division by zero. + num_records: int or None + Only process first num_records datapoints. Set to + None to process all records. Returns ------- - dSdt : ndarray + dSdt : xarray Dataset Fourth-order numerical derivative S'(t). norm_deriv: xarray Dataset Normalized derivative S'(t) / S(t). """ - spectra = ds.isel( + + if num_records is None: + num_records = S.dims['event_index'] + + time = S['time'].values + if time.ndim == 2: + time = time[0] # Take the first event's time axis if 2D + dt = 1 + n_time = 1 + for i_records in range(num_records): + S_i = S.isel(event_index=i_records) + S_i = S_i[['Data_ch0', 'Data_ch4']] + #S_i['Data_ch0'] = S_i['Data_ch0'].clip(min=min_signal) + #S_i['Data_ch4'] = S_i['Data_ch4'].clip(min=min_signal) + + dSdt = xr.full_like(S_i, fill_value=np.nan) + + # --- Forward difference (i = 0, 1) --- + dSdt[0] = ( + -25*S_i[0] + 48*S_i[1] - 36*S_i[2] + 16*S_i[3] - 3*S_i[4] + ) / (12*dt) + + dSdt[1] = ( + -25*S_i[1] + 48*S_i[2] - 36*S_i[3] + 16*S_i[4] - 3*S_i[5] + ) / (12*dt) + + # --- Central difference (i = 2 ... n-3) --- + for i in range(2, n_time-2): + dSdt[i] = ( + -S_i[i+2] + 8*S_i[i+1] + - 8*S_i[i-1] + S_i[i-2] + ) / (12*dt) + + # --- Backward difference (i = n-2, n-1) --- + dSdt[n_time-2] = ( + 25*S_i[n_time-2] - 48*S_i[n_time-3] + + 36*S_i[n_time-4] - 16*S_i[n_time-5] + + 3*S_i[n_time-6] + ) / (12*dt) + + dSdt[n_time-1] = ( + 25*S_i[n_time-1] - 48*S_i[n_time-2] + + 36*S_i[n_time-3] - 16*S_i[n_time-4] + + 3*S_i[n_time-5] + ) / (12*dt) + + # Compute the normalized derivative + norm_deriv = dSdt / S_i + + dSdt[i_records, :] = dSdt + norm_deriv[i_records, :] = norm_deriv + print(f'Processed record {i_records+1}/{num_records}') + print(f'Min S(t) = {S_i.min().values}, Max S(t) = {S_i.max().values}') + print(f'Min dS/dt = {dSdt.min().values}, Max dS/dt = {dSdt.max().values}') + print(f'Min normalized derivative = {norm_deriv.min().values}, Max normalized derivative = {norm_deriv.max().values}') + + return dSdt, norm_deriv - return dSdt, norm_deriv \ No newline at end of file diff --git a/tests/test_ndm.py b/tests/test_ndm.py new file mode 100644 index 0000000..9555bbc --- /dev/null +++ b/tests/test_ndm.py @@ -0,0 +1,14 @@ +import pysp2 +import numpy as np + +def test_central_difference(): + my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B) + my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI) + my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False) + dSdt, norm_deriv = pysp2.util.central_difference(my_binary) + + # Check that the outputs have the expected dimensions and contain finite values + assert dSdt.dims == my_binary.dims + assert norm_deriv.dims == my_binary.dims + assert np.isfinite(dSdt).all() + assert np.isfinite(norm_deriv).all() \ No newline at end of file From 906fc74e7f8ee63adf853ae734e1b9b4608b200d Mon Sep 17 00:00:00 2001 From: jtgasparik Date: Thu, 12 Feb 2026 10:43:31 -0500 Subject: [PATCH 3/8] passing now, need to extend to all events --- pysp2/util/normalized_derivative_method.py | 80 ++++++++-------------- tests/test_ndm.py | 11 +-- 2 files changed, 35 insertions(+), 56 deletions(-) diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py index 6f22082..769194e 100644 --- a/pysp2/util/normalized_derivative_method.py +++ b/pysp2/util/normalized_derivative_method.py @@ -2,7 +2,7 @@ from scipy.optimize import curve_fit import xarray as xr -def central_difference(S, min_signal=1e-12, num_records=None): +def central_difference(S, min_signal=1e-13, num_records=None): """ Compute normalized derivative S'(t) / S(t) using the central difference scheme (Moteki & Kondo, Eq. A.2). @@ -31,58 +31,36 @@ def central_difference(S, min_signal=1e-12, num_records=None): if num_records is None: num_records = S.dims['event_index'] + print(f'Processing {num_records} records...') time = S['time'].values + print(f'Time:' f' shape: {time.shape}, dtype: {time.dtype}, min: {time.min()}, max: {time.max()}') if time.ndim == 2: time = time[0] # Take the first event's time axis if 2D - dt = 1 - n_time = 1 - for i_records in range(num_records): - S_i = S.isel(event_index=i_records) - S_i = S_i[['Data_ch0', 'Data_ch4']] - #S_i['Data_ch0'] = S_i['Data_ch0'].clip(min=min_signal) - #S_i['Data_ch4'] = S_i['Data_ch4'].clip(min=min_signal) - - dSdt = xr.full_like(S_i, fill_value=np.nan) - - # --- Forward difference (i = 0, 1) --- - dSdt[0] = ( - -25*S_i[0] + 48*S_i[1] - 36*S_i[2] + 16*S_i[3] - 3*S_i[4] - ) / (12*dt) - - dSdt[1] = ( - -25*S_i[1] + 48*S_i[2] - 36*S_i[3] + 16*S_i[4] - 3*S_i[5] - ) / (12*dt) - - # --- Central difference (i = 2 ... n-3) --- - for i in range(2, n_time-2): - dSdt[i] = ( - -S_i[i+2] + 8*S_i[i+1] - - 8*S_i[i-1] + S_i[i-2] - ) / (12*dt) - - # --- Backward difference (i = n-2, n-1) --- - dSdt[n_time-2] = ( - 25*S_i[n_time-2] - 48*S_i[n_time-3] - + 36*S_i[n_time-4] - 16*S_i[n_time-5] - + 3*S_i[n_time-6] - ) / (12*dt) - - dSdt[n_time-1] = ( - 25*S_i[n_time-1] - 48*S_i[n_time-2] - + 36*S_i[n_time-3] - 16*S_i[n_time-4] - + 3*S_i[n_time-5] - ) / (12*dt) - - # Compute the normalized derivative - norm_deriv = dSdt / S_i - - dSdt[i_records, :] = dSdt - norm_deriv[i_records, :] = norm_deriv - print(f'Processed record {i_records+1}/{num_records}') - print(f'Min S(t) = {S_i.min().values}, Max S(t) = {S_i.max().values}') - print(f'Min dS/dt = {dSdt.min().values}, Max dS/dt = {dSdt.max().values}') - print(f'Min normalized derivative = {norm_deriv.min().values}, Max normalized derivative = {norm_deriv.max().values}') - - return dSdt, norm_deriv + dt = np.diff(time).mean() + if np.issubdtype(dt.dtype, np.timedelta64): + dt = dt / np.timedelta64(1, 's') # convert to seconds as float + print(f'Calculated dt: {dt}') + if dt == 0 or np.isnan(dt): + print('Warning: dt is zero or NaN!') + n_time = 100 + print(f'Number of time points: {n_time}') + # Only process the first record and 'Data_ch0' for now + S_i = S.isel(event_index=0) + ch = 'Data_ch0' + y = S_i[ch].values + print(f'Input signal y: shape: {y.shape}, dtype: {y.dtype}, min: {y.min()}, max: {y.max()}') + dSdt = np.full_like(y, np.nan, dtype=np.float64) + # --- Forward difference (i = 0, 1) --- + dSdt[0] = (-25*y[0] + 48*y[1] - 36*y[2] + 16*y[3] - 3*y[4]) / (12*dt) + dSdt[1] = (-25*y[1] + 48*y[2] - 36*y[3] + 16*y[4] - 3*y[5]) / (12*dt) + # --- Central difference (i = 2 ... n-3) --- + for i in range(2, n_time-2): + dSdt[i] = (-y[i+2] + 8*y[i+1] - 8*y[i-1] + y[i-2]) / (12*dt) + # --- Backward difference (i = n-2, n-1) --- + dSdt[n_time-2] = (25*y[n_time-2] - 48*y[n_time-3] + 36*y[n_time-4] - 16*y[n_time-5] + 3*y[n_time-6]) / (12*dt) + dSdt[n_time-1] = (25*y[n_time-1] - 48*y[n_time-2] + 36*y[n_time-3] - 16*y[n_time-4] + 3*y[n_time-5]) / (12*dt) + # Wrap as DataArray with same coords as input + dSdt_da = xr.DataArray(dSdt, dims=S_i[ch].dims, coords=S_i[ch].coords, name=f'd{ch}_dt') + return dSdt_da diff --git a/tests/test_ndm.py b/tests/test_ndm.py index 9555bbc..d4042ac 100644 --- a/tests/test_ndm.py +++ b/tests/test_ndm.py @@ -5,10 +5,11 @@ def test_central_difference(): my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B) my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI) my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False) - dSdt, norm_deriv = pysp2.util.central_difference(my_binary) + dSdt = pysp2.util.central_difference(my_binary) # Check that the outputs have the expected dimensions and contain finite values - assert dSdt.dims == my_binary.dims - assert norm_deriv.dims == my_binary.dims - assert np.isfinite(dSdt).all() - assert np.isfinite(norm_deriv).all() \ No newline at end of file + print(f'dSdt shape: {dSdt.dims}') + print(f'my_binary shape: {my_binary.dims}') + + #assert dSdt.dims == my_binary.dims + assert np.isfinite(dSdt).all() \ No newline at end of file From 2a7b765308b52bdba8d02bddd2dde204c204dcf0 Mon Sep 17 00:00:00 2001 From: jtgasparik Date: Thu, 12 Feb 2026 15:59:10 -0500 Subject: [PATCH 4/8] i think the function is working as intended, need to add tests to check --- pysp2/util/normalized_derivative_method.py | 63 +++++++++++----------- 1 file changed, 32 insertions(+), 31 deletions(-) diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py index 769194e..4b87151 100644 --- a/pysp2/util/normalized_derivative_method.py +++ b/pysp2/util/normalized_derivative_method.py @@ -25,42 +25,43 @@ def central_difference(S, min_signal=1e-13, num_records=None): ------- dSdt : xarray Dataset Fourth-order numerical derivative S'(t). - norm_deriv: xarray Dataset - Normalized derivative S'(t) / S(t). """ if num_records is None: num_records = S.dims['event_index'] print(f'Processing {num_records} records...') - time = S['time'].values - print(f'Time:' f' shape: {time.shape}, dtype: {time.dtype}, min: {time.min()}, max: {time.max()}') - if time.ndim == 2: - time = time[0] # Take the first event's time axis if 2D - dt = np.diff(time).mean() - if np.issubdtype(dt.dtype, np.timedelta64): - dt = dt / np.timedelta64(1, 's') # convert to seconds as float - print(f'Calculated dt: {dt}') - if dt == 0 or np.isnan(dt): - print('Warning: dt is zero or NaN!') - n_time = 100 - print(f'Number of time points: {n_time}') + dt = 200e-9 # Time step in seconds + n_time = 100 # Number of time bins (assuming 100 bins from 0 to 100) + # Only process the first record and 'Data_ch0' for now - S_i = S.isel(event_index=0) - ch = 'Data_ch0' - y = S_i[ch].values - print(f'Input signal y: shape: {y.shape}, dtype: {y.dtype}, min: {y.min()}, max: {y.max()}') - dSdt = np.full_like(y, np.nan, dtype=np.float64) - # --- Forward difference (i = 0, 1) --- - dSdt[0] = (-25*y[0] + 48*y[1] - 36*y[2] + 16*y[3] - 3*y[4]) / (12*dt) - dSdt[1] = (-25*y[1] + 48*y[2] - 36*y[3] + 16*y[4] - 3*y[5]) / (12*dt) - # --- Central difference (i = 2 ... n-3) --- - for i in range(2, n_time-2): - dSdt[i] = (-y[i+2] + 8*y[i+1] - 8*y[i-1] + y[i-2]) / (12*dt) - # --- Backward difference (i = n-2, n-1) --- - dSdt[n_time-2] = (25*y[n_time-2] - 48*y[n_time-3] + 36*y[n_time-4] - 16*y[n_time-5] + 3*y[n_time-6]) / (12*dt) - dSdt[n_time-1] = (25*y[n_time-1] - 48*y[n_time-2] + 36*y[n_time-3] - 16*y[n_time-4] + 3*y[n_time-5]) / (12*dt) - # Wrap as DataArray with same coords as input - dSdt_da = xr.DataArray(dSdt, dims=S_i[ch].dims, coords=S_i[ch].coords, name=f'd{ch}_dt') - return dSdt_da + channels = ['Data_ch0', 'Data_ch4'] + dSdt = {} + for ch in channels: + dSdt_ch = np.full((num_records, n_time), np.nan, dtype=np.float64) + for i_event in range(num_records): + S_i = S.isel(event_index=i_event) + y = S_i[ch].values + dSdt_i = np.full_like(y, np.nan, dtype=np.float64) + # --- Forward difference (i = 0, 1) --- + dSdt_i[0] = (-25*y[0] + 48*y[1] - 36*y[2] + 16*y[3] - 3*y[4]) / (12*dt) + dSdt_i[1] = (-25*y[1] + 48*y[2] - 36*y[3] + 16*y[4] - 3*y[5]) / (12*dt) + # --- Central difference (i = 2 ... n-3) --- + for i in range(2, n_time-2): + dSdt_i[i] = (-y[i+2] + 8*y[i+1] - 8*y[i-1] + y[i-2]) / (12*dt) + # --- Backward difference (i = n-2, n-1) --- + dSdt_i[n_time-2] = (25*y[n_time-2] - 48*y[n_time-3] + 36*y[n_time-4] - 16*y[n_time-5] + 3*y[n_time-6]) / (12*dt) + dSdt_i[n_time-1] = (25*y[n_time-1] - 48*y[n_time-2] + 36*y[n_time-3] - 16*y[n_time-4] + 3*y[n_time-5]) / (12*dt) + dSdt_ch[i_event, :] = dSdt_i + # Wrap as DataArray with event_index and time dims + dSdt_ch0_ch4 = xr.DataArray( + dSdt_ch, + dims=('event_index', 'time'), + coords={'event_index': S['event_index'], 'time': S['time']}, + name=f'd{ch}_dt' + ) + dSdt[ch] = dSdt_ch0_ch4 + # Return as a Dataset with both channels + print('dSdt_dict', dSdt) + return xr.Dataset(dSdt) From e68e08c53dc3aca88862efa2392d76f679504c9e Mon Sep 17 00:00:00 2001 From: jtgasparik Date: Fri, 13 Feb 2026 10:37:38 -0500 Subject: [PATCH 5/8] centralized_difference function done and tests written --- pysp2/util/normalized_derivative_method.py | 16 ++++++++-------- tests/test_ndm.py | 11 ++++++----- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py index 4b87151..1b888fa 100644 --- a/pysp2/util/normalized_derivative_method.py +++ b/pysp2/util/normalized_derivative_method.py @@ -2,7 +2,7 @@ from scipy.optimize import curve_fit import xarray as xr -def central_difference(S, min_signal=1e-13, num_records=None): +def central_difference(S, num_records=None): """ Compute normalized derivative S'(t) / S(t) using the central difference scheme (Moteki & Kondo, Eq. A.2). @@ -15,8 +15,6 @@ def central_difference(S, min_signal=1e-13, num_records=None): ---------- S: xarray Dataset The scattering signal dataset. - min_signal: float - Minimum S(t) signal value to avoid division by zero. num_records: int or None Only process first num_records datapoints. Set to None to process all records. @@ -29,12 +27,11 @@ def central_difference(S, min_signal=1e-13, num_records=None): if num_records is None: num_records = S.dims['event_index'] - print(f'Processing {num_records} records...') dt = 200e-9 # Time step in seconds n_time = 100 # Number of time bins (assuming 100 bins from 0 to 100) - # Only process the first record and 'Data_ch0' for now + # Process 'Data_ch0' and 'Data_ch4' only channels = ['Data_ch0', 'Data_ch4'] dSdt = {} for ch in channels: @@ -43,15 +40,19 @@ def central_difference(S, min_signal=1e-13, num_records=None): S_i = S.isel(event_index=i_event) y = S_i[ch].values dSdt_i = np.full_like(y, np.nan, dtype=np.float64) + # --- Forward difference (i = 0, 1) --- dSdt_i[0] = (-25*y[0] + 48*y[1] - 36*y[2] + 16*y[3] - 3*y[4]) / (12*dt) dSdt_i[1] = (-25*y[1] + 48*y[2] - 36*y[3] + 16*y[4] - 3*y[5]) / (12*dt) - # --- Central difference (i = 2 ... n-3) --- + + # --- Central difference (i = 2 ... n-2) --- for i in range(2, n_time-2): dSdt_i[i] = (-y[i+2] + 8*y[i+1] - 8*y[i-1] + y[i-2]) / (12*dt) + # --- Backward difference (i = n-2, n-1) --- dSdt_i[n_time-2] = (25*y[n_time-2] - 48*y[n_time-3] + 36*y[n_time-4] - 16*y[n_time-5] + 3*y[n_time-6]) / (12*dt) dSdt_i[n_time-1] = (25*y[n_time-1] - 48*y[n_time-2] + 36*y[n_time-3] - 16*y[n_time-4] + 3*y[n_time-5]) / (12*dt) + dSdt_ch[i_event, :] = dSdt_i # Wrap as DataArray with event_index and time dims dSdt_ch0_ch4 = xr.DataArray( @@ -61,7 +62,6 @@ def central_difference(S, min_signal=1e-13, num_records=None): name=f'd{ch}_dt' ) dSdt[ch] = dSdt_ch0_ch4 - # Return as a Dataset with both channels - print('dSdt_dict', dSdt) + return xr.Dataset(dSdt) diff --git a/tests/test_ndm.py b/tests/test_ndm.py index d4042ac..8867890 100644 --- a/tests/test_ndm.py +++ b/tests/test_ndm.py @@ -7,9 +7,10 @@ def test_central_difference(): my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False) dSdt = pysp2.util.central_difference(my_binary) - # Check that the outputs have the expected dimensions and contain finite values - print(f'dSdt shape: {dSdt.dims}') - print(f'my_binary shape: {my_binary.dims}') - - #assert dSdt.dims == my_binary.dims + np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=0).item(), + 8.3333333333e6, decimal=2) + np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=99).item(), + 7.166666666e7, decimal=2) + np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=19).item(), + 1.5e7, decimal=2) assert np.isfinite(dSdt).all() \ No newline at end of file From 2e4a93d4fe6f1ce7db2b32b2a2d6001c2c9622db Mon Sep 17 00:00:00 2001 From: jtgasparik Date: Fri, 13 Feb 2026 15:06:32 -0500 Subject: [PATCH 6/8] removing loop for computational efficiency --- pysp2/util/normalized_derivative_method.py | 70 +++++++--------------- 1 file changed, 23 insertions(+), 47 deletions(-) diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py index 1b888fa..a312fbe 100644 --- a/pysp2/util/normalized_derivative_method.py +++ b/pysp2/util/normalized_derivative_method.py @@ -3,65 +3,41 @@ import xarray as xr def central_difference(S, num_records=None): - """ - Compute normalized derivative S'(t) / S(t) using the - central difference scheme (Moteki & Kondo, Eq. A.2). - - Fourth-order central difference - Edge cases: - - Fourth-order forward difference at the beginning - - Fourth-order backward difference at the end - Parameters - ---------- - S: xarray Dataset - The scattering signal dataset. - num_records: int or None - Only process first num_records datapoints. Set to - None to process all records. - - Returns - ------- - dSdt : xarray Dataset - Fourth-order numerical derivative S'(t). - """ - if num_records is None: - num_records = S.dims['event_index'] - - dt = 200e-9 # Time step in seconds - n_time = 100 # Number of time bins (assuming 100 bins from 0 to 100) + num_records = S.sizes['event_index'] - # Process 'Data_ch0' and 'Data_ch4' only + dt = 200e-9 channels = ['Data_ch0', 'Data_ch4'] dSdt = {} + for ch in channels: - dSdt_ch = np.full((num_records, n_time), np.nan, dtype=np.float64) - for i_event in range(num_records): - S_i = S.isel(event_index=i_event) - y = S_i[ch].values - dSdt_i = np.full_like(y, np.nan, dtype=np.float64) + y = S[ch].isel(event_index=slice(0, num_records)).values + d = np.full_like(y, np.nan, dtype=np.float64) + + # Interior points (vectorized) + d[:, 2:-2] = ( + -y[:, 4:] + 8*y[:, 3:-1] - 8*y[:, 1:-3] + y[:, 0:-4] + ) / (12 * dt) - # --- Forward difference (i = 0, 1) --- - dSdt_i[0] = (-25*y[0] + 48*y[1] - 36*y[2] + 16*y[3] - 3*y[4]) / (12*dt) - dSdt_i[1] = (-25*y[1] + 48*y[2] - 36*y[3] + 16*y[4] - 3*y[5]) / (12*dt) + # Forward edges + d[:, 0] = (-25*y[:, 0] + 48*y[:, 1] - 36*y[:, 2] + 16*y[:, 3] - 3*y[:, 4]) / (12*dt) + d[:, 1] = (-25*y[:, 1] + 48*y[:, 2] - 36*y[:, 3] + 16*y[:, 4] - 3*y[:, 5]) / (12*dt) - # --- Central difference (i = 2 ... n-2) --- - for i in range(2, n_time-2): - dSdt_i[i] = (-y[i+2] + 8*y[i+1] - 8*y[i-1] + y[i-2]) / (12*dt) + # Backward edges + n = y.shape[1] + d[:, n-2] = (25*y[:, n-2] - 48*y[:, n-3] + 36*y[:, n-4] - 16*y[:, n-5] + 3*y[:, n-6]) / (12*dt) + d[:, n-1] = (25*y[:, n-1] - 48*y[:, n-2] + 36*y[:, n-3] - 16*y[:, n-4] + 3*y[:, n-5]) / (12*dt) - # --- Backward difference (i = n-2, n-1) --- - dSdt_i[n_time-2] = (25*y[n_time-2] - 48*y[n_time-3] + 36*y[n_time-4] - 16*y[n_time-5] + 3*y[n_time-6]) / (12*dt) - dSdt_i[n_time-1] = (25*y[n_time-1] - 48*y[n_time-2] + 36*y[n_time-3] - 16*y[n_time-4] + 3*y[n_time-5]) / (12*dt) - - dSdt_ch[i_event, :] = dSdt_i - # Wrap as DataArray with event_index and time dims - dSdt_ch0_ch4 = xr.DataArray( - dSdt_ch, + dSdt[ch] = xr.DataArray( + d, dims=('event_index', 'time'), - coords={'event_index': S['event_index'], 'time': S['time']}, + coords={ + 'event_index': S['event_index'][:num_records], + 'time': S['time'] + }, name=f'd{ch}_dt' ) - dSdt[ch] = dSdt_ch0_ch4 return xr.Dataset(dSdt) From 62c126bee32760ddf0288909cc23c5a571d5b5fe Mon Sep 17 00:00:00 2001 From: jtgasparik Date: Fri, 13 Feb 2026 15:10:56 -0500 Subject: [PATCH 7/8] added description back --- pysp2/util/normalized_derivative_method.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py index a312fbe..0995f2d 100644 --- a/pysp2/util/normalized_derivative_method.py +++ b/pysp2/util/normalized_derivative_method.py @@ -4,6 +4,28 @@ def central_difference(S, num_records=None): + """ Compute normalized derivative S'(t) / S(t) using the + central difference scheme (Moteki & Kondo, Eq. A.2). + Interior points: + - Fourth-order central difference + Edge cases: + - First two points: fourth-order forward difference + - Last two points: fourth-order backward difference + + Parameters + ---------- + S: xarray Dataset + The scattering signal dataset. + num_records: int or None + Only process first num_records datapoints. + Set to None to process all records. + + Returns + ------- + dSdt : xarray Dataset + Fourth-order numerical derivative S'(t). + """ + if num_records is None: num_records = S.sizes['event_index'] From 709cec60287471b25939b037affe481e7a5808c2 Mon Sep 17 00:00:00 2001 From: jtgasparik Date: Fri, 13 Feb 2026 15:15:55 -0500 Subject: [PATCH 8/8] some formatting stuff --- pysp2/util/normalized_derivative_method.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py index 0995f2d..b89ae66 100644 --- a/pysp2/util/normalized_derivative_method.py +++ b/pysp2/util/normalized_derivative_method.py @@ -4,7 +4,8 @@ def central_difference(S, num_records=None): - """ Compute normalized derivative S'(t) / S(t) using the + """ + Compute normalized derivative S'(t) / S(t) using the central difference scheme (Moteki & Kondo, Eq. A.2). Interior points: - Fourth-order central difference