diff --git a/examples/cdip_example.html b/examples/cdip_example.html index 26853489a..bce98ba94 100644 --- a/examples/cdip_example.html +++ b/examples/cdip_example.html @@ -1,12 +1,12 @@ -MHKiT CDIP IO

MHKiT CDIP IO

MHKiT includes functions to pull data directly from the Coastal Data Information Program (CDIP), an extensive network for monitoring waves and beaches along the coastlines of the United States.

Quick Start

The cdip_request_parse_workflow function is your one-stop shop for pulling and parsing data. If a station number is passed the function will pull the buoy data from the CDIP servers and then proceed to process the data for the user. cdip_request_parse_workflow allows the user to slice data between a start and/or end date. Alternatively, the user may specify a single year or multiple years of interest. To reduce processing time cdip_request_parse_workflow will allow the user to specify parameters of interest and only return data for the specified parameters. By default, the cdip_request_parse_workflow will pull historic data, but the user may also specify real-time data if interested. Lastly, due to the long processing time, two-dimensional (2D) data is not returned by default. If all 2D data is desired the boolean all_2D_variables may be specified as true. It is recommended, however, that if 2D data is needed that the user passes that parameter in the parameters field. For a full list of 1D, 2D, and metadata variables see the CDIP reference here (see BUOY VARIABLES).

1. Request Data from CDIP

In the above example, cdip_request_parse_workflow is used to both request and parse the CDIP data using multiple MHKiT functions in a prespecfied workflow for the user.
MHKiT can be used to request historical or realtime data from the CDIP buoys. A station table and map can be found here to determine a buoy of interest <http://cdip.ucsd.edu/m/stn_table/>. To get started we will call historic data from station number 100. The function will return a structure with all historic 1D data and metadata. The 'data' and 'metadata' each hold different types of data based based on a prefix used on the variable returned by CDIP. Example prefixes are 'wave', 'sst', and 'gps', which are described in further detail below.
station_number = '100';
data_type = 'historic';
buoy_data = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'data_type', data_type);
disp("Returned data: " + string(strjoin(fieldnames(buoy_data))))
Returned data: data metadata

2. View the Returned Data

2.a The 'data' field As can be seen above the function has returned a structure with two field names:
1. data - field with subfields referring to the data type. Possible subfields include:
2. metadata - any other data not of length time or a 2D variable. This has the same fields as the 'data' field above.
By selecting the subfield data.wave we can see the associated array for the wave data. The array is summarized below where one can see the values of each metric. By selecting the subfield data.waveTime we can see the associated timestamps.
buoy_data.data.wave
ans = struct with fields:
waveDp: [379560×1 double] - waveFlagPrimary: [379560×1 double] - waveFlagSecondary: [379560×1 double] - waveHs: [379560×1 double] - wavePeakPSD: [379560×1 double] - waveSourceIndex: [379560×1 double] - waveTa: [379560×1 double] - waveTime: [379560×1 datetime] - waveTp: [379560×1 double] - waveTz: [379560×1 double] -
2.b The metadata field The metadata field stores anything not of length time for each of the previous fields. This field also has the meta key for any data returned that started with the meta prefix.
fieldnames(buoy_data.metadata)
ans = 5×1 cell
'meta'
'other2D'
'wave'
'wave2D'
'name'
 
buoy_data.metadata.meta
ans = struct with fields:
metaDeclination: [18×1 single] - metaDeployLatitude: [18×1 single] - metaDeployLongitude: [18×1 single] +.S11 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

MHKiT CDIP IO

MHKiT includes functions to pull data directly from the Coastal Data Information Program (CDIP), an extensive network for monitoring waves and beaches along the coastlines of the United States.

Quick Start

The cdip_request_parse_workflow function is your one-stop shop for pulling and parsing data. If a station number is passed the function will pull the buoy data from the CDIP servers and then proceed to process the data for the user. cdip_request_parse_workflow allows the user to slice data between a start and/or end date. Alternatively, the user may specify a single year or multiple years of interest. To reduce processing time cdip_request_parse_workflow will allow the user to specify parameters of interest and only return data for the specified parameters. By default, the cdip_request_parse_workflow will pull historic data, but the user may also specify real-time data if interested. Lastly, due to the long processing time, two-dimensional (2D) data is not returned by default. If all 2D data is desired the boolean all_2D_variables may be specified as true. It is recommended, however, that if 2D data is needed that the user passes that parameter in the parameters field. For a full list of 1D, 2D, and metadata variables see the CDIP reference here (see BUOY VARIABLES).

1. Request Data from CDIP

In this example, cdip_request_parse_workflow is used to both request and parse the CDIP data using multiple MHKiT functions in a prespecified workflow for the user.
MHKiT can be used to request historical or realtime data from the CDIP buoys. A station table and map can be used to determine a buoy of interest. To get started we will call historic data from station number 100. The function will return a structure with all historic 1D data and metadata. The 'data' and 'metadata' each hold different types of data based on a prefix used on the variable returned by CDIP. Example prefixes are 'wave', 'sst', and 'gps', which are described in further detail below.
station_number = '100';
data_type = 'historic';
buoy_data = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'data_type', data_type);
disp("Returned data: " + string(strjoin(fieldnames(buoy_data))))
Returned data: data metadata

2. View the Returned Data

2.a The 'data' field As can be seen above the function has returned a structure with two field names:
1. data - field with subfields referring to the data type. Possible subfields include:
  • wave - includes processed variables such as Hs, and Tp
  • sst - timeseries of sea surface temperature
  • gps - timeseries of buoy latitude and longitude
  • dwr - directional waverider contains information about accelerometer and buoy battery level
  • wave2D - arrays which are of height time and have columns of wave frequency.
2. metadata - any other data not of length time or a 2D variable. This has the same fields as the 'data' field above.
  • meta - processed variables which started with the meta prefix and are not part of the variable prefixes listed below
  • wave - includes processed variables such as Hs, and Tp
  • sst - timeseries of sea surface temperature
  • gps - timeseries of buoy latitude and longitude
  • dwr - directional Waverider contains information about accelerometer and buoy battery level
  • wave2D - arrays which are of height time and have columns of wave frequency.
By selecting the subfield data.wave we can see the associated array for the wave data. The array is summarized below where one can see the values of each metric. By selecting the subfield data.waveTime we can see the associated timestamps.
buoy_data.data.wave
ans = struct with fields:
waveDp: [436469×1 double] + waveFlagPrimary: [436469×1 double] + waveFlagSecondary: [436469×1 double] + waveHs: [436469×1 double] + wavePeakPSD: [436469×1 double] + waveSourceIndex: [436469×1 double] + waveTa: [436469×1 double] + waveTime: [436469×1 datetime] + waveTp: [436469×1 double] + waveTz: [436469×1 double] +
2.b The metadata field The metadata field stores anything not of length time for each of the previous fields. This field also has the meta key for any data returned that started with the meta prefix.
fieldnames(buoy_data.metadata)
ans = 5×1 cell
'meta'
'other2D'
'wave'
'wave2D'
'name'
 
buoy_data.metadata.meta
ans = struct with fields:
metaDeclination: [19×1 single] + metaDeployLatitude: [19×1 single] + metaDeployLongitude: [19×1 single] metaGridMapping: [64×1 char] metaInstrumentation: [64×1 char] metaPlatform: [64×1 char] - metaStationLatitude: 32.9314 - metaStationLongitude: -117.3918 + metaStationLatitude: 32.9312 + metaStationLongitude: -117.3919 metaStationName: [64×1 char] - metaWaterDepth: [18×1 single] -

3. Plotting

The wave graphics package has two CDIP stylized plots, boxplot and compendium. The CDIP scripts these are based on are hyperlinked.
3.a The Box Plot The CDIP stylized box plot is designed to group a passed parameter by months and display descriptive statistics. Below the plot is a general box plot legend showing the outliers as red plus signs, and the 25-percentile and 75-percentile mark the edges of the box. The median is shown as a red line in the box and the mean is shown as a green line with the mean annotated in green above the line. The x-axis has the month name at the bottom and the number of data points in each month at the top.
buoy_data_2011 = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'years', 2011, ...
'data_type', data_type);
Hs_2011_data = buoy_data_2011.data.wave.waveHs;
time_2011 = buoy_data_2011.data.wave.waveTime;
buoy_name = buoy_data_2011.metadata.name;
plot_boxplot(Hs_2011_data, time_2011, 'buoy_title', buoy_name);
3.b CDIP Compendium Plot The compendium plot is a time series summary plot of three variables (Hs, Tp, and Dp). To return the plot simply pass these variable from the data.wave structure.
wave_data_May_2011 = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'start_date', '2011-05-01', ...
'end_date', '2011-05-31', ...
'data_type', data_type);
Hs = wave_data_May_2011.data.wave.waveHs;
Tp = wave_data_May_2011.data.wave.waveTp;
Dp = wave_data_May_2011.data.wave.waveDp;
time = wave_data_May_2011.data.wave.waveTime;
plot_compendium(Hs, Tp, Dp, time, "buoy_title", buoy_name);

4. Advanced Queries

In this section we will show more on how to slice on time ranges of interest. 4.a start_date, end_date and parameters The function cdip_request_parse_workflow can iterate over and extract variables from CDIP buoy netCDF data. It takes inputs of start_date, end_date, parameters, and all_2D_variables. The start and end date parameters allow us to specify string dates of interest in which to include data after the start date and before the end date. These are specified as a string with precision down to the day in 'YYYY-MM-DD' format. Specifying particular parameters of interest allows for faster processing times and data handling. For a complete list of parameters see Buoy Variables here. Lastly, we can specify to return all 2D variables, although this is not recommended unless a user truly plans to use all the data because of the significant processing time.
For an example we will create a box plot of Hs for the year 2011.
buoy_data = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'start_date', '2011-01-01', ...
'end_date', '2011-12-31', ...
'parameters', {'waveHs'}, ...
'data_type', data_type);
buoy_name = buoy_data.metadata.name;
plot_boxplot( ...
buoy_data.data.wave.waveHs, ...
buoy_data.data.wave.waveTime, ...
"buoy_title", buoy_name);
4.b years and 2D data The function cdip_request_parse_workflow is especially useful for processing 2D data, which is only processed when specifically requested due to the amount of time it takes. Next, cdip_request_parse_workflow will slice on time by years, start_date, or end date. Years can be a single integer or an array of integers and is not required to be consecutive. If specified, the start date will remove any data prior to the specified string (e.g. '2011-01-01') and end_date will remove any data after the specified date. start_date and end_date may be used together, separately or not at all. Years works independently of start and end date. Next, the data_type defaults to historic but specifying this as realtime will return realtime data from the buoy. Lastly, there is the boolean all_2D_variables. If set to true, the function will return all of the wave 2D variables. It is not recommended to do this due to the computational expense, instead it is recommended to specify 2D quantities of interest using the parameters parameter.
For an example we will create a compendium of Hs, Tp and Dp for the year 2011. In this case we can use the years parameter instead of start and end dates.
buoy_data = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'years', 2011, ...
'parameters', {'waveHs', 'waveTp', 'waveDp'});
Hs = buoy_data.data.wave.waveHs;
Tp = buoy_data.data.wave.waveTp;
Dp = buoy_data.data.wave.waveDp;
time = buoy_data.data.wave.waveTime;
buoy_name = buoy_data.metadata.name;
plot_compendium(Hs, Tp, Dp, time, "buoy_title", buoy_name);
+ metaWaterDepth: [19×1 single] +

3. Plotting

The wave graphics package has two CDIP stylized plots, boxplot and compendium. The CDIP scripts these are based on are hyperlinked.
3.a The Box Plot The CDIP stylized box plot is designed to group a passed parameter by months and display descriptive statistics. Below the plot is a general box plot legend showing the outliers as red plus signs, and the 25-percentile and 75-percentile mark the edges of the box. The median is shown as a red line in the box and the mean is shown as a green line with the mean annotated in green above the line. The x-axis has the month name at the bottom and the number of data points in each month at the top.
buoy_data_2011 = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'years', 2011, ...
'data_type', data_type);
Hs_2011_data = buoy_data_2011.data.wave.waveHs;
time_2011 = buoy_data_2011.data.wave.waveTime;
buoy_name = buoy_data_2011.metadata.name;
plot_boxplot(Hs_2011_data, time_2011, 'buoy_title', buoy_name);
3.b CDIP Compendium Plot The compendium plot is a time series summary plot of three variables (Hs, Tp, and Dp). To return the plot simply pass these variable from the data.wave structure.
wave_data_May_2011 = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'start_date', '2011-05-01', ...
'end_date', '2011-05-31', ...
'data_type', data_type);
Hs = wave_data_May_2011.data.wave.waveHs;
Tp = wave_data_May_2011.data.wave.waveTp;
Dp = wave_data_May_2011.data.wave.waveDp;
time = wave_data_May_2011.data.wave.waveTime;
plot_compendium(Hs, Tp, Dp, time, "buoy_title", buoy_name);

4. Advanced Queries

In this section we will show more on how to slice on time ranges of interest. 4.a start_date, end_date and parameters The function cdip_request_parse_workflow can iterate over and extract variables from CDIP buoy netCDF data. It takes inputs of start_date, end_date, parameters, and all_2D_variables. The start and end date parameters allow us to specify string dates of interest in which to include data after the start date and before the end date. These are specified as a string with precision down to the day in 'YYYY-MM-DD' format. Specifying particular parameters of interest allows for faster processing times and data handling. For a complete list of parameters see Buoy Variables here. Lastly, we can specify to return all 2D variables, although this is not recommended unless a user truly plans to use all the data because of the significant processing time.
For an example we will create a box plot of Hs for the year 2011.
buoy_data = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'start_date', '2011-01-01', ...
'end_date', '2011-12-31', ...
'parameters', {'waveHs'}, ...
'data_type', data_type);
buoy_name = buoy_data.metadata.name;
plot_boxplot( ...
buoy_data.data.wave.waveHs, ...
buoy_data.data.wave.waveTime, ...
"buoy_title", buoy_name);
4.b years and 2D data The function cdip_request_parse_workflow is especially useful for processing 2D data, which is only processed when specifically requested due to the amount of time it takes. Next, cdip_request_parse_workflow will slice on time by years, start_date, or end date. Years can be a single integer or an array of integers and is not required to be consecutive. If specified, the start date will remove any data prior to the specified string (e.g. '2011-01-01') and end_date will remove any data after the specified date. start_date and end_date may be used together, separately or not at all. Years works independently of start and end date. Next, the data_type defaults to historic but specifying this as realtime will return realtime data from the buoy. Lastly, there is the boolean all_2D_variables. If set to true, the function will return all of the wave 2D variables. It is not recommended to do this due to the computational expense, instead it is recommended to specify 2D quantities of interest using the parameters parameter.
For an example we will create a compendium of Hs, Tp and Dp for the year 2011. In this case we can use the years parameter instead of start and end dates.
buoy_data = cdip_request_parse_workflow( ...
'station_number', station_number, ...
'years', 2011, ...
'parameters', {'waveHs', 'waveTp', 'waveDp'});
Hs = buoy_data.data.wave.waveHs;
Tp = buoy_data.data.wave.waveTp;
Dp = buoy_data.data.wave.waveDp;
time = buoy_data.data.wave.waveTime;
buoy_name = buoy_data.metadata.name;
plot_compendium(Hs, Tp, Dp, time, "buoy_title", buoy_name);

