Merge pull request #176 from simmsa/ndbc_examples_wave_module_only

Split PR #167 into wave module changes

Author: simmsa

changelog.md

@@ -1,5 +1,17 @@
 # Current development
 
+## PR 176 - Wave Module Native MATLAB Implementation
+
+- Authors: @MShabara, @simmsa
+- Convert wave module functions to native MATLAB code:
+  - wave/resource/standardize_wave_spectra_frequency.m:
+  - wave/resource/frequency_moment.m:
+  - wave/resource/energy_period.m:
+  - wave/resource/significant_wave_height.m:
+  - wave/resource/jonswap_spectrum.m:
+  - wave/resource/pierson_moskowitz_spectrum.m:
+  - wave/resource/surface_elevation.m:
+
 ## PR 153 - Mooring Module
 
 - Author: @hivanov-nrel

examples/SWAN_example.html

@@ -17,7 +17,7 @@
 % # An ASCII block file ('SWANOUTBlock.DAT')
 % # A binary block file ('SWANOUT.mat')
 
-swan_path = "./examples/data/wave/SWAN/";
+swan_path = "./data/wave/SWAN/";
 swan_table_file = append(swan_path,"SWANOUT.DAT");
 swan_block_file = append(swan_path,"SWANOUTBlock.DAT");
 swan_block_mat_file = append(swan_path,"SWANOUT.mat") ;

examples/SWAN_example.m

@@ -6,7 +6,7 @@
 % We can use MHKiT to load data downloaded from <https://www.ndbc.noaa.gov/
 % https://www.ndbc.noaa.gov>.
 
-relative_file_name = './examples/data/wave/data.txt';
+relative_file_name = './data/wave/data.txt';
 current_dir = fileparts(matlab.desktop.editor.getActiveFilename);
 full_file_name = fullfile(current_dir, relative_file_name);
 ndbc_data = read_NDBC_file(full_file_name);

examples/SWAN_example.mlx

@@ -1,49 +1,54 @@
-function Te=energy_period(S,varargin)
+function Te = energy_period(S, varargin)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
+% Calculate wave energy period (Te) seconds from power spectral density (PSD)
 %
 % 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 spectral density array (vector or 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
 % ---------
-%    Te: float
-%        Wave energy Period (s)
+%     Te: double
+%         Wave energy period [seconds] indexed by S.columns
 %
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-% 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);
+arguments
+    S
 end
 
-S_py = typecast_spectra_to_mhkit_python(S);
+arguments (Repeating)
+    varargin
+end
 
-Te = py.mhkit.wave.resource.energy_period(S_py, pyargs('frequency_bins',freq_bins));
+    % Calculate moments using frequency_moment function
+    if isstruct(S)
+        % Pass struct directly
+        m0 = frequency_moment(S, 0, varargin{:});
+        m_neg1 = frequency_moment(S, -1, varargin{:});
+    elseif isnumeric(S)
+        % Pass numeric spectrum with frequency vector
+        if nargin < 2
+            error('When S is numeric, frequency vector must be provided as second argument');
+        end
+        m0 = frequency_moment(S, 0, varargin{:});
+        m_neg1 = frequency_moment(S, -1, varargin{:});
+    else
+        error('Input S must be a struct or numeric array');
+    end
 
-Te = typecast_from_mhkit_python(Te).data;
+    % Calculate energy period
+    Te = m_neg1 ./ m0;
+
+end

examples/wave_example.html

@@ -1,48 +1,68 @@
-function m=frequency_moment(S,N,varargin)
+function m = frequency_moment(S, N, varargin)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%   Calculates the Nth frequency moment of the spectrum
 %
-% 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.
+% Calculates the Nth frequency moment of the spectrum
 %
-%           S.frequency: frequency (Hz)
+% Parameters
+% -----------
+%     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
 %
-%    N: int
-%       Moment (0 for 0th, 1 for 1st ....)
+%     N: int
+%         Moment (0 for 0th, 1 for 1st ....)
 %
-%    frequency_bins: vector (optional)
-%       Bin widths for frequency of S. Required for unevenly sized bins
+%     frequency_bins: vector (optional)
+%         Bin widths for frequency of S. Required for unevenly sized bins
 %
 % Returns
-% ---------
-%    m: double
-%
-%
+% -------
+%     m: double
+%         Nth Frequency Moment indexed by S.columns
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-S_py = typecast_spectra_to_mhkit_python(S);
+arguments
+    S
+    N {mustBeNumeric, mustBeFinite, mustBeInteger}
+end
 
