""" Tools for fatigue analysis Main functions: - equivalent_load: calculate damage equivalent load for a given signal - find_range_count: returns range and number of cycles for a given signal Subfunctions: - eq_load: calculate equivalent loads using one of the two rain flow counting methods - cycle_matrix: calculates a matrix of cycles (binned on amplitude and mean value) - eq_load_and_cycles: calculate eq_loads of multiple time series (e.g. life time equivalent load) Main aglorithms for rain flow counting: - rainflow_windap: taken from [2], based on [3] - rainflow_astm: taken from [2], based [4] - fatpack: using Fatpack - Python package References: [1] Hayman (2012) MLife theory manual for Version 1.00 [2] Wind energy toolbox, wetb.fatigue_tools, DTU wind energy, Denmark [3] "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A [4] Adam Nieslony - Rainflow Counting Algorithm, MATLAB Central File Exchange """ import warnings import numpy as np import pandas as pd fast_out_file = "WT_Sample_Loads.out" rainflow_func_dict = ["rainflow_windap", "rainflow_astm"] __all__ = ['equivalent_load', 'find_range_count'] __all__ += ['rainflow_astm', 'rainflow_windap','eq_load','eq_load_and_cycles', 'cycle_matrix','cycle_matrix2'] class SignalConstantError(Exception): pass def equivalent_load(time, signal, m=3, Teq=1, bins=100, method='rainflow_astm', meanBin=True, binStartAt0=False, outputMore=False, debug=False): """Equivalent load calculation Calculate the damage equivalent load for a given signal and a given Wohler exponent INPUTS - time : array-like, the time values corresponding to the signal (s) - signals : array-like, the load signal - m : Wohler exponent (default is 3) - Teq : The equivalent period (Default 1, for 1 Hz) - bins : Number of bins in rainflow count histogram - method: rain flow counting algorithm: 'rainflow_windap', 'rainflow_astm' or 'fatpack' - meanBin: if True, use the mean of the ranges within a bin (recommended) otherwise use the middle of the bin (not recommended). - binStartAt0: if True bins start at zero. Otherwise, start a lowest range - outputMore: if True, returns range, cycles and bins as well OUTPUTS - Leq : the equivalent load for given m and Teq or (if outputMore is True ) - Leq, S, N, bins, DELi: - S: ranges - N: cycles - bins: bin edges - DELi: component 'i' of the DEL (for cycle i) """ time = np.asarray(time) signal = np.asarray(signal) # Remove nan, might not be the cleanest b = ~np.isnan(signal) signal = signal[b] time = time[b] try: if len(time)<=1: raise Exception("Length of time array <= 1\n") if type(time[0]) is np.datetime64: T = (time[-1] - time[0]) / np.timedelta64(1, 's') else: # time length of signal (s). Will fail for signal of length 1 T = time[-1]-time[0] if T==0: raise Exception("Invalid time-step!\n") neq = T/Teq # number of equivalent periods, see Eq. (26) of [1] # --- Range (S) and counts (N) N, S, bins = find_range_count(signal, bins=bins, method=method, meanBin=meanBin, binStartAt0=binStartAt0 ) print("\n---Printing counts N ---") print(N) print("\n---Printing range S ---") print(S) # --- get DEL DELi = S**m * N / neq Leq = DELi.sum() ** (1/m) # See e.g. eq. (30) of [1] except Exception as e: print(f"Error in equivalent_load(): {type(e).__name__}: {e}") #print("Error:", e) raise ''' except: if outputMore: return np.nan, np.nan, np.nan, np.nan, np.nan else: print("Something went wrong!\n") return np.nan ''' if debug: for i,(b,n,s,DEL) in enumerate(zip(bins, N, S, DELi)): if n>0: print('Bin {:3d}: [{:6.1f}-{:6.1f}] Mid:{:6.1f} - \ Mean:{:6.1f} Counts:{:4.1f} DEL:{:8.1f} Fraction:{:3.0f}%'\ .format(i,b,bins[i+1],(b+bins[i+1])/2,s,n,DEL,DEL/Leq**m*100)) if outputMore: return Leq, S, N, bins, DELi else: return Leq def find_range_count(signal, bins, method='rainflow_astm', meanBin=True, binStartAt0=True): """ Returns number of cycles `N` for each range range `S` Equidistant bins are setup based on the min/max of the signal. INPUTS: - signal: array - bins : 1d-array, int If bins is a sequence, left edges (and the rightmost edge) of the bins. If bins is an int, a sequence is created dividing the range `min`--`max` of signal into `bins` number of equally sized bins. OUTPUTS: - N: number of cycles for each bin - S: Ranges for each bin S is either the center of the bin (meanBin=False) or S is the mean of the ranges within this bin (meanBin=True) - S_bin_edges: edges of the bins """ if method in rainflow_func_dict: rainflow_func = globals()[method] try: print(f"Calculating Markow load cycle matrix. \n") #cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges N, S, S_bin_edges, _, _ = cycle_matrix( signal, ampl_bins=bins, mean_bins=1, rainflow_func=rainflow_func, binStartAt0=binStartAt0 ) except SignalConstantError: print("Markow load cycle matrix calculation failed! \n") return np.nan, np.nan, np.nan S_bin_edges = S_bin_edges.flatten() N = N.flatten() S = S.flatten() S_mid = (S_bin_edges[:-1] + S_bin_edges[1:]) / 2 if not meanBin: S=S_mid elif method=='fatpack': print("Specified method {method} not supported, using [fatpack]\n") import fatpack # find rainflow ranges try: ranges = fatpack.find_rainflow_ranges(signal) except IndexError: # Currently fails for constant signal return np.nan, np.nan, np.nan S_bin_edges = create_bins(ranges, bins, binStartAt0=binStartAt0) # --- Using bin_count to get value at center of bins N, S = bin_count(ranges, S_bin_edges, meanBin=meanBin) else: raise NotImplementedError('Rain flow algorithm {}'.format(method)) # Remove NaN b = np.isnan(S) S[b] = 0 N[b] = 0 return N, S, S_bin_edges def create_bins(x, bins, binStartAt0=False): """ Equidistant bins are setup based on the min/max of the x, unless the user provided the bins as a sequence. INPUTS: - x: array - bins : 1d-array, int If bins is a sequence, left edges (and the rightmost edge) of the bins. If bins is an int, a sequence is created dividing the range `min`--`max` of x into `bins` number of equally sized bins. OUTPUTS: - bins: """ if isinstance(bins, int): print("Creating bins!\n") xmax = np.max(x) xmin, xmax = np.min(x), np.max(x) if binStartAt0: xmin = 0 else: print("Calculating min / max of bins.\n") xmin = np.min(x) if xmin == xmax: # This is what is done by histogram - double check xmin = xmin-0.5 xmax = xmax+0.5 bins = np.linspace(xmin, xmax, num=bins + 1) else: print("Non-integer values of bins specified!\n") return bins def bin_count(x, bins, meanBin=True): """ Return counts of x within bins """ if not meanBin: # Use the middle of the bin N, bns = np.histogram(x, bins=bins) S = bns[:-1] + np.diff(bns) / 2. else: bins = create_bins(x, bins, binStartAt0=False) import pandas as pd df = pd.DataFrame(data=x, columns=['x']) xmid = (bins[:-1]+bins[1:])/2 # Adding a column that has bin attribute df['x_mid']= pd.cut(df['x'], bins= bins, labels = xmid ) # Average by bin df2 = df.groupby('x_mid').mean() df['N'] = 1 dfCount = df[['N','x_mid']].groupby('x_mid').sum() df2['N'] = dfCount['N'] # Just in case some bins are missing (will be nan) df2 = df2.reindex(xmid) df2 = df2.fillna(0) S = df2['x'].values N = df2['N'].values return N, S def rainflowcount(sig): #cpdef rainflowcount(np.ndarray[double,ndim=1] sig): """Cython compilable rain ampl_mean count without time analysis This implemementation is based on the c-implementation by Adam Nieslony found at the MATLAB Central File Exchange References ---------- Adam Nieslony, "Determination of fragments of multiaxial service loading strongly influencing the fatigue of machine components," Mechanical Systems and Signal Processing 23, no. 8 (2009): 2712-2721. And is based on the following standard: ASTM E 1049-85 (Reapproved 1997), Standard practices for cycle counting in fatigue analysis, in: Annual Book of ASTM Standards, vol. 03.01, ASTM, Philadelphia, 1999, pp. 710-718. Copyright (c) 1999-2002 by Adam Nieslony Ported to Cython compilable Python by Mads M Pedersen. In addition peak amplitude is changed to peak to peak amplitude """ a = [] sig_ptr = 0 ampl_mean = [] for _ in range(len(sig)): a.append(sig[sig_ptr]) sig_ptr += 1 while len(a) > 2 and abs(a[-3] - a[-2]) <= abs(a[-2] - a[-1]): ampl = abs(a[-3] - a[-2]) mean = (a[-3] + a[-2]) / 2; if len(a) == 3: del a[0] if ampl > 0: ampl_mean.append((ampl, mean)) elif len(a) > 3: del a[-3:-1] if ampl > 0: ampl_mean.append((ampl, mean)) ampl_mean.append((ampl, mean)) for index in range(len(a) - 1): ampl = abs(a[index] - a[index + 1]) mean = (a[index] + a[index + 1]) / 2; if ampl > 0: ampl_mean.append((ampl, mean)) return ampl_mean def find_extremes(signal): """return indexes of local minima and maxima plus first and last element of signal""" # sign of gradient sign_grad = np.int8(np.sign(np.diff(signal))) # remove plateaus(sign_grad==0) by # sign_grad[plateau_index]=sign_grad[plateau_index-1] plateau_indexes, = np.where(sign_grad == 0) if len(plateau_indexes) > 0 and plateau_indexes[0] == 0: # first element is a plateau if len(plateau_indexes) == len(sign_grad): # All values are equal to crossing level! return np.array([0]) # set first element = first element which is not a # plateau and delete plateau index i = 0 while sign_grad[i] == 0: i += 1 sign_grad[0] = sign_grad[i] plateau_indexes = np.delete(plateau_indexes, 0) for pi in plateau_indexes.tolist(): sign_grad[pi] = sign_grad[pi - 1] extremes, = np.where(np.r_[1, (sign_grad[1:] * sign_grad[:-1] < 0), 1]) return signal[extremes] def peak_trough(x, R): """ Returns list of local maxima/minima. x: 1-dimensional numpy array containing signal R: Thresshold (minimum difference between succeeding min and max This routine is implemented directly as described in "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A """ BEGIN, MINZO, MAXZO, ENDZO = 0, 1, 2, 3 S = np.zeros(x.shape[0] + 1, dtype=int) L = x.shape[0] goto = BEGIN while 1: if goto == BEGIN: trough = x[0] peak = x[0] i = 0 p = 1 f = 0 while goto == BEGIN: i += 1 if i == L: goto = ENDZO continue else: if x[i] > peak: peak = x[i] if peak - trough >= R: S[p] = trough goto = MAXZO continue elif x[i] < trough: trough = x[i] if peak - trough >= R: S[p] = peak goto = MINZO continue elif goto == MINZO: f = -1 while goto == MINZO: i += 1 if i == L: goto = ENDZO continue else: if x[i] < trough: trough = x[i] else: if x[i] - trough >= R: p += 1 S[p] = trough peak = x[i] goto = MAXZO continue elif goto == MAXZO: f = 1 while goto == MAXZO: i += 1 if i == L: goto = ENDZO continue else: if x[i] > peak: peak = x[i] else: if peak - x[i] >= R: p += 1 S[p] = peak trough = x[i] goto = MINZO continue elif goto == ENDZO: n = p + 1 if abs(f) == 1: if f == 1: S[n] = peak else: S[n] = trough else: S[n] = (trough + peak) / 2 S = S[1:n + 1] return S def pair_range_amplitude_mean(x): # cpdef pair_range(np.ndarray[long,ndim=1] x): """ Returns a list of half-cycle-amplitudes x: Peak-Trough sequence (integer list of local minima and maxima) This routine is implemented according to "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A except that a list of half-cycle-amplitudes are returned instead of a from_level-to_level-matrix """ x = x - np.min(x) k = np.max(x) n = x.shape[0] S = np.zeros(n + 1) ampl_mean = [] A = np.zeros((k + 1, k + 1)) S[1] = x[0] ptr = 1 p = 1 q = 1 f = 0 # phase 1 while True: p += 1 q += 1 # read S[p] = x[ptr] ptr += 1 if q == n: f = 1 while p >= 4: if (S[p - 2] > S[p - 3] and S[p - 1] >= S[p - 3] and S[p] >= S[p - 2]) \ or\ (S[p - 2] < S[p - 3] and S[p - 1] <= S[p - 3] and S[p] <= S[p - 2]): # Extract two intermediate half cycles ampl = abs(S[p - 2] - S[p - 1]) mean = (S[p - 2] + S[p - 1]) / 2 ampl_mean.append((ampl, mean)) ampl_mean.append((ampl, mean)) S[p - 2] = S[p] p -= 2 else: break if f == 0: pass else: break # phase 2 q = 0 while True: q += 1 if p == q: break else: ampl = abs(S[q + 1] - S[q]) mean = (S[q + 1] + S[q]) / 2 ampl_mean.append((ampl, mean)) return ampl_mean def check_signal(signal): # check input data validity if not type(signal).__name__ == 'ndarray': raise TypeError('signal must be ndarray, not: ' + type(signal).__name__) elif len(signal.shape) not in (1, 2): raise TypeError('signal must be 1D or 2D, not: ' + str(len(signal.shape))) if len(signal.shape) == 2: if signal.shape[1] > 1: raise TypeError('signal must have one column only, not: ' + str(signal.shape[1])) if np.min(signal) == np.max(signal): err_str = "Signal is constant, cannot compute DLC and range" raise SignalConstantError(err_str) def rainflow_windap(signal, levels=255, thresshold=(255 / 50)): """Windap equivalent rainflow counting Calculate the amplitude and mean values of half cycles in signal. This algorithms used by this routine is implemented directly as described in "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads", 2. edition 1990, Appendix A Parameters ---------- Signal : array-like The raw signal levels : int, optional The signal is discretize into this number of levels. 255 is equivalent to the implementation in Windap thresshold : int, optional Cycles smaller than this thresshold are ignored 255/50 is equivalent to the implementation in Windap Returns ------- ampl : array-like Peak to peak amplitudes of the half cycles mean : array-like Mean values of the half cycles Examples -------- >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0]) >>> ampl, mean = rainflow_windap(signal) """ check_signal(signal) #type is required by and signal = signal.astype(np.double) if np.all(np.isnan(signal)): return None offset = np.nanmin(signal) signal -= offset if np.nanmax(signal) > 0: gain = np.nanmax(signal) / levels signal = signal / gain signal = np.round(signal).astype(int) # Convert to list of local minima/maxima where difference > threshold sig_ext = peak_trough(signal, thresshold) # rainflow count ampl_mean = pair_range_amplitude_mean(sig_ext) ampl_mean = np.array(ampl_mean) ampl_mean = np.round(ampl_mean / thresshold) * gain * thresshold ampl_mean[:, 1] += offset return ampl_mean.T def rainflow_astm(signal): """Matlab equivalent rainflow counting Calculate the amplitude and mean values of half cycles in signal This implementation is based on [3] Parameters ---------- Signal : array-like The raw signal Returns ------- ampl : array-like peak to peak amplitudes of the half cycles (note that the matlab implementation uses peak amplitude instead of peak to peak) mean : array-like Mean values of the half cycles Examples -------- >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0]) >>> ampl, mean = rainflow_astm(signal) """ check_signal(signal) # type is reuqired by and signal = signal.astype(np.double) # Remove points which is not local minimum/maximum sig_ext = find_extremes(signal) # rainflow count ampl_mean = np.array(rainflowcount(sig_ext)) return np.array(ampl_mean).T def eq_load(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1, rainflow_func=rainflow_windap): """Equivalent load calculation Calculate the equivalent loads for a list of Wohler exponent and number of equivalent loads Parameters ---------- signals : list of tuples or array_like - if list of tuples: list must have format [(sig1_weight, sig1),(sig2_weight, sig1),...] where\n - sigx_weight is the weight of signal x\n - sigx is signal x\n - if array_like: The signal no_bins : int, optional Number of bins in rainflow count histogram m : int, float or array-like, optional Wohler exponent (default is [3, 4, 6, 8, 10, 12]) neq : int, float or array-like, optional The equivalent number of load cycles (default is 1, but normally the time duration in seconds is used) rainflow_func : {rainflow_windap, rainflow_astm}, optional The rainflow counting function to use (default is rainflow_windap) Returns ------- eq_loads : array-like List of lists of equivalent loads for the corresponding equivalent number(s) and Wohler exponents Examples -------- >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0]) >>> eq_load(signal, no_bins=50, neq=[1, 17], m=[3, 4, 6], rainflow_func=rainflow_windap) [[10.31, 9.594, 9.079], # neq = 1, m=[3,4,6] [ 4.010, 4.725, 5.662]], # neq = 17, m=[3,4,6] eq_load( [(.4, signal), (.6, signal)], no_bins=50, neq=[1, 17], m=[3, 4, 6], rainflow_func=rainflow_windap ) [[10.31, 9.594, 9.079], # neq = 1, m=[3,4,6] [ 4.010, 4.725, 5.662]], # neq = 17, m=[3,4,6] """ try: return eq_load_and_cycles(signals, no_bins, m, neq, rainflow_func)[0] except TypeError: return [[np.nan] * len(np.atleast_1d(m))] * len(np.atleast_1d(neq)) def eq_load_and_cycles(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=[10 ** 6, 10 ** 7, 10 ** 8], rainflow_func=rainflow_windap): """Calculate combined fatigue equivalent load Parameters ---------- signals : list of tuples or array_like - if list of tuples: list must have format [(sig1_weight, sig1),(sig2_weight, sig1),...] where\n - sigx_weight is the weight of signal x\n - sigx is signal x\n - if array_like: The signal no_bins : int, optional Number of bins for rainflow counting m : int, float or array-like, optional Wohler exponent (default is [3, 4, 6, 8, 10, 12]) neq : int or array-like, optional Equivalent number, default is [10^6, 10^7, 10^8] rainflow_func : {rainflow_windap, rainflow_astm}, optional The rainflow counting function to use (default is rainflow_windap) Returns ------- eq_loads : array-like List of lists of equivalent loads for the corresponding equivalent number(s) and Wohler exponents cycles : array_like 2d array with shape = (no_ampl_bins, 1) ampl_bin_mean : array_like mean amplitude of the bins ampl_bin_edges Edges of the amplitude bins """ cycles, ampl_bin_mean, ampl_bin_edges, _, _ = cycle_matrix( signals, no_bins, 1, rainflow_func ) if 0: #to be similar to windap ampl_bin_mean = (ampl_bin_edges[:-1] + ampl_bin_edges[1:]) / 2 cycles, ampl_bin_mean = cycles.flatten(), ampl_bin_mean.flatten() with warnings.catch_warnings(): warnings.simplefilter("ignore") eq_loads = [ [((np.nansum(cycles * ampl_bin_mean ** _m) / _neq) ** (1. / _m)) for _m in np.atleast_1d(m)] for _neq in np.atleast_1d(neq)] return eq_loads, cycles, ampl_bin_mean, ampl_bin_edges def cycle_matrix( signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_astm, binStartAt0=True ): """Markow load cycle matrix Calculate the Markow load cycle matrix Parameters ---------- Signals : array-like or list of tuples - if array-like, the raw signal\n - if list of tuples, list of (weight, signal), e.g. [(0.1,sig1), (0.8,sig2), (0.1,sig3)]\n ampl_bins : int or array-like, optional if int, Number of amplitude value bins (default is 10) if array-like, the bin edges for amplitude mean_bins : int or array-like, optional if int, Number of mean value bins (default is 10) if array-like, the bin edges for mea rainflow_func : {rainflow_windap, rainflow_astm}, optional The rainflow counting function to use (default is rainflow_astm) binStartAt0 : boolean Start the bins at 0. Otherwise, start at the min of ranges Returns ------- cycles : ndarray, shape(ampl_bins, mean_bins) A bi-dimensional histogram of load cycles(full cycles). Amplitudes are histogrammed along the first dimension and mean values are histogrammed along the second dimension. ampl_bin_mean : ndarray, shape(ampl_bins,) The average cycle amplitude of the bins ampl_edges : ndarray, shape(ampl_bins+1,) The amplitude bin edges mean_bin_mean : ndarray, shape(ampl_bins,) The average cycle mean of the bins mean_edges : ndarray, shape(mean_bins+1,) The mean bin edges Examples -------- >>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0]) >>> cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges = cycle_matrix(signal) >>> cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges = cycle_matrix([(.4, signal), (.6,signal)]) """ if isinstance(signals[0], tuple): print("Calculation-1: cycle-matrix.\n") weights, ampls, means = np.array( [(np.zeros_like(ampl)+weight, ampl, mean) for weight, signal in signals for ampl,mean in rainflow_func(signal[:]).T], dtype=np.float64).T else: print("Calculation-2: cycle-matrix.\n") ampls, means = rainflow_func(signals[:]) weights = np.ones_like(ampls) if isinstance(ampl_bins, int): ampl_bins = create_bins( ampls[weights>0], ampl_bins, binStartAt0=binStartAt0 ) print("\n---The amplitudes of specified channel are:") print(ampls) cycles, ampl_edges, mean_edges = np.histogram2d( ampls, means, [ampl_bins, mean_bins], weights=weights ) with warnings.catch_warnings(): warnings.simplefilter("ignore") ampl_bin_sum = np.histogram2d( ampls, means, [ampl_bins, mean_bins], weights=weights * ampls )[0] ampl_bin_mean = np.nanmean( ampl_bin_sum / np.where(cycles,cycles,np.nan), 1 ) mean_bin_sum = np.histogram2d( ampls, means, [ampl_bins, mean_bins], weights=weights * means )[0] mean_bin_mean = np.nanmean( mean_bin_sum / np.where(cycles, cycles, np.nan), 1 ) cycles = cycles / 2 # to get full cycles return cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges def cycle_matrix2(signal, nrb_amp, nrb_mean, rainflow_func=rainflow_windap): """ Same as wetb.fatigue.cycle_matrix but bin from min_amp to max_amp instead of 0 to max_amp. Parameters ---------- Signal : ndarray(n) 1D Raw signal array nrb_amp : int Number of bins for the amplitudes nrb_mean : int Number of bins for the means rainflow_func : {rainflow_windap, rainflow_astm}, optional The rainflow counting function to use (default is rainflow_windap) Returns ------- cycles : ndarray, shape(ampl_bins, mean_bins) A bi-dimensional histogram of load cycles(full cycles). Amplitudes are\ histogrammed along the first dimension and mean values are histogrammed along the second dimension. ampl_edges : ndarray, shape(no_bins+1,n) The amplitude bin edges mean_edges : ndarray, shape(no_bins+1,n) The mean bin edges """ bins = [nrb_amp, nrb_mean] ampls, means = rainflow_func(signal) weights = np.ones_like(ampls) cycles, ampl_edges, mean_edges = np.histogram2d(ampls, means, bins, weights=weights) cycles = cycles / 2 # to get full cycles return cycles, ampl_edges, mean_edges # Example usage: df = pd.read_csv(fast_out_file) # Compute equivalent load for one signal and Wohler slope m = 5 # Wohler slope Leq = equivalent_load(df['Time'], df['RootMyc1'], Teq=1, bins=50, m=m, method='rainflow_astm') print('\n---Equivalent Load = ', round(Leq, 4), '\n')