\ No newline at end of file diff --git a/examples/wave_example.m b/examples/wave_example.m index 3921da543..2bde1c373 100644 --- a/examples/wave_example.m +++ b/examples/wave_example.m @@ -1,108 +1,193 @@ %% Example: MHKiT-MATLAB Wave Module % The following example runs an application of the to 1) generate a capture length matrix, 2) calculate MAEP, +% MHKiT wave module> to 1) generate a capture width matrix, 2) calculate MAEP, % and 3) plot the scatter diagrams. -%% Load NDBC Data -% We can use MHKiT to load data downloaded from . - -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); -disp(ndbc_data) - +%% Load CDIP Wave Measurement Data PacWave North, CDIP 277 +% This example uses one year (2025) of wave spectra from the PacWave +% North buoy off Newport, Oregon, published by CDIP as +% . The spectra are +% requested from the CDIP THREDDS server with |cdip_request_parse_workflow|. +% +% A saved copy of that request is shipped with the examples so this script +% runs without network access. A working request is shown commented out +% below. + +station_number = '277'; + +% cdip = cdip_request_parse_workflow('station_number', station_number, ... +% 'data_type', 'historic', ... +% 'start_date', '2025-01-01', 'end_date', '2025-12-31', ... +% 'parameters', {'waveEnergyDensity', 'metaWaterDepth'}); +% % waveEnergyDensity is stored [time x frequency]; MHKiT uses frequency +% % down the rows and one spectrum per column +% S.frequency = double(cdip.metadata.wave.waveFrequency(:)); +% S.spectrum = double(cdip.data.wave2D.waveEnergyDensity'); +% S.time = cdip.data.wave.waveTime(:); +% S.time.TimeZone = 'UTC'; +% % Water depth at the buoy [m], a scalar metadata variable in the NetCDF +% h = double(cdip.metadata.meta.metaWaterDepth); + +% Saved copy of the request above +saved = load('./data/wave/cdip_277_2025.mat'); +S.frequency = saved.cdip_277.frequency; +S.spectrum = double(saved.cdip_277.spectrum); % stored as single, as CDIP publishes it +S.time = saved.cdip_277.time; +h = saved.cdip_277.water_depth; % [m], from metaWaterDepth +fprintf('Using saved CDIP 277 data downloaded %s UTC\n', saved.cdip_277.downloaded); +disp(S) +time = S.time; + +% Shared plot title suffix +title_suffix = sprintf('CDIP %s, %s', station_number, datetime(time(1), 'Format', 'yyyy')); %% Compute Wave Metrics % We will now use MHKiT to compute the significant wave height, energy period, -% and energy flux. - -% Compute the energy periods from the NDBC spectra data -Te = energy_period(ndbc_data); -disp(Te) -% Compute the significant wave height from the NDBC spectra data -Hm0 = significant_wave_height(ndbc_data); -disp(Hm0) -% Set water depth to 60 m -h = 60; - -% Compute the energy flux from the NDBC spectra data and water depth -J = energy_flux(ndbc_data,h); -disp(J) -%% Generate Random Power Data -% For demonstration purposes, this example uses synthetic power data generated -% from statistical distributions. In a real application, the user would provide -% power values from a WEC. - -% generating 1,000,000 random power values -Power = randi([40,200], 743,1); -disp(Power) -%% Capture Length Matrices -% The following operations create capture length matrices, as specified by the -% IEC/TS 62600-100. But first, we need to calculate capture length and define -% bin centers. The mean capture length matrix is printed below. Keep in mind that -% this data has been artificially generated, so it may not be representative of -% what a real-world scatter diagram would look like. - -% calculating capture length with power and wave flux in vectors -L = capture_length(Power,J) - -% Need to set our Hm0 and Te bins for the capture length matrix -Hm0_bins = -0.5:0.5:max(fix(Hm0))+0.5; % Input is min, max, and n indices for vector -Hm0_bins = Hm0_bins+0.25 ; -Te_bins = 0:1:max(fix(Te)); -Te_bins = Te_bins+0.5; - -% Calculate the necessary capture length matrices for each statistic based +% and energy flux from each wave spectrum. + +Hm0 = significant_wave_height(S); % [m] +Te = energy_period(S); % [s] +J = energy_flux(S, h); % [W/m], uses the water depth + +figure('Position', [100 100 900 700]); +tiledlayout(3, 1, 'TileSpacing', 'compact'); +sgtitle(sprintf('Wave Metrics | %s', title_suffix)) +nexttile; plot(time, Hm0); ylabel('H_{m0} [m]'); title('Significant Wave Height'); grid on +nexttile; plot(time, Te); ylabel('T_e [sec]'); title('Energy Period'); grid on +nexttile; plot(time, J / 1000); ylabel('J [kW/m]'); title('Wave Energy Flux'); grid on; xlabel('Time [UTC]') +%% Joint Probability Distribution +% The joint probability distribution (JPD) is the fraction of records that +% fall in each Hm0 and Te bin, and it depends only on the wave resource. +% |plot_wave_joint_probability_distribution| bins the sea states on the +% IEC/TS 62600-100 grid (0.5 m Hm0 and 1 s Te bins with edges on 0, 0.5, 1, +% ... m and 0, 1, 2, ... s), plots the JPD as a percentage with the empty bins +% left blank, and returns the JPD matrix struct. + +figure('Position', [100 100 1000 800]); +[~, jpd] = plot_wave_joint_probability_distribution(Hm0, Te, ... + 'title', sprintf('Joint Probability Distribution | %s', title_suffix)); + +% The same bins are used for the capture width and power matrices below +Hm0_bins = jpd.y_bins; +Te_bins = jpd.x_bins; +%% Model Device Power Using DOE Reference Models +% In a real application, the user would provide power values measured from a +% WEC. This example instead looks up each sea state in the power matrix of a +% device from the , a U.S. Department of Energy effort that developed +% open-source marine energy point designs to benchmark technology performance +% and costs. MHKiT ships the matrices for Reference Model 3 (wave point +% absorber), 5 (oscillating surge flap), and 6 (oscillating water column) from +% the ; change |device| to switch. The devices are +% described in Neary et al. (2014), , +% SAND2014-9040, Sandia National Laboratories. + +device = "RM3"; % "RM3", "RM5", or "RM6" +reference_model = reference_model_power_matrix(device); +disp(reference_model.description) +disp(reference_model.data_source) + +reference_model_kW = reference_model; +reference_model_kW.values = reference_model.values / 1000; +figure('Position', [100 100 1000 800]); +ax_device = mhkit_plot_matrix(reference_model_kW, ... + 'xlabel', 'Energy Period, T_e [sec]', 'ylabel', 'Significant Wave Height, H_{m0} [m]', ... + 'zlabel', 'Power [kW]', 'trim_to_data', true, 'value_format', '%.0f', 'font_size', 8); +title(ax_device, sprintf('%s Power Matrix | Neary et al. (2014)', device)); + +% Interpolate between bin centers; sea states outside the matrix produce no power +Power = interp2(reference_model.x_bins, reference_model.y_bins, reference_model.values, ... + Te, Hm0, 'linear', 0); % [W] + +figure('Position', [100 100 900 350]); +plot(time, Power / 1000); ylabel('Modeled Power [kW]'); xlabel('Time [UTC]'); grid on +title(sprintf('Modeled %s Power | %s', device, title_suffix)) +%% Capture Width +% The following operations create capture width matrices, as specified by the +% IEC/TS 62600-100, on the same Hm0 and Te bins as the JPD. But first, we need +% to calculate capture width. Keep in mind that the power comes from a modeled +% power matrix, not measurements, so the capture width varies smoothly with +% the sea state. + +% calculating capture width with power and wave flux in vectors +CW = capture_width(Power, J); + +% Calculate the necessary capture width matrices for each statistic based % on IEC/TS 62600-100 -clmat.mean = capture_length_matrix(Hm0,Te,L,"mean",Hm0_bins,Te_bins); -clmat.std = capture_length_matrix(Hm0,Te,L,"std",Hm0_bins,Te_bins); -clmat.count = capture_length_matrix(Hm0,Te,L,"count",Hm0_bins,Te_bins); -clmat.min = capture_length_matrix(Hm0,Te,L,"min",Hm0_bins,Te_bins); -clmat.max = capture_length_matrix(Hm0,Te,L,"max",Hm0_bins,Te_bins); - -% Calculate the frequency matrix for convenience -clmat.freq = capture_length_matrix(Hm0,Te,L,"frequency",Hm0_bins,Te_bins); +cwmat.mean = capture_width_matrix(Hm0, Te, CW, "mean", Hm0_bins, Te_bins); +cwmat.std = capture_width_matrix(Hm0, Te, CW, "std", Hm0_bins, Te_bins); +cwmat.count = capture_width_matrix(Hm0, Te, CW, "count", Hm0_bins, Te_bins); +cwmat.min = capture_width_matrix(Hm0, Te, CW, "min", Hm0_bins, Te_bins); +cwmat.max = capture_width_matrix(Hm0, Te, CW, "max", Hm0_bins, Te_bins); %% -% Let's see what the data in the mean matrix looks like. Keep in mind that this -% data has been artificially generated, so it may not be representative of what -% a real-world scatter diagram would look like. - -disp(clmat.mean.values) -%% Power Matrices +% The graphics function |plot_matrix| visualizes a matrix struct. Each bin is +% drawn between its edges with the tick marks at the bin edges, and +% |trim_to_data| limits the axes to the bins that hold data plus one empty +% bin around them. + +figure('Position', [100 100 1000 800]); +ax_cw = plot_matrix(cwmat.mean, "Capture Width", "zlabel", "Capture Width [m]", ... + "trim_to_data", true, "font_size", 8); +title(ax_cw, sprintf('Modeled %s Mean Capture Width Matrix | %s', device, title_suffix)); +%% Power Matrix % As specified in IEC/TS 62600-100, the power matrix is generated from the capture -% length matrix and wave energy flux matrix, as shown below +% width matrix and wave energy flux matrix, as shown below % Create wave energy flux matrix using mean -jmat = wave_energy_flux_matrix(Hm0,Te,J,"mean",Hm0_bins,Te_bins); +jmat = wave_energy_flux_matrix(Hm0, Te, J, "mean", Hm0_bins, Te_bins); % Create power matrix using mean -avg_power_mat = power_matrix(clmat.mean, jmat); +avg_power_mat = power_matrix(cwmat.mean, jmat); % Create power matrix using standard deviation -std_power_mat = power_matrix(clmat.std, jmat); +std_power_mat = power_matrix(cwmat.std, jmat); %% -% The |capture_length_matrix| function can also be used as an arbitrary scatter +% The |capture_width_matrix| function can also be used as an arbitrary scatter % plot generator. To do this, simply pass a different array in the place of capture -% length (L). For example, while not specified by the IEC standards, if the user +% width (CW). For example, while not specified by the IEC standards, if the user % doesn't have the omnidirectional wave flux, the average power matrix could hypothetically % be generated in the following manner: -avgpowmat_not_standard = capture_length_matrix(Hm0,Te,Power,'mean',Hm0_bins,Te_bins); -%% MAEP +avgpowmat_not_standard = capture_width_matrix(Hm0, Te, Power, 'mean', Hm0_bins, Te_bins); +%% +% The underlying |mhkit_plot_matrix| function accepts any matrix struct and +% custom axis labels. Here the mean power matrix is plotted in kW. + +avg_power_mat_kW = avg_power_mat; +avg_power_mat_kW.values = avg_power_mat.values / 1000; +figure('Position', [100 100 1000 800]); +ax_power = mhkit_plot_matrix(avg_power_mat_kW, 'xlabel', 'Energy Period, T_e [sec]', 'ylabel', 'Significant Wave Height, H_{m0} [m]', ... + 'zlabel', 'Mean Power [kW]', 'trim_to_data', true, 'value_format', '%.1f kW', 'font_size', 8); +title(ax_power, sprintf('Modeled %s Mean Power Matrix | %s', device, title_suffix)); +%% Energy Production +% The modeled power series gives the energy production of the device at this +% site, first month by month and then as the mean annual energy production +% (MAEP) defined by IEC/TS 62600-100. +%% Monthly Energy Production +% Each record is a 30 minute average, so the energy per record is the power +% times that duration. Summing the records by month shows how the energy +% production is distributed through the year. + +record_hours = hours(median(diff(time))); +energy_MWh = Power * record_hours / 1e6; % [MWh] per record +monthly_MWh = accumarray(month(time), energy_MWh, [12 1]); + +figure('Position', [100 100 900 400]); +bar(1:12, monthly_MWh); +xticks(1:12); xticklabels(month(datetime(2025, 1:12, 1), 'shortname')); +ylabel('Energy [MWh]'); grid on +yline(mean(monthly_MWh), '--', 'Monthly mean', 'LabelHorizontalAlignment', 'left'); +title(sprintf('Modeled %s Monthly Energy Production | %s', device, title_suffix)) +%% Mean Annual Energy Production % There are two ways to calculate mean annual energy production (MAEP). One -% is from capture length and wave energy flux matrices, the other is from time +% is from capture width and wave energy flux matrices, the other is from time % series data, as shown below. -% -% % Calculate maep from timeseries -maep_timeseries = mean_annual_energy_production_timeseries(L,J) +maep_timeseries = mean_annual_energy_production_timeseries(CW, J); % [W*h] % Calculate maep from matrix -maep_matrix = mean_annual_energy_production_matrix(clmat.mean, jmat , clmat.freq) -%% Graphics -% The graphics function |plot_matrix| can be used to visualize results. It is -% important to note that the plotting function assumes the step size between bins -% to be linear. - -% Plot the capture length matrix -p1 = plot_matrix(clmat.mean,"Capture Length", "annotate", false); +maep_matrix = mean_annual_energy_production_matrix(cwmat.mean, jmat, jpd); % [W*h] + +% MHKiT returns energy in W*h; report in MWh +fprintf('Modeled %s MAEP from time series: %.0f MWh\n', device, maep_timeseries / 1e6); +fprintf('Modeled %s MAEP from matrices: %.0f MWh\n', device, maep_matrix / 1e6); diff --git a/examples/wave_example.mlx b/examples/wave_example.mlx index 1020bbc25..d73197373 100644 Binary files a/examples/wave_example.mlx and b/examples/wave_example.mlx differ diff --git a/examples/wecsim_power_performance_example.html b/examples/wecsim_power_performance_example.html index 1ef065882..169374bc6 100644 --- a/examples/wecsim_power_performance_example.html +++ b/examples/wecsim_power_performance_example.html @@ -1,11 +1,11 @@ -MHKiT WEC-Sim Power Performance Example

MHKiT WEC-Sim Power Performance Example

This notebook demonstrates using the MHKiT wave module and WEC-Sim together to perform a resource characterization study in MHKiT, simulate representative cases with WEC-Sim, and visualize the results in MHKiT to estimate MAEP (Mean Annual Energy Production).
1. Characterize the available resource at a location
- See the PacWave example notebook
2. Write a WEC-Sim batch file for the given clusters
3. Simulate the device in WEC-Sim.
- Ensure that the spectra used in WEC-Sim is identical to the one used in MHKiT.
4. Load WEC-Sim batch results
5. Assess results and visualize quantities of interest
This example uses WEC-Sim to simulate the [Oscillating Surge Wave Energy Converter (OSWEC)](https://wec-sim.github.io/WEC-Sim/main/user/tutorials.html#oscillating-surge-wec-oswec), a flap-type device.

1. Characterize the available resource at a location

This example will use an abbreviated version of `PacWave_resource_characterization_example.ipynb`.
For full details on downloading, calculating, and visualizing the k-means clusters representation of the site's wave resouce, see that example.
We select the N=32 cluster as it's total energy flux is closet to the total energy flux of the site considering all wave conditions. We will load the PacWave example output, which can be easily saved after running the example with the command `results[32].to_csv("pacwave_cluster_32.csv")`. We will start this example by reading in that csv output and formatting it for WEC-Sim.
% Relative location and filename of simulated WEC-Sim data (run with mooring)
filename = './data/wave/pacwave_cluster_32.csv'
filename = './data/wave/pacwave_cluster_32.csv'
results = readtable(filename)
results = 32×6 table
 Var1TeHm0weightsTpJ
107.97451.25400.05899.29436.0311e+03
2110.79452.64140.035212.58103.7004e+04
326.90201.95310.05208.04431.2516e+04
4312.66767.31010.005114.76413.6745e+05
5412.89372.26230.016015.02763.6455e+04
6510.55764.75430.017312.30491.1678e+05
768.76672.73940.043610.21763.1826e+04
876.53741.30560.05077.61945.2724e+03
989.66631.34070.037111.26618.4436e+03
10912.78733.92040.016514.90361.0768e+05
111011.60591.82100.022313.52671.9446e+04
12117.58411.87870.05468.83931.2827e+04
131210.17546.13390.008811.85951.8819e+05
14139.31944.58740.018810.86179.5155e+04
15147.22901.25670.05778.42545.4490e+03
16155.64701.33910.03216.58154.7508e+03
17167.61602.63460.03658.87642.5340e+04
18179.46043.38110.033811.02615.2509e+04
191816.00043.04420.004318.64858.7447e+04
201910.55031.56360.030312.29641.2622e+04
212011.81742.98290.021713.77325.3748e+04
222112.10115.30570.014514.10391.7731e+05
232210.40003.58830.032312.12136.5405e+04
24239.13252.01160.048110.64401.7913e+04
25249.88032.46290.050611.51552.9157e+04
26258.32182.00310.05429.69911.6105e+04
27266.13141.79440.03547.14619.2960e+03
282711.43073.97990.025313.32259.0735e+04
292814.26382.78170.009416.62456.6183e+04
302913.74425.46520.007516.01882.4048e+05
31308.73531.27060.047110.18096.8213e+03
32318.30173.67680.02219.67575.4123e+04

2. Write a WEC-Sim batch file for the given clusters

WEC-Sim MCR (multiple condition run) files should contain a structure `mcr` that contains two variables: `header` and `cases`. Each column of `header` and `cases` denotes a variable and it's value respectively. Each row is distinct simulation. WEC-Sim defines waves using the significant wave height and peak period. We will isolate these values from the results of the cluster analysis and create a dictionary that is written to the `.mat` file.
mcr = struct();
mcr.header = {'waves.height','waves.period'};
mcr.cases = [results.Hm0 results.Tp]
mcr = struct with fields:
header: {'waves.height' 'waves.period'} +.S8 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 0px none rgb(33, 33, 33); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S9 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 1px solid rgb(217, 217, 217); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; } +.S10 { border-left: 1px solid rgb(217, 217, 217); border-right: 1px solid rgb(217, 217, 217); border-top: 0px none rgb(33, 33, 33); border-bottom: 1px solid rgb(217, 217, 217); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 18.004px; min-height: 0px; white-space: nowrap; color: rgb(33, 33, 33); font-family: Menlo, Monaco, Consolas, "Courier New", monospace, Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; }

MHKiT WEC-Sim Power Performance Example

This notebook demonstrates using the MHKiT wave module and WEC-Sim together to perform a resource characterization study in MHKiT, simulate representative cases with WEC-Sim, and visualize the results in MHKiT to estimate MAEP (Mean Annual Energy Production).
1. Characterize the available resource at a location
- See the PacWave example notebook
2. Write a WEC-Sim batch file for the given clusters
3. Simulate the device in WEC-Sim.
- Ensure that the spectra used in WEC-Sim is identical to the one used in MHKiT.
4. Load WEC-Sim batch results
5. Assess results and visualize quantities of interest
This example uses WEC-Sim to simulate the Oscillating Surge Wave Energy Converter (OSWEC), a flap-type device.

1. Characterize the available resource at a location

This example will use an abbreviated version of PacWave_resource_characterization_example.ipynb.
For full details on downloading, calculating, and visualizing the k-means clusters representation of the site's wave resouce, see that example.
We select the N=32 cluster as it's total energy flux is closet to the total energy flux of the site considering all wave conditions. We will load the PacWave example output, which can be easily saved after running the example with the command results[32].to_csv("pacwave_cluster_32.csv"). We will start this example by reading in that csv output and formatting it for WEC-Sim.
% Relative location and filename of simulated WEC-Sim data (run with mooring)
filename = './data/wave/pacwave_cluster_32.csv'
filename = './data/wave/pacwave_cluster_32.csv'
results = readtable(filename)
results = 32×6 table
Var1TeHm0weightsTpJ
107.97451.25400.05899.29436.0311e+03
2110.79452.64140.035212.58103.7004e+04
326.90201.95310.05208.04431.2516e+04
4312.66767.31010.005114.76413.6745e+05
5412.89372.26230.016015.02763.6455e+04
6510.55764.75430.017312.30491.1678e+05
768.76672.73940.043610.21763.1826e+04
876.53741.30560.05077.61945.2724e+03
989.66631.34070.037111.26618.4436e+03
10912.78733.92040.016514.90361.0768e+05
111011.60591.82100.022313.52671.9446e+04
12117.58411.87870.05468.83931.2827e+04
131210.17546.13390.008811.85951.8819e+05
14139.31944.58740.018810.86179.5155e+04
⋮

2. Write a WEC-Sim batch file for the given clusters

WEC-Sim MCR (multiple condition run) files should contain a structure mcr that contains two variables: header and cases. Each column of header and cases denotes a variable and it's value respectively. Each row is distinct simulation. WEC-Sim defines waves using the significant wave height and peak period. We will isolate these values from the results of the cluster analysis and create a dictionary that is written to the .mat file.
mcr = struct();
mcr.header = {'waves.height','waves.period'};
mcr.cases = [results.Hm0 results.Tp]
mcr = struct with fields:
header: {'waves.height' 'waves.period'} cases: [32×2 double] -
save('mcr_mhkit.mat', 'mcr');

3. Simulate the device in WEC-Sim

Now that the MCR file is created, we need to go simulate WEC performance in these wave conditions using WEC-Sim. To recreate the data used in the next step, use the created MCR file with WEC-Sim's OSWEC example. For an accurate comparison to the power calculated in the resource characterization, we should ensure that the WEC-Sim cases use irregular JONSWAP wave spectra as in the PacWave example.
For convenience in this demonstration, we enforce OSWEC model stability in the extreme wave conditions by arbitrarily applying a large PTO stiffness and damping in the wecSimInputFile.m:
% pto(1).stiffness = 1e5;
% pto(1).damping = 5e7;
To reduce the amount of extranenous data saved for this example, we limit the WEC-Sim output to the PTO's power output in the `userDefinedFunctions.m` script:
% if exist('imcr','var')
% if imcr == 1
% nmcr = size(mcr.cases,1);
% power = nan(1, nmcr);
% end
%
% iRampEnd = simu.rampTime./simu.dtOut + 1;
% power(imcr) = -mean(output.ptos(1).powerInternalMechanics(iRampEnd:end,5));
%
% if imcr == nmcr
% % % Save output
% save('mcr_mhkit_power.mat', 'power');
% end
% end
% bodies = output.bodies;

4. Load WEC-Sim batch results

Note that in this example we do not save the entire WEC-Sim `output` structure for each case. See the `wecsim_example.ipynb` for information on loading the WEC-Sim responseClass. Here, the output is one array of average power output that we will load and compare to the resource characterization.
Note that the power output [W] is significantly larger than the energy flux [W/m] due to the large width of the OSWEC.
% Relative location and filename of simulated WEC-Sim data (run with mooring)
filename = './data/wave/mcr_mhkit_power.mat';
% Load the WEC-Sim output data which contains the variable `power`.
load(filename)
results.Power = power'
results = 32×7 table
 Var1TeHm0weightsTpJPower
107.97451.25400.05899.29436.0311e+036.8613e+04
2110.79452.64140.035212.58103.7004e+043.8735e+05
326.90201.95310.05208.04431.2516e+041.9234e+05
4312.66767.31010.005114.76413.6745e+051.9512e+06
5412.89372.26230.016015.02763.6455e+043.1159e+05
6510.55764.75430.017312.30491.1678e+058.4103e+05
768.76672.73940.043610.21763.1826e+043.3986e+05
876.53741.30560.05077.61945.2724e+038.3326e+04
989.66631.34070.037111.26618.4436e+036.9516e+04
10912.78733.92040.016514.90361.0768e+055.0688e+05
111011.60591.82100.022313.52671.9446e+041.6020e+05
12117.58411.87870.05468.83931.2827e+042.1125e+05
131210.17546.13390.008811.85951.8819e+052.4020e+06
14139.31944.58740.018810.86179.5155e+041.0677e+06
15147.22901.25670.05778.42545.4490e+036.4143e+04
16155.64701.33910.03216.58154.7508e+037.3289e+04
17167.61602.63460.03658.87642.5340e+043.6720e+05
18179.46043.38110.033811.02615.2509e+044.0202e+05
191816.00043.04420.004318.64858.7447e+042.4145e+05
201910.55031.56360.030312.29641.2622e+041.3351e+05
212011.81742.98290.021713.77325.3748e+043.5522e+05
222112.10115.30570.014514.10391.7731e+051.3432e+06
232210.40003.58830.032312.12136.5405e+047.3683e+05
24239.13252.01160.048110.64401.7913e+041.9988e+05
25249.88032.46290.050611.51552.9157e+042.7157e+05
26258.32182.00310.05429.69911.6105e+042.3197e+05
27266.13141.79440.03547.14619.2960e+031.4780e+05
282711.43073.97990.025313.32259.0735e+049.3911e+05
292814.26382.78170.009416.62456.6183e+042.0003e+05
302913.74425.46520.007516.01882.4048e+058.6862e+05
31308.73531.27060.047110.18096.8213e+038.5504e+04
32318.30173.67680.02219.67575.4123e+047.5789e+05

5. Assess results and visualize quantities of interest

Now that we have loaded the OSWEC's modeled power, we can assess it's performance relative to the incoming wave and calculate the mean annual energy production (MAEP) using MHKiT.
results.CW = capture_length(results.Power, results.J)';
 
oswec_width = 18;
results.CWR = results.CW / oswec_width
results = 32×9 table
 Var1TeHm0weightsTpJPowerCWCWR
107.97451.25400.05899.29436.0311e+036.8613e+0411.37650.6320
2110.79452.64140.035212.58103.7004e+043.8735e+0510.46770.5815
326.90201.95310.05208.04431.2516e+041.9234e+0515.36730.8537
4312.66767.31010.005114.76413.6745e+051.9512e+065.31000.2950
5412.89372.26230.016015.02763.6455e+043.1159e+058.54710.4748
6510.55764.75430.017312.30491.1678e+058.4103e+057.20150.4001
768.76672.73940.043610.21763.1826e+043.3986e+0510.67870.5933
876.53741.30560.05077.61945.2724e+038.3326e+0415.80410.8780
989.66631.34070.037111.26618.4436e+036.9516e+048.23290.4574
10912.78733.92040.016514.90361.0768e+055.0688e+054.70730.2615
111011.60591.82100.022313.52671.9446e+041.6020e+058.23800.4577
12117.58411.87870.05468.83931.2827e+042.1125e+0516.46880.9149
131210.17546.13390.008811.85951.8819e+052.4020e+0612.76380.7091
14139.31944.58740.018810.86179.5155e+041.0677e+0611.22080.6234
15147.22901.25670.05778.42545.4490e+036.4143e+0411.77140.6540
16155.64701.33910.03216.58154.7508e+037.3289e+0415.42660.8570
17167.61602.63460.03658.87642.5340e+043.6720e+0514.49100.8051
18179.46043.38110.033811.02615.2509e+044.0202e+057.65630.4253
191816.00043.04420.004318.64858.7447e+042.4145e+052.76110.1534
201910.55031.56360.030312.29641.2622e+041.3351e+0510.57700.5876
212011.81742.98290.021713.77325.3748e+043.5522e+056.60900.3672
222112.10115.30570.014514.10391.7731e+051.3432e+067.57580.4209
232210.40003.58830.032312.12136.5405e+047.3683e+0511.26560.6259
24239.13252.01160.048110.64401.7913e+041.9988e+0511.15810.6199
25249.88032.46290.050611.51552.9157e+042.7157e+059.31410.5174
26258.32182.00310.05429.69911.6105e+042.3197e+0514.40370.8002
27266.13141.79440.03547.14619.2960e+031.4780e+0515.89990.8833
282711.43073.97990.025313.32259.0735e+049.3911e+0510.35010.5750
292814.26382.78170.009416.62456.6183e+042.0003e+053.02240.1679
302913.74425.46520.007516.01882.4048e+058.6862e+053.61210.2007
31308.73531.27060.047110.18096.8213e+038.5504e+0412.53480.6964
32318.30173.67680.02219.67575.4123e+047.5789e+0514.00310.7779
Calculate and display the mean annual energy production.
CW.values = results.CW;
J.values = results.J;
weights.values = results.weights;
MAEP = mean_annual_energy_production_matrix(CW, J, weights) / 1000 % kWh
MAEP = 2.8651e+06
+
save('mcr_mhkit.mat', 'mcr');

3. Simulate the device in WEC-Sim

Now that the MCR file is created, we need to go simulate WEC performance in these wave conditions using WEC-Sim. To recreate the data used in the next step, use the created MCR file with WEC-Sim's OSWEC example. For an accurate comparison to the power calculated in the resource characterization, we should ensure that the WEC-Sim cases use irregular JONSWAP wave spectra as in the PacWave example.
For convenience in this demonstration, we enforce OSWEC model stability in the extreme wave conditions by arbitrarily applying a large PTO stiffness and damping in the wecSimInputFile.m:
% pto(1).stiffness = 1e5;
% pto(1).damping = 5e7;
To reduce the amount of extranenous data saved for this example, we limit the WEC-Sim output to the PTO's power output in the userDefinedFunctions.m script:
% if exist('imcr','var')
% if imcr == 1
% nmcr = size(mcr.cases,1);
% power = nan(1, nmcr);
% end
%
% iRampEnd = simu.rampTime./simu.dtOut + 1;
% power(imcr) = -mean(output.ptos(1).powerInternalMechanics(iRampEnd:end,5));
%
% if imcr == nmcr
% % % Save output
% save('mcr_mhkit_power.mat', 'power');
% end
% end
% bodies = output.bodies;

4. Load WEC-Sim batch results

Note that in this example we do not save the entire WEC-Sim output structure for each case. See the wecsim_example.ipynb for information on loading the WEC-Sim responseClass. Here, the output is one array of average power output that we will load and compare to the resource characterization.
Note that the power output [W] is significantly larger than the energy flux [W/m] due to the large width of the OSWEC.
% Relative location and filename of simulated WEC-Sim data (run with mooring)
filename = './data/wave/mcr_mhkit_power.mat';
% Load the WEC-Sim output data which contains the variable |power|.
load(filename)
results.Power = power'
results = 32×7 table
Var1TeHm0weightsTpJPower
107.97451.25400.05899.29436.0311e+036.8613e+04
2110.79452.64140.035212.58103.7004e+043.8735e+05
326.90201.95310.05208.04431.2516e+041.9234e+05
4312.66767.31010.005114.76413.6745e+051.9512e+06
5412.89372.26230.016015.02763.6455e+043.1159e+05
6510.55764.75430.017312.30491.1678e+058.4103e+05
768.76672.73940.043610.21763.1826e+043.3986e+05
876.53741.30560.05077.61945.2724e+038.3326e+04
989.66631.34070.037111.26618.4436e+036.9516e+04
10912.78733.92040.016514.90361.0768e+055.0688e+05
111011.60591.82100.022313.52671.9446e+041.6020e+05
12117.58411.87870.05468.83931.2827e+042.1125e+05
131210.17546.13390.008811.85951.8819e+052.4020e+06
14139.31944.58740.018810.86179.5155e+041.0677e+06
⋮

5. Assess results and visualize quantities of interest

Now that we have loaded the OSWEC's modeled power, we can assess it's performance relative to the incoming wave and calculate the mean annual energy production (MAEP) using MHKiT.
results.CW = capture_width(results.Power, results.J);
 
oswec_width = 18;
results.CWR = results.CW / oswec_width
results = 32×9 table
Var1TeHm0weightsTpJPowerCWCWR
107.97451.25400.05899.29436.0311e+036.8613e+0411.37650.6320
2110.79452.64140.035212.58103.7004e+043.8735e+0510.46770.5815
326.90201.95310.05208.04431.2516e+041.9234e+0515.36730.8537
4312.66767.31010.005114.76413.6745e+051.9512e+065.31000.2950
5412.89372.26230.016015.02763.6455e+043.1159e+058.54710.4748
6510.55764.75430.017312.30491.1678e+058.4103e+057.20150.4001
768.76672.73940.043610.21763.1826e+043.3986e+0510.67870.5933
876.53741.30560.05077.61945.2724e+038.3326e+0415.80410.8780
989.66631.34070.037111.26618.4436e+036.9516e+048.23290.4574
10912.78733.92040.016514.90361.0768e+055.0688e+054.70730.2615
111011.60591.82100.022313.52671.9446e+041.6020e+058.23800.4577
12117.58411.87870.05468.83931.2827e+042.1125e+0516.46880.9149
131210.17546.13390.008811.85951.8819e+052.4020e+0612.76380.7091
14139.31944.58740.018810.86179.5155e+041.0677e+0611.22080.6234
⋮
Calculate and display the mean annual energy production.
CW.values = results.CW;
J.values = results.J;
weights.values = results.weights;
MAEP = mean_annual_energy_production_matrix(CW, J, weights) / 1000 % kWh
MAEP = 2.8651e+06

') + % Column headers are 2 lines after this + if i+2 <= length(lines) + col_line_idx = i + 2; + end + break; + end + end + + % Parse column headers if found + columns = []; + if col_line_idx > 0 && col_line_idx <= length(lines) + col_str = regexprep(lines{col_line_idx}, '^%\s*', ''); + col_parts = strsplit(strtrim(col_str)); + columns = str2double(col_parts); + columns = columns(~isnan(columns)); + end + + % Find data start (first non-% line after headers) + data_start = 0; + for i = start_idx+1:end_idx + if i <= length(lines) && ~startsWith(lines{i}, '%') && ~isempty(strtrim(lines{i})) + data_start = i; + break; + end + end + + if data_start == 0 + return; + end + + % Parse data rows + row_data = {}; + y_indices = []; + + for i = data_start:end_idx + if i > length(lines) + break; + end + line = lines{i}; + + if startsWith(line, '%') + continue; + end + + trimmed = strtrim(line); + if isempty(trimmed) + continue; + end + + % Split on whitespace and periods + % SWAN block format: " 100 100.100.100.100..." + % First split by whitespace to get Y index and data string + parts = regexp(trimmed, '\s+', 'split'); + if length(parts) < 2 + continue; + end + + y_idx = str2double(parts{1}); + if isnan(y_idx) + continue; + end + + % Join remaining parts and split by period + data_str = strjoin(parts(2:end), ''); + + % Split by period (but handle **** for NaN) + % Replace **** patterns first + data_str = regexprep(data_str, '\*+', 'NaN.'); + + % Split by period + val_parts = strsplit(data_str, '.'); + val_parts = val_parts(~cellfun(@isempty, val_parts)); + + % Convert to numbers + vals = zeros(1, length(val_parts)); + for k = 1:length(val_parts) + if strcmpi(val_parts{k}, 'NaN') + vals(k) = NaN; + else + vals(k) = str2double(val_parts{k}); + end + end + + y_indices(end+1) = y_idx; %#ok + row_data{end+1} = vals; %#ok + end + + if isempty(row_data) + return; + end + + % Determine matrix size + n_rows = length(row_data); + n_cols = max(cellfun(@length, row_data)); + + % Build matrix (Y indices go from high to low in file) + matrix = NaN(n_rows, n_cols); + for r = 1:n_rows + vals = row_data{r}; + matrix(r, 1:length(vals)) = vals; + end + + % Apply unit multiplier + matrix = matrix * unit_multiplier; +end + + +function metaDict = parse_line_metadata(line) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Parse the key:value pairs of a SWAN "% Run" header line, e.g. +% "% Run:TEST Frame: COMPGRID ** Significant wave height, Unit: 0.1000E-01 m" +% +% Parameters +% ------------ +% line : char +% One "% Run" header line +% +% Returns +% --------- +% metaDict : structure +% One field per key (Run, Frame, vars, Unit, etc.) holding the +% trimmed value text. "**" is read as the key "vars" +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + metaDict = struct(); + + % Remove leading % and whitespace + line = regexprep(line, '^%\s*', ''); + + % Replace ** with vars: + line = strrep(line, '**', 'vars:'); + + % Replace commas and extra whitespace + line = regexprep(line, ',', ' '); + line = regexprep(line, '\s+', ' '); + + % Split by colon + parts = strsplit(line, ':'); + + % Parse key:value pairs + for i = 1:length(parts)-1 + % Current part ends with key name + current_words = strsplit(strtrim(parts{i})); + key = current_words{end}; + + % Next part starts with value (everything except last word which is next key) + next_words = strsplit(strtrim(parts{i+1})); + if i < length(parts)-1 + % Value is all words except the last (which is the next key) + val = strjoin(next_words(1:end-1), ' '); + else + % Last value - take all words + val = strjoin(next_words, ' '); + end + + metaDict.(key) = strtrim(val); + end +end diff --git a/mhkit/wave/graphics/plot_matrix.m b/mhkit/wave/graphics/plot_matrix.m index 949619697..3ab9da4ab 100644 --- a/mhkit/wave/graphics/plot_matrix.m +++ b/mhkit/wave/graphics/plot_matrix.m @@ -1,67 +1,66 @@ -function figure=plot_matrix(M,Mtype, options) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Plots the matrix with Hm0 and Te on the y and x axis -% -% Parameters -% ---------- -% M: structure -% -% M.values: matrix -% -% M.Hm0_bins -% -% M.Te_bins -% -% M.stat +function ax = plot_matrix(M, Mtype, options) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Mtype: string -% type of matrix (i.e. power, capture length, etc.) to be used -% in plot title +% Plots a wave performance matrix with Hm0 on the y axis and Te on the x axis % -% options: name-value pairs -% savepath: string (optional) -% path and filename to save figure. -% annotate: logical (optional) -% toggle text annotations on/off (default: true) -% to call: plot_matrix(M, Mtype,"savepath",savepath,"annotate",false) +% Parameters +% ------------ +% M : struct +% Matrix from capture_width_matrix, wave_energy_flux_matrix, or power_matrix +% M.values : matrix +% One value per bin, Hm0 bins down the rows and Te bins across the columns +% M.stat : string +% Statistic used, shown in the title +% M.x_edges : vector [s] +% Te bin edges +% M.y_edges : vector [m] +% Hm0 bin edges +% Mtype : string +% Type of matrix (e.g. "Capture Width", "Power") used in the plot title +% savepath : string (optional) +% Path and filename to save the figure. Default none +% annotate : logical (optional) +% Print each bin value in its cell. Default true +% value_format : string (optional) +% sprintf format for the bin labels, e.g. '%.2f m'. Default '%.2f' +% font_size : double (optional) +% Font size of the bin labels in points. Default 7 +% trim_to_data : logical (optional) +% Limit the axes to the bins with data plus one empty bin around them. Default false +% zlabel : string (optional) +% Colorbar label, e.g. "Capture Width [m]". Default none % % Returns % --------- -% figure: plot of the matrix -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -arguments - M - Mtype - options.savepath = ""; - options.annotate = true; +% ax : matlab.graphics.axis.Axes +% Axes containing the plot +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + M (1,1) struct + Mtype {mustBeTextScalar} + options.savepath {mustBeTextScalar} = "" + options.annotate (1,1) logical = true + options.zlabel {mustBeTextScalar} = "" + options.trim_to_data (1,1) logical = false + options.value_format {mustBeTextScalar} = '%.2f' + options.font_size (1,1) {mustBeNumeric, mustBePositive} = 7 +end + +arguments (Output) + ax (1,1) matlab.graphics.axis.Axes end -figure=pcolor(M.Te_bins,M.Hm0_bins,M.values); -colormap(flipud(hot(256))) -ylabel('Hm0 [m]','FontSize',20) -xlabel('Te [s]','FontSize',20) -x=strcat(Mtype,': ',M.stat); -title(x) -colorbar -pos=get(gca,'position'); -[rows,cols]=size(M.values); -width=pos(3)/(cols-1); -height =pos(4)/(rows-1); +ax = mhkit_plot_matrix(M, 'xlabel', 'Energy Period, T_e [sec]', 'ylabel', 'Significant Wave Height, H_{m0} [m]', ... + 'show_values', options.annotate, 'zlabel', options.zlabel, ... + 'trim_to_data', options.trim_to_data, 'value_format', options.value_format, ... + 'font_size', options.font_size); +title(ax, string(Mtype) + ": " + string(M.stat)); -% create textbox annotations if annotate is true -if options.annotate - for i=1:cols-1 - for j=rows-1:-1:1 - if ~isnan(M.values(j,i)) - annotation('textbox',[pos(1)+width*(i-1),pos(2)+height*(j-1),width,height], ... - 'string',num2str(M.values(j,i)),'LineStyle','none','HorizontalAlignment','center',... - 'VerticalAlignment','middle'); - end - end - end +if strlength(options.savepath) > 0 + saveas(ancestor(ax, 'figure'), options.savepath); end -len = strlength(options.savepath); -if len > 1 - saveas(figure, options.savepath); end diff --git a/mhkit/wave/graphics/plot_wave_joint_probability_distribution.m b/mhkit/wave/graphics/plot_wave_joint_probability_distribution.m new file mode 100644 index 000000000..863e3ace4 --- /dev/null +++ b/mhkit/wave/graphics/plot_wave_joint_probability_distribution.m @@ -0,0 +1,121 @@ +function [ax, jpd] = plot_wave_joint_probability_distribution(Hm0, Te, options) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Plot the joint probability distribution (JPD) of sea states, the +% fraction of records in each Hm0 and Te bin, as a percentage +% +% By default the bins follow IEC TS 62600-100: 0.5 m wide for Hm0 and +% 1 s wide for Te with edges on 0, 0.5, 1, ... m and 0, 1, 2, ... s, +% extending to the largest sea state. The returned matrix struct has +% the same layout as capture_width_matrix so its bins can be reused. +% +% Parameters +% ------------ +% Hm0 : vector +% Significant wave height of each record [m] +% Te : vector +% Energy period of each record [s], same length as Hm0 +% Hm0_bins : vector (optional) +% Hm0 bin centers [m]. Default IEC 0.5 m bins +% Te_bins : vector (optional) +% Te bin centers [s]. Default IEC 1 s bins +% xlabel : string (optional) +% Default "Energy Period, T_e [sec]" +% ylabel : string (optional) +% Default "Significant Wave Height, H_{m0} [m]" +% title : string (optional) +% Default "Joint Probability Distribution" +% colormap : string or matrix (optional) +% Colormap name ("viridis", any cmocean or MATLAB colormap name, +% "-" or "_r" to flip) or an N x 3 RGB matrix, see mhkit_colormap. +% Default "viridis" +% annotate : logical (optional) +% Print the percentage in each occupied bin. Default true +% font_size : double (optional) +% Font size of the bin labels in points. Default 8 +% trim_to_data : logical (optional) +% Limit the axes to the occupied bins plus one empty bin around +% them. Default true +% ax : axes handle (optional) +% Axes to plot into. Default is the current axes +% savepath : string (optional) +% Path and filename to save the figure. Default none +% +% Returns +% --------- +% ax : axes handle +% Axes the JPD was plotted into +% jpd : structure +% jpd.values : matrix +% Fraction of records in each bin, Hm0 rows by Te columns, +% summing to one +% jpd.stat : char +% 'probability' +% jpd.x_bins, x_edges : vector +% Te bin centers and edges [s] +% jpd.y_bins, y_edges : vector +% Hm0 bin centers and edges [m] +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + Hm0 {mustBeNumeric, mustBeVector} + Te {mustBeNumeric, mustBeVector} + options.Hm0_bins {mustBeNumeric} = [] + options.Te_bins {mustBeNumeric} = [] + options.xlabel {mustBeTextScalar} = "Energy Period, T_e [sec]" + options.ylabel {mustBeTextScalar} = "Significant Wave Height, H_{m0} [m]" + options.title {mustBeTextScalar} = "Joint Probability Distribution" + options.colormap = "viridis" + options.annotate (1,1) logical = true + options.font_size (1,1) {mustBeNumeric, mustBePositive} = 8 + options.trim_to_data (1,1) logical = true + options.ax = [] + options.savepath {mustBeTextScalar} = "" +end + +arguments (Output) + ax matlab.graphics.axis.Axes + jpd struct +end + +if numel(Hm0) ~= numel(Te) + error('MHKiT:plot_wave_joint_probability_distribution:InvalidInput', ... + 'Hm0 and Te must have the same number of elements'); +end + +% IEC TS 62600-100 bin centers, so the edges fall on multiples of the width +if isempty(options.Hm0_bins) + options.Hm0_bins = 0.25:0.5:ceil(max(Hm0) / 0.5) * 0.5; +end +if isempty(options.Te_bins) + options.Te_bins = 0.5:1:ceil(max(Te)); +end + +bin_spec = struct('x', struct('centers', options.Te_bins), ... + 'y', struct('centers', options.Hm0_bins)); +% 'probability' ignores the values argument, so Hm0 is passed as a placeholder +jpd = mhkit_binned_statistic_2d(Te, Hm0, Hm0, 'probability', bin_spec, ... + 'function_name', mfilename); + +% Plot as a percentage with the empty bins left blank +jpd_percent = jpd; +jpd_percent.values = 100 * jpd.values; +jpd_percent.values(jpd_percent.values == 0) = NaN; + +if isempty(options.ax) + options.ax = gca; +end +ax = mhkit_plot_matrix(jpd_percent, 'xlabel', options.xlabel, 'ylabel', options.ylabel, ... + 'zlabel', 'Occurrence [%]', 'colormap', options.colormap, ... + 'show_values', options.annotate, 'value_format', '%.2f %%', ... + 'font_size', options.font_size, 'trim_to_data', options.trim_to_data, ... + 'ax', options.ax); +title(ax, options.title); + +if strlength(options.savepath) > 0 + saveas(ancestor(ax, 'figure'), options.savepath); +end + +end diff --git a/mhkit/wave/performance/capture_length.m b/mhkit/wave/performance/capture_length.m index b8f8743e9..6e121c154 100644 --- a/mhkit/wave/performance/capture_length.m +++ b/mhkit/wave/performance/capture_length.m @@ -1,28 +1,41 @@ -function L = capture_length(Power, J) +function L = capture_length(P, J) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the capture length (often called capture width). +% +% Deprecated alias for capture_width. +% +% IEC TS 62600-100 Ed. 2.0 replaces "capture length" with "capture +% width". This function will be removed in MHKiT-MATLAB v1.3; use +% capture_width instead. % % Parameters % ------------ -% P: array or vector -% Power [W] -% -% J: array or vector -% Omnidirectional wave energy flux [W/m] +% P : numeric [W] +% Power (scalar, vector, or array) +% J : numeric [W/m] +% Omnidirectional wave energy flux (same shape as P) % % Returns % --------- -% L: vector -% Capture length [m] +% L : numeric [m] +% Capture length (same shape as P and J) % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +arguments (Input) + P {mustBeNumeric} + J {mustBeNumeric} +end +arguments (Output) + L {mustBeNumeric} +end -Power = py.numpy.array(Power); -J = py.numpy.array(J); +warning('MHKiT:capture_length:DeprecatedFunction', ... + ['IEC TS 62600-100 Ed. 2.0 replaces "capture length" with "capture ' ... + 'width". capture_length will be removed in MHKiT-MATLAB v1.3. ' ... + 'Use capture_width instead.']); -L = py.mhkit.wave.performance.capture_length(Power, J); +L = capture_width(P, J); -L = typecast_from_mhkit_python(L).data; +end diff --git a/mhkit/wave/performance/capture_length_matrix.m b/mhkit/wave/performance/capture_length_matrix.m index 25cf9b16f..00946478c 100644 --- a/mhkit/wave/performance/capture_length_matrix.m +++ b/mhkit/wave/performance/capture_length_matrix.m @@ -1,62 +1,56 @@ -function clm=capture_length_matrix(Hm0,Te,L,statistic,Hm0_bins,Te_bins) +function clm = capture_length_matrix(Hm0, Te, L, statistic, Hm0_bins, Te_bins) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Generates a capture length matrix for a given statistic +% Deprecated alias for capture_width_matrix. % -% Note that IEC/TS 62600-100 requires capture length matrices for -% the mean, std, count, min, and max. +% IEC TS 62600-100 Ed. 2.0 replaces "capture length" with "capture +% width". This function will be removed in MHKiT-MATLAB v1.3; use +% capture_width_matrix instead. % % Parameters % ------------ -% Hm0: numpy array or vector -% Significant wave height from spectra [m] -% -% Te: numpy array or vector -% Energy period from spectra [s] -% -% L : numpy array or vector -% Capture length [m] -% -% statistic: string -% Statistic for each bin, options include: 'mean', 'std', 'median', -% 'count', 'sum', 'min', 'max', and 'frequency'. Note that 'std' uses -% a degree of freedom of 1 in accordance with IEC/TS 62600-100. or -% a callable function of python type -% -% Hm0_bins: numpy array or vector -% Bin centers for Hm0 [m] -% -% Te_bins: numpy array or vector -% Bin centers for Te [s] +% Hm0 : vector [m] +% Significant wave height from spectra +% Te : vector [s] +% Energy period from spectra +% L : vector [m] +% Capture length +% statistic : char or string +% Statistic for each bin. Options: 'mean', 'std', 'median', +% 'count', 'sum', 'min', 'max', 'probability', or 'frequency'. +% 'probability' and 'frequency' are the same statistic. +% Hm0_bins : numeric vector [m] +% Hm0 bin centers, see capture_width_matrix +% Te_bins : numeric vector [s] +% Te bin centers, see capture_width_matrix % % Returns % --------- -% clm: structure -% +% clm : struct +% Same layout as capture_width_matrix: values (Hm0 rows x Te columns), +% stat, x_bins, y_bins, x_edges, y_edges with Te on x and Hm0 on y % -% clm.values -% -% clm.stat -% -% clm.Hm0_bins -% -% clm.Te_bins -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Hm0 = py.numpy.array(Hm0); -Te = py.numpy.array(Te); -L = py.numpy.array(L); +arguments (Input) + Hm0 {mustBeNumeric} + Te {mustBeNumeric} + L {mustBeNumeric} + statistic {mustBeTextScalar} + Hm0_bins {mustBeNumeric, mustBeVector} + Te_bins {mustBeNumeric, mustBeVector} +end -py_Hm0_bins = py.numpy.array(Hm0_bins); -py_Te_bins = py.numpy.array(Te_bins); +arguments (Output) + clm struct +end -LM = py.mhkit.wave.performance.capture_length_matrix(Hm0, Te, L, statistic, py_Hm0_bins, py_Te_bins); +warning('MHKiT:capture_length_matrix:DeprecatedFunction', ... + ['IEC TS 62600-100 Ed. 2.0 replaces "capture length" with "capture ' ... + 'width". capture_length_matrix will be removed in MHKiT-MATLAB v1.3. ' ... + 'Use capture_width_matrix instead.']); -LM = typecast_from_mhkit_python(LM); +clm = capture_width_matrix(Hm0, Te, L, statistic, Hm0_bins, Te_bins); -clm.values = LM.capture_length_matrix; -clm.stat = statistic; -clm.Hm0_bins = Hm0_bins; -clm.Te_bins = Te_bins; +end diff --git a/mhkit/wave/performance/capture_width.m b/mhkit/wave/performance/capture_width.m new file mode 100644 index 000000000..7498df5a5 --- /dev/null +++ b/mhkit/wave/performance/capture_width.m @@ -0,0 +1,34 @@ +function CW = capture_width(P, J) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Calculates the capture width (sometimes called capture length) +% +% CW = P / J +% +% Parameters +% ------------ +% P : numeric [W] +% Power (scalar, vector, or array) +% J : numeric [W/m] +% Omnidirectional wave energy flux (same shape as P) +% +% Returns +% --------- +% CW : numeric [m] +% Capture width (same shape as P and J) +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + P {mustBeNumeric} + J {mustBeNumeric} +end + +arguments (Output) + CW {mustBeNumeric} +end + +CW = P ./ J; + +end diff --git a/mhkit/wave/performance/capture_width_matrix.m b/mhkit/wave/performance/capture_width_matrix.m new file mode 100644 index 000000000..6382b3d4d --- /dev/null +++ b/mhkit/wave/performance/capture_width_matrix.m @@ -0,0 +1,76 @@ +function cwm = capture_width_matrix(Hm0, Te, CW, statistic, Hm0_bins, Te_bins) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Generates a capture width matrix for a given statistic +% +% Note that IEC TS 62600-100 Ed. 2.0 en 2024 section 9.2.4 requires +% capture width matrices for the mean, std, count, min, and max. +% +% Parameters +% ------------ +% Hm0 : vector [m] +% Significant wave height from spectra +% Te : vector [s] +% Energy period from spectra +% CW : vector [m] +% Capture width +% statistic : char or string +% Statistic for each bin. Options: 'mean', 'std', 'median', +% 'count', 'sum', 'min', 'max', 'probability', or 'frequency'. +% 'probability' and 'frequency' are the same statistic. +% 'std' is the population standard deviation (1/N), following +% MHKiT-Python convention. +% Hm0_bins : numeric vector [m] +% Hm0 bin centers, strictly increasing. Edges are the midpoints between +% centers, extended half a spacing beyond the end centers, following +% MHKiT-Python convention. Bins may be non-uniform. +% Te_bins : numeric vector [s] +% Te bin centers, same rules as Hm0_bins +% +% Returns +% --------- +% cwm : struct +% cwm.values : matrix +% Capture width matrix, Hm0 bins down the rows and Te bins across the columns +% cwm.stat : string +% Statistic used +% cwm.x_bins : row vector [s] +% Te bin centers +% cwm.y_bins : row vector [m] +% Hm0 bin centers +% cwm.x_edges : row vector [s] +% Te bin edges, one more than the number of columns +% cwm.y_edges : row vector [m] +% Hm0 bin edges, one more than the number of rows +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + Hm0 {mustBeNumeric} + Te {mustBeNumeric} + CW {mustBeNumeric} + statistic {mustBeTextScalar, mustBeMember(statistic, {'mean', 'std', 'median', 'count', 'sum', 'min', 'max', 'probability', 'frequency'})} + Hm0_bins {mustBeNumeric, mustBeVector} + Te_bins {mustBeNumeric, mustBeVector} +end + +arguments (Output) + cwm struct +end + +% Te is x and Hm0 is y, following the IEC TS 62600-100 scatter diagram +bin_spec = struct('x', struct('centers', Te_bins), 'y', struct('centers', Hm0_bins)); + +cwm = mhkit_binned_statistic_2d(Te, Hm0, CW, statistic, bin_spec, 'function_name', mfilename); + +if any(diff(cwm.x_edges) > 1.0) + warning('MHKiT:capture_width_matrix:BinSpacing', ... + 'Energy period bins are greater than the IEC TS 62600-100 limit of 1.0 seconds.'); +end +if any(diff(cwm.y_edges) > 0.5) + warning('MHKiT:capture_width_matrix:BinSpacing', ... + 'Significant wave height bins are greater than the IEC TS 62600-100 limit of 0.5 meters.'); +end + +end diff --git a/mhkit/wave/performance/mean_annual_energy_production_matrix.m b/mhkit/wave/performance/mean_annual_energy_production_matrix.m index 506203932..3fb6fcc8c 100644 --- a/mhkit/wave/performance/mean_annual_energy_production_matrix.m +++ b/mhkit/wave/performance/mean_annual_energy_production_matrix.m @@ -1,70 +1,72 @@ -function maep=mean_annual_energy_production_matrix(LM,JM,frequency) +function maep = mean_annual_energy_production_matrix(CWM, JM, frequency) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Calculates mean annual energy production (MAEP) from matrix data -% along with data frequency in each bin +% Calculates mean annual energy production (MAEP) from matrix data +% +% MAEP = T * nansum(CWM * JM * frequency) +% where T = 8766 hours (average length of a year) % % Parameters % ------------ -% LM: Capture Length -% 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(Hm0_bins,L) -% -% OR -% -% structure of form: -% -% LM.values -% -% LM.stat -% -% LM.Hm0_bins -% -% LM.Te_bins -% -% -% JM: Wave Energy Flux -% 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(Hm0_bins,J) -% -% OR -% -% structure of form: -% -% JM.values: Wave energy flux matrix -% -% JM.Hm0_bins -% -% JM.Te_bins -% -% frequency: Data frequency for each bin. -% 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(Hm0_bins,frequency) -% -% OR -% -% structure of form: -% -% frequency.values -% -% frequency.Hm0_bins -% -% frequency.Te_bins +% CWM : struct or matrix +% Capture width matrix. If struct, uses CWM.values +% JM : struct or matrix +% Wave energy flux matrix. If struct, uses JM.values +% frequency : struct or matrix +% Data frequency for each bin (must sum to 1). If struct, uses frequency.values % % Returns % --------- -% maep: float -% Mean annual energy production +% maep : double [W*h] +% Mean annual energy production % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + CWM + JM + frequency +end + +arguments (Output) + maep {mustBeNumeric} +end + +% Extract values from structs if needed +if isstruct(CWM) + CWM_vals = CWM.values; +else + CWM_vals = CWM; +end + +if isstruct(JM) + JM_vals = JM.values; +else + JM_vals = JM; +end + +if isstruct(frequency) + freq_vals = frequency.values; +else + freq_vals = frequency; +end + +% Validate dimensions match +if ~isequal(size(CWM_vals), size(JM_vals)) || ~isequal(size(CWM_vals), size(freq_vals)) + error('MHKiT:mean_annual_energy_production_matrix:DimensionMismatch', ... + 'CWM, JM, and frequency must have the same dimensions'); +end + +% Validate frequency sums to 1 +freq_sum = sum(freq_vals(:), 'omitnan'); +if abs(freq_sum - 1) > 1e-6 + error('MHKiT:mean_annual_energy_production_matrix:InvalidFrequency', ... + 'Frequency components must sum to one. Got: %g', freq_sum); +end + +T = 8766; % Average length of a year in hours -LM_py = py.pandas.DataFrame(py.numpy.array(LM.values)); -JM_py = py.pandas.DataFrame(py.numpy.array(JM.values)); -freq_py = py.pandas.DataFrame(py.numpy.array(frequency.values)); +maep = T * sum(CWM_vals(:) .* JM_vals(:) .* freq_vals(:), 'omitnan'); -maep = py.mhkit.wave.performance.mean_annual_energy_production_matrix(LM_py,JM_py,freq_py); -maep = double(maep); +end diff --git a/mhkit/wave/performance/mean_annual_energy_production_timeseries.m b/mhkit/wave/performance/mean_annual_energy_production_timeseries.m index 3f49eec52..1f85d3b8d 100644 --- a/mhkit/wave/performance/mean_annual_energy_production_timeseries.m +++ b/mhkit/wave/performance/mean_annual_energy_production_timeseries.m @@ -1,29 +1,106 @@ -function maep=mean_annual_energy_production_timeseries(L,J) +function maep = mean_annual_energy_production_timeseries(CW, J, options) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Calculates mean annual energy production (MAEP) from timeseries +% Calculates mean annual energy production (MAEP) from time-series +% +% MAEP = (T/n) * sum(CW * J) +% where T = 8766 hours (average length of a year) +% +% Each sample is weighted equally, so the time series is assumed to be +% representative of a full year of sea states. A warning is issued when +% a time axis is available and the record spans less than one year. % % Parameters % ------------ -% L: numpy array or vector -% Capture length -% -% J: numpy array or vector -% Wave energy flux +% CW : vector or timetable [m] +% Capture width. A timetable must have one variable, and its row +% times are used as the time axis. +% J : vector or timetable [W/m] +% Wave energy flux. A timetable must have one variable, and its row +% times are used as the time axis. +% time : datetime, duration, or numeric vector [s] (optional) +% Time axis for CW and J, used only to check that the record spans at +% least one year. Ignored if CW or J is a timetable. % % Returns % --------- -% maep: float -% Mean annual energy production +% maep : double [W*h] +% Mean annual energy production % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + CW {mustBeA(CW, {'numeric', 'timetable'})} + J {mustBeA(J, {'numeric', 'timetable'})} + options.time {mustBeA(options.time, {'numeric', 'datetime', 'duration'})} = [] +end + +arguments (Output) + maep {mustBeNumeric} +end + +[CW_vals, CW_time] = extract_values_and_time(CW, 'CW'); +[J_vals, J_time] = extract_values_and_time(J, 'J'); + +if length(CW_vals) ~= length(J_vals) + error('MHKiT:mean_annual_energy_production_timeseries:InvalidInput', ... + 'CW length (%d) must match J length (%d)', ... + length(CW_vals), length(J_vals)); +end + +% Prefer timetable row times over the optional time argument +if ~isempty(CW_time) + time = CW_time; +elseif ~isempty(J_time) + time = J_time; +else + time = options.time(:); +end + +T = 8766; % Average length of a year in hours per IEC 62600-101 Ed 2.0 Section A.2 +n = length(CW_vals); +if ~isempty(time) + if length(time) ~= n + error('MHKiT:mean_annual_energy_production_timeseries:InvalidInput', ... + 'time length (%d) must match CW length (%d)', length(time), n); + end + span = max(time) - min(time); + if isnumeric(span) + span_hours = span / 3600; + else + span_hours = hours(span); + end + if span_hours < T + warning('MHKiT:mean_annual_energy_production_timeseries:RecordShorterThanOneYear', ... + ['Time series spans %.1f days, less than one year (%.2f days). ' ... + 'MAEP weights every sample equally and assumes the record is ' ... + 'representative of a full year.'], span_hours / 24, T / 24); + end +end -py.importlib.import_module('mhkit'); +maep = (T / n) * sum(CW_vals(:) .* J_vals(:)); -J=py.numpy.array(J); -L=py.numpy.array(L); +end -maep=double(py.mhkit.wave.performance.mean_annual_energy_production_timeseries(L,J)); +function [vals, time] = extract_values_and_time(x, name) +% Returns the data as a column vector plus row times for timetable input. +if istimetable(x) + if width(x) ~= 1 + error('MHKiT:mean_annual_energy_production_timeseries:InvalidInput', ... + '%s timetable must have exactly one variable, got %d', name, width(x)); + end + vals = x{:, 1}; + time = x.Properties.RowTimes; +else + vals = x; + time = []; +end +if ~isvector(vals) + error('MHKiT:mean_annual_energy_production_timeseries:InvalidInput', ... + '%s must be a vector', name); +end +vals = vals(:); +end diff --git a/mhkit/wave/performance/power_matrix.m b/mhkit/wave/performance/power_matrix.m index 51a2f4d3b..fda0bd74b 100644 --- a/mhkit/wave/performance/power_matrix.m +++ b/mhkit/wave/performance/power_matrix.m @@ -1,106 +1,75 @@ -function PM=power_matrix(LM,JM) +function PM = power_matrix(CWM, JM) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Generates a power matrix from a capture length matrix and wave energy -% flux matrix +% Generates a power matrix from a capture width matrix and wave energy flux matrix +% +% PM = CWM * JM (element-wise multiplication) % % Parameters % ------------ -% LM: Capture Length -% 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(Hm0_bins,L) -% -% OR -% -% structure of form: -% -% LM.values -% -% LM.stat -% -% LM.Hm0_bins -% -% LM.Te_bins -% -% -% JM: Wave Energy Flux -% 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(Hm0_bins,J) -% -% OR -% -% structure of form: -% -% JM.values: Wave energy flux matrix -% -% JM.Hm0_bins -% -% JM.Te_bins +% CWM : struct +% Capture width matrix from capture_width_matrix: +% CWM.values : matrix +% Capture width values +% CWM.stat : string +% Statistic used (e.g., 'mean') +% CWM.x_bins, CWM.y_bins : row vectors +% Te and Hm0 bin centers +% CWM.x_edges, CWM.y_edges : row vectors +% Te and Hm0 bin edges, must match JM +% JM : struct +% Wave energy flux matrix structure: +% JM.values : matrix +% Wave energy flux values +% JM.x_bins, JM.y_bins : row vectors +% Te and Hm0 bin centers +% JM.x_edges, JM.y_edges : row vectors +% Te and Hm0 bin edges, must match CWM % % Returns % --------- -% PM: Structure -% -% -% PM.values: Power matrix +% PM : struct +% Power matrix structure: +% PM.values : matrix +% Power matrix values +% PM.stat : string +% Statistic from CWM +% PM.x_bins, PM.y_bins : row vectors +% Te and Hm0 bin centers from CWM +% PM.x_edges, PM.y_edges : row vectors +% Te and Hm0 bin edges from CWM % -% PM.stat: statistic of the matrix (i.e. "mean", "max", etc.) -% (string) -% -% PM.Hm0_bins -% -% PM.Te_bins -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -py.importlib.import_module('mhkit'); -py.importlib.import_module('mhkit_python_utils'); - -% need to add asserts for pandas -if (isa(LM,'py.pandas.core.frame.DataFrame')~=1) - x=size(LM.values); - - li=py.list(); - if x(2)>1 - for i = 1:x(2) - app=py.list(LM.values(:,i)); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); - - end - end - - LMpan=py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(li,py.list(LM.Hm0_bins),int32(x(2))); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +arguments (Input) + CWM struct + JM struct end -if (isa(JM,'py.pandas.core.frame.DataFrame')~=1) - x=size(JM.values); - li=py.list(); - if x(2)>1 - for i = 1:x(2) - app=py.list(JM.values(:,i)); - li=py.mhkit_python_utils.pandas_dataframe.lis(li,app); +arguments (Output) + PM struct +end - end +% Both inputs must be matrix structs from mhkit_binned_statistic_2d on the same grid +for f = {'values', 'stat', 'x_bins', 'y_bins', 'x_edges', 'y_edges'} + if ~isfield(CWM, f{1}) || ~isfield(JM, f{1}) + error('MHKiT:power_matrix:InvalidInput', ... + 'CWM and JM must both be matrix structs from capture_width_matrix and wave_energy_flux_matrix, missing field %s', f{1}); end - JMpan=py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(li,py.list(JM.Hm0_bins),int32(x(2))); +end +if ~isequal(CWM.x_edges, JM.x_edges) || ~isequal(CWM.y_edges, JM.y_edges) + error('MHKiT:power_matrix:InvalidInput', ... + 'CWM and JM must be binned on the same x and y edges'); end -PMpan=py.mhkit.wave.performance.power_matrix(LMpan,JMpan); - -% convert vals1d to original shape -vals=double(py.array.array('d',py.numpy.nditer(PMpan.values))); -sha=cell(PMpan.values.shape); -x=int64(sha{1,1}); -y=int64(sha{1,2}); -vals=reshape(vals,[y,x]); -vals=transpose(vals); - -PM.values=vals; -PM.stat=LM.stat; -PM.Hm0_bins=double(LM.Hm0_bins); -PM.Te_bins=double(LM.Te_bins); +% Build output structure +PM = struct(); +PM.values = CWM.values .* JM.values; +PM.stat = CWM.stat; +PM.x_bins = CWM.x_bins; +PM.y_bins = CWM.y_bins; +PM.x_edges = CWM.x_edges; +PM.y_edges = CWM.y_edges; +end diff --git a/mhkit/wave/performance/power_performance_workflow.m b/mhkit/wave/performance/power_performance_workflow.m index ef1d12500..e6a4c37a8 100644 --- a/mhkit/wave/performance/power_performance_workflow.m +++ b/mhkit/wave/performance/power_performance_workflow.m @@ -1,8 +1,9 @@ -function [clmat,maep_matrix] = power_performance_workflow(S, h, P, statistic, options) +function [cwmat, maep_matrix] = power_performance_workflow(S, h, P, statistic, options) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % High-level function to compute power performance quantities of -% interest following IEC TS 62600-100 for given wave spectra. +% interest following IEC TS 62600-100 Ed. 2.0 en 2024 for given wave +% spectra. % % Parameters % ------------ @@ -17,14 +18,13 @@ % P: array or vector % Power [W] % -% statistic: string or array of strings -% Capture length statistics for plotting -% options include: "mean", "std", "median", -% "count", "sum", "min", "max", and "frequency". -% Note that "std" -% uses a degree of freedom of 1 in accordance with IEC/TS 62600-100. -% To output capture length matrices for multiple binning parameters, -% define as a string array: statistic = ["", "", ""]; +% statistic: string or string array +% Capture width statistics to compute and plot, one figure each. +% Options: "mean", "std", "median", "count", "sum", "min", "max", +% "probability", "frequency". "probability" and "frequency" are the +% same statistic. "std" is the population standard deviation, +% matching MHKiT-Python. "mean" and "probability" are always +% computed because MAEP needs them. % % savepath: string (optional) % Path to save figure. @@ -43,18 +43,63 @@ % % Returns % --------- -% cl_matrix: figure -% Capture length matrix +% cwmat: struct +% One capture width matrix struct per computed statistic, keyed by +% statistic name, e.g. cwmat.mean +% Capture width matrix % % maep_matrix: float % Mean annual energy production +% +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> % cast to double: fzero (inside wave_number, called by energy_flux +% >> % below) requires double, but CDIP returns frequency as single. +% >> S.frequency = double(data.metadata.wave.waveFrequency); % [Hz], 64x1 double +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> S.spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% >> S.time = data.data.wave.waveTime; % 17520x1 datetime +% +% >> % Water depth at CDIP 225 (WETS Berth 3), per PacIOOS buoy +% >> % metadata: "Moored in water 80 meters deep". +% >> h = 80; % [m] +% +% >> % SYNTHETIC power output - NOT measured. CDIP 225 is a wave buoy +% >> % only, with no WEC device attached, so no real matching power +% >> % record exists for this site. A real capture width matrix needs +% >> % power telemetry from an actual device. +% >> rng(1); +% >> P = 50000 + 20000*randn(size(S.time)); % [W], illustrative only, 17520x1 double +% +% >> [cwmat, maep_matrix] = power_performance_workflow(S, h, P, "mean"); +% maep_matrix = +% 456328960.94 % [W*h] +% +% >> cwmat.mean +% ans = +% struct with fields: +% values: [8×14 double] +% stat: 'mean' +% x_bins: [0.5000 1.5000 2.5000 3.5000 4.5000 5.5000 … ] (1×14 double) +% y_bins: [0.2500 0.7500 1.2500 1.7500 2.2500 2.7500 3.2500 3.7500] +% x_edges: [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14] +% y_edges: [0 0.5000 1 1.5000 2 2.5000 3 3.5000 4] +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% arguments S h P - statistic + statistic (1,:) string {mustBeMember(statistic, ["mean", "std", "median", "count", "sum", "min", "max", "probability", "frequency"])} options.rho = 1025; options.g = 9.80665; options.frequency_bins = ""; @@ -62,85 +107,56 @@ end if ~isstruct(S) - ME = MException('MATLAB:power_performance_wave','S must be a structure. See API header for formatting'); - throw(ME); + error('MHKiT:power_performance_workflow:InvalidInput', ... + 'S must be a structure with spectrum and frequency fields.'); end if any([~isnumeric(h), ~isnumeric(P)]) - ME = MException('MATLAB:power_performance_wave','h and P must be numbers'); - throw(ME); + error('MHKiT:power_performance_workflow:InvalidInput', 'h and P must be numeric.'); end -% Check dimensions of the spectrum data -[num_frequencies, num_spectra] = size(S.spectrum); - -% Initialize arrays to store results -Te = zeros(num_spectra, 1); -Hm0 = zeros(num_spectra, 1); -J = zeros(num_spectra, 1); - -% Process each spectrum separately -for i = 1:num_spectra - % Create a temporary structure for the current spectrum - S_temp = struct(); - S_temp.spectrum = S.spectrum(:,i); - S_temp.frequency = S.frequency; - - % Compute the energy periods from spectra data - Te(i) = energy_period(S_temp); - - % Compute the significant wave height from spectra data - Hm0(i) = significant_wave_height(S_temp); - - % Compute the energy flux from spectra data and water depth - J(i) = energy_flux(S_temp, h); -end; clear S_temp - -% calculating capture length with power and wave flux in vectors -L = capture_length(P,J); - -% Need to set our Hm0 and Te bins for the capture length matrix -Hm0_bins = -0.5:0.5:max(fix(Hm0))+0.5; % Input is min, max, and n indecies for vector -Hm0_bins = Hm0_bins+0.25 ; -Te_bins = 0:1:max(fix(Te)); -Te_bins = Te_bins+0.5; - -% Calculate the necessary capture length matrices for each statistic based -% on IEC/TS 62600-100 -clmat.mean = capture_length_matrix(Hm0, Te, L, "mean", Hm0_bins, Te_bins); -clmat.std = capture_length_matrix(Hm0, Te, L, "std", Hm0_bins, Te_bins); -clmat.median = capture_length_matrix(Hm0 ,Te, L, "median", Hm0_bins, Te_bins); -clmat.count = capture_length_matrix(Hm0 ,Te, L, "count", Hm0_bins, Te_bins); -clmat.sum = capture_length_matrix(Hm0 ,Te, L, "sum", Hm0_bins, Te_bins); -clmat.min = capture_length_matrix(Hm0 ,Te, L, "min", Hm0_bins, Te_bins); -clmat.max = capture_length_matrix(Hm0 ,Te, L, "max", Hm0_bins, Te_bins); -clmat.freq = capture_length_matrix(Hm0 ,Te, L, "frequency", Hm0_bins, Te_bins); +% energy_period/significant_wave_height/energy_flux accept S.spectrum as +% a full NxK matrix directly, one column per spectrum - no per-spectrum +% loop needed. +if isnumeric(options.frequency_bins) + Te = energy_period(S, options.frequency_bins); + Hm0 = significant_wave_height(S, options.frequency_bins); +else + Te = energy_period(S); + Hm0 = significant_wave_height(S); +end +J = energy_flux(S, h, 'rho', options.rho, 'g', options.g); + +% calculating capture width with power and wave flux in vectors +% Ensure P is a column vector to match J +P = P(:); +CW = capture_width(P, J); + +% Bin centers with IEC TS 62600-100 widths of 0.5 m for Hm0 and 1 s for Te, +% so the edges fall on 0, 0.5, 1, ... m and 0, 1, 2, ... s. +Hm0_bins = 0.25:0.5:ceil(max(Hm0) / 0.5) * 0.5; +Te_bins = 0.5:1:ceil(max(Te)); + +% mean and probability are always needed for MAEP. Any other requested +% statistic is computed once and stored under its own name. +cwmat = struct(); +for stat = unique(["mean", "probability", statistic]) + cwmat.(stat) = capture_width_matrix(Hm0, Te, CW, stat, Hm0_bins, Te_bins); +end -% Create wave energy flux matrix using statistic jmat = wave_energy_flux_matrix(Hm0, Te, J, "mean", Hm0_bins, Te_bins); -% Calcaulte MAEP from matrix -maep_matrix = mean_annual_energy_production_matrix(clmat.mean, jmat, clmat.freq); -stats_cell = {'mean', 'std', 'median','count', 'sum', 'min', 'max','frequency'}; +maep_matrix = mean_annual_energy_production_matrix(cwmat.mean, jmat, cwmat.probability); -% Capture Length Matrix using statistic -cl_matrix = []; -len = strlength(options.savepath); -for i = 1:length(statistic) - if any(strcmp(stats_cell,statistic(i))) - figure('Name',sprintf('Capture Length Matrix %s', statistic(i)),'NumberTitle','off') - cl_matrix(i) = plot_matrix(clmat.(statistic(i)),"Capture Length"); - name = [options.savepath, filesep, sprintf('Capture Length Matrix %s', statistic(i)), '.png']; - - if len > 1 - saveas(cl_matrix(i), name); - end - else - ME = MException('MATLAB:power_performance_wave',... - 'statistic must be a string or string array defined', ... - 'by one or multiple of the following: "mean", "std", "median","count", "sum", "min", "max", "frequency"'); - throw(ME); +% Plot each requested statistic +cw_matrix = gobjects(1, numel(statistic)); +for i = 1:numel(statistic) + figure('Name', sprintf('Capture Width Matrix %s', statistic(i)), 'NumberTitle', 'off') + cw_matrix(i) = plot_matrix(cwmat.(statistic(i)), "Capture Width"); + if strlength(options.savepath) > 0 + name = fullfile(options.savepath, sprintf('Capture Width Matrix %s.png', statistic(i))); + saveas(cw_matrix(i), name); end end -end \ No newline at end of file +end diff --git a/mhkit/wave/performance/wave_energy_flux_matrix.m b/mhkit/wave/performance/wave_energy_flux_matrix.m index 880b71b1a..eff8dc435 100644 --- a/mhkit/wave/performance/wave_energy_flux_matrix.m +++ b/mhkit/wave/performance/wave_energy_flux_matrix.m @@ -1,68 +1,73 @@ -function WEFM=wave_energy_flux_matrix(Hm0,Te,J,statistic,Hm0_bins,Te_bins) +function WEFM = wave_energy_flux_matrix(Hm0, Te, J, statistic, Hm0_bins, Te_bins) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Generates a wave eneergy flux matrix for a given statistic -% -% Note that IEC/TS 62600-100 requires capture length matrices for -% the mean, std, count, min, and max. +% Generates a wave energy flux matrix for a given statistic % % Parameters % ------------ -% Hm0: numpy array or vector -% Significant wave height from spectra [m] -% -% Te: numpy array or vector -% Energy period from spectra [s] -% -% J : numpy array or vector -% wave energy flux from spectra [W/m] -% -% statistic: string -% Statistic for each bin, options include: 'mean', 'std', 'median', -% 'count', 'sum', 'min', 'max', and 'frequency'. Note that 'std' uses -% a degree of freedom of 1 in accordance with IEC/TS 62600-100. -% -% Hm0_bins: numpy array or vector -% Bin centers for Hm0 [m] -% -% Te_bins: numpy array or vector -% Bin centers for Te [s] +% Hm0 : vector [m] +% Significant wave height from spectra +% Te : vector [s] +% Energy period from spectra +% J : vector [W/m] +% Wave energy flux from spectra +% statistic : char or string +% Statistic for each bin. Options: 'mean', 'std', 'median', +% 'count', 'sum', 'min', 'max', 'probability', or 'frequency'. +% 'probability' and 'frequency' are the same statistic. +% 'std' is the population standard deviation (1/N), following +% MHKiT-Python convention. +% Hm0_bins : numeric vector [m] +% Hm0 bin centers, strictly increasing. Edges are the midpoints between +% centers, extended half a spacing beyond the end centers, following +% MHKiT-Python convention. Bins may be non-uniform. +% Te_bins : numeric vector [s] +% Te bin centers, same rules as Hm0_bins % % Returns % --------- -% WEFM: Structure +% WEFM : struct +% WEFM.values : matrix +% Wave energy flux matrix, Hm0 bins down the rows and Te bins across the columns +% WEFM.stat : string +% Statistic used +% WEFM.x_bins : row vector [s] +% Te bin centers +% WEFM.y_bins : row vector [m] +% Hm0 bin centers +% WEFM.x_edges : row vector [s] +% Te bin edges, one more than the number of columns +% WEFM.y_edges : row vector [m] +% Hm0 bin edges, one more than the number of rows % -% -% WEFM.values -% -% WEFM.stat -% -% WEFM.Hm0_bins -% -% WEFM.Te_bins -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +arguments (Input) + Hm0 {mustBeNumeric} + Te {mustBeNumeric} + J {mustBeNumeric} + statistic {mustBeTextScalar, mustBeMember(statistic, {'mean', 'std', 'median', 'count', 'sum', 'min', 'max', 'probability', 'frequency'})} + Hm0_bins {mustBeNumeric, mustBeVector} + Te_bins {mustBeNumeric, mustBeVector} +end -py.importlib.import_module('mhkit'); +arguments (Output) + WEFM struct +end -Hm0=py.numpy.array(Hm0); -Te=py.numpy.array(Te); -J=py.numpy.array(J); -Hm0_bins=py.numpy.array(Hm0_bins); -Te_bins=py.numpy.array(Te_bins); +% Te is x and Hm0 is y, following the IEC TS 62600-100 scatter diagram +bin_spec = struct('x', struct('centers', Te_bins), 'y', struct('centers', Hm0_bins)); -JM=py.mhkit.wave.performance.wave_energy_flux_matrix(Hm0,Te,J,statistic,Hm0_bins,Te_bins); -vals=double(py.array.array('d',py.numpy.nditer(JM.values))); -sha=cell(JM.values.shape); -x=int64(sha{1,1}); -y=int64(sha{1,2}); -vals=reshape(vals,[y,x]); -vals=transpose(vals); +WEFM = mhkit_binned_statistic_2d(Te, Hm0, J, statistic, bin_spec, 'function_name', mfilename); -WEFM.values=vals; -WEFM.stat=statistic; -WEFM.Hm0_bins=double(Hm0_bins); -WEFM.Te_bins=double(Te_bins); +if any(diff(WEFM.x_edges) > 1.0) + warning('MHKiT:wave_energy_flux_matrix:BinSpacing', ... + 'Energy period bins are greater than the IEC TS 62600-100 limit of 1.0 seconds.'); +end +if any(diff(WEFM.y_edges) > 0.5) + warning('MHKiT:wave_energy_flux_matrix:BinSpacing', ... + 'Significant wave height bins are greater than the IEC TS 62600-100 limit of 0.5 meters.'); +end +end diff --git a/mhkit/wave/performance/wec_reference_power_matrices/RM3.csv b/mhkit/wave/performance/wec_reference_power_matrices/RM3.csv new file mode 100644 index 000000000..927bb7964 --- /dev/null +++ b/mhkit/wave/performance/wec_reference_power_matrices/RM3.csv @@ -0,0 +1,24 @@ +# RM3: Heaving Buoy, Hydraulic Drivetrain, characteristic diameter 20 m, unballasted structural mass 687 t, floating foundation, slack mooring, A36 Steel +# Source: National Laboratory of the Rockies (NLR) System Advisor Model (SAM) 2025.4.16, libraries/Wave Energy Converters.csv, BSD-3-Clause, https://github.com/NatLabRockies/SAM +# First row: Te bin centers [s]; first column: Hs bin centers [m]; body: device power [kW] +0,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5 +0.25,0,0,0,0,0.4,0.6,0.8,1,1.1,1.1,1,0.8,0.7,0.6,0.5,0.4,0.3,0.3,0.2,0.2,0 +0.75,0,0,0,0,3.2,5.3,7.4,9.1,9.8,9.5,8.6,7.4,6.2,5.1,4.1,3.4,2.8,2.3,1.9,1.6,0 +1.25,0,0,0,0,9,14.8,20.5,25,26.8,25.9,23.3,20,16.8,13.8,11.3,9.2,7.6,6.3,5.2,4.3,0 +1.75,0,0,0,0,17.6,28.9,39.9,48.3,51.6,49.7,44.7,38.4,32.2,26.5,21.7,17.8,14.6,12.1,10,8.4,0 +2.25,0,0,0,0,29,47.5,65.4,78.8,83.8,80.6,72.4,62.3,52.2,43,35.3,28.9,23.8,19.7,16.3,13.7,0 +2.75,0,0,0,0,43.2,70.7,97,116.3,123.1,118.1,106.1,91.3,76.5,63.2,51.9,42.5,35,28.9,24.1,20.1,0 +3.25,0,0,0,0,60.2,98.3,134.5,160.5,169.3,162.1,145.5,125.2,105,86.8,71.3,58.5,48.2,39.9,33.2,27.8,0 +3.75,0,0,0,0,79.9,130.4,177.8,211.2,222,212.2,190.4,164,137.6,113.8,93.6,76.9,63.3,52.5,43.7,36.6,0 +4.25,0,0,0,0,102.4,166.7,226.7,268.3,281.1,268.2,240.5,207.2,174.1,144.1,118.5,97.4,80.3,66.6,55.5,46.5,0 +4.75,0,0,0,0,127.6,207.4,281.2,286,286,286,286,255,214.3,177.5,146.1,120.2,99.2,82.2,68.6,57.6,0 +5.25,0,0,0,0,155.4,252.4,286,286,286,286,286,286,258.2,214,176.3,145.1,119.8,99.4,83,69.7,0 +5.75,0,0,0,0,186,286,286,286,286,286,286,286,286,253.6,209,172.2,142.2,118.1,98.6,82.8,0 +6.25,0,0,0,0,219.2,286,286,286,286,286,286,286,286,286,244.2,201.2,166.4,138.2,115.5,97.1,0 +6.75,0,0,0,0,255,286,286,286,286,286,286,286,286,286,281.9,232.4,192.2,159.7,133.5,112.3,0 +7.25,0,0,0,0,286,286,286,286,286,286,286,286,286,286,286,265.6,219.8,182.8,152.9,128.7,0 +7.75,0,0,0,0,286,286,286,286,286,286,286,286,286,286,286,286,249,207.2,173.4,146,0 +8.25,0,0,0,0,286,286,286,286,286,286,286,286,286,286,286,286,279.9,233,195.1,164.4,0 +8.75,0,0,0,0,286,286,286,286,286,286,286,286,286,286,286,286,286,260.2,218,183.8,0 +9.25,0,0,0,0,286,286,286,286,286,286,286,286,286,286,286,286,286,286,242.1,204.1,0 +9.75,0,0,0,0,286,286,286,286,286,286,286,286,286,286,286,286,286,286,267.4,225.6,0 diff --git a/mhkit/wave/performance/wec_reference_power_matrices/RM5.csv b/mhkit/wave/performance/wec_reference_power_matrices/RM5.csv new file mode 100644 index 000000000..2c7554347 --- /dev/null +++ b/mhkit/wave/performance/wec_reference_power_matrices/RM5.csv @@ -0,0 +1,24 @@ +# RM5: Oscillating Surge Wave Converter, Hydraulic Drivetrain, characteristic diameter 25 m, unballasted structural mass 800 t, floating foundation, taut mooring, Fiberglass and A36 Steel +# Source: National Laboratory of the Rockies (NLR) System Advisor Model (SAM) 2025.4.16, libraries/Wave Energy Converters.csv, BSD-3-Clause, https://github.com/NatLabRockies/SAM +# First row: Te bin centers [s]; first column: Hs bin centers [m]; body: device power [kW] +0,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5 +0.25,0,0,0,0,4.8,6.7,7.9,9.3,10.2,10.1,9.7,9,8.8,7.6,7.3,6.4,5.6,5,4.5,4,0 +0.75,0,0,0,0,12.3,16.5,18.8,21.2,22.9,22.2,20.9,19.4,18.7,16.5,16,14.2,12.8,11.5,10.4,9.4,0 +1.25,0,0,0,0,31.8,40.7,44.6,48.5,51.7,48.8,45.1,41.8,40.1,36.2,35.1,31.9,29.2,26.5,24.3,22,0 +1.75,0,0,0,0,58.3,72.3,77.1,81.7,86.5,80.8,74,69.7,66.7,59.7,57.6,52.7,48.7,44.5,41.1,37.6,0 +2.25,0,0,0,0,91.3,110.4,115.7,119.3,126.5,117.3,107.9,102,97.1,86.4,82.6,75.6,70.5,64.7,60.3,55.3,0 +2.75,0,0,0,0,130.5,154.9,160,162.7,171.7,158.5,145.4,137.5,130.4,115.6,109.7,101.4,94.6,86.6,80.8,74,0 +3.25,0,0,0,0,174.9,204.4,208.9,210.4,220.5,202.7,185.4,175.4,165.9,148,140.3,129.7,120.5,110.1,102.2,93.4,0 +3.75,0,0,0,0,223.9,258.5,261.9,261.6,272.4,249.5,227.7,215.3,204.5,183.2,173,159.8,147.9,134.8,124.8,113.7,0 +4.25,0,0,0,0,277.2,316.8,318.5,316,327,298.4,271.6,257.2,245.5,220.2,207.3,191.5,177.1,161.8,149.7,136.8,0 +4.75,0,0,0,0,334.5,360,360,360,360,349.4,317.2,302.2,288.2,258.7,243.1,225.4,208.6,190.3,176.1,160.7,0 +5.25,0,0,0,0,360,360,360,360,360,360,360,348.9,332.4,298.6,280.1,261.3,241.4,220,203.3,185.2,0 +5.75,0,0,0,0,360,360,360,360,360,360,360,360,360,339.7,319.1,298.4,275.5,250.8,231.5,210.7,0 +6.25,0,0,0,0,360,360,360,360,360,360,360,360,360,360,360,340.8,314.3,285.8,263.5,239.7,0 +6.75,0,0,0,0,360,360,360,360,360,360,360,360,360,360,360,360,358.6,325.8,300.1,272.6,0 +7.25,0,0,0,0,360,360,360,360,360,360,360,360,360,360,360,360,360,360,341.6,310.1,0 +7.75,0,0,0,0,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,352.8,0 +8.25,0,0,0,0,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,0 +8.75,0,0,0,0,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,0 +9.25,0,0,0,0,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,0 +9.75,0,0,0,0,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,360,0 diff --git a/mhkit/wave/performance/wec_reference_power_matrices/RM6.csv b/mhkit/wave/performance/wec_reference_power_matrices/RM6.csv new file mode 100644 index 000000000..ac360c7d8 --- /dev/null +++ b/mhkit/wave/performance/wec_reference_power_matrices/RM6.csv @@ -0,0 +1,24 @@ +# RM6: Backward Bent Duct Buoy Floating Oscillating Water Column, Wells Air Turbine, characteristic diameter 27 m, unballasted structural mass 1809 t, floating foundation, slack mooring, A36 Steel +# Source: National Laboratory of the Rockies (NLR) System Advisor Model (SAM) 2025.4.16, libraries/Wave Energy Converters.csv, BSD-3-Clause, https://github.com/NatLabRockies/SAM +# First row: Te bin centers [s]; first column: Hs bin centers [m]; body: device power [kW] +0,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5 +0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 +0.75,0,0,0,0,0,0,0,0,0,0,0,0,0,0.2,0.2,0.2,0.2,0.1,0.1,0,0 +1.25,0,0,0,0,0,0,0,0,4.3,14.7,19,21.3,18.8,16.2,13,10.1,7.6,5.7,4.2,3.1,0 +1.75,0,0,0,0,0,0,0.2,35.2,68.4,75.1,70.9,64.1,55.3,46.6,38.3,30.8,24.1,18.9,14.8,11.6,0 +2.25,0,0,0,0,0,0.1,26.6,128.3,145.6,142,132.2,120,106.1,91.4,76.1,61.2,48.8,39.2,31.4,25.1,0 +2.75,0,0,0,0,1.6,0.6,121.3,215.7,229,219.9,202.5,182,159.3,137.2,117,98.4,81.5,65.3,52.4,42,0 +3.25,0,0,0,0,19.8,36.5,212.4,288.1,321.7,305,278.8,250.2,220.9,188.3,159.8,133.6,110.3,93.2,77.3,64.1,0 +3.75,0,0,0,0,46.3,120.2,212.4,288.1,338.7,350.5,348.6,319.8,278.8,237.8,202.8,170.4,140.7,119.9,101.5,86,0 +4.25,0,0,0,0,84.7,142.3,212.4,288.1,338.7,350.5,344.8,350.3,337.8,284.8,233.4,186.4,149.2,129.5,112.8,98.3,0 +4.75,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,349.9,289.4,241.1,196,160.2,137.8,118,101.1,0 +5.25,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,309.3,261.9,206.1,172.1,146.6,123.5,104,0 +5.75,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,216.8,184.8,156,129.2,107,0 +6.25,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,205.4,189.1,169.8,153.4,145.6,0 +6.75,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,205.9,193.4,181.1,171.3,166.9,0 +7.25,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,207.9,198.3,189.8,183.5,180.7,0 +7.75,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,209.9,202.5,196.3,191.9,190,0 +8.25,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,211.4,205.7,200.9,197.6,196.2,0 +8.75,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,212.5,207.8,204,201.4,200.3,0 +9.25,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,213.2,209.1,205.8,203.5,202.6,0 +9.75,0,0,0,0,114,142.3,212.4,288.1,338.7,350.5,344.8,339.7,350.1,276.1,217.7,213.4,209.6,206.5,204.4,203.5,0 diff --git a/mhkit/wave/performance/wec_reference_power_matrix.m b/mhkit/wave/performance/wec_reference_power_matrix.m new file mode 100644 index 000000000..5ec99637d --- /dev/null +++ b/mhkit/wave/performance/wec_reference_power_matrix.m @@ -0,0 +1,90 @@ +function M = wec_reference_power_matrix(device) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Power matrix of a DOE Reference Model Wave Energy Converter +% +% The Reference Model Project, sponsored by the U.S. Department of Energy, +% developed open-source marine energy point designs as reference models to +% benchmark technology performance and costs: +% https://openei.org/wiki/PRIMRE/Signature_Projects/Reference_Model +% Returns the device power [kW] over significant wave height and energy +% period bins for Reference Model 3 (wave point absorber), 5 (oscillating +% surge flap), or 6 (oscillating water column), as distributed with the +% National Laboratory of the Rockies (NLR) System Advisor Model (SAM, +% BSD-3-Clause, https://github.com/NatLabRockies/SAM). The struct has the same layout +% as the matrices from capture_width_matrix and power_matrix, so it can be +% plotted with mhkit_plot_matrix and used with interp2 to look up power for +% measured sea states. +% +% Parameters +% ------------ +% device : string +% Reference model name: "RM3" (wave point absorber), "RM5" (oscillating +% surge flap), or "RM6" (oscillating water column) +% +% Returns +% --------- +% M : struct +% M.values : matrix [kW] +% Device power, Hs bins down the rows and Te bins across the columns +% M.stat : char +% 'power' +% M.x_bins : row vector [s] +% Te bin centers, 0.5 to 20.5 s by 1 s +% M.y_bins : row vector [m] +% Hs bin centers, 0.25 to 9.75 m by 0.5 m +% M.x_edges : row vector [s] +% Te bin edges +% M.y_edges : row vector [m] +% Hs bin edges +% M.device : char +% Device name +% M.description : char +% Device description from the SAM library +% M.source : char +% Data source and license +% +% Examples +% -------- +% >> M = wec_reference_power_matrix("RM3"); +% >> size(M.values) +% ans = +% 20 21 +% +% >> max(M.values(:)) +% ans = +% 286 +% +% >> % Power [W] for measured sea states, zero outside the matrix +% >> P = 1000 * interp2(M.x_bins, M.y_bins, M.values, Te, Hm0, 'linear', 0); +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + device {mustBeTextScalar, mustBeMember(device, ["RM3", "RM5", "RM6"])} +end + +arguments (Output) + M (1,1) struct +end + +file = fullfile(fileparts(mfilename('fullpath')), 'wec_reference_power_matrices', char(device) + ".csv"); + +% Header comment lines hold the device description and the data source +header = strsplit(strtrim(fileread(file)), newline); +header = header(startsWith(header, '#')); +grid = readmatrix(file, 'CommentStyle', '#'); + +M = struct(); +M.values = grid(2:end, 2:end); +M.stat = 'power'; +M.x_bins = grid(1, 2:end); +M.y_bins = grid(2:end, 1)'; +M.x_edges = mhkit_define_bins_1d(struct('centers', M.x_bins)).edges'; +M.y_edges = mhkit_define_bins_1d(struct('centers', M.y_bins)).edges'; +M.device = char(device); +M.description = strtrim(erase(header{1}, '#')); +M.source = strtrim(erase(header{2}, '# Source:')); + +end diff --git a/mhkit/wave/reference_models/RM3.mat b/mhkit/wave/reference_models/RM3.mat new file mode 100644 index 000000000..16e5819c6 Binary files /dev/null and b/mhkit/wave/reference_models/RM3.mat differ diff --git a/mhkit/wave/reference_models/RM5.mat b/mhkit/wave/reference_models/RM5.mat new file mode 100644 index 000000000..8f25f1857 Binary files /dev/null and b/mhkit/wave/reference_models/RM5.mat differ diff --git a/mhkit/wave/reference_models/RM6.mat b/mhkit/wave/reference_models/RM6.mat new file mode 100644 index 000000000..415483d77 Binary files /dev/null and b/mhkit/wave/reference_models/RM6.mat differ diff --git a/mhkit/wave/reference_models/reference_model_power_matrix.m b/mhkit/wave/reference_models/reference_model_power_matrix.m new file mode 100644 index 000000000..96dbfa73d --- /dev/null +++ b/mhkit/wave/reference_models/reference_model_power_matrix.m @@ -0,0 +1,57 @@ +function reference_model = reference_model_power_matrix(device) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Load the power matrix of a Reference Model wave energy converter +% +% The Reference Model Project, sponsored by the U.S. Department of +% Energy, developed open-source marine energy point designs as reference +% models to benchmark technology performance and costs, see +% https://openei.org/wiki/PRIMRE/Signature_Projects/Reference_Model. +% The power matrices are from the National Laboratory of the Rockies +% (NLR) System Advisor Model (SAM) wave energy converter library, +% https://github.com/NatLabRockies/SAM, and are stored in +% mhkit/wave/reference_models/.mat by +% scripts/generate_reference_model_power_matrices.m. +% +% Parameters +% ------------ +% device : string +% "RM3" (wave point absorber), "RM5" (oscillating surge flap), or +% "RM6" (oscillating water column) +% +% Returns +% --------- +% reference_model : structure +% reference_model.values : matrix +% Device power [W], Hs rows by Te columns +% reference_model.stat : char +% 'power' +% reference_model.x_bins, x_edges : vector +% Te bin centers and edges [s] +% reference_model.y_bins, y_edges : vector +% Hs bin centers and edges [m] +% reference_model.device, technology_type, pto_type, description : char +% Device metadata from the SAM library +% reference_model.characteristic_diameter : double +% Characteristic diameter [m] +% reference_model.mass : double +% Unballasted structural mass [kg] +% reference_model.data_source : char +% Where the matrix came from +% reference_model.reference : char +% Report describing the devices, Neary et al. (2014), SAND2014-9040 +% reference_model.last_modified : datetime +% When the .mat file was generated [UTC] +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments + device {mustBeTextScalar, mustBeMember(device, ["RM3", "RM5", "RM6"])} +end + +mat_file = fullfile(fileparts(mfilename('fullpath')), char(device) + ".mat"); +loaded = load(mat_file, 'reference_model'); +reference_model = loaded.reference_model; + +end diff --git a/mhkit/wave/resource/average_crest_period.m b/mhkit/wave/resource/average_crest_period.m index e642b8fab..b62138579 100644 --- a/mhkit/wave/resource/average_crest_period.m +++ b/mhkit/wave/resource/average_crest_period.m @@ -1,47 +1,199 @@ -function Tavg=average_crest_period(S,varargin) +function Tc = average_crest_period(S, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the average crest period % -% Parameters -% ------------ -% S: Spectral Density (m^2/Hz) +% Calculates wave average crest period from spectra % -% 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) +% Computed as sqrt(m2/m4) where m2 and m4 are spectral moments. % -% OR +% Parameters +% ------------ +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as second argument +% table: one variable per spectrum (rows), frequency required as second +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as second argument +% frequency_bins : vector (optional) +% Bin widths for frequency of S. Required for unevenly sized bins. % -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) +% Returns +% --------- +% Tc : double, column vector, table, or timetable [s] +% Average crest period, one value per spectrum. Matches the container +% style of S: numeric for struct/numeric input, table/timetable for +% table/timetable input. % -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 % -% S.frequency: frequency (Hz) +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 % -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 % -% Returns -% --------- -% Tavg: double -% Average crest Period (s) +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Tc = average_crest_period(S); +% Tc = +% 3.5816 % [s] % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> Tc = average_crest_period(spectrum, frequency); +% Tc = % [s] +% 3.5438 +% 3.5760 +% 3.5199 +% : +% 3.5226 +% 3.7637 +% 3.5816 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> Tc_table = average_crest_period(T, frequency); +% Tc_table = % [s] +% 17520×1 table +% average_crest_period +% ____________________ +% 3.5438 +% 3.576 +% 3.5199 +% : +% 3.5226 +% 3.7637 +% 3.5816 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> Tc_tt = average_crest_period(TT, frequency); +% Tc_tt = % [s] +% 17520×1 timetable +% time average_crest_period +% ____________________ ____________________ +% 01-Jan-2025 00:00:00 3.5438 +% 01-Jan-2025 00:30:00 3.576 +% 01-Jan-2025 01:00:00 3.5199 +% : +% 31-Dec-2025 22:30:00 3.5226 +% 31-Dec-2025 23:00:00 3.7637 +% 31-Dec-2025 23:30:00 3.5816 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 +% +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 +% +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Tc = average_crest_period(S); +% Tc = +% 2.0289 % [s] +% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> Tc = average_crest_period(spectrum, frequency); +% Tc = +% 2.0289 % [s] +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -if nargin == 2 - freq_bins = py.numpy.array(varargin{1}); -elseif nargin == 1 - freq_bins = py.None; +arguments (Input) + S +end + +arguments (Repeating) + varargin +end + +arguments (Output) + Tc +end + +[spectrum, frequency, time, input_style, remaining] = mhkit_standardize_spectrum_input(S, 'average_crest_period', varargin{:}); + +% Compute on primitive numeric arrays so the formula never touches a +% table/timetable object directly. +if isempty(remaining) + m2 = frequency_moment(spectrum, 2, frequency); + m4 = frequency_moment(spectrum, 4, frequency); else - ME = MException('MATLAB:average_crest_period','Incorrect number of input arguments'); - throw(ME); + m2 = frequency_moment(spectrum, 2, frequency, remaining{1}); + m4 = frequency_moment(spectrum, 4, frequency, remaining{1}); end -S_py = typecast_spectra_to_mhkit_python(S); +Tc = sqrt(m2 ./ m4); +Tc = Tc(:); +mhkit_verify_is_column_vector(Tc, 'function_name', mfilename); -Tavg = py.mhkit.wave.resource.average_crest_period(S_py, pyargs('frequency_bins', freq_bins)); +Tc = mhkit_restore_spectrum_output(Tc, input_style, 'average_crest_period', time); -Tavg = typecast_from_mhkit_python(Tavg).data; +end diff --git a/mhkit/wave/resource/average_wave_period.m b/mhkit/wave/resource/average_wave_period.m index 027f4348d..71d6c29bf 100644 --- a/mhkit/wave/resource/average_wave_period.m +++ b/mhkit/wave/resource/average_wave_period.m @@ -1,47 +1,197 @@ -function Tavg=average_wave_period(S,varargin) +function Tm = average_wave_period(S, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the average wave period +% +% Calculates mean wave period from spectra +% +% Computed as m0/m1 where m0 and m1 are spectral moments. % % 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) +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as second argument +% table: one variable per spectrum (rows), frequency required as second +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as second argument +% frequency_bins : vector (optional) +% Bin widths for frequency of S. Required for unevenly sized bins. % -% OR +% Returns +% --------- +% Tm : double, column vector, table, or timetable [s] +% Mean wave period, one value per spectrum. Matches the container style +% of S: numeric for struct/numeric input, table/timetable for +% table/timetable input. % -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 % -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 % -% S.frequency: frequency (Hz) +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 % -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Tm = average_wave_period(S); +% Tm = +% 7.0229 % [s] % -% Returns -% -------- -% Tavg: float -% Mean wave period (s) +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> Tm = average_wave_period(spectrum, frequency); +% Tm = % [s] +% 5.9716 +% 5.8869 +% 5.7310 +% : +% 6.9671 +% 7.4852 +% 7.0229 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> Tm_table = average_wave_period(T, frequency); +% Tm_table = % [s] +% 17520×1 table +% average_wave_period +% ___________________ +% 5.9716 +% 5.8869 +% 5.731 +% : +% 6.9671 +% 7.4852 +% 7.0229 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> Tm_tt = average_wave_period(TT, frequency); +% Tm_tt = % [s] +% 17520×1 timetable +% time average_wave_period +% ____________________ ___________________ +% 01-Jan-2025 00:00:00 5.9716 +% 01-Jan-2025 00:30:00 5.8869 +% 01-Jan-2025 01:00:00 5.731 +% : +% 31-Dec-2025 22:30:00 6.9671 +% 31-Dec-2025 23:00:00 7.4852 +% 31-Dec-2025 23:30:00 7.0229 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 +% +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Tm = average_wave_period(S); +% Tm = +% 5.6826 % [s] +% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> Tm = average_wave_period(spectrum, frequency); +% Tm = +% 5.6826 % [s] +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + S +end +arguments (Repeating) + varargin +end + +arguments (Output) + Tm +end -if nargin == 2 - freq_bins = py.numpy.array(varargin{1}); -elseif nargin == 1 - freq_bins = py.None; +[spectrum, frequency, time, input_style, remaining] = mhkit_standardize_spectrum_input(S, 'average_wave_period', varargin{:}); + +if isempty(remaining) + m0 = frequency_moment(spectrum, 0, frequency); + m1 = frequency_moment(spectrum, 1, frequency); else - ME = MException('MATLAB:average_wave_period','Incorrect number of input arguments'); - throw(ME); + m0 = frequency_moment(spectrum, 0, frequency, remaining{1}); + m1 = frequency_moment(spectrum, 1, frequency, remaining{1}); end -S_py = typecast_spectra_to_mhkit_python(S); +Tm = m0 ./ m1; +Tm = Tm(:); +mhkit_verify_is_column_vector(Tm, 'function_name', mfilename); -Tavg = py.mhkit.wave.resource.average_wave_period(S_py, pyargs('frequency_bins',freq_bins)); +Tm = mhkit_restore_spectrum_output(Tm, input_style, 'average_wave_period', time); -Tavg = typecast_from_mhkit_python(Tavg).data; +end diff --git a/mhkit/wave/resource/average_zero_crossing_period.m b/mhkit/wave/resource/average_zero_crossing_period.m index 8c44aeff2..6eafc8d8a 100644 --- a/mhkit/wave/resource/average_zero_crossing_period.m +++ b/mhkit/wave/resource/average_zero_crossing_period.m @@ -1,45 +1,199 @@ -function Tz=average_zero_crossing_period(S,varargin) +function Tz = average_zero_crossing_period(S, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates wave average zero crossing period from spectra +% +% Calculates wave average zero crossing period, T_z, [seconds] from spectra +% +% Computed as sqrt(m0/m2) where m0 and m2 are spectral moments +% (Eq 15 in IEC 62600-101 Ed. 2.0 en 2024). % % 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,x) +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as second argument +% table: one variable per spectrum (rows), frequency required as second +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as second argument +% frequency_bins : vector (optional) +% Bin widths for frequency of S. Required for unevenly sized bins. % -% OR +% Returns +% --------- +% Tz : double, column vector, table, or timetable [s] +% Average zero crossing period, one value per spectrum. Matches the +% container style of S: numeric for struct/numeric input, +% table/timetable for table/timetable input. % -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 % -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 % -% S.frequency: frequency (Hz) +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 % -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Tz = average_zero_crossing_period(S); +% Tz = +% 6.1209 % [s] % -% Returns -% --------- -% Tz: double -% Average Zero Crossing Period (s) +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> Tz = average_zero_crossing_period(spectrum, frequency); +% Tz = % [s] +% 5.4396 +% 5.3685 +% 5.2301 +% : +% 6.0647 +% 6.5259 +% 6.1209 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> Tz_table = average_zero_crossing_period(T, frequency); +% Tz_table = % [s] +% 17520×1 table +% average_zero_crossing_period +% ____________________________ +% 5.4396 +% 5.3685 +% 5.2301 +% : +% 6.0647 +% 6.5259 +% 6.1209 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> Tz_tt = average_zero_crossing_period(TT, frequency); +% Tz_tt = % [s] +% 17520×1 timetable +% time average_zero_crossing_period +% ____________________ ____________________________ +% 01-Jan-2025 00:00:00 5.4396 +% 01-Jan-2025 00:30:00 5.3685 +% 01-Jan-2025 01:00:00 5.2301 +% : +% 31-Dec-2025 22:30:00 6.0647 +% 31-Dec-2025 23:00:00 6.5259 +% 31-Dec-2025 23:30:00 6.1209 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 +% +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 +% +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Tz = average_zero_crossing_period(S); +% Tz = +% 5.0958 % [s] % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> Tz = average_zero_crossing_period(spectrum, frequency); +% Tz = +% 5.0958 % [s] +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + S +end + +arguments (Repeating) + varargin +end -if nargin == 2 - freq_bins = py.numpy.array(varargin{1}); -elseif nargin == 1 - freq_bins = py.None; +arguments (Output) + Tz +end + +[spectrum, frequency, time, input_style, remaining] = mhkit_standardize_spectrum_input(S, 'average_zero_crossing_period', varargin{:}); + +% Eq 15 in IEC 62600-101 Ed. 2.0 en 2024 +if isempty(remaining) + m0 = frequency_moment(spectrum, 0, frequency); + m2 = frequency_moment(spectrum, 2, frequency); else - ME = MException('MATLAB:average_zero_crossing_period','Incorrect number of input arguments'); - throw(ME); + m0 = frequency_moment(spectrum, 0, frequency, remaining{1}); + m2 = frequency_moment(spectrum, 2, frequency, remaining{1}); end -S_py = typecast_spectra_to_mhkit_python(S); +Tz = sqrt(m0 ./ m2); +Tz = Tz(:); +mhkit_verify_is_column_vector(Tz, 'function_name', mfilename); -Tz = py.mhkit.wave.resource.average_zero_crossing_period(S_py ,pyargs('frequency_bins',freq_bins)); -Tz = typecast_from_mhkit_python(Tz).data; +Tz = mhkit_restore_spectrum_output(Tz, input_style, 'average_zero_crossing_period', time); + +end diff --git a/mhkit/wave/resource/depth_regime.m b/mhkit/wave/resource/depth_regime.m index 91e96ab55..daca98c49 100644 --- a/mhkit/wave/resource/depth_regime.m +++ b/mhkit/wave/resource/depth_regime.m @@ -1,45 +1,45 @@ -function depth_reg=depth_regime(l,h,options) +function depth_reg = depth_regime(l, h, options) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates the depth regime based on wavelength and height -% Deep water: h/l > ratio -% This function exists so sinh in wave celerity doesn't blow -% up to infinity. % -% P.K. Kundu, I.M. Cohen (2000) suggest h/l >> 1 for deep water (pg 209) -% Same citation as above, they also suggest for 3% accuracy, h/l > 0.28 (pg 210) -% However, since this function allows multiple wavelengths, higher ratio -% numbers are more accurate across varying wavelengths. +% Calculates the depth regime based on wavelength and depth % -% Parameters -% ------------ -% l: vector -% wave length (m) -% -% h: double -% Water depth (m) +% Deep water: h/l > ratio +% This function exists so sinh in wave celerity doesn't blow +% up to infinity. % -% ratio: double or int (optional) -% If h/l > ratio, -% water depth will be set to deep. Default ratio = 2 -% to call: energy_flux(k,h,"ratio",ratio) +% P.K. Kundu, I.M. Cohen (2000) suggest h/l >> 1 for deep water (pg 209) +% Same citation, they also suggest for 3% accuracy, h/l > 0.28 (pg 210) +% However, since this function allows multiple wavelengths, higher ratio +% numbers are more accurate across varying wavelengths. % +% Parameters +% ------------ +% l : numeric [m] +% Wavelength (scalar, vector, or array) +% h : double [m] +% Water column depth +% ratio : double (optional) +% If h/l > ratio, water depth is classified as deep. Default = 2 % % Returns -% ------- -% depth_reg: boolean or boolean array -% Boolean True if deep water, False otherwise +% --------- +% depth_reg : logical +% Boolean True if deep water, False otherwise. +% Same shape as input l. % -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -arguments - l - h - options.ratio = 2; +arguments (Input) + l {mustBeNumeric} + h (1,1) {mustBeNumeric, mustBePositive} + options.ratio (1,1) {mustBeNumeric, mustBePositive} = 2 +end +arguments (Output) + depth_reg logical end -depth_reg=double(py.array.array('d',py.numpy.nditer... - (py.mhkit.wave.resource.depth_regime(py.numpy.array(l),h,pyargs('ratio',options.ratio))))); +depth_reg = (h ./ l) > options.ratio; +end diff --git a/mhkit/wave/resource/energy_flux.m b/mhkit/wave/resource/energy_flux.m index 542047891..7564af046 100644 --- a/mhkit/wave/resource/energy_flux.m +++ b/mhkit/wave/resource/energy_flux.m @@ -1,62 +1,181 @@ -function J=energy_flux(S,h,options) +function J = energy_flux(S, h, options) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % +% Calculates the omnidirectional wave energy flux of the spectra +% +% For deep water (deep=true), uses simplified formula (Eq 8 in +% IEC 62600-100 Ed. 2.0 en 2024): +% J = (rho * g^2 / (64*pi)) * Hm0^2 * Te +% +% For general depth (deep=false), uses Eq 9 in IEC 62600-101 Ed. 2.0 en 2024: +% J = rho * g * sum(Cg * S * df) % % Parameters % ------------ -% S: Spectral Density (m^2/Hz) -% 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 : struct +% Wave spectrum structure: +% S.spectrum : vector or matrix [m^2/Hz] +% Spectral density +% S.frequency : vector [Hz] +% Frequency +% S.time : datetime or duration (optional) +% One timestamp per spectrum/column +% Numeric, table, and timetable input are not supported here: h (water +% depth) occupies the positional slot the multi-container convention +% uses for an explicit frequency argument, so only struct input (which +% carries frequency inline) is unambiguous. +% h : double [m] +% Water depth +% deep : logical (optional) +% If true, use the deep water approximation. Default = false. +% When false, a depth check is run to check for shallow water. +% rho : double [kg/m^3] (optional) +% Water density. Default = 1025 kg/m^3 +% g : double [m/s^2] (optional) +% Gravitational acceleration. Default = 9.80665 m/s^2 +% ratio : double (optional) +% Only applied if deep=false. If h/l > ratio, +% water depth is set to deep. Default = 2 % -% S.frequency: frequency (Hz) +% Returns +% --------- +% J : double or column vector [W/m] +% Omni-directional wave energy flux, one value per spectrum. A table +% or timetable if S.time is present, matching mhkit_restore_spectrum_output. % -% h: double -% Water depth (m) +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity', 'metaWaterDepth'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> % cast to double: fzero (inside wave_number, called by energy_flux +% >> % below) requires double, but CDIP returns these as single. +% >> frequency = double(data.metadata.wave.waveFrequency); % [Hz], 64x1 double +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% >> % CDIP's own buoy depth reading, not a literature value. +% >> h = double(data.metadata.meta.metaWaterDepth); % [m], 84 % -% deep: logical (optional) -% If True use the deep water approximation. Default False. When -% False a depth check is run to check for shallow water. The ratio -% of the shallow water regime can be changed using the ratio -% keyword. -% to call: energy_flux(S,h,"deep",py.True) +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> J = energy_flux(S, h); +% J = +% 14495.52 % [W/m] % -% rho: double (optional) -% water density (kg/m^3) -% to call: energy_flux(S,h,"rho",rho) +% Struct (CDIP real-world data): full dataset, one value per spectrum +% >> S.spectrum = spectrum; +% >> J = energy_flux(S, h); +% J = % [W/m] +% 10152.76 +% 9548.60 +% 9148.13 +% : +% 16898.06 +% 20193.00 +% 14495.52 % -% g: double (optional) -% gravitational acceleration (m/s^2) -% to call: energy_flux(S,h,"g",g) +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 % -% ratio: double or int (optional) -% Only applied if depth=False. If h/l > ratio, -% water depth will be set to deep. Default ratio = 2. -% to call: energy_flux(S,h,"ratio",1.5) -% Returns -% ------- -% J: double -% Omni-directional wave energy flux (W/m) +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 % +% >> % The official WEC-Sim RM3 tutorial's hydrodynamic data was +% >> % computed by WAMIT at infinite depth (rm3.out: "Water depth: +% >> % infinite"), so this example uses the deep-water approximation. +% >> % h is unused when deep=true; Inf documents that assumption +% >> % rather than standing in for a real finite depth. +% Struct (WEC-Sim output): the single spectrum, deep-water approximation +% >> S.frequency = frequency; +% >> S.spectrum = spectrum; +% >> J = energy_flux(S, Inf, 'deep', true); +% J = +% 10741.91 % [W/m] % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + S struct + h (1,1) {mustBeNumeric} + options.deep (1,1) logical = false + options.rho (1,1) {mustBeNumeric, mustBePositive} = 1025 + options.g (1,1) {mustBeNumeric, mustBePositive} = 9.80665 + options.ratio (1,1) {mustBeNumeric, mustBePositive} = 2 +end + +arguments (Output) + J +end + +[spectrum, frequency, time, input_style] = mhkit_standardize_spectrum_input(S, 'energy_flux'); -arguments - S - h - options.deep = py.False; - options.rho = 1025; - options.g = 9.80665; - options.ratio = 2; +rho = options.rho; +g = options.g; +% Validate depth is positive when using general depth calculation +if ~options.deep && h <= 0 + error('MHKiT:energy_flux:InvalidInput', ... + 'h must be positive when deep=false'); end -S_py = typecast_spectra_to_mhkit_python(S); +if options.deep + % Eq 8 in IEC 62600-100 Ed. 2.0 en 2024 (deep water simplification) + Te = energy_period(spectrum, frequency); + Hm0 = significant_wave_height(spectrum, frequency); + + coeff = rho * (g^2) / (64 * pi); + J = coeff * (Hm0.^2) .* Te; +else + % Calculate wave number + k = wave_number(frequency, h, 'rho', rho, 'g', g); + + % Calculate wave celerity (group velocity) + Cg = wave_celerity(k, h, 'g', g, 'depth_check', true, 'ratio', options.ratio); + + % Calculate frequency bin widths + delta_f = diff(frequency); + delta_f = [frequency(2) - frequency(1); delta_f]; % Prepend first bin width -J = py.mhkit.wave.resource.energy_flux(S_py,h,pyargs('rho',options.rho,'g',options.g,... - 'deep',options.deep, 'ratio',options.ratio)); + % Eq 9 in IEC 62600-101 Ed. 2.0 en 2024 + Cg_values = Cg.values(:); -J = typecast_from_mhkit_python(J).data; + if isvector(spectrum) + J = rho * g * sum(Cg_values .* spectrum .* delta_f); + else + % Matrix case: sum along frequency dimension (first dimension) + J = rho * g * sum(Cg_values .* spectrum .* delta_f, 1); + end +end + +J = J(:); +mhkit_verify_is_column_vector(J, 'function_name', mfilename); +J = mhkit_restore_spectrum_output(J, input_style, 'energy_flux', time); + +end diff --git a/mhkit/wave/resource/energy_period.m b/mhkit/wave/resource/energy_period.m index c83f4cb2c..d61b19de9 100644 --- a/mhkit/wave/resource/energy_period.m +++ b/mhkit/wave/resource/energy_period.m @@ -2,29 +2,167 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Calculate wave energy period (Te) seconds from power spectral density (PSD) +% Calculates wave energy period (Te) [seconds] from power spectral density (PSD) % % 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 -% -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as second argument +% table: one variable per spectrum (rows), frequency required as second +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as second argument % % Returns % --------- -% Te: double -% Wave energy period [seconds] indexed by S.columns +% Te : double, column vector, table, or timetable [s] +% Wave energy period, one value per spectrum. Matches the container +% style of S: numeric for struct/numeric input, table/timetable for +% table/timetable input. +% +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 +% +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 +% +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 +% +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Te = energy_period(S); +% Te = +% 8.8095 % [s] +% +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> Te = energy_period(spectrum, frequency); +% Te = % [s] +% 6.9209 +% 6.8977 +% 6.7150 +% : +% 8.7120 +% 9.2423 +% 8.8095 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> Te_table = energy_period(T, frequency); +% Te_table = % [s] +% 17520×1 table +% energy_period +% _____________ +% 6.9209 +% 6.8977 +% 6.715 +% : +% 8.712 +% 9.2423 +% 8.8095 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> Te_tt = energy_period(TT, frequency); +% Te_tt = % [s] +% 17520×1 timetable +% time energy_period +% ____________________ _____________ +% 01-Jan-2025 00:00:00 6.9209 +% 01-Jan-2025 00:30:00 6.8977 +% 01-Jan-2025 01:00:00 6.715 +% : +% 31-Dec-2025 22:30:00 8.712 +% 31-Dec-2025 23:00:00 9.2423 +% 31-Dec-2025 23:30:00 8.8095 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 +% +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 +% +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Te = energy_period(S); +% Te = +% 6.8700 % [s] +% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> Te = energy_period(spectrum, frequency); +% Te = +% 6.8700 % [s] % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -arguments +arguments (Input) S end @@ -32,23 +170,25 @@ varargin end - % 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 - - % Calculate energy period - Te = m_neg1 ./ m0; +arguments (Output) + Te +end + +[spectrum, frequency, time, input_style, remaining] = mhkit_standardize_spectrum_input(S, 'energy_period', varargin{:}); + +% Compute on primitive numeric arrays only, following MHKiT-Python convention. +if isempty(remaining) + m0 = frequency_moment(spectrum, 0, frequency); + m_neg1 = frequency_moment(spectrum, -1, frequency); +else + m0 = frequency_moment(spectrum, 0, frequency, remaining{1}); + m_neg1 = frequency_moment(spectrum, -1, frequency, remaining{1}); +end + +Te = m_neg1 ./ m0; +Te = Te(:); +mhkit_verify_is_column_vector(Te, 'function_name', mfilename); + +Te = mhkit_restore_spectrum_output(Te, input_style, 'energy_period', time); end diff --git a/mhkit/wave/resource/energy_period_to_peak_period.m b/mhkit/wave/resource/energy_period_to_peak_period.m new file mode 100644 index 000000000..8881eb10d --- /dev/null +++ b/mhkit/wave/resource/energy_period_to_peak_period.m @@ -0,0 +1,39 @@ +function Tp = energy_period_to_peak_period(Te, gamma) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Convert spectral energy period (Te) to peak period (Tp) +% +% Uses ITTC approximation for JONSWAP Spectrum from "The Specialist +% Committee on Waves, Final Report and Recommendations to the 23rd ITTC", +% Proceedings of the 23rd ITTC - Volume 2, Table A4. +% +% Parameters +% ------------ +% Te : numeric [s] +% Spectral energy period (scalar, vector, or array) +% gamma : double +% Peak enhancement factor for JONSWAP spectrum +% +% Returns +% --------- +% Tp : numeric [s] +% Spectral peak period (same shape as Te) +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + Te {mustBeNumeric} + gamma (1,1) {mustBeNumeric, mustBePositive} +end + +arguments (Output) + Tp {mustBeNumeric} +end + +% ITTC approximation factor +factor = 0.8255 + 0.03852 * gamma - 0.005537 * gamma^2 + 0.0003154 * gamma^3; + +Tp = Te / factor; + +end diff --git a/mhkit/wave/resource/frequency_moment.m b/mhkit/wave/resource/frequency_moment.m index 85345a59e..b97bbd3a0 100644 --- a/mhkit/wave/resource/frequency_moment.m +++ b/mhkit/wave/resource/frequency_moment.m @@ -4,26 +4,169 @@ % % Calculates the Nth frequency moment of the spectrum % +% Computes: m_N = sum(f^N * S(f) * df) per Eq 8 in IEC 62600-101 Ed. 2.0 en 2024 +% % 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 +% ------------ +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as next argument +% table: one variable per spectrum (rows), frequency required as next +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as next argument +% N : integer +% Moment order (0 for 0th, 1 for 1st, -1 for inverse, etc.) +% frequency_bins : vector (optional) +% Bin widths for frequency of S. Required for unevenly sized bins. +% Passed as the next argument after frequency (numeric/table/timetable) +% or as the only extra argument (struct). % -% N: int -% Moment (0 for 0th, 1 for 1st ....) +% Returns +% --------- +% m : double, column vector, table, or timetable +% Nth frequency moment, one value per spectrum. Matches the container +% style of S: numeric for struct/numeric input, table/timetable for +% table/timetable input. % -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 % -% Returns -% ------- -% m: double -% Nth Frequency Moment indexed by S.columns +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 +% +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 +% +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> m = frequency_moment(S, 2); +% m = +% 0.0053 +% +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> m = frequency_moment(spectrum, 2, frequency); +% m = +% 0.0062 +% 0.0060 +% 0.0062 +% : +% 0.0064 +% 0.0062 +% 0.0053 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> m_table = frequency_moment(T, 2, frequency); +% m_table = +% 17520×1 table +% m2 +% _________ +% 0.0062259 +% 0.0060092 +% 0.0062446 +% : +% 0.0064063 +% 0.0061927 +% 0.0053146 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> m_tt = frequency_moment(TT, 2, frequency); +% m_tt = +% 17520×1 timetable +% time m2 +% ____________________ _________ +% 01-Jan-2025 00:00:00 0.0062259 +% 01-Jan-2025 00:30:00 0.0060092 +% 01-Jan-2025 01:00:00 0.0062446 +% : +% 31-Dec-2025 22:30:00 0.0064063 +% 31-Dec-2025 23:00:00 0.0061927 +% 31-Dec-2025 23:30:00 0.0053146 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 +% +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 +% +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> m = frequency_moment(S, 2); +% m = +% 0.0077 +% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> m = frequency_moment(spectrum, 2, frequency); +% m = +% 0.0077 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -36,24 +179,16 @@ varargin end - % 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 +arguments (Output) + m +end + + [spectrum, frequency, time, input_style, remaining] = ... + mhkit_standardize_spectrum_input(S, 'frequency_moment', varargin{:}); % Standardize frequency, spectrum, and frequency bins - if ~isempty(varargin) - [frequency, spectrum, freq_bins] = standardize_wave_spectra_frequency(frequency, spectrum, varargin{1}); + if ~isempty(remaining) + [frequency, spectrum, freq_bins] = standardize_wave_spectra_frequency(frequency, spectrum, remaining{1}); else [frequency, spectrum, freq_bins] = standardize_wave_spectra_frequency(frequency, spectrum); end @@ -64,5 +199,16 @@ % Calculate Nth moment: m_N = sum(f^N * S * df) m = sum((frequency.^N) .* spectrum .* freq_bins, 1); + m = m(:); + mhkit_verify_is_column_vector(m, 'function_name', mfilename); + + % Name follows MHKiT-Python's m.name = "m" + str(N), but stays a valid MATLAB + % identifier for negative N (dot-indexing t.m-1 would parse as subtraction). + if N < 0 + statistic_name = sprintf('m_neg%d', abs(N)); + else + statistic_name = sprintf('m%d', N); + end + m = mhkit_restore_spectrum_output(m, input_style, statistic_name, time); end diff --git a/mhkit/wave/resource/peak_period.m b/mhkit/wave/resource/peak_period.m index 1ab94c1fc..74dd61e27 100644 --- a/mhkit/wave/resource/peak_period.m +++ b/mhkit/wave/resource/peak_period.m @@ -1,34 +1,206 @@ -function Tp=peak_period(S) +function Tp = peak_period(S, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates wave energy period from spectra +% +% Calculates wave peak period [seconds] from spectra per Eq 14 in IEC 62600-101 Ed. 2.0 en 2024 % % 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) +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as second argument +% table: one variable per spectrum (rows), frequency required as second +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as second argument % -% OR +% Returns +% --------- +% Tp : double, column vector, table, or timetable [s] +% Wave peak period, one value per spectrum. Matches the container style +% of S: numeric for struct/numeric input, table/timetable for +% table/timetable input. % -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 % -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 % -% S.frequency: frequency (Hz) +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 % -% Returns -% --------- -% Tp float -% Wave Peak Period (s) +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Tp = peak_period(S); +% Tp = +% 14.2857 % [s] +% +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> Tp = peak_period(spectrum, frequency); +% Tp = % [s] +% 7.1429 +% 7.1429 +% 6.6667 +% : +% 13.3333 +% 12.5000 +% 14.2857 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> Tp_table = peak_period(T, frequency); +% Tp_table = % [s] +% 17520×1 table +% peak_period +% ___________ +% 7.1429 +% 7.1429 +% 6.6667 +% : +% 13.333 +% 12.5 +% 14.286 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> Tp_tt = peak_period(TT, frequency); +% Tp_tt = % [s] +% 17520×1 timetable +% time peak_period +% ____________________ ___________ +% 01-Jan-2025 00:00:00 7.1429 +% 01-Jan-2025 00:30:00 7.1429 +% 01-Jan-2025 01:00:00 6.6667 +% : +% 31-Dec-2025 22:30:00 13.333 +% 31-Dec-2025 23:00:00 12.5 +% 31-Dec-2025 23:30:00 14.286 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 % +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Tp = peak_period(S); +% Tp = +% 5.0000 % [s] +% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> Tp = peak_period(spectrum, frequency); +% Tp = +% 5.0000 % [s] +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + S +end + +arguments (Repeating) + varargin +end + +arguments (Output) + Tp +end + +[spectrum, frequency, time, input_style] = mhkit_standardize_spectrum_input(S, 'peak_period', varargin{:}); + +% Handle both vector and matrix spectra +if isvector(spectrum) + spectrum = spectrum(:); + if length(frequency) ~= length(spectrum) + error('MHKiT:peak_period:InvalidInput', ... + 'frequency length (%d) must match spectrum length (%d)', ... + length(frequency), length(spectrum)); + end + [~, idx] = max(spectrum); + fp = frequency(idx); +else + if length(frequency) ~= size(spectrum, 1) + error('MHKiT:peak_period:InvalidInput', ... + 'frequency length (%d) must match number of rows in spectrum (%d)', ... + length(frequency), size(spectrum, 1)); + end + % Matrix case: find max along first dimension (frequency) + [~, idx] = max(spectrum, [], 1); + fp = frequency(idx); +end + +Tp = 1 ./ fp; +Tp = Tp(:); +mhkit_verify_is_column_vector(Tp, 'function_name', mfilename); -S_py = typecast_spectra_to_mhkit_python(S); +Tp = mhkit_restore_spectrum_output(Tp, input_style, 'peak_period', time); -Tp = py.mhkit.wave.resource.peak_period(S_py); -Tp = typecast_from_mhkit_python(Tp).data; +end diff --git a/mhkit/wave/resource/significant_wave_height.m b/mhkit/wave/resource/significant_wave_height.m index 59d9734f6..785c321c0 100644 --- a/mhkit/wave/resource/significant_wave_height.m +++ b/mhkit/wave/resource/significant_wave_height.m @@ -2,29 +2,170 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Calculates significant wave height Hm0 from spectra +% Calculates significant wave height Hm0 from spectra per Eq 12 in +% IEC 62600-101 Ed. 2.0 en 2024 % % Parameters % ------------ -% 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 +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as second argument +% table: one variable per spectrum (rows), frequency required as second +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as second argument +% frequency_bins : vector (optional) +% Bin widths for frequency of S. Required for unevenly sized bins. % % Returns % --------- -% Hm0: double -% Significant wave height [m] index by S.columns +% Hm0 : double, column vector, table, or timetable [m] +% Significant wave height, one value per spectrum. Matches the +% container style of S: numeric for struct/numeric input, +% table/timetable for table/timetable input. +% +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 +% +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 +% +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 +% +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Hm0 = significant_wave_height(S); +% Hm0 = +% 1.7849 % [m] +% +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> Hm0 = significant_wave_height(spectrum, frequency); +% Hm0 = % [m] +% 1.7168 +% 1.6646 +% 1.6532 +% : +% 1.9416 +% 2.0542 +% 1.7849 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> Hm0_table = significant_wave_height(T, frequency); +% Hm0_table = % [m] +% 17520×1 table +% significant_wave_height +% _______________________ +% 1.7168 +% 1.6646 +% 1.6532 +% : +% 1.9416 +% 2.0542 +% 1.7849 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> Hm0_tt = significant_wave_height(TT, frequency); +% Hm0_tt = % [m] +% 17520×1 timetable +% time significant_wave_height +% ____________________ _______________________ +% 01-Jan-2025 00:00:00 1.7168 +% 01-Jan-2025 00:30:00 1.6646 +% 01-Jan-2025 01:00:00 1.6532 +% : +% 31-Dec-2025 22:30:00 1.9416 +% 31-Dec-2025 23:00:00 2.0542 +% 31-Dec-2025 23:30:00 1.7849 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 +% +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 +% +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> Hm0 = significant_wave_height(S); +% Hm0 = +% 1.7859 % [m] +% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> Hm0 = significant_wave_height(spectrum, frequency); +% Hm0 = +% 1.7859 % [m] % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -arguments +arguments (Input) S end @@ -32,21 +173,22 @@ varargin end - % Calculate zeroth moment m0 using frequency_moment function - if isstruct(S) - % Pass struct directly - m0 = frequency_moment(S, 0, 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{:}); - else - error('Input S must be either a struct with fields .spectrum and .frequency, or a numeric array'); - end - - % Calculate significant wave height Hm0 - Hm0 = 4 * sqrt(m0); +arguments (Output) + Hm0 +end + +[spectrum, frequency, time, input_style, frequency_bins] = mhkit_standardize_spectrum_input(S, 'significant_wave_height', varargin{:}); + +if isempty(frequency_bins) + m0 = frequency_moment(spectrum, 0, frequency); +else + m0 = frequency_moment(spectrum, 0, frequency, frequency_bins{1}); +end + +Hm0 = 4 * sqrt(m0); +Hm0 = Hm0(:); +mhkit_verify_is_column_vector(Hm0, 'function_name', mfilename); + +Hm0 = mhkit_restore_spectrum_output(Hm0, input_style, 'significant_wave_height', time); end diff --git a/mhkit/wave/resource/spectral_bandwidth.m b/mhkit/wave/resource/spectral_bandwidth.m index e8753ce48..8b05c901c 100644 --- a/mhkit/wave/resource/spectral_bandwidth.m +++ b/mhkit/wave/resource/spectral_bandwidth.m @@ -1,47 +1,204 @@ -function e=spectral_bandwidth(S,varargin) +function e = spectral_bandwidth(S, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates bandwidth from spectra +% +% Calculates spectral bandwidth from spectra +% +% Computed as sqrt(1 - m2^2 / (m0 * m4)) where m0, m2, and m4 are +% spectral moments % % 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) +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as second argument +% table: one variable per spectrum (rows), frequency required as second +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as second argument +% frequency_bins : vector (optional) +% Bin widths for frequency of S. Required for unevenly sized bins. % -% OR +% Returns +% --------- +% e : double, column vector, table, or timetable [dimensionless] +% Spectral bandwidth, one value per spectrum. Matches the container +% style of S: numeric for struct/numeric input, table/timetable for +% table/timetable input. % -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 % -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 % -% S.frequency: frequency (Hz) +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 % -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> e = spectral_bandwidth(S); +% e = +% 0.8109 % [dimensionless] % -% Returns -% --------- -% e: double -% Spectral BandWidth +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> e = spectral_bandwidth(spectrum, frequency); +% e = % [dimensionless] +% 0.7587 +% 0.7459 +% 0.7396 +% : +% 0.8140 +% 0.8169 +% 0.8109 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> e_table = spectral_bandwidth(T, frequency); +% e_table = % [dimensionless] +% 17520×1 table +% spectral_bandwidth +% __________________ +% 0.75866 +% 0.74585 +% 0.73964 +% : +% 0.81402 +% 0.81694 +% 0.81093 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> e_tt = spectral_bandwidth(TT, frequency); +% e_tt = % [dimensionless] +% 17520×1 timetable +% time spectral_bandwidth +% ____________________ __________________ +% 01-Jan-2025 00:00:00 0.75866 +% 01-Jan-2025 00:30:00 0.74585 +% 01-Jan-2025 01:00:00 0.73964 +% : +% 31-Dec-2025 22:30:00 0.81402 +% 31-Dec-2025 23:00:00 0.81694 +% 31-Dec-2025 23:30:00 0.81093 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 +% +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 % +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> e = spectral_bandwidth(S); +% e = +% 0.9173 % [dimensionless] % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> e = spectral_bandwidth(spectrum, frequency); +% e = +% 0.9173 % [dimensionless] +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + S +end + +arguments (Repeating) + varargin +end + +arguments (Output) + e +end -if nargin == 2 - freq_bins = py.numpy.array(varargin{1}); -elseif nargin == 1 - freq_bins = py.None; +[spectrum, frequency, time, input_style, remaining] = mhkit_standardize_spectrum_input(S, 'spectral_bandwidth', varargin{:}); + +% Moments are computed on primitive numeric arrays, never on S directly. +if isempty(remaining) + m0 = frequency_moment(spectrum, 0, frequency); + m2 = frequency_moment(spectrum, 2, frequency); + m4 = frequency_moment(spectrum, 4, frequency); else - ME = MException('MATLAB:spectral_bandwidth','Incorrect number of input arguments'); - throw(ME); + m0 = frequency_moment(spectrum, 0, frequency, remaining{1}); + m2 = frequency_moment(spectrum, 2, frequency, remaining{1}); + m4 = frequency_moment(spectrum, 4, frequency, remaining{1}); end -S_py = typecast_spectra_to_mhkit_python(S); +% Equation 28 from: https://coastalmonitoring.org/ccoresources/waveparameterhandbook/ +% MHKiT-Python possibly has a bug here and divides m0 by m4, which returns ~1 for every spectrum +% The equation cited above. multiplies m0 by m4 +e = sqrt(1 - (m2.^2) ./ (m0 .* m4)); +e = e(:); +mhkit_verify_is_column_vector(e, 'function_name', mfilename); -e0 = py.mhkit.wave.resource.spectral_bandwidth(S_py,pyargs('frequency_bins',freq_bins)); +e = mhkit_restore_spectrum_output(e, input_style, 'spectral_bandwidth', time); -e = typecast_from_mhkit_python(e0).data; +end diff --git a/mhkit/wave/resource/spectral_width.m b/mhkit/wave/resource/spectral_width.m index fb20d6fd6..36c2f7a00 100644 --- a/mhkit/wave/resource/spectral_width.m +++ b/mhkit/wave/resource/spectral_width.m @@ -1,45 +1,201 @@ -function v = spectral_width(S,varargin) +function v = spectral_width(S, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % Calculates wave spectral width from spectra % +% Computed as sqrt(m0 * m_{-2} / m_{-1}^2 - 1) where m0, m_{-1}, and m_{-2} +% are spectral moments (Eq 16 in IEC 62600-101 Ed. 2.0 en 2024). +% % 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) +% S : struct, numeric, table, or timetable +% struct: S.spectrum (vector or matrix [m^2/Hz]), S.frequency (vector +% [Hz]), optional S.time +% numeric: spectral density array, frequency required as second argument +% table: one variable per spectrum (rows), frequency required as second +% argument, optional 'time' variable +% timetable: one variable per spectrum (rows = RowTimes), frequency +% required as second argument +% frequency_bins : vector (optional) +% Bin widths for frequency of S. Required for unevenly sized bins. % -% OR +% Returns +% --------- +% v : double, column vector, table, or timetable [dimensionless] +% Spectral width, one value per spectrum. Matches the container style +% of S: numeric for struct/numeric input, table/timetable for +% table/timetable input. % -% structure of form: -% S.spectrum: Spectral Density (m^2/Hz) +% Examples +% -------- +% CDIP Example: Station Number 225: Kaneohe Bay, WETS, Oahu, HI +% https://cdip.ucsd.edu/themes/cdip?pb=1&u2=s:225:st:1&d2=p70 +% >> station_number = '225'; +% >> data_type = 'realtime'; +% >> years = 2025; +% >> parameters = {'waveEnergyDensity'}; +% >> data = cdip_request_parse_workflow('station_number', station_number, 'data_type', data_type, 'years', years, 'parameters', parameters); +% >> frequency = data.metadata.wave.waveFrequency; % [Hz], 64x1 single +% frequency = +% 0.0250 +% 0.0300 +% 0.0350 +% : +% 0.5600 +% 0.5700 +% 0.5800 % -% S.type: String of the spectra type, i.e. Bretschneider, -% time series, date stamp etc. +% >> % waveEnergyDensity is stored [time x frequency]; transpose to +% >> % match MHKiT's frequency-as-rows, spectra-as-columns convention. +% >> spectrum = data.data.wave2D.waveEnergyDensity'; % [m^2/Hz], 64x17520 double +% spectrum = +% 0.0002 0.0001 0.0002 ... 0.0006 0.0003 0.0004 +% 0.0007 0.0003 0.0005 ... 0.0021 0.0007 0.0010 +% 0.0032 0.0030 0.0020 ... 0.0047 0.0029 0.0038 +% : +% 0.0165 0.0122 0.0138 ... 0.0186 0.0079 0.0128 +% 0.0124 0.0116 0.0080 ... 0.0128 0.0120 0.0088 +% 0.0184 0.0097 0.0116 ... 0.0095 0.0084 0.0080 % -% S.frequency: frequency (Hz) +% >> time = data.data.wave.waveTime; % 17520x1 datetime +% time = +% 01-Jan-2025 00:00:00 +% 01-Jan-2025 00:30:00 +% 01-Jan-2025 01:00:00 +% : +% 31-Dec-2025 22:30:00 +% 31-Dec-2025 23:00:00 +% 31-Dec-2025 23:30:00 % -% frequency_bins: vector (optional) -% Bin widths for frequency of S. Required for unevenly sized bins +% Struct (CDIP real-world data): single (most recent) spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> v = spectral_width(S); +% v = +% 0.4207 % [dimensionless] % -% Returns -% --------- -% e0: float -% Spectral Width +% Numeric (CDIP real-world data): matrix, one spectrum per column +% >> v = spectral_width(spectrum, frequency); +% v = % [dimensionless] +% 0.3669 +% 0.3937 +% 0.3947 +% : +% 0.4155 +% 0.3892 +% 0.4207 +% +% Table (CDIP real-world data): one row per spectrum, one variable per frequency bin +% >> col_names = mhkit_frequency_to_column_names(frequency); +% >> T = array2table(spectrum', 'VariableNames', col_names); +% >> v_table = spectral_width(T, frequency); +% v_table = % [dimensionless] +% 17520×1 table +% spectral_width +% ______________ +% 0.36695 +% 0.39367 +% 0.3947 +% : +% 0.41548 +% 0.3892 +% 0.42073 +% +% Timetable (CDIP real-world data): RowTimes carried through to the output +% >> TT = array2timetable(spectrum', 'RowTimes', time, 'VariableNames', col_names); +% >> v_tt = spectral_width(TT, frequency); +% v_tt = % [dimensionless] +% 17520×1 timetable +% time spectral_width +% ____________________ ______________ +% 01-Jan-2025 00:00:00 0.36695 +% 01-Jan-2025 00:30:00 0.39367 +% 01-Jan-2025 01:00:00 0.3947 +% : +% 31-Dec-2025 22:30:00 0.41548 +% 31-Dec-2025 23:00:00 0.3892 +% 31-Dec-2025 23:30:00 0.42073 +% +% WEC-Sim Output Example +% >> S = load('examples/data/RM3MooringMatrix_matlabWorkspace.mat', 'output'); +% >> elevation = S.output.wave.elevation; % [m], RM3 float, 40001x1 double +% >> raw_time = S.output.wave.time; % [s], 40001x1 double +% >> sample_rate = 1 / (raw_time(2) - raw_time(1)); % [Hz], 100 +% >> % Note: IEC 62600-101 Ed. 2.0 en 2024, "Wave energy resource +% >> % assessment and characterization", section 6.5.2 specifies a +% >> % wave record length of minimum 1200 s (20 min), up to 3600 s +% >> % (60 min) for better spectral resolution/precision. This +% >> % WEC-Sim record is only ~6.7 min, so treat this as illustrative, +% >> % not a fully representative sea-state estimate. +% >> % +% >> % Note: spectral statistics assume an irregular wave record. +% >> % They are not meant for regular, single-frequency (sine) wave +% >> % tests, which WEC-Sim is often used to run. +% >> Sxx = elevation_spectrum(elevation, sample_rate, 1000, raw_time); +% >> frequency = Sxx.frequency; % [Hz], 501x1 double +% frequency = +% 0.0000 +% 0.1000 +% 0.2000 +% : +% 49.8000 +% 49.9000 +% 50.0000 +% +% >> spectrum = Sxx.spectrum; % [m^2/Hz], 501x1 double +% spectrum = +% 0.1439 +% 0.8485 +% 0.8739 +% : +% 0.0000 +% 0.0000 +% 0.0000 +% +% Struct (WEC-Sim output): the single spectrum +% >> S.frequency = frequency; +% >> S.spectrum = spectrum(:,end); +% >> v = spectral_width(S); +% v = +% 0.4031 % [dimensionless] % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Numeric (WEC-Sim output): matrix, one spectrum per column +% >> v = spectral_width(spectrum, frequency); +% v = +% 0.4031 % [dimensionless] +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + S +end + +arguments (Repeating) + varargin +end -if nargin == 2 - freq_bins = py.numpy.array(varargin{1}); -elseif nargin == 1 - freq_bins = py.None; +arguments (Output) + v +end + +[spectrum, frequency, time, input_style, remaining] = mhkit_standardize_spectrum_input(S, 'spectral_width', varargin{:}); + +% Eq 16 in IEC 62600-101 Ed. 2.0 en 2024 +if isempty(remaining) + m0 = frequency_moment(spectrum, 0, frequency); + mn1 = frequency_moment(spectrum, -1, frequency); + mn2 = frequency_moment(spectrum, -2, frequency); else - ME = MException('MATLAB:spectral_width','Incorrect number of input arguments'); - throw(ME); + m0 = frequency_moment(spectrum, 0, frequency, remaining{1}); + mn1 = frequency_moment(spectrum, -1, frequency, remaining{1}); + mn2 = frequency_moment(spectrum, -2, frequency, remaining{1}); end -S_py = typecast_spectra_to_mhkit_python(S); +v = sqrt(((m0 .* mn2) ./ (mn1.^2)) - 1); +v = v(:); +mhkit_verify_is_column_vector(v, 'function_name', mfilename); -v = py.mhkit.wave.resource.spectral_width(S_py,pyargs('frequency_bins',freq_bins)); -v = typecast_from_mhkit_python(v).data; +v = mhkit_restore_spectrum_output(v, input_style, 'spectral_width', time); + +end diff --git a/mhkit/wave/resource/standardize_wave_spectra_frequency.m b/mhkit/wave/resource/standardize_wave_spectra_frequency.m index e5e3236b9..10a0ba28e 100644 --- a/mhkit/wave/resource/standardize_wave_spectra_frequency.m +++ b/mhkit/wave/resource/standardize_wave_spectra_frequency.m @@ -1,7 +1,7 @@ function [frequency, spectrum, freq_bins] = standardize_wave_spectra_frequency(frequency, spectrum, varargin) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Standardize frequency and spectrum data for wave calculations +% Standardize frequency and spectrum data for spectral calculations % % This utility function standardizes frequency vectors, spectra, and frequency bins % following MHKiT-Python conventions for consistent numerical integration. @@ -35,19 +35,6 @@ varargin end - % Handle optional frequency bins input - if ~isempty(varargin) - freq_bins = varargin{1}; - else - % Calculate individual frequency bin widths to match MHKiT-Python implementation - % MHKiT-Python uses: delta_f = f.diff() then prepends first difference - % This creates a vector where each frequency has its own bin width, - % which is critical for accurate numerical integration when frequencies - % are not perfectly uniform (using mean(df) introduces systematic error) - df = diff(frequency); - freq_bins = [df(1); df(:)]; % Prepend first difference, ensure column vector - end - % Ensure proper shapes frequency = frequency(:); if isvector(spectrum) @@ -56,25 +43,36 @@ % Check size consistency if length(frequency) ~= size(spectrum,1) - error('Length of frequency vector must match number of rows in spectrum'); + error('MHKiT:standardize_wave_spectra_frequency:InvalidInput', ... + 'Length of frequency vector must match number of rows in spectrum'); end - % Filter out zero and near-zero frequencies to avoid division by zero - % Following MHKiT-Python implementation: omit frequencies <= 1e-12 + % Omit near-zero frequencies, following MHKiT-Python convention: + % https://github.com/MHKiT-Software/MHKiT-Python/blob/6bad8fe4f2bd8a9bff66fb9607ed0900f09d0258/mhkit/wave/resource.py#L455-L457 valid_idx = frequency > 1e-12; frequency = frequency(valid_idx); spectrum = spectrum(valid_idx, :); - - % Handle frequency bins - filter if vector, keep scalar as-is - if ~isscalar(freq_bins) - freq_bins = freq_bins(:); - if length(freq_bins) ~= length(valid_idx) - error('Length of freq_bins must match original frequency vector'); + + % Handle optional frequency bins input + if ~isempty(varargin) + freq_bins = varargin{1}; + if ~isscalar(freq_bins) + freq_bins = freq_bins(:); + if length(freq_bins) ~= length(valid_idx) + error('MHKiT:standardize_wave_spectra_frequency:InvalidInput', ... + 'Length of freq_bins must match original frequency vector'); + end + freq_bins = freq_bins(valid_idx); + else + % Ensure scalar freq_bins is column vector for consistency + freq_bins = freq_bins(:); end - freq_bins = freq_bins(valid_idx); else - % Ensure scalar freq_bins is column vector for consistency - freq_bins = freq_bins(:); + % Bin widths computed after filtering, following MHKiT-Python convention: + % https://github.com/MHKiT-Software/MHKiT-Python/blob/6bad8fe4f2bd8a9bff66fb9607ed0900f09d0258/mhkit/wave/resource.py#L460-L464 + % Computing them before filtering gives the wrong first bin width for non-uniform spacing. + df = diff(frequency); + freq_bins = [df(1); df(:)]; % Prepend first difference, ensure column vector end end diff --git a/mhkit/wave/resource/surface_elevation.m b/mhkit/wave/resource/surface_elevation.m index eb78060da..6f57a6299 100644 --- a/mhkit/wave/resource/surface_elevation.m +++ b/mhkit/wave/resource/surface_elevation.m @@ -6,75 +6,89 @@ % % Parameters % ------------ -% S: structure -% Spectral data with fields: -% S.spectrum: (n_freq x 1) spectral density [m^2/Hz] -% S.frequency: (n_freq x 1) frequency vector [Hz] -% -% time_index: vector -% Time used to create the wave elevation time-series [s], -% for example, time = 0:0.01:100 -% -% options: structure (optional) -% Optional fields: -% seed: Random seed (default = 123) -% frequency_bins: Bin widths for frequency of S. Required for unevenly sized bins -% phases: Explicit phases for frequency components (overrides seed) -% for example, phases = rand(length(S.frequency),1) * 2 * pi -% method: Method used to calculate the surface elevation. 'ifft' -% (Inverse Fast Fourier Transform) used by default if the -% given frequency_bins==[] or is evenly spaced. -% 'sum_of_sines' explicitly sums each frequency component -% and used by default if uneven frequency_bins are provided. -% The 'ifft' method is significantly faster. +% S : struct +% Wave spectrum structure: +% S.spectrum : vector [m^2/Hz] +% Spectral density +% S.frequency : vector [Hz] +% Frequency +% time_index : vector [s] +% Time used to create the wave elevation time-series, +% for example, time_index = 0:0.01:100 +% seed : double (optional) +% Random seed. Default = [] (unseeded, non-reproducible) +% frequency_bins : vector [Hz] (optional) +% Bin widths for frequency of S. Required for unevenly sized bins. +% phases : vector [rad] (optional) +% Explicit phases for frequency components (overrides seed), +% for example, phases = rand(length(S.frequency), 1) * 2 * pi +% method : char (optional) +% Method used to calculate the surface elevation. 'ifft' (Inverse +% Fast Fourier Transform) used by default if the given frequency_bins +% is empty or evenly spaced. 'sum_of_sines' explicitly sums each +% frequency component and is used by default if uneven frequency_bins +% are provided. The 'ifft' method is significantly faster. +% Default = 'ifft' % % Returns % --------- -% wave_elevation: structure -% Generated wave elevation with fields: -% elevation: Wave surface elevation [m] -% time: time vector [s] -% type: description string +% wave_elevation : struct +% wave_elevation.elevation : vector [m] +% Wave surface elevation +% wave_elevation.time : vector [s] +% Time vector +% wave_elevation.type : char +% 'Time Series from Spectra' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -arguments +arguments (Input) S struct time_index (:,1) double - options.seed (1,1) double = 123 + options.seed {mustBeNumeric} = [] options.frequency_bins = [] options.phases = [] options.method char {mustBeMember(options.method, {'ifft','sum_of_sines'})} = 'ifft' end +arguments (Output) + wave_elevation struct +end + +if ~isempty(options.seed) && ~isscalar(options.seed) + error('MHKiT:surface_elevation:InvalidInput', 'seed must be a scalar or empty (unseeded).'); +end + % 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) - % Calculate individual frequency bin widths to match MHKiT-Python implementation - % MHKiT-Python uses: delta_f = f.diff() then prepends first difference - % This creates a vector where each frequency has its own bin width, - % which is critical for accurate numerical integration when frequencies - % are not perfectly uniform (using mean(df) introduces systematic error) + delta_f = f(2) - f(1); df = diff(f); - delta_f = [df(1); df(:)]; % Prepend first difference, ensure column vector df_uniform = all(abs(df - df(1)) < 1e-8); else - delta_f = options.frequency_bins(:); - if length(delta_f) ~= length(f) + freq_bins = options.frequency_bins(:); + if length(freq_bins) ~= 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); + df_uniform = all(abs(freq_bins - freq_bins(1)) < 1e-8); + if df_uniform + delta_f = freq_bins(1); + else + delta_f = freq_bins; + end end -% Handle phases +% An empty seed is left unseeded (non-reproducible), +% following MHKiT-Python convention. +% https://github.com/MHKiT-Software/MHKiT-Python/blob/6bad8fe4f2bd8a9bff66fb9607ed0900f09d0258/mhkit/wave/resource.py#L372-L373 if isempty(options.phases) - rng(options.seed); + if ~isempty(options.seed) + rng(options.seed); + end phase = 2*pi*rand(Nf, 1); else phase = options.phases(:); @@ -83,31 +97,37 @@ end end -% Choose method +method = options.method; +if strcmp(method, 'ifft') + if f(1) ~= 0 + warning('MHKiT:surface_elevation:MethodFallback', ... + ['ifft method must have zero frequency defined. Setting ' ... + 'method to less efficient sum_of_sines method.']); + method = 'sum_of_sines'; + end + if ~df_uniform + warning('MHKiT:surface_elevation:MethodFallback', ... + ['ifft method must have evenly spaced frequency bins. ' ... + 'Setting method to less efficient sum_of_sines method.']); + method = 'sum_of_sines'; + end +end + 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 - % Match MHKiT-Python implementation: A_cmplx = A * (cos(phase) + 1j * sin(phase)) - A_complex = A .* (cos(phase) + 1i * sin(phase)); - % Match MHKiT-Python scaling: 0.5 * A_cmplx * time_index.size - A_scaled = 0.5 * A_complex * length(time_index); - % Use MATLAB's equivalent of irfft - eta = real(ifft(A_scaled, length(time_index), 'symmetric')); - 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; +if strcmp(method, 'ifft') + A_complex = A .* (cos(phase) + 1i * sin(phase)); + A_scaled = 0.5 * A_complex * length(time_index); + % Use MATLAB's equivalent of irfft + eta = real(ifft(A_scaled, length(time_index), 'symmetric')); + wave_elevation.elevation = eta; +elseif strcmp(method, 'sum_of_sines') + B = omega .* time_index'; + B = B'; % (Nt x Nf) + C = cos(B + phase'); + eta = C * A; + wave_elevation.elevation = eta; end wave_elevation.time = time_index; diff --git a/mhkit/wave/resource/wave_celerity.m b/mhkit/wave/resource/wave_celerity.m index d689da6ad..cddda6c89 100644 --- a/mhkit/wave/resource/wave_celerity.m +++ b/mhkit/wave/resource/wave_celerity.m @@ -1,71 +1,102 @@ -function Cg=wave_celerity(k,h,options) +function Cg = wave_celerity(k, h, options) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates wave celerity (group velocity) -% -% Parameters -% ------------ -% k: wave number (1/m) -% structure of form: -% k.values= wave number -% -% k.frequency= frequency (Hz) % -% h: double -% Water depth (m) +% Calculates wave celerity (group velocity) % -% g: double (optional) -% gravitational acceleration (m/s^2) -% to call: energy_flux(k,h,"g",g) +% Uses the formula from Eq 10 in IEC 62600-101 Ed. 2.0 en 2024: +% Cg = (pi * f / k) * (1 + (2*h*k) / sinh(2*h*k)) % -% depth_check: bool (optional) -% If True check depth regime. Default False. -% to call: energy_flux(k,h,"depth_check",py.True) +% For deep water (when depth_check=true and h/l > ratio), uses the +% simplified formula: Cg = pi * f / k % -% ratio: double or int (optional) -% Only applied if depth_check=True. If h/l > ratio, -% water depth will be set to deep. Default ratio = 2 -% to call: energy_flux(k,h,"ratio",ratio) +% Parameters +% ------------ +% k : struct +% Wave number structure: +% k.values : vector [1/m] +% Wave number values +% k.frequency : vector [Hz] +% Frequency +% h : double [m] +% Water depth +% g : double [m/s^2] (optional) +% Gravitational acceleration. Default = 9.80665 m/s^2 +% depth_check : logical (optional) +% If true, check depth regime and use deep water approximation +% where applicable. Default = false +% ratio : double (optional) +% Only applied if depth_check=true. If h/l > ratio, +% water depth is set to deep. Default = 2 % % Returns -% ------- -% Cg: structure -% -% -% Cg.values: water celerity -% -% Cg.frequency [Hz] -% -% Cg.h: height [m] +% --------- +% Cg : struct +% Cg.values : vector [m/s] +% Wave celerity (group velocity) +% Cg.frequency : vector [Hz] +% Frequency +% Cg.h : double [m] +% Water depth % -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + k struct + h (1,1) {mustBeNumeric, mustBePositive} + options.g (1,1) {mustBeNumeric, mustBePositive} = 9.80665 + options.depth_check (1,1) logical = false + options.ratio (1,1) {mustBeNumeric, mustBePositive} = 2 +end -arguments - k - h - options.g = 9.80665; - options.depth_check = py.False; - options.ratio = 2; +arguments (Output) + Cg struct +end +% Validate input structure +if ~isfield(k, 'values') || ~isfield(k, 'frequency') + error('MHKiT:wave_celerity:InvalidInput', ... + 'k must be a structure with values and frequency fields'); end -if (isa(k,'py.pandas.core.frame.DataFrame')~=1) - if (isstruct(k)==1) +k_values = k.values(:); +f = k.frequency(:); + +if options.depth_check + % Calculate wavelength + l = wave_length(k_values); - k=py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas(k.frequency,py.numpy.array(k.values),1); + % Get depth regime (true = deep water) + dr = depth_regime(l, h, 'ratio', options.ratio); - else - ME = MException('MATLAB:wave_celerity','k needs to be a structure or Pandas dataframe, use py.mhkit_python_utils.pandas_dataframe.spectra_to_pandas to create one'); - throw(ME); + % Initialize output + Cg_values = zeros(size(k_values)); + + % Deep water approximation for deep frequencies + if any(dr) + Cg_values(dr) = pi * f(dr) ./ k_values(dr); end -end + % Full formula for shallow/intermediate frequencies + if any(~dr) + sf = f(~dr); + sk = k_values(~dr); + Cg_values(~dr) = (pi * sf ./ sk) .* (1 + (2 * h * sk) ./ sinh(2 * h * sk)); + end +else + % Eq 10 in IEC 62600-101 Ed. 2.0 en 2024 + Cg_values = (pi * f ./ k_values) .* (1 + (2 * h * k_values) ./ sinh(2 * h * k_values)); +end -Cgdf=py.mhkit.wave.resource.wave_celerity(k,h,pyargs('g',options.g,... - 'depth_check',options.depth_check,'ratio',options.ratio)); +% Handle shape to match input +if isrow(k.values) + Cg_values = Cg_values'; +end -Cg.values=double(py.array.array('d',py.numpy.nditer(Cgdf.values))); -Cg.frequency=double(py.array.array('d',py.numpy.nditer(Cgdf.index))); -Cg.h=h; +% Build output structure +Cg = struct(); +Cg.values = Cg_values; +Cg.frequency = k.frequency; +Cg.h = h; +end diff --git a/mhkit/wave/resource/wave_length.m b/mhkit/wave/resource/wave_length.m index b2320b0e1..619297eaf 100644 --- a/mhkit/wave/resource/wave_length.m +++ b/mhkit/wave/resource/wave_length.m @@ -1,21 +1,48 @@ -function l=wave_length(k) +function l = wave_length(k) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % Calculates wave length from wave number -% To compute: 2*pi/wavenumber +% +% Computes: l = 2*pi/k % % Parameters % ------------ -% k: wave number (1/m) -% intiger, double, or vector +% k : numeric (scalar, vector, matrix, or struct) +% Wave number [1/m]. If struct, uses k.values field. % % Returns -% ------- -% l: double or array -% Wave length [m] indexed by frequency +% --------- +% l : numeric [m] +% Wave length. Same type/shape as input k (or k.values if struct). % -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + k +end + +arguments (Output) + l {mustBeNumeric} +end + +% Handle struct input (k.values) +if isstruct(k) + if ~isfield(k, 'values') + error('MHKiT:wave_length:InvalidInput', ... + 'Structure input must have a values field'); + end + k_values = k.values; +else + k_values = k; +end + +% Validate input +if ~isnumeric(k_values) + error('MHKiT:wave_length:InvalidInput', ... + 'k must be numeric or a struct with numeric values field'); +end -l=double(py.array.array('d',py.numpy.nditer(py.mhkit.wave.resource.wave_length(py.numpy.array(k))))); +l = 2 * pi ./ k_values; +end diff --git a/mhkit/wave/resource/wave_number.m b/mhkit/wave/resource/wave_number.m index 536468e03..011d302d1 100644 --- a/mhkit/wave/resource/wave_number.m +++ b/mhkit/wave/resource/wave_number.m @@ -1,53 +1,89 @@ -function k=wave_number(f,h,options) +function k = wave_number(f, h, options) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Calculates wave number % -% Parameters -% ------------ -% f: frequency (Hz) -% vector or numpy array -% h: float -% Water depth (m) -% -% rho: float (optional) -% water density (kg/m^3) -% to call: wave_number(f,h,"rho",rho) -% -% g: float (optional) -% gravitational acceleration (m/s^2) -% to call: wave_number(f,h,"g",g) +% Calculates wave number from frequency and water depth % +% Solves the linear dispersion relation (Eq 11 in IEC 62600-101 Ed. 2.0 en 2024) % +% To compute wave number from angular frequency (w), convert w to f before +% using this function (f = w / (2*pi)) % +% Parameters +% ------------ +% f : vector or scalar [Hz] +% Frequency +% h : double [m] +% Water depth +% rho : double [kg/m^3] (optional) +% Water density. Default = 1025 kg/m^3 +% g : double [m/s^2] (optional) +% Gravitational acceleration. Default = 9.80665 m/s^2 % % Returns -% ------- -% k: structure -% -% -% k.values: wave number -% -% k.frequency: frequency [Hz] -% +% --------- +% k : struct +% k.values : vector [1/m] +% Wave number +% k.frequency : vector [Hz] +% Frequency % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +arguments (Input) + f {mustBeNumeric} + h (1,1) {mustBeNumeric, mustBePositive} + options.rho (1,1) {mustBeNumeric, mustBePositive} = 1025 + options.g (1,1) {mustBeNumeric, mustBePositive} = 9.80665 +end -arguments - f - h - options.rho = 1025; - options.g = 9.80665; +arguments (Output) + k struct end -f_py = py.numpy.array(f); +g = options.g; + +% Ensure f is a column vector for consistent processing +f_input = f; +f = f(:); + +% Angular frequency +w = 2 * pi * f; +% note: = h*wa/sqrt(h*g/h) +xi = w / sqrt(g / h); +yi = xi.^2 ./ (1 - exp(-xi.^2.4908)).^0.4015; +k0 = yi / h; -k_py = py.mhkit.wave.resource.wave_number(f_py, h, pyargs('rho', options.rho, 'g', options.rho)); +% Solve dispersion relation: w^2 = g*k*tanh(k*h) +% Rearranged as: w^2 - g*k*tanh(k*h) = 0 +k_values = k0; -k_py = typecast_from_mhkit_python(k_py); +% Set solver options to suppress output +fzero_options = optimset('Display', 'off', 'TolX', 1e-12); +% Only solve for points where initial guess isn't accurate enough +for i = 1:length(f) + % Eq 11 in IEC 62600-101 Ed. 2.0 en 2024 using initial guess from Guo (2002) + residual = w(i)^2 - g * k0(i) * tanh(k0(i) * h); + if abs(residual) > 1e-9 + func = @(kk) w(i)^2 - g * kk * tanh(kk * h); + try + k_values(i) = fzero(func, k0(i), fzero_options); + catch ME + error('MHKiT:wave_number:SolverFailed', ... + 'Wave number solver failed for f=%g Hz: %s', f(i), ME.message); + end + end +end + +% Reshape output to match input shape +if isrow(f_input) + k_values = k_values'; +end + +% Build output structure k = struct(); +k.values = k_values; +k.frequency = f_input; -k.values = k_py.data; -% MHKiT Python does not output frequency -k.frequency = f; +end