-if nargin == 3
-    m = py.mhkit.wave.resource.frequency_moment(S_py, int32(N),pyargs('frequency_bins',py.numpy.array(varargin{1})));
-elseif nargin == 2
-    m = py.mhkit.wave.resource.frequency_moment(S_py, int32(N));
-else
-    ME = MException('MATLAB:frequency_moment','Incorrect number of arguments');
-        throw(ME);
+arguments (Repeating)
+    varargin
 end
 
-m = typecast_from_mhkit_python(m).data;
+    % Extract spectrum and frequency
+    if isstruct(S)
+        spectrum = S.spectrum;
+        frequency = S.frequency;
+    elseif isnumeric(S)
+        if nargin < 3
+            error('When S is numeric, frequency vector must be provided as third argument');
+        end
+        spectrum = S;
+        frequency = varargin{1};
+        varargin(1) = [];
+    else
+        error('Input S must be a struct or numeric array');
+    end
+
+    % Standardize frequency, spectrum, and frequency bins
+    if ~isempty(varargin)
+        [frequency, spectrum, freq_bins] = standardize_wave_spectra_frequency(frequency, spectrum, varargin{1});
+    else
+        [frequency, spectrum, freq_bins] = standardize_wave_spectra_frequency(frequency, spectrum);
+    end
+
+    if isscalar(freq_bins)
+        freq_bins = freq_bins(:);
+    end
+
+    % Calculate Nth moment: m_N = sum(f^N * S * df)
+    m = sum((frequency.^N) .* spectrum .* freq_bins, 1);
+
+end

examples/wave_example.m

@@ -1,55 +1,90 @@
-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 from IEC TS 62600-2 ED2 Annex C.2 (2019)
 %
 % Parameters
 % ------------
-%
-%     Frequency: float
-%         Wave frequency (Hz)
+%     frequency: vector
+%         Frequency [Hz]
 %
 %     Tp: float
-%         Peak Period (s)
+%         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: structure with fields
+%         .spectrum  = Spectral density [m^2/Hz]
+%         .frequency = Frequency [Hz]
+%         .type      = 'JONSWAP (Hs,Tp)'
 %
-%         S.frequency= frequency (Hz)
-%
-%         S.Te
-%
-%         S.Hm0
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-if (isa(frequency,'py.numpy.ndarray') ~= 1)
-    frequency = py.numpy.array(frequency);
+arguments
+    frequency {mustBeNumeric, mustBeFinite}
+    Tp {mustBeNumeric, mustBeFinite, mustBePositive}
+    Hs {mustBeNumeric, mustBeFinite, mustBePositive}
+    gamma {mustBeNumeric, mustBeFinite, mustBePositive} = []
 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}));
+% Ensure column vector and sort (match MHKiT-Python behavior)
+frequency = frequency(:);
+f = sort(frequency);
+
+% Constants (same as MHKiT-Python)
+B_PM = (5 / 4) * (1 / Tp)^4;
+A_PM = B_PM * (Hs / 2)^2;
+
+% Avoid divide by zero if the 0 frequency is provided
+% The zero frequency should always have 0 amplitude, otherwise
+% we end up with a mean offset when computing the surface elevation.
+S_f = zeros(size(f));
+if f(1) == 0.0
+    inds = 2:length(f);  % Skip first element like MHKiT-Python
+else
+    inds = 1:length(f);  % Use all elements
+end
+
+S_f(inds) = A_PM * f(inds).^(-5) .* exp(-B_PM * f(inds).^(-4));
+
+% Gamma computation (match MHKiT-Python logic)
+if 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);
+% Cutoff frequencies for gamma function (match MHKiT-Python)
+siga = 0.07;
+sigb = 0.09;
+
+% Peak frequency
+fp = 1 / Tp;
+lind = f <= fp;
+hind = f > fp;
+Gf = zeros(size(f));
+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));
 
-S = struct();
+C = 1 - 0.287 * log(gamma);
+Sf = C * S_f .* Gf;
 
-S.frequency = S_py.index.data;
-S.spectrum = S_py.data;
-S.type = S_py.columns{1};
+% Output structure (match MHKiT-Python naming style)
+name = sprintf('JONSWAP (%.1fm,%.1fs)', Hs, Tp);
+S.frequency = f;
+S.spectrum = Sf;
+S.type = name;
+
+end