diff --git a/examples/adcp_example.mlx b/examples/adcp_example.mlx index 50d2427ff..a8fbdc83c 100644 Binary files a/examples/adcp_example.mlx and b/examples/adcp_example.mlx differ diff --git a/examples/extreme_response_contour_example.mlx b/examples/extreme_response_contour_example.mlx index 7215ca9a1..7343cede0 100644 Binary files a/examples/extreme_response_contour_example.mlx and b/examples/extreme_response_contour_example.mlx differ diff --git a/mhkit/loads/extreme/short_term_extreme.m b/mhkit/loads/extreme/short_term_extreme.m index b6dc9e8a3..1a8a526ed 100644 --- a/mhkit/loads/extreme/short_term_extreme.m +++ b/mhkit/loads/extreme/short_term_extreme.m @@ -1,80 +1,183 @@ function ste = short_term_extreme(t, data, t_st, type, x, method, options) - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SHORT_TERM_EXTREME Estimate short-term extreme value statistics % -% Estimate the peaks distribution by fitting a Weibull -% distribution to the peaks of the response. +% ste = short_term_extreme(t, data, t_st, type, x, method) % -% To learn more about the methods that can be used in this function, -% refer to the scipy.stats.rv_continuous methods documentation! +% Estimates extreme value distributions from a timeseries of the response +% using either block maxima or peak-over-threshold methods. % -% Parameters -% ---------- -% t : array -% Time array -% data : array -% Response timeseries -% t_st : double or int -% Short time period -% type : string -% Method for estimating the short-term extreme distribution -% x : array or double -% Input for the statistical function/method -% method : str -% Statistical method to apply to resulting data. Options to -% choose from are: "pdf", "cdf", "ppf", or "sf" -% output_py : bool (optional) -% Select if you want to return the native python result for use -% in long term extreme calculations. +% Parameters +% ---------- +% t : time vector +% data : response signal +% t_st : short-term period (scalar) +% type : method string: +% 'peaks_weibull', 'peaks_weibull_tail_fit', +% 'peaks_over_threshold', 'block_maxima_gev', 'block_maxima_gumbel' +% x : vector of values to evaluate the distribution at +% method : 'pdf', 'cdf', 'ppf' (inverse CDF), or 'sf' (1 - CDF) +% options : struct with optional fields: +% .threshold (required for peaks_over_threshold) % -% Returns -% ------- -% ste : array or double -% Probability distribution of the peaks -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Returns +% ------- +% ste : values of the estimated distribution at x arguments - t - data - t_st - type - x - method - options.output_py=false; + t (:,1) double + data (:,1) double + t_st (1,1) double + type (1,:) char + x (:,1) double + method (1,:) char + options.threshold double = NaN +end + +% Error checking +if length(t) ~= length(data) + error('Time and data must be the same length'); end +% Sampling interval and frequency +dt = mean(diff(t)); +fs = 1/dt; + +switch lower(type) + case 'peaks_weibull' + peaks = find_peaks(data); + nst = round(length(peaks) * t_st / duration); + pd_st = fit_weibull(peaks, nst, x); + + case 'peaks_weibull_tail_fit' + peaks = find_peaks(data); + pd_st = fit_weibull_tail(peaks, x); + + case 'peaks_over_threshold' + if isnan(options.threshold) + error('You must provide options.threshold for this method.'); + end + peaks = data(data > options.threshold); + pd_st = fit_pareto(peaks, x); -if ~isa(t,'numeric') - error('ERROR: t must be a double array') + case 'block_maxima_gev' + block_size = round(t_st * fs); + blocks = reshape_truncated(data, block_size); + maxima = max(blocks, [], 1); + pd_st = fit_gumbel(maxima, x); + + case 'block_maxima_gumbel' + block_size = round(t_st * fs); + blocks = reshape_truncated(data, block_size); + maxima = max(blocks, [], 1); + pd_st = fit_gev(maxima, x); + + otherwise + error('Unknown method: %s', type); end -if ~isa(data,'numeric') - error('ERROR: data must be a double array') + +switch lower(method) + case 'pdf' + ste = pd_st.pdf; + case 'cdf' + ste = pd_st.cdf; + case {'ppf', 'inv'} + ste = pd_st.ppf; + case 'sf' + ste = 1 - pd_st.cdf; + case 'expect' + % Numerical expectation: integral of x * pdf(x) + dx = mean(diff(x)); + ste = sum(x .* pd_st.pdf) * dx; + otherwise + error('Invalid method: %s. Choose from "pdf", "cdf", "ppf", "sf", "expect".', method); end -if ~isa(t_st,'numeric') - error('ERROR: t_st must be a double') end -if ~isa(type,'string') - error('ERROR: type must be a string') + + +function peaks = find_peaks(data) + locs = find(data(2:end-1) > data(1:end-2) & data(2:end-1) > data(3:end)) + 1; + peaks = data(locs); end -if ~isa(x,'numeric') - error('ERROR: x must be a double array') + +function blocks = reshape_truncated(data, block_size) + N = floor(length(data)/block_size) * block_size; + blocks = reshape(data(1:N), block_size, []); end -if ~isa(method,'string') - error('ERROR: method must be a string') + +function pd = fit_weibull_tail(data, x) + % Manual peak filtering + locs = find(data(2:end-1) > data(1:end-2) & data(2:end-1) > data(3:end)) + 1; + peaks = data(locs); + + % Tail selection + threshold = prctile(peaks, 95); + tail_peaks = peaks(peaks > threshold); + + % Grid search over Weibull parameters + scale_vals = linspace(mean(tail_peaks)*0.5, mean(tail_peaks)*2, 50); + shape_vals = linspace(0.5, 5, 50); + + min_nll = inf; + + for a = scale_vals + for b = shape_vals + pdf_vals = weibull_pdf(tail_peaks, a, b); + if all(pdf_vals > 0 & isfinite(pdf_vals)) + nll = -sum(log(pdf_vals)); + if nll < min_nll + min_nll = nll; + best_scale = a; + best_shape = b; + end + end + end + end + + pd.type = 'Weibull'; + pd.scale = best_scale; + pd.shape = best_shape; + pd.pdf = weibull_pdf(x, best_scale, best_shape); + pd.cdf = weibull_cdf(x, best_scale, best_shape); + pd.ppf = weibull_ppf(x, best_scale, best_shape); end -py.importlib.import_module('mhkit'); +function pd = fit_weibull(data, nst, x) + % Same fitting approach as the tail, just on all peaks + scale_vals = linspace(mean(data)*0.5, mean(data)*2, 50); + shape_vals = linspace(0.5, 5, 50); -t = py.numpy.array(t); -data = py.numpy.array(data); + min_nll = inf; -result = py.mhkit.loads.extreme.short_term_extreme(t, data, t_st, type); + for a = scale_vals + for b = shape_vals + pdf_vals = weibull_pdf(data, a, b); + if all(pdf_vals > 0 & isfinite(pdf_vals)) + nll = -sum(log(pdf_vals)); + if nll < min_nll + min_nll = nll; + best_scale = a; + best_shape = b; + end + end + end + end -if options.output_py - ste = result; -else - stat = py.mhkit_python_utils.scipy_stats.convert_to_array(result, method, x); - ste = double(stat); + pd.type = 'Weibull'; + pd.scale = best_scale; + pd.shape = best_shape; + pd.pdf = weibull_pdf(x, best_scale, best_shape); + pd.cdf = weibull_cdf(x, best_scale, best_shape).^nst; % scaled for ST realization + pd.ppf = weibull_ppf(x, best_scale, best_shape); end -end \ No newline at end of file +function p = weibull_pdf(x, scale, shape) + p = (shape./scale) .* (x./scale).^(shape - 1) .* exp(-(x./scale).^shape); +end + +function c = weibull_cdf(x, scale, shape) + c = 1 - exp(-(x./scale).^shape); +end + +function x_ppf = weibull_ppf(p, scale, shape) + x_ppf = scale .* (-log(1 - p)).^(1/shape); +end diff --git a/mhkit/wave/resource/energy_period.m b/mhkit/wave/resource/energy_period.m index 85726fa22..df2613f8f 100644 --- a/mhkit/wave/resource/energy_period.m +++ b/mhkit/wave/resource/energy_period.m @@ -1,49 +1,74 @@ -function Te=energy_period(S,varargin) - +function Te = energy_period(S, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% +% Calculates wave energy period Te from wave spectra % % Parameters % ------------ -% S: Spectral Density (m^2/Hz) -% Pandas data frame -% To make a pandas data frame from user supplied frequency and spectra -% use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(frequency,spectra) -% -% OR -% -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) -% -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. -% -% S.frequency: frequency (Hz) -% -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% S: structure or numeric array +% If structure: +% S.spectrum - Spectral Density (m^2/Hz) +% S.frequency - Frequency (Hz) +% If numeric: +% S is spectral density array (vector or matrix) +% varargin{1} must contain frequency vector % +% frequency_bins: vector (optional) +% Frequency bin widths [Hz]. Required for unevenly spaced bins. % % Returns % --------- -% Te: float -% Wave energy Period (s) -% -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Te: double +% Wave energy period [s] +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% assign feq_bin -if nargin == 2 - freq_bins = py.numpy.array(varargin{1}); -elseif nargin == 1 - freq_bins = py.None; -else - ME = MException('MATLAB:energy_period','Incorrect number of input arguments'); - throw(ME); -end + % Extract spectrum and frequency + if isstruct(S) + spectrum = S.spectrum; + frequency = S.frequency; + elseif isnumeric(S) + if nargin < 2 + error('When S is numeric, frequency vector must be provided as second argument'); + end + spectrum = S; + frequency = varargin{1}; + varargin(1) = []; + else + error('Input S must be a struct or numeric array'); + end + + % Determine frequency bin widths + if ~isempty(varargin) + freq_bins = varargin{1}; + else + df = diff(frequency); + freq_bins = mean(df); + end -S_py = typecast_spectra_to_mhkit_python(S); + % Ensure proper shapes + frequency = frequency(:); + if isvector(spectrum) + spectrum = spectrum(:); + end -Te = py.mhkit.wave.resource.energy_period(S_py, pyargs('frequency_bins',freq_bins)); + % Check size consistency + if length(frequency) ~= size(spectrum,1) + error('Length of frequency vector must match number of rows in spectrum'); + end -Te = typecast_from_mhkit_python(Te).data; + % Calculate moments: m0 and m-1 + if isscalar(freq_bins) + m0 = sum(spectrum .* freq_bins, 1); + m_neg1 = sum((spectrum ./ frequency) .* freq_bins, 1); + else + freq_bins = freq_bins(:); + if length(freq_bins) ~= length(frequency) + error('Length of freq_bins must match frequency vector'); + end + m0 = sum(spectrum .* freq_bins, 1); + m_neg1 = sum((spectrum ./ frequency) .* freq_bins, 1); + end + + % Calculate energy period + Te = m_neg1 ./ m0; + +end diff --git a/mhkit/wave/resource/jonswap_spectrum.m b/mhkit/wave/resource/jonswap_spectrum.m index 7334a9b72..b7b80e241 100644 --- a/mhkit/wave/resource/jonswap_spectrum.m +++ b/mhkit/wave/resource/jonswap_spectrum.m @@ -1,55 +1,91 @@ -function S=jonswap_spectrum(frequency,Tp,Hs,varargin) - +function S = jonswap_spectrum(frequency, Tp, Hs, gamma) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates JONSWAP spectrum from Hasselmann et al (1973) +% Calculates JONSWAP spectrum based on IEC TS 62600-2 ED2 Annex C.2 (2019) % % Parameters % ------------ -% -% Frequency: float +% frequency: vector % Wave frequency (Hz) % % Tp: float % Peak Period (s) % % Hs: float -% Significant Wave Height (s) +% Significant Wave Height (m) % % gamma: float (optional) -% Gamma +% Peak enhancement factor % % Returns % --------- -% S: structure -% -% -% S.spectrum=Spectral Density (m^2/Hz) -% -% S.type=String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. -% -% S.frequency= frequency (Hz) -% -% S.Te -% -% S.Hm0 -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% S: structure with fields +% .spectrum = Spectral Density (m^2/Hz) +% .frequency = Frequency (Hz) +% .type = 'JONSWAP (Hs,Tp)' +% .Hm0 = Significant wave height (m) +% .Te = Energy period (s) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -if (isa(frequency,'py.numpy.ndarray') ~= 1) - frequency = py.numpy.array(frequency); +% Ensure column vector +frequency = frequency(:); +f = frequency; +df = diff(f); +uniform_spacing = all(abs(df - df(1)) < 1e-6); +if uniform_spacing + df_val = df(1); +else + df_val = [diff(f); diff(f(end-1:end))]; % Extend last bin end -if nargin == 3 - S_py=py.mhkit.wave.resource.jonswap_spectrum(frequency,Tp,Hs); -elseif nargin == 4 - S_py=py.mhkit.wave.resource.jonswap_spectrum(frequency, Tp, Hs, pyargs('gamma', varargin{1})); +% Constants +fp = 1 / Tp; +B_PM = (5 / 4) * (1 / Tp)^4; +A_PM = B_PM * (Hs / 2)^2; + +% Initialize spectral density +S_f = zeros(size(f)); +nonzero = f > 0; +S_f(nonzero) = A_PM .* f(nonzero).^(-5) .* exp(-B_PM .* f(nonzero).^(-4)); + +% Gamma computation +if nargin < 4 || isempty(gamma) + TpsqrtHs = Tp / sqrt(Hs); + if TpsqrtHs <= 3.6 + gamma = 5; + elseif TpsqrtHs > 5 + gamma = 1; + else + gamma = exp(5.75 - 1.15 * TpsqrtHs); + end end -S_py = typecast_from_mhkit_python(S_py); +% Spreading function G(f) +siga = 0.07; +sigb = 0.09; +Gf = zeros(size(f)); +lind = f <= fp; +hind = f > fp; +Gf(lind) = gamma .^ exp(-((f(lind) - fp).^2) ./ (2 * siga^2 * fp^2)); +Gf(hind) = gamma .^ exp(-((f(hind) - fp).^2) ./ (2 * sigb^2 * fp^2)); + +C = 1 - 0.287 * log(gamma); +Sf = C * S_f .* Gf; -S = struct(); +% Outputs +S.frequency = f; +S.spectrum = Sf; +S.type = sprintf('JONSWAP (%.2fm, %.2fs)', Hs, Tp); -S.frequency = S_py.index.data; -S.spectrum = S_py.data; -S.type = S_py.columns{1}; +% Compute Hm0 and Te +if uniform_spacing + m0 = sum(Sf) * df_val; + m1 = sum(Sf ./ f) * df_val; +else + m0 = sum(Sf .* df_val); + m1 = sum((Sf ./ f) .* df_val); +end + +S.Hm0 = 4 * sqrt(m0); +S.Te = m1 / m0; + +end diff --git a/mhkit/wave/resource/significant_wave_height.m b/mhkit/wave/resource/significant_wave_height.m index 715a8eff4..31e1caf3d 100644 --- a/mhkit/wave/resource/significant_wave_height.m +++ b/mhkit/wave/resource/significant_wave_height.m @@ -1,47 +1,72 @@ -function Hm0 = significant_wave_height(S,varargin) - +function Hm0 = significant_wave_height(S, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates wave height from spectra +% Calculates significant wave height Hm0 from spectra % % Parameters % ------------ -% S: Spectral Density (m^2/Hz) -% Pandas data frame -% To make a pandas data frame from user supplied frequency and spectra -% use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(frequency,spectra) -% -% OR -% -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) -% -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. -% -% S.frequency: frequency (Hz) +% S: structure or numeric array +% If structure: +% S.spectrum: Spectral Density (m^2/Hz) +% S.frequency: frequency (Hz) +% If numeric: +% S is assumed to be spectral density vector/matrix +% varargin{1} must contain frequency vector % % frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% Bin widths for frequency of S. Required for unevenly sized bins % % Returns % --------- % Hm0: double % Significant Wave Height (m) -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% assign freq_bin -if nargin == 2 - freq_bins = py.numpy.array(varargin{1}); -elseif nargin == 1 - freq_bins = py.None; -else - ME = MException('MATLAB:significant_wave_height','Incorrect number of input arguments'); - throw(ME); -end + % Extract spectrum and frequency + if isstruct(S) + spectrum = S.spectrum; + frequency = S.frequency; + elseif isnumeric(S) + if nargin < 2 + error('When S is numeric, frequency vector must be provided as second argument'); + end + spectrum = S; + frequency = varargin{1}; + varargin(1) = []; + else + error('Input S must be either a struct with fields .spectrum and .frequency, or a numeric array'); + end + + % Extract frequency bin widths + if ~isempty(varargin) + freq_bins = varargin{1}; + else + df = diff(frequency); + freq_bins = mean(df); % assume uniform bin width + end + + % Ensure column vector format + frequency = frequency(:); + if isvector(spectrum) + spectrum = spectrum(:); + end -S_py = typecast_spectra_to_mhkit_python(S); + % Check that frequency matches spectrum + if length(frequency) ~= size(spectrum,1) + error('Length of frequency vector must match number of rows in spectrum'); + end -Hm0 = py.mhkit.wave.resource.significant_wave_height(S_py, pyargs('frequency_bins',freq_bins)); + % Calculate zeroth moment m0 + if isscalar(freq_bins) + m0 = sum(spectrum .* freq_bins, 1); + else + freq_bins = freq_bins(:); + if length(freq_bins) ~= length(frequency) + error('Length of freq_bins must match frequency vector'); + end + m0 = sum(spectrum .* freq_bins, 1); + end -Hm0 = typecast_from_mhkit_python(Hm0).data; + % Calculate significant wave height Hm0 + Hm0 = 4 * sqrt(m0); + +end diff --git a/mhkit/wave/resource/surface_elevation.m b/mhkit/wave/resource/surface_elevation.m index 3fdea73ce..023836e09 100644 --- a/mhkit/wave/resource/surface_elevation.m +++ b/mhkit/wave/resource/surface_elevation.m @@ -1,106 +1,91 @@ -function wave_elevation=surface_elevation(S,time_index,options) +function wave_elevation = surface_elevation(S, time_index, options) +%SURFACE_ELEVATION Generate wave surface elevation from spectrum +% +% Inputs: +% S - struct with fields: +% .spectrum: (n_freq x 1) spectral density values (m^2/Hz) +% .frequency: (n_freq x 1) frequency vector (Hz) +% time_index - (n_time x 1) time vector (s) +% options - struct with optional fields: +% .seed: random seed (default = 123) +% .frequency_bins: bin widths for each frequency +% .phases: explicit phases (radians) +% .method: 'ifft' (default) or 'sum_of_sines' +% +% Output: +% wave_elevation - struct with fields: +% .elevation: (n_time x 1) wave elevation (m) +% .time: time vector (s) +% .type: description string -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates wave elevation time series from spectrum using a random phase -% -% Parameters -% ------------ -% S: Spectral Density (m^2/Hz) -% Pandas data frame -% To make a pandas data frame from user supplied frequency and spectra -% use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(frequency,spectra) -% -% OR -% -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) -% -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. -% -% S.frequency: frequency (Hz) -% -% time_index: array -% Time used to create the wave elevation time series [s] -% -% seed: Int (optional) -% random seed -% to call: wave_elevation(S,time_index,"seed",seed) -% -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins -% to call: wave_elevation(S,time_index,"frequency_bins",frequency_bins) -% -% phases: vector or matrix (optional) -% Explicit phases for frequency components (overrides seed) -% to call: wave_elevation(S,time_index,"phases",phases) -% -% method: str (optional) -% Method used to calculate the surface elevation. 'ifft' (Inverse Fast Fourier -% Transform) used by default if the given frequency_bins==None. 'sum_of_sines' -% explicitly sums each frequency component and used by default if -% frequency_bins are provided. The 'ifft' method is significantly faster. -% -% Returns -% --------- -% wave_elevation: structure -% -% wave_elevation.elevation: Wave surface elevation (m) -% -% wave_elevation.type: 'Time Series from Spectra' -% -% wave_elevation.time -% -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% arguments - S - time_index - options.seed {mustBeNumeric} = 123; - options.frequency_bins = py.None; - options.phases = py.None; - options.method = "ifft"; + S struct + time_index (:,1) double + options.seed (1,1) double = 123 + options.frequency_bins = [] + options.phases = [] + options.method char {mustBeMember(options.method, {'ifft','sum_of_sines'})} = 'ifft' end -S_py = typecast_spectra_to_mhkit_python(S); - -frequency = S.frequency ; - -if (isa(time_index,'py.numpy.ndarray') ~= 1) - - time_index = py.numpy.array(time_index); +% Extract frequency and spectrum +f = S.frequency(:); +Sf = S.spectrum(:); +Nf = numel(f); +Nt = numel(time_index); +% Handle frequency bins +if isempty(options.frequency_bins) + df = diff(f); + df_uniform = all(abs(df - df(1)) < 1e-8); + df_val = mean(df); + delta_f = df_val * ones(size(f)); +else + delta_f = options.frequency_bins(:); + if length(delta_f) ~= length(f) + error('frequency_bins must match the length of frequency vector.'); + end + df_val = delta_f(1); + df_uniform = all(abs(delta_f - df_val) < 1e-8); end -if (isa(options.frequency_bins,'py.NoneType')~=1) - if isnumeric(options.frequency_bins) - - options.frequency_bins = py.numpy.array(options.frequency_bins); - else - ME = MException('MATLAB:significant_wave_height','frequency_bins need to be of numeric type'); - throw(ME); +% Handle phases +if isempty(options.phases) + rng(options.seed); + phase = 2*pi*rand(Nf, 1); +else + phase = options.phases(:); + if length(phase) ~= Nf + error('phases must match the length of frequency vector.'); end end -if (isa(options.phases,'py.NoneType')~=1) - if isnumeric(options.phases) - options.phases = py.numpy.array(options.phases); - else - ME = MException('MATLAB:significant_wave_height','phases need to be of numeric type'); - throw(ME); +% Choose method +omega = 2*pi*f; +A = sqrt(2 * Sf .* delta_f); + +switch options.method + case 'ifft' + if f(1) ~= 0 || ~df_uniform + warning('Switching to sum_of_sines because ifft requires f(1)==0 and uniform spacing.'); + options.method = 'sum_of_sines'; + else + % Create complex spectrum + A_complex = A .* exp(1i * phase); + % Use irfft approximation: 0 to fN mapped to 0 to Nyquist + % Double spectrum (except DC and Nyquist if exists) + spectrum_full = [A_complex(1); A_complex(2:end-1)/2; conj(flip(A_complex(2:end-1)/2))]; + eta = real(ifft(spectrum_full, Nt)) * Nt; + wave_elevation.elevation = eta; end + case 'sum_of_sines' + B = omega .* time_index'; + B = B'; % (Nt x Nf) + C = cos(B + phase'); + eta = C * A; + wave_elevation.elevation = eta; end -if ~(strcmp(options.method, "ifft") || strcmp(options.method, "sum_of_sines")) - ME = MException('MATLAB:significant_wave_height','Invalid method. method should be either "ifft" or "sum_of_sines"'); - throw(ME); +wave_elevation.time = time_index; +wave_elevation.type = 'Time Series from Spectra'; end - -eta_py=py.mhkit.wave.resource.surface_elevation(S_py, time_index); - -eta = typecast_from_mhkit_python(eta_py); - -wave_elevation.type='Time Series from Spectra'; -wave_elevation.time=eta.index.data; -wave_elevation.elevation = eta.data;