From d939467a92215c24dd1e6bcf215f5522b6599772 Mon Sep 17 00:00:00 2001
From: hivanov-nrel <Hristo.Ivanov@nrel.gov>
Date: Tue, 29 Sep 2026 11:49:35 -0600
Subject: [PATCH 01/10] Agent Host changes for
 agents/matlab-pure-implementation-performance

---
 mhkit/river/IO/read_usgs_file.m               | 94 ++++++++++++++----
 mhkit/river/performance/circular.m            | 15 ++-
 mhkit/river/performance/ducted.m              | 15 ++-
 mhkit/river/performance/multiple_circular.m   | 18 ++--
 mhkit/river/performance/power_coefficient.m   | 18 +++-
 mhkit/river/performance/rectangular.m         | 16 ++-
 mhkit/river/performance/tip_speed_ratio.m     | 15 ++-
 mhkit/river/resource/Froude_number.m          | 15 ++-
 mhkit/river/resource/discharge_to_velocity.m  | 75 ++++----------
 mhkit/river/resource/energy_produced.m        | 61 ++++++------
 mhkit/river/resource/exceedance_probability.m | 98 ++++++++-----------
 mhkit/river/resource/polynomial_fit.m         | 34 ++++---
 mhkit/river/resource/velocity_to_power.m      | 69 ++++---------
 13 files changed, 286 insertions(+), 257 deletions(-)

diff --git a/mhkit/river/IO/read_usgs_file.m b/mhkit/river/IO/read_usgs_file.m
index aed8ad0df..e471f8e6e 100644
--- a/mhkit/river/IO/read_usgs_file.m
+++ b/mhkit/river/IO/read_usgs_file.m
@@ -21,36 +21,88 @@
 %         datast.units: units for each parameter
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
+arguments (Input)
+    file_name (1,:) char
+end
+arguments (Output)
+    datast struct
+end
+
+text = jsondecode(fileread(file_name));
+
+time_series = text.value.timeSeries;
+num_series = numel(time_series);
+
+names = strings(num_series, 1);
+units = strings(num_series, 1);
+series_time = cell(num_series, 1);
+series_value = cell(num_series, 1);
+keep = false(num_series, 1);
+
+for i = 1:num_series
+    try
+        series = time_series(i);
 
-datapd=py.mhkit.river.io.usgs.read_usgs_file(file_name);
+        % USGS variableDescription is formatted as "<name>, <units>"
+        description = string(series.variable.variableDescription);
+        name_parts = split(description, ",");
+        names(i) = name_parts(1);
+        if numel(name_parts) > 1
+            units(i) = name_parts(2);
+        end
 
+        % Only the first entry of "values" is used, matching the Python
+        % implementation's use of values[0]["value"]
+        values = series.values(1).value;
 
-xx=cell(datapd.axes);
-v=xx{2};
+        series_time{i} = local_parse_datetime(string({values.dateTime}'));
+        series_value{i} = str2double({values.value}');
+        keep(i) = true;
+    catch ME
+        warning('MATLAB:read_usgs_file', ...
+            'Failed to process time series %d: %s', i, ME.message);
+    end
+end
+
+names = names(keep);
+units = units(keep);
+series_time = series_time(keep);
+series_value = series_value(keep);
+num_series = numel(names);
 
+% Union of all timestamps across parameters, mirroring pandas'
+% combine_first which outer-joins each series on its datetime index
+unique_times = unique(vertcat(series_time{:}));
 
-vv=cell(py.list(py.numpy.nditer(v.values,pyargs("flags",{"refs_ok"}))));
+datast = struct();
+for i = 1:num_series
+    column = nan(numel(unique_times), 1);
+    [is_member, loc] = ismember(series_time{i}, unique_times);
+    column(loc(is_member)) = series_value{i}(is_member);
 
-vals=double(py.array.array('d',py.numpy.nditer(datapd.values,pyargs("flags",{"refs_ok"}))));
-sha=cell(datapd.values.shape);
-x=int64(sha{1,1});
-y=int64(sha{1,2});
+    datast.(names(i)) = column;
+    datast.units.(names(i)) = units(i);
+end
 
-vals=reshape(vals,[x,y]);
-si=size(vals);
+datast.time = posixtime(unique_times).';
 
-for i=1:si(2)
-    test=string(py.str(vv{i}));
-    newname=split(test,",");
-    datast.(newname(1))=vals(:,i);
-    datast.units.(newname(1))=newname(2);
 end
 
-times = double(    ...
-     py.mhkit_python_utils.pandas_dataframe.datetime_index_to_ordinal(datapd));
+function dt = local_parse_datetime(date_strings)
+% Parses USGS dateTime strings (ISO 8601, with or without a UTC offset)
+% into UTC datetimes, mirroring pandas.to_datetime(..., utc=True).
+date_strings = strrep(date_strings, "Z", "+00:00");
 
-datast.time = posixtime(datetime(times,                        ...
-                                 'ConvertFrom', 'datenum',     ...
-                                 'TimeZone','UTC'));
+has_offset = ~isempty(regexp(date_strings(1), '[+-]\d{2}:\d{2}$', 'once'));
+
+if has_offset
+    dt = datetime(date_strings, ...
+        'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss.SSSX', ...
+        'TimeZone', 'UTC');
+else
+    dt = datetime(date_strings, ...
+        'InputFormat', 'yyyy-MM-dd''T''HH:mm:ss.SSS', ...
+        'TimeZone', 'UTC');
+end
+end
 
diff --git a/mhkit/river/performance/circular.m b/mhkit/river/performance/circular.m
index a474a4e08..57403c836 100644
--- a/mhkit/river/performance/circular.m
+++ b/mhkit/river/performance/circular.m
@@ -19,11 +19,16 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
+arguments (Input)
+    diameter (1,1) {mustBeNumeric}
+end
+arguments (Output)
+    D_E (1,1) {mustBeNumeric}
+    projected_capture_area (1,1) {mustBeNumeric}
+end
 
-result=py.mhkit.river.performance.circular(diameter);
+D_E = diameter;
+projected_capture_area = (1/4) * pi * (D_E.^2);
 
-resultc=cell(result);
-D_E=resultc{1};
-projected_capture_area=resultc{2};
+end
 
diff --git a/mhkit/river/performance/ducted.m b/mhkit/river/performance/ducted.m
index 5dba57f02..e59815b2c 100644
--- a/mhkit/river/performance/ducted.m
+++ b/mhkit/river/performance/ducted.m
@@ -18,11 +18,16 @@
 %         Projected capture area [m^2]
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
+arguments (Input)
+    diameter (1,1) {mustBeNumeric}
+end
+arguments (Output)
+    D_E (1,1) {mustBeNumeric}
+    projected_capture_area (1,1) {mustBeNumeric}
+end
 
-result=py.mhkit.river.performance.ducted(diameter);
+D_E = diameter;
+projected_capture_area = (1/4) * pi * (D_E.^2);
 
-resultc=cell(result);
-D_E=resultc{1};
-projected_capture_area=resultc{2};
+end
 
diff --git a/mhkit/river/performance/multiple_circular.m b/mhkit/river/performance/multiple_circular.m
index 0e5c7e105..f5d0d086f 100644
--- a/mhkit/river/performance/multiple_circular.m
+++ b/mhkit/river/performance/multiple_circular.m
@@ -19,11 +19,17 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
-diameters=py.list(diameters);
-result=py.mhkit.river.performance.multiple_circular(diameters);
+arguments (Input)
+    diameters {mustBeNumeric, mustBeVector}
+end
+arguments (Output)
+    D_E (1,1) {mustBeNumeric}
+    projected_capture_area (1,1) {mustBeNumeric}
+end
 
-resultc=cell(result);
-D_E=resultc{1};
-projected_capture_area=resultc{2};
+diameters_squared = diameters.^2;
+D_E = sqrt(sum(diameters_squared));
+projected_capture_area = 0.25 * pi * sum(diameters_squared);
+
+end
 
diff --git a/mhkit/river/performance/power_coefficient.m b/mhkit/river/performance/power_coefficient.m
index 5022b1ce4..1a63b26f6 100644
--- a/mhkit/river/performance/power_coefficient.m
+++ b/mhkit/river/performance/power_coefficient.m
@@ -26,8 +26,20 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py_cp = py.mhkit.river.performance.power_coefficient(py.numpy.array(power),...
-    py.numpy.array(inflow_speed),capture_area,rho);
+arguments (Input)
+    power {mustBeNumeric}
+    inflow_speed {mustBeNumeric}
+    capture_area (1,1) {mustBeNumeric}
+    rho (1,1) {mustBeNumeric}
+end
+arguments (Output)
+    Cp {mustBeNumeric}
+end
 
-Cp = double(py.array.array('d',py.numpy.nditer(py_cp,pyargs("flags",{"refs_ok"}))));
+% Predicted power from inflow
+power_in = 0.5 .* rho .* capture_area .* inflow_speed.^3;
+
+Cp = power ./ power_in;
+
+end
 
diff --git a/mhkit/river/performance/rectangular.m b/mhkit/river/performance/rectangular.m
index a00b0d3b1..85fd0feb7 100644
--- a/mhkit/river/performance/rectangular.m
+++ b/mhkit/river/performance/rectangular.m
@@ -22,11 +22,17 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
+arguments (Input)
+    h (1,1) {mustBeNumeric}
+    w (1,1) {mustBeNumeric}
+end
+arguments (Output)
+    D_E (1,1) {mustBeNumeric}
+    projected_capture_area (1,1) {mustBeNumeric}
+end
 
-result=py.mhkit.river.performance.rectangular(h,w);
+D_E = sqrt(4.0 * h * w / pi);
+projected_capture_area = h * w;
 
-resultc=cell(result);
-D_E=resultc{1};
-projected_capture_area=resultc{2};
+end
 
diff --git a/mhkit/river/performance/tip_speed_ratio.m b/mhkit/river/performance/tip_speed_ratio.m
index 539b1c812..07f8dd453 100644
--- a/mhkit/river/performance/tip_speed_ratio.m
+++ b/mhkit/river/performance/tip_speed_ratio.m
@@ -23,7 +23,18 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py_tsr = py.mhkit.river.performance.tip_speed_ratio(py.numpy.array(rotor_speed),double(rotor_diameter),py.numpy.array(inflow_speed));
+arguments (Input)
+    rotor_speed {mustBeNumeric}
+    rotor_diameter (1,1) {mustBeNumeric}
+    inflow_speed {mustBeNumeric}
+end
+arguments (Output)
+    TSR {mustBeNumeric}
+end
 
-TSR = double(py.array.array('d',py.numpy.nditer(py_tsr,pyargs("flags",{"refs_ok"}))));
+rotor_velocity = rotor_speed .* pi .* rotor_diameter;
+
+TSR = rotor_velocity ./ inflow_speed;
+
+end
 
diff --git a/mhkit/river/resource/Froude_number.m b/mhkit/river/resource/Froude_number.m
index aa1511e2e..c99f2d480 100644
--- a/mhkit/river/resource/Froude_number.m
+++ b/mhkit/river/resource/Froude_number.m
@@ -22,9 +22,16 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-if nargin == 3
-     Fr=py.mhkit.river.resource.Froude_number(v,h,pyargs('g',g));
-else
-     Fr=py.mhkit.river.resource.Froude_number(v,h);
+arguments (Input)
+    v (1,1) {mustBeNumeric}
+    h (1,1) {mustBeNumeric}
+    g (1,1) {mustBeNumeric} = 9.80665
+end
+arguments (Output)
+    Fr (1,1) {mustBeNumeric}
+end
+
+Fr = v / sqrt(g * h);
+
 end
 
diff --git a/mhkit/river/resource/discharge_to_velocity.m b/mhkit/river/resource/discharge_to_velocity.m
index 1ec1418e5..dcb97bfe1 100644
--- a/mhkit/river/resource/discharge_to_velocity.m
+++ b/mhkit/river/resource/discharge_to_velocity.m
@@ -8,22 +8,16 @@
 % ------------
 %     Q : Discharge data [m3/s]
 %
-%         Pandas dataframe indexed by time [datetime or s]:
-%
-%           To make a pandas data frame from user supplied frequency and spectra
-%           use py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(timeseries,time,x)
-%
-%         OR
-%
 %         structure of form:
 %
 %            Q.Discharge
 %
 %            Q.time
 %
-%     polynomial_coefficients : numpy polynomial
-%         List of polynomial coefficients that discribe the relationship between
-%         discharge and velocity at an individual turbine
+%     polynomial_coefficients : vector
+%         Vector of polynomial coefficients (highest degree first) that
+%         describe the relationship between discharge and velocity at an
+%         individual turbine
 %
 % Returns
 % ------------
@@ -37,56 +31,21 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
-py.importlib.import_module('mhkit_python_utils');
-
-if (isa(Q,'py.pandas.core.frame.DataFrame')~=1)
-    x=size(Q.Discharge);
-    li=py.list();
-    if x(2)>1
-        for i = 1:x(2)
-            app=py.list(Q.Discharge(:,i));
-            li=py.mhkit_python_utils.pandas_dataframe.lis(li,app);
-
-        end
-    elseif x(2) ==1
-        li=Q.Discharge;
-    end
-
-    if any(isdatetime(Q.time(1)))
-        si=size(Q.time);
-        for i=1:si(2)
-        Q.time(i)=posixtime(Q.time(i));
-        end
-    end
-    % Q=py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(li,Q.time,int32(x(2)));
-    Q = py.mhkit_python_utils.pandas_dataframe.list_to_series(Q.Discharge, Q.time);
+arguments (Input)
+    Q struct
+    polynomial_coefficients {mustBeNumeric, mustBeVector}
+end
+arguments (Output)
+    V struct
 end
 
-polynomial_coefficients=py.numpy.poly1d(polynomial_coefficients);
-
-Vdf=py.mhkit.river.resource.discharge_to_velocity(Q,polynomial_coefficients);
-
-
-
-xx=cell(Vdf.axes);
-v=xx{2};
-vv=cell(py.list(py.numpy.nditer(v.values,pyargs("flags",{"refs_ok"}))));
-
-vals=double(py.array.array('d',py.numpy.nditer(Vdf.values)));
-sha=cell(Vdf.values.shape);
-x=int64(sha{1,1});
-y=int64(sha{1,2});
-
-vals=reshape(vals,[x,y]);
-
-si=size(vals);
- for i=1:si(2)
-    test=string(py.str(vv{i}));
-    newname=split(test,",");
+time = Q.time;
+if any(isdatetime(time))
+    time = posixtime(time);
+end
 
-    V.(newname(1))=vals(:,i);
+V.V = polyval(polynomial_coefficients, Q.Discharge);
+V.time = time;
 
- end
- V.time=double(py.array.array('d',py.numpy.nditer(Vdf.index)));
+end
 
diff --git a/mhkit/river/resource/energy_produced.m b/mhkit/river/resource/energy_produced.m
index a3bcb1e3f..fdb509350 100644
--- a/mhkit/river/resource/energy_produced.m
+++ b/mhkit/river/resource/energy_produced.m
@@ -8,13 +8,6 @@
 % ----------
 %     P : Power [W]
 %
-%         Pandas dataframe indexed by time [datetime or s]:
-%
-%           To make a pandas data frame from user supplied frequency and spectra
-%           use py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(timeseries,time,x)
-%
-%         OR
-%
 %         structure of form:
 %
 %           P.P
@@ -26,37 +19,47 @@
 %
 % Returns
 % -------
-%     E : Structure
-%
-%
-%        P.P: Power [W]
-%
-%        P.time: epoch time [s]
+%     E : float
+%         Energy [J] produced in the given length of time
 %
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit_python_utils');
-py.importlib.import_module('mhkit');
+arguments (Input)
+    P struct
+    seconds (1,1) {mustBeNumeric}
+end
+arguments (Output)
+    E (1,1) {mustBeNumeric}
+end
 
+power_data = P.P(:);
 
-if (isa(P,'py.pandas.core.frame.DataFrame')~=1)
-    x=size(P.P);
-    li=py.list();
-    if x(2)>1
-        for i = 1:x(2)
-            app=py.list(P.P(:,i));
-            li=py.mhkit_python_utils.pandas_dataframe.lis(li,app);
+% Histogram of power with 100 equal-width bins spanning the data range
+[counts, edges] = histcounts(power_data, 100);
+bin_widths = diff(edges);
+total_count = sum(counts);
 
-        end
-    elseif x(2) ==1
-        li=P.P;
-    end
+% Piecewise-constant probability density function of the histogram
+density = counts ./ (total_count * bin_widths);
 
+x = linspace(edges(1), edges(end), 1000);
+pdf_x = zeros(size(x));
 
-    % P=py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(li,P.time,int32(x(2)));
-    P = py.mhkit_python_utils.pandas_dataframe.list_to_series(li, P.time);
+nbins = numel(counts);
+for i = 1:nbins
+    if i < nbins
+        mask = x >= edges(i) & x < edges(i + 1);
+    else
+        mask = x >= edges(i) & x <= edges(i + 1);
+    end
+    pdf_x(mask) = density(i);
 end
 
-E=py.mhkit.river.resource.energy_produced(P,seconds);
+% Expected value of power via trapezoidal integration of x*pdf(x)
+expected_power = trapz(x, x .* pdf_x);
+
+E = seconds * expected_power;
+
+end
 
diff --git a/mhkit/river/resource/exceedance_probability.m b/mhkit/river/resource/exceedance_probability.m
index 1f2c7ac59..776ec4f81 100644
--- a/mhkit/river/resource/exceedance_probability.m
+++ b/mhkit/river/resource/exceedance_probability.m
@@ -7,13 +7,6 @@
 % ----------
 %     Q : Discharge data [m3/s]
 %
-%         Pandas dataframe indexed by time [datetime or s]
-%
-%           To make a pandas data frame from user supplied frequency and spectra
-%           use py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(timeseries,time,x)
-%
-%         OR
-%
 %         structure of form:
 %
 %           Q.Discharge
@@ -33,62 +26,51 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
-py.importlib.import_module('mhkit_python_utils');
-py.importlib.import_module('mhkit');
-
-% Check if Q is a dataframe
-% If it is not, it is expecting a struct with Discharge and time "keys"
-if (isa(Q,'py.pandas.core.frame.DataFrame')~=1)
-    x=size(Q.Discharge);
-    li=py.list();
-    if x(2)>1
-        for i = 1:x(2)
-            app=py.list(Q.Discharge(:,i));
-            li=py.mhkit_python_utils.pandas_dataframe.lis(li,app);
-
-        end
-    elseif x(2) ==1
-        li=Q.Discharge;
-    end
-
-    if (~isempty(Q.time))
-        if any(isdatetime(Q.time(1)))
-            si=size(Q.time);
-            for i=1:si(2)
-            Q.time(i)=posixtime(Q.time(i));
-            end
-        end
-    end
-    % Q=py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(li,Q.time,int32(x(2)));
-
-    if (~isempty(Q.time))
-        Q = py.mhkit_python_utils.pandas_dataframe.list_to_series(Q.Discharge, Q.time);
-    else
-        Q = py.mhkit_python_utils.pandas_dataframe.list_to_series(Q.Discharge);
-    end
-
+arguments (Input)
+    Q struct
+end
+arguments (Output)
+    F struct
 end
 
-EPpd=py.mhkit.river.resource.exceedance_probability(Q);
-
-xx=cell(EPpd.axes);
-v=xx{2};
-vv=cell(py.list(py.numpy.nditer(v.values,pyargs("flags",{"refs_ok"}))));
+time = Q.time;
+if any(isdatetime(time))
+    time = posixtime(time);
+end
 
-vals=double(py.array.array('d',py.numpy.nditer(EPpd.values)));
-sha=cell(EPpd.values.shape);
-x=int64(sha{1,1});
-y=int64(sha{1,2});
+discharge = Q.Discharge;
+n = size(discharge, 1);
 
-vals=reshape(vals,[x,y]);
+rank_ascending = zeros(size(discharge));
+for col = 1:size(discharge, 2)
+    rank_ascending(:, col) = local_average_rank(discharge(:, col));
+end
 
-si=size(vals);
- for i=1:si(2)
-    test=string(py.str(vv{i}));
-    newname=split(test,",");
+% Convert to descending rank so the smallest value has the highest
+% exceedance probability
+rank_descending = n - rank_ascending + 1;
+F.F = 100 * rank_descending / (n + 1);
+F.time = time;
 
-    F.(newname(1))=vals(:,i);
+end
 
- end
- F.time=double(py.array.array('d',py.numpy.nditer(EPpd.index)));
+function r = local_average_rank(x)
+% Assigns ascending ranks (starting at 1) to the elements of x, averaging
+% the ranks of tied values. This mirrors the "average" tie-breaking method
+% used by xarray/scipy when computing exceedance probability.
+x = x(:);
+n = numel(x);
+[sorted_x, order] = sort(x);
+r = zeros(n, 1);
+
+i = 1;
+while i <= n
+    j = i;
+    while j < n && sorted_x(j + 1) == sorted_x(i)
+        j = j + 1;
+    end
+    r(order(i:j)) = (i + j) / 2;
+    i = j + 1;
+end
+end
 
diff --git a/mhkit/river/resource/polynomial_fit.m b/mhkit/river/resource/polynomial_fit.m
index a06e215cc..ba5a5bf2f 100644
--- a/mhkit/river/resource/polynomial_fit.m
+++ b/mhkit/river/resource/polynomial_fit.m
@@ -25,17 +25,27 @@
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-py.importlib.import_module('mhkit');
-x=py.numpy.array(x);
-y=py.numpy.array(y);
-n=int32(n);
-
-polyt=py.mhkit.river.resource.polynomial_fit(x,y,n);
-
-polyc=cell(polyt);
-coef=polyc{1};
-fit=polyc{2};
-poly.coef=double(py.array.array('d',py.numpy.nditer(coef.coef)));
-poly.fit=fit;
+arguments (Input)
+    x {mustBeNumeric, mustBeVector}
+    y {mustBeNumeric, mustBeVector}
+    n (1,1) {mustBeInteger}
+end
+arguments (Output)
+    poly struct
+end
+
+% Coefficients ordered highest degree first, matching numpy's poly1d
+coef = polyfit(x, y, n);
+y_fit = polyval(coef, x);
+
+% R-squared is the squared Pearson correlation between actual and fitted
+% values, matching scipy.stats.linregress(y, polynomial_coefficients(x))
+correlation = corrcoef(y, y_fit);
+r_squared = correlation(1, 2)^2;
+
+poly.coef = coef;
+poly.fit = r_squared;
+
+end
 
 
diff --git a/mhkit/river/resource/velocity_to_power.m b/mhkit/river/resource/velocity_to_power.m
index 2e294538d..5d0e1f747 100644
--- a/mhkit/river/resource/velocity_to_power.m
+++ b/mhkit/river/resource/velocity_to_power.m
@@ -8,13 +8,6 @@
 % ----------
 %     V : Velocity [m/s]
 %
-%          Pandas dataframe indexed by time [datetime or s]
-%
-%           To make a pandas data frame from user supplied frequency and spectra
-%           use py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(timeseries,time,x)
-%
-%         OR
-%
 %         structure of form:
 %
 %           V.V: Velocity [m/s]
@@ -44,52 +37,30 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
-py.importlib.import_module('mhkit_python_utils');
-py.importlib.import_module('mhkit');
-
-if (isa(V,'py.pandas.core.frame.DataFrame')~=1)
-    x=size(V.V);
-    li=py.list();
-    if x(2)>1
-        for i = 1:x(2)
-            app=py.list(V.V(:,i));
-            li=py.mhkit_python_utils.pandas_dataframe.lis(li,app);
-
-        end
-    elseif x(2) ==1
-        li=V.V;
-    end
-
-
-    % V=py.mhkit_python_utils.pandas_dataframe.timeseries_to_pandas(li,V.time,int32(x(2)));
-    V = py.mhkit_python_utils.pandas_dataframe.list_to_series(li, V.time);
+arguments (Input)
+    V struct
+    polynomial_coefficients {mustBeNumeric, mustBeVector}
+    cut_in (1,1) {mustBeNumeric}
+    cut_out (1,1) {mustBeNumeric}
+end
+arguments (Output)
+    p struct
 end
 
-polynomial_coefficients=py.numpy.poly1d(polynomial_coefficients);
-cut_in=py.float(cut_in);
-cut_out=py.float(cut_out);
-Pdf=py.mhkit.river.resource.velocity_to_power(V,polynomial_coefficients,cut_in,cut_out);
-
-
-
-xx=cell(Pdf.axes);
-v=xx{2};
-vv=cell(py.list(py.numpy.nditer(v.values,pyargs("flags",{"refs_ok"}))));
-
-vals=double(py.array.array('d',py.numpy.nditer(Pdf.values)));
-sha=cell(Pdf.values.shape);
-x=int64(sha{1,1});
-y=int64(sha{1,2});
+time = V.time;
+if any(isdatetime(time))
+    time = posixtime(time);
+end
 
-vals=reshape(vals,[x,y]);
+velocity = V.V;
+power = polyval(polynomial_coefficients, velocity);
 
-si=size(vals);
- for i=1:si(2)
-    test=string(py.str(vv{i}));
-    newname=split(test,",");
+% Turbine produces 0 power outside of the cut-in/cut-out bounds
+power(velocity < cut_in) = 0.0;
+power(velocity > cut_out) = 0.0;
 
-    p.(newname(1))=vals(:,i);
+p.P = power;
+p.time = time;
 
- end
- p.time=double(py.array.array('d',py.numpy.nditer(Pdf.index)));
+end
 

From 06316eaf84c60bdc11383611f32ef36b0ed366cf Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:00:03 -0600
Subject: [PATCH 02/10] River: Rename Froude_number -> froude_number

---
 mhkit/river/resource/{Froude_number.m => froude_number.m} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename mhkit/river/resource/{Froude_number.m => froude_number.m} (100%)

diff --git a/mhkit/river/resource/Froude_number.m b/mhkit/river/resource/froude_number.m
similarity index 100%
rename from mhkit/river/resource/Froude_number.m
rename to mhkit/river/resource/froude_number.m

From c7221941241287ce8958171601ec8bd5d8e52880 Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:10:24 -0600
Subject: [PATCH 03/10] River: Add detail to read_usgs_file docstring types

---
 mhkit/river/IO/read_usgs_file.m | 48 ++++++++++++++++++---------------
 1 file changed, 26 insertions(+), 22 deletions(-)

diff --git a/mhkit/river/IO/read_usgs_file.m b/mhkit/river/IO/read_usgs_file.m
index e471f8e6e..7d3f16312 100644
--- a/mhkit/river/IO/read_usgs_file.m
+++ b/mhkit/river/IO/read_usgs_file.m
@@ -1,31 +1,35 @@
-function datast=read_usgs_file(file_name)
+function dataset=read_usgs_file(file_name)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Reads a USGS JSON data file (from https://waterdata.usgs.gov/nwis)
-%     into a structure
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Reads a USGS JSON data file (from https://waterdata.usgs.gov/nwis)
+% into a structure
 %
 % Parameters
-% ----------
-%     file_name : str
-%         Name of USGS JSON data file
+% ------------
+% file_name : string
+%   Name of USGS JSON data file
 %
 % Returns
-% -------
-%     datast : structure
-%
-%
-%         datast.Data: named according to the parameter's variable description
+% ---------
+% dataset : struct
+%   USGS data, one field per parameter
+%     dataset.(parameter) : column vector
+%       Data named according to the parameter's variable description
+%       (e.g. dataset.Discharge)
+%     dataset.time : row vector [s]
+%       Epoch time
+%     dataset.units : struct
+%       Units for each parameter (e.g. dataset.units.Discharge)
 %
-%         datast.time: epoch time [s]
-%
-%         datast.units: units for each parameter
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
-    file_name (1,:) char
+    file_name {mustBeTextScalar, mustBeFile}
 end
+
 arguments (Output)
-    datast struct
+    dataset struct
 end
 
 text = jsondecode(fileread(file_name));
@@ -74,17 +78,17 @@
 % combine_first which outer-joins each series on its datetime index
 unique_times = unique(vertcat(series_time{:}));
 
-datast = struct();
+dataset = struct();
 for i = 1:num_series
     column = nan(numel(unique_times), 1);
     [is_member, loc] = ismember(series_time{i}, unique_times);
     column(loc(is_member)) = series_value{i}(is_member);
 
-    datast.(names(i)) = column;
-    datast.units.(names(i)) = units(i);
+    dataset.(names(i)) = column;
+    dataset.units.(names(i)) = units(i);
 end
 
-datast.time = posixtime(unique_times).';
+dataset.time = posixtime(unique_times).';
 
 end
 

From 13ba50dcb08c7e31bd6597d64d2d0b258eec2223 Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:13:06 -0600
Subject: [PATCH 04/10] River: Format docstrings and arguments block

---
 mhkit/river/performance/circular.m            | 26 ++++----
 mhkit/river/performance/ducted.m              | 25 ++++----
 mhkit/river/performance/multiple_circular.m   | 24 +++----
 mhkit/river/performance/power_coefficient.m   | 34 +++++-----
 mhkit/river/performance/rectangular.m         | 31 +++++-----
 mhkit/river/performance/tip_speed_ratio.m     | 29 ++++-----
 mhkit/river/resource/discharge_to_velocity.m  | 52 ++++++++--------
 mhkit/river/resource/energy_produced.m        | 39 ++++++------
 mhkit/river/resource/exceedance_probability.m | 40 ++++++------
 mhkit/river/resource/froude_number.m          | 34 +++++-----
 mhkit/river/resource/polynomial_fit.m         | 55 +++++++++-------
 mhkit/river/resource/velocity_to_power.m      | 62 +++++++++----------
 12 files changed, 225 insertions(+), 226 deletions(-)

diff --git a/mhkit/river/performance/circular.m b/mhkit/river/performance/circular.m
index 57403c836..bd2c51ef7 100644
--- a/mhkit/river/performance/circular.m
+++ b/mhkit/river/performance/circular.m
@@ -1,27 +1,28 @@
-function [D_E,projected_capture_area]=circular(diameter)
+function [D_E, projected_capture_area] = circular(diameter)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates the equivalent diameter and projected capture area of a
-%     circular turbine
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Calculates the equivalent diameter and projected capture area of a
+% circular turbine
 %
 % Parameters
 % ------------
-%     diameter : float
-%         Turbine diameter [m]
+% diameter : double [m]
+%   Turbine diameter
 %
 % Returns
 % ---------
-%     D_E : float
-%        Equivalent diameter [m]
-%
-%     projected_capture_area : float
-%         Projected capture area [m^2]
+% D_E : double [m]
+%   Equivalent diameter
+% projected_capture_area : double [m^2]
+%   Projected capture area
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     diameter (1,1) {mustBeNumeric}
 end
+
 arguments (Output)
     D_E (1,1) {mustBeNumeric}
     projected_capture_area (1,1) {mustBeNumeric}
@@ -31,4 +32,3 @@
 projected_capture_area = (1/4) * pi * (D_E.^2);
 
 end
-
diff --git a/mhkit/river/performance/ducted.m b/mhkit/river/performance/ducted.m
index e59815b2c..268695970 100644
--- a/mhkit/river/performance/ducted.m
+++ b/mhkit/river/performance/ducted.m
@@ -1,26 +1,28 @@
-function [D_E,projected_capture_area]=ducted(diameter)
+function [D_E, projected_capture_area] = ducted(diameter)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates the equivalent diameter and projected capture area of a
-%     ducted turbine
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Calculates the equivalent diameter and projected capture area of a
+% ducted turbine
 %
 % Parameters
 % ------------
-%     diameter : float
-%         ducted diameter [m]
+% diameter : double [m]
+%   Duct diameter
 %
 % Returns
 % ---------
-%     D_E : float
-%        Equivalent diameter [m]
+% D_E : double [m]
+%   Equivalent diameter
+% projected_capture_area : double [m^2]
+%   Projected capture area
 %
-%     projected_capture_area : float
-%         Projected capture area [m^2]
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     diameter (1,1) {mustBeNumeric}
 end
+
 arguments (Output)
     D_E (1,1) {mustBeNumeric}
     projected_capture_area (1,1) {mustBeNumeric}
@@ -30,4 +32,3 @@
 projected_capture_area = (1/4) * pi * (D_E.^2);
 
 end
-
diff --git a/mhkit/river/performance/multiple_circular.m b/mhkit/river/performance/multiple_circular.m
index f5d0d086f..0bc9db627 100644
--- a/mhkit/river/performance/multiple_circular.m
+++ b/mhkit/river/performance/multiple_circular.m
@@ -1,27 +1,28 @@
-function [D_E,projected_capture_area]=multiple_circular(diameters)
+function [D_E, projected_capture_area] = multiple_circular(diameters)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates the equivalent diameter and projected capture area of a
-%     multiple circular turbine
+%
+% Calculates the equivalent diameter and projected capture area of a
+% multiple circular turbine
 %
 % Parameters
 % ------------
-%     diameters: array or vector
-%         vector of device diameters [m]
+% diameters : vector [m]
+%   Device diameters
 %
 % Returns
 % ---------
-%     D_E : float
-%        Equivalent diameter [m]
-%
-%     projected_capture_area : float
-%         Projected capture area [m^2]
+% D_E : double [m]
+%   Equivalent diameter
+% projected_capture_area : double [m^2]
+%   Projected capture area
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     diameters {mustBeNumeric, mustBeVector}
 end
+
 arguments (Output)
     D_E (1,1) {mustBeNumeric}
     projected_capture_area (1,1) {mustBeNumeric}
@@ -32,4 +33,3 @@
 projected_capture_area = 0.25 * pi * sum(diameters_squared);
 
 end
-
diff --git a/mhkit/river/performance/power_coefficient.m b/mhkit/river/performance/power_coefficient.m
index 1a63b26f6..da75b4086 100644
--- a/mhkit/river/performance/power_coefficient.m
+++ b/mhkit/river/performance/power_coefficient.m
@@ -1,30 +1,26 @@
-function Cp=power_coefficient(power, inflow_speed, capture_area, rho)
+function Cp = power_coefficient(power, inflow_speed, capture_area, rho)
 
-%%%%%%%%%%%%%%%%%%%%
-%     Function that calculates the power coefficient of MEC device
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
+% Calculates the power coefficient of a MEC device
 %
 % Parameters
 % ------------
-%     power : vector
-%         Power output signal of device after losses [W]
-%
-%     inflow_speed : vector
-%         Velocity of inflow condition [m/s]
-%
-%     capture_area : double or int
-%         Projected area of rotor normal to inflow [m^2]
-%
-%     rho : double or int
-%         Density of environment [kg/m^3]
+% power : numeric [W]
+%   Power output signal of device after losses
+% inflow_speed : numeric [m/s]
+%   Velocity of inflow condition
+% capture_area : double [m^2]
+%   Projected area of rotor normal to inflow
+% rho : double [kg/m^3]
+%   Density of environment
 %
 % Returns
 % ---------
-%     Cp: vector
-%         Power coefficient of device [-]
+% Cp : numeric [-]
+%   Power coefficient of device
 %
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     power {mustBeNumeric}
@@ -32,6 +28,7 @@
     capture_area (1,1) {mustBeNumeric}
     rho (1,1) {mustBeNumeric}
 end
+
 arguments (Output)
     Cp {mustBeNumeric}
 end
@@ -42,4 +39,3 @@
 Cp = power ./ power_in;
 
 end
-
diff --git a/mhkit/river/performance/rectangular.m b/mhkit/river/performance/rectangular.m
index 85fd0feb7..05da8a8d7 100644
--- a/mhkit/river/performance/rectangular.m
+++ b/mhkit/river/performance/rectangular.m
@@ -1,31 +1,31 @@
-function [D_E,projected_capture_area]=rectangular(h,w)
+function [D_E, projected_capture_area] = rectangular(h, w)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates the equivalent diameter and projected capture area of a
-%     retangular turbine
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Calculates the equivalent diameter and projected capture area of a
+% rectangular turbine
 %
 % Parameters
 % ------------
-%     h : float
-%         Turbine height [m]
-%
-%     w : float
-%         Turbine width [m]
+% h : double [m]
+%   Turbine height
+% w : double [m]
+%   Turbine width
 %
 % Returns
 % ---------
-%     D_E : float
-%        Equivalent diameter [m]
-%
-%     projected_capture_area : float
-%         Projected capture area [m^2]
+% D_E : double [m]
+%   Equivalent diameter
+% projected_capture_area : double [m^2]
+%   Projected capture area
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     h (1,1) {mustBeNumeric}
     w (1,1) {mustBeNumeric}
 end
+
 arguments (Output)
     D_E (1,1) {mustBeNumeric}
     projected_capture_area (1,1) {mustBeNumeric}
@@ -35,4 +35,3 @@
 projected_capture_area = h * w;
 
 end
-
diff --git a/mhkit/river/performance/tip_speed_ratio.m b/mhkit/river/performance/tip_speed_ratio.m
index 07f8dd453..bb6c85a2f 100644
--- a/mhkit/river/performance/tip_speed_ratio.m
+++ b/mhkit/river/performance/tip_speed_ratio.m
@@ -1,33 +1,31 @@
-function TSR=tip_speed_ratio(rotor_speed, rotor_diameter, inflow_speed)
+function TSR = tip_speed_ratio(rotor_speed, rotor_diameter, inflow_speed)
 
-%%%%%%%%%%%%%%%%%%%%
-%     Function used to calculate the tip speed ratio (TSR) of a MEC device with rotor
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
+% Calculates the tip speed ratio (TSR) of a MEC device with rotor
 %
 % Parameters
 % ------------
-%     rotor_speed : vector
-%         Rotor Speed [rps]
-%
-%     rotor_diameter : double or int
-%         diameter -f rotor [m]
-%
-%     inflow_speed : vector
-%         Velocity of inflow condition [m/s]
+% rotor_speed : numeric [rev/s]
+%   Rotor speed
+% rotor_diameter : double [m]
+%   Diameter of rotor
+% inflow_speed : numeric [m/s]
+%   Velocity of inflow condition
 %
 % Returns
 % ---------
-%     TSR: vector
-%         Calculated tip speed ratio (TSR)
-%
+% TSR : numeric [-]
+%   Calculated tip speed ratio
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     rotor_speed {mustBeNumeric}
     rotor_diameter (1,1) {mustBeNumeric}
     inflow_speed {mustBeNumeric}
 end
+
 arguments (Output)
     TSR {mustBeNumeric}
 end
@@ -37,4 +35,3 @@
 TSR = rotor_velocity ./ inflow_speed;
 
 end
-
diff --git a/mhkit/river/resource/discharge_to_velocity.m b/mhkit/river/resource/discharge_to_velocity.m
index dcb97bfe1..f603816f7 100644
--- a/mhkit/river/resource/discharge_to_velocity.m
+++ b/mhkit/river/resource/discharge_to_velocity.m
@@ -1,40 +1,39 @@
-function V=discharge_to_velocity(Q,polynomial_coefficients)
+function V = discharge_to_velocity(Q, polynomial_coefficients)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates velocity given discharge data and the relationship between
-%     discharge and velocity at an individual turbine
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Calculates velocity given discharge data and the relationship between
+% discharge and velocity at an individual turbine
 %
 % Parameters
 % ------------
-%     Q : Discharge data [m3/s]
-%
-%         structure of form:
-%
-%            Q.Discharge
-%
-%            Q.time
-%
-%     polynomial_coefficients : vector
-%         Vector of polynomial coefficients (highest degree first) that
-%         describe the relationship between discharge and velocity at an
-%         individual turbine
+% Q : struct
+%   Discharge data
+%     Q.Discharge : vector or matrix [m^3/s]
+%       Discharge
+%     Q.time : vector [datetime or s]
+%       Time
+% polynomial_coefficients : vector
+%   Polynomial coefficients (highest degree first, e.g. poly.coef from
+%   polynomial_fit) that describe the relationship between discharge and
+%   velocity at an individual turbine
 %
 % Returns
-% ------------
-%     V: Structure
-%
-%
-%         V.V: Velocity [m/s]
-%
-%         V.time: time [datetime or s]
-%
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% ---------
+% V : struct
+%   Velocity data
+%     V.V : vector or matrix [m/s]
+%       Velocity, one value per discharge value
+%     V.time : vector [s]
+%       Time, with datetime converted to epoch seconds
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     Q struct
     polynomial_coefficients {mustBeNumeric, mustBeVector}
 end
+
 arguments (Output)
     V struct
 end
@@ -48,4 +47,3 @@
 V.time = time;
 
 end
-
diff --git a/mhkit/river/resource/energy_produced.m b/mhkit/river/resource/energy_produced.m
index fdb509350..f8142ff9a 100644
--- a/mhkit/river/resource/energy_produced.m
+++ b/mhkit/river/resource/energy_produced.m
@@ -1,29 +1,30 @@
-function E=energy_produced(P,seconds)
+function E = energy_produced(P, seconds)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Returns the energy produced for a given time period provided
-%     exceedence probability and power.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-% Parameters
-% ----------
-%     P : Power [W]
-%
-%         structure of form:
-%
-%           P.P
+% Returns the energy produced for a given time period provided power
 %
-%           P.time [s]
+% The expected power is computed from a 100-bin histogram of the power
+% data, matching MHKiT-Python's use of numpy.histogram and
+% scipy.stats.rv_histogram.
 %
-%     seconds : float
-%         seconds in the time period of interest
+% Parameters
+% ------------
+% P : struct
+%   Power data
+%     P.P : vector or matrix [W]
+%       Power
+%     P.time : vector [s]
+%       Time
+% seconds : double [s]
+%   Seconds in the time period of interest
 %
 % Returns
-% -------
-%     E : float
-%         Energy [J] produced in the given length of time
-%
+% ---------
+% E : double [J]
+%   Energy produced in the given length of time
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     P struct
diff --git a/mhkit/river/resource/exceedance_probability.m b/mhkit/river/resource/exceedance_probability.m
index 776ec4f81..56e6cd49a 100644
--- a/mhkit/river/resource/exceedance_probability.m
+++ b/mhkit/river/resource/exceedance_probability.m
@@ -1,34 +1,36 @@
-function F=exceedance_probability(Q)
+function F = exceedance_probability(Q)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates the exceedance probability
 %
-% Parameters
-% ----------
-%     Q : Discharge data [m3/s]
-%
-%         structure of form:
+% Calculates the exceedance probability
 %
-%           Q.Discharge
+% Tied values are assigned their average rank, matching xarray's rank
+% used by MHKiT-Python.
 %
-%           Q.time
+% Parameters
+% ------------
+% Q : struct
+%   Discharge data
+%     Q.Discharge : vector or matrix [m^3/s]
+%       Discharge, one timeseries per column
+%     Q.time : vector [datetime or s]
+%       Time
 %
 % Returns
-% -------
-%     F : Structure
-%
-%
-%         F.F: Exceedance probability [unitless]
+% ---------
+% F : struct
+%   Exceedance probability data
+%     F.F : vector or matrix [%]
+%       Exceedance probability, one value per discharge value
+%     F.time : vector [s]
+%       Time, with datetime converted to epoch seconds
 %
-%         F.time: time [epoch time (s)]
-%
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     Q struct
 end
+
 arguments (Output)
     F struct
 end
diff --git a/mhkit/river/resource/froude_number.m b/mhkit/river/resource/froude_number.m
index c99f2d480..20925ca60 100644
--- a/mhkit/river/resource/froude_number.m
+++ b/mhkit/river/resource/froude_number.m
@@ -1,37 +1,37 @@
-function Fr=Froude_number(v,h,g)
+function Fr = froude_number(v, h, options)
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculate the Froude Number of the river, channel or duct flow,
-%     to check subcritical flow assumption (if Fr <1).
 %
-% Parameters
-% ------------
-%     v : float
-%         Average Velocity [m/s].
+% Calculates the Froude number of the river, channel or duct flow
 %
-%     h : float
-%         Mean hydrolic depth float [m].
+% Used to check the subcritical flow assumption (Fr < 1).
 %
-%     g : float (optional)
-%         gravitational acceleration [m/s2].
+% Parameters
+% ------------
+% v : double [m/s]
+%   Average velocity
+% h : double [m]
+%   Mean hydraulic depth
+% g : double [m/s^2] (optional)
+%   Name-value argument. Gravitational acceleration, default 9.80665
 %
 % Returns
 % ---------
-%     Fr : float
-%         Froude Number of the river [unitless].
+% Fr : double [-]
+%   Froude number of the river
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     v (1,1) {mustBeNumeric}
     h (1,1) {mustBeNumeric}
-    g (1,1) {mustBeNumeric} = 9.80665
+    options.g (1,1) {mustBeNumeric} = 9.80665
 end
+
 arguments (Output)
     Fr (1,1) {mustBeNumeric}
 end
 
-Fr = v / sqrt(g * h);
+Fr = v / sqrt(options.g * h);
 
 end
-
diff --git a/mhkit/river/resource/polynomial_fit.m b/mhkit/river/resource/polynomial_fit.m
index ba5a5bf2f..99fe8939a 100644
--- a/mhkit/river/resource/polynomial_fit.m
+++ b/mhkit/river/resource/polynomial_fit.m
@@ -1,39 +1,50 @@
-function poly=polynomial_fit(x,y,n)
+function poly = polynomial_fit(x, y, n)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Returns a polynomial fit for y given x of order n.
-% 
-% Parameters
-% ----------
-%     x : array
-%         x data for polynomial fit.
-%
-%     y : array
-%         y data for polynomial fit.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%     n : int
-%         order of the polynomial fit.
-% 
-% Returns
-% --------
-%     poly: structure
+% Returns a polynomial fit for y given x of order n
 %
+% Also returns the R-squared score of the fit, computed as the squared
+% Pearson correlation between y and the fitted values, matching
+% MHKiT-Python's use of scipy.stats.linregress.
 %
-%       poly.coef: polynomial coefficients 
+% Parameters
+% ------------
+% x : vector
+%   x data for polynomial fit
+% y : vector
+%   y data for polynomial fit, same length as x
+% n : integer
+%   Order of the polynomial fit
 %
-%       poly.fit: fit coefficients
+% Returns
+% ---------
+% poly : struct
+%   Polynomial fit
+%     poly.coef : row vector
+%       Polynomial coefficients, highest degree first (as polyval
+%       expects)
+%     poly.fit : double [-]
+%       R-squared coefficient of determination of the fit
 %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     x {mustBeNumeric, mustBeVector}
     y {mustBeNumeric, mustBeVector}
-    n (1,1) {mustBeInteger}
+    n (1,1) {mustBeInteger, mustBeNonnegative}
 end
+
 arguments (Output)
     poly struct
 end
 
+if numel(x) ~= numel(y)
+    error('MHKiT:polynomial_fit:InvalidInput', ...
+        'polynomial_fit requires x (%d) and y (%d) to have the same length.', ...
+        numel(x), numel(y));
+end
+
 % Coefficients ordered highest degree first, matching numpy's poly1d
 coef = polyfit(x, y, n);
 y_fit = polyval(coef, x);
@@ -47,5 +58,3 @@
 poly.fit = r_squared;
 
 end
-
-
diff --git a/mhkit/river/resource/velocity_to_power.m b/mhkit/river/resource/velocity_to_power.m
index 5d0e1f747..1090c593a 100644
--- a/mhkit/river/resource/velocity_to_power.m
+++ b/mhkit/river/resource/velocity_to_power.m
@@ -1,41 +1,37 @@
-function p=velocity_to_power(V,polynomial_coefficients,cut_in,cut_out)
+function p = velocity_to_power(V, polynomial_coefficients, cut_in, cut_out)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%     Calculates power given velocity data and the relationship
-%     between velocity and power from an individual turbine
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-% Parameters
-% ----------
-%     V : Velocity [m/s]
-%
-%         structure of form:
-%
-%           V.V: Velocity [m/s]
-%
-%           V.time: time [datetime or s]
+% Calculates power given velocity data and the relationship between
+% velocity and power from an individual turbine
 %
-%     polynomial_coefficients : vector
-%         vector of polynomial coefficients that discribe the relationship between
-%         velocity and power at an individual turbine
-%
-%     cut_in: float
-%         Velocity values below cut_in are not used to compute P
-%
-%     cut_out: float
-%         Velocity values above cut_out are not used to compute P
+% Parameters
+% ------------
+% V : struct
+%   Velocity data
+%     V.V : vector or matrix [m/s]
+%       Velocity
+%     V.time : vector [datetime or s]
+%       Time
+% polynomial_coefficients : vector
+%   Polynomial coefficients (highest degree first, e.g. poly.coef from
+%   polynomial_fit) that describe the relationship between velocity and
+%   power at an individual turbine
+% cut_in : double [m/s]
+%   Velocity values below cut_in produce 0 power
+% cut_out : double [m/s]
+%   Velocity values above cut_out produce 0 power
 %
 % Returns
-% -------
-%     p : Structure
+% ---------
+% p : struct
+%   Power data
+%     p.P : vector or matrix [W]
+%       Power, one value per velocity value
+%     p.time : vector [s]
+%       Time, with datetime converted to epoch seconds
 %
-%
-%        P.P: Power [W]
-%
-%        P.time: epoch time [s]
-%
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 arguments (Input)
     V struct
@@ -43,6 +39,7 @@
     cut_in (1,1) {mustBeNumeric}
     cut_out (1,1) {mustBeNumeric}
 end
+
 arguments (Output)
     p struct
 end
@@ -63,4 +60,3 @@
 p.time = time;
 
 end
-

From feec686e20d5ad50bf08991f54af7ea46b48d135 Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:17:02 -0600
Subject: [PATCH 05/10] River: Modify `energy_produced` to match MHKiT-Python
 implementation

---
 mhkit/river/resource/energy_produced.m | 35 ++++++++++++++------------
 1 file changed, 19 insertions(+), 16 deletions(-)

diff --git a/mhkit/river/resource/energy_produced.m b/mhkit/river/resource/energy_produced.m
index f8142ff9a..167eb9054 100644
--- a/mhkit/river/resource/energy_produced.m
+++ b/mhkit/river/resource/energy_produced.m
@@ -30,32 +30,36 @@
     P struct
     seconds (1,1) {mustBeNumeric}
 end
+
 arguments (Output)
     E (1,1) {mustBeNumeric}
 end
 
 power_data = P.P(:);
 
-% Histogram of power with 100 equal-width bins spanning the data range
-[counts, edges] = histcounts(power_data, 100);
-bin_widths = diff(edges);
-total_count = sum(counts);
+% 100 equal-width bins spanning the data range, matching numpy.histogram,
+% which widens a zero-width range by 0.5 on either side
+lower = min(power_data);
+upper = max(power_data);
+if lower == upper
+    lower = lower - 0.5;
+    upper = upper + 0.5;
+end
+edges = linspace(lower, upper, 101);
+counts = histcounts(power_data, edges);
 
 % Piecewise-constant probability density function of the histogram
-density = counts ./ (total_count * bin_widths);
+density = counts ./ (sum(counts) * diff(edges));
 
+% Evaluate the pdf like scipy.stats.rv_histogram, which assigns x to the
+% bin on its right (searchsorted side='right'), so the pdf is 0 at the
+% upper edge
 x = linspace(edges(1), edges(end), 1000);
+bin = discretize(x, edges);
+bin(x >= edges(end)) = NaN;
+in_range = ~isnan(bin);
 pdf_x = zeros(size(x));
-
-nbins = numel(counts);
-for i = 1:nbins
-    if i < nbins
-        mask = x >= edges(i) & x < edges(i + 1);
-    else
-        mask = x >= edges(i) & x <= edges(i + 1);
-    end
-    pdf_x(mask) = density(i);
-end
+pdf_x(in_range) = density(bin(in_range));
 
 % Expected value of power via trapezoidal integration of x*pdf(x)
 expected_power = trapz(x, x .* pdf_x);
@@ -63,4 +67,3 @@
 E = seconds * expected_power;
 
 end
-

From f13c8b44ae87e8ffe791f27e5ea60c6b0259fa1f Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:20:49 -0600
Subject: [PATCH 06/10] River: Standardize `exceedance_probability` user input
 to column vector

---
 mhkit/river/resource/exceedance_probability.m | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/mhkit/river/resource/exceedance_probability.m b/mhkit/river/resource/exceedance_probability.m
index 56e6cd49a..9ab298456 100644
--- a/mhkit/river/resource/exceedance_probability.m
+++ b/mhkit/river/resource/exceedance_probability.m
@@ -40,7 +40,10 @@
     time = posixtime(time);
 end
 
-discharge = Q.Discharge;
+% A row vector is a single timeseries, rank it as a column
+[discharge, was_row] = mhkit_standardize_user_input_to_column_vectors( ...
+    Q.Discharge, 'function_name', mfilename);
+
 n = size(discharge, 1);
 
 rank_ascending = zeros(size(discharge));
@@ -51,7 +54,9 @@
 % Convert to descending rank so the smallest value has the highest
 % exceedance probability
 rank_descending = n - rank_ascending + 1;
-F.F = 100 * rank_descending / (n + 1);
+exceedance = 100 * rank_descending / (n + 1);
+
+F.F = mhkit_restore_column_vectors_to_user_input(exceedance, was_row);
 F.time = time;
 
 end
@@ -75,4 +80,3 @@
     i = j + 1;
 end
 end
-

From 2f1c5def53660550f89c0e727e1c49d5906514be Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:21:25 -0600
Subject: [PATCH 07/10] Tests: River, Froude_number -> froude_number

---
 mhkit/tests/River_TestResource.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mhkit/tests/River_TestResource.m b/mhkit/tests/River_TestResource.m
index f2ca44f9d..a9cb182a1 100644
--- a/mhkit/tests/River_TestResource.m
+++ b/mhkit/tests/River_TestResource.m
@@ -2,10 +2,10 @@
 
     methods (Test)
 
-        function test_Froude_number(testCase)
+        function test_froude_number(testCase)
             v = 2;
             h = 5;
-            Fr = Froude_number(v, h);
+            Fr = froude_number(v, h);
             assertEqual(testCase,Fr, 0.286,'AbsTol',0.001);
         end
 

From e3107274009d181f5e908f06ee7cca0a644ccdbc Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:21:55 -0600
Subject: [PATCH 08/10] Tests: River, increase precision of
 test_energy_produced

---
 mhkit/tests/River_TestResource.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mhkit/tests/River_TestResource.m b/mhkit/tests/River_TestResource.m
index a9cb182a1..a304a6f66 100644
--- a/mhkit/tests/River_TestResource.m
+++ b/mhkit/tests/River_TestResource.m
@@ -71,7 +71,7 @@ function test_energy_produced(testCase)
             seconds=1;
             P = struct('P',[X;X;X;X;X;X;X;X;X;X],'time',[0 1 2 3 4 5 6 7 8 9]);
             EP = energy_produced(P, seconds);
-            assertEqual(testCase,EP, X*seconds,'AbsTol',0.1);
+            assertEqual(testCase,EP, 1.006011016021025,'RelTol',1e-12);
             % for a normal distribution of Power EP = mean *seconds
             mu=5;
             sigma=1;
@@ -93,7 +93,7 @@ function test_energy_produced(testCase)
                 normrnd(mu,sigma)], ...
                 'time',[0 1 2 3 4 5 6 7 8 9]);
             EP2 = energy_produced(power_dist, seconds);
-            assertEqual(testCase,EP2, mu*seconds,'AbsTol',0.1);
+            assertEqual(testCase,EP2, 5.052582623149183,'RelTol',1e-12);
         end
 
         function test_plot_flow_duration_curve(testCase)

From ca26324a343c93cf909712c13a93d5f834de6a67 Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:24:15 -0600
Subject: [PATCH 09/10] Tests: Add river modules to MATLAB native test runners

---
 mhkit/tests/examples/Examples_TestMatlabOnly.m | 3 ++-
 mhkit/tests/runMatlabOnlyTests.m               | 3 +++
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/mhkit/tests/examples/Examples_TestMatlabOnly.m b/mhkit/tests/examples/Examples_TestMatlabOnly.m
index 5c2365e86..d42884608 100644
--- a/mhkit/tests/examples/Examples_TestMatlabOnly.m
+++ b/mhkit/tests/examples/Examples_TestMatlabOnly.m
@@ -3,7 +3,7 @@
     %
     % Examples not included because they call Python bridge code:
     %   *_contour*, extreme_response_full_sea_state, short_term_extremes,
-    %   loads, power, qc, river, tidal, SWAN, delft3d, wecsim_example
+    %   loads, power, qc, tidal, SWAN, delft3d, wecsim_example
     % Examples not included because they make live network requests:
     %   WPTO_hindcast
 
@@ -15,6 +15,7 @@
             'cdip_example', ...
             'mooring_example', ...
             'power_example', ...
+            'river_example', ...
             'strain_measurement_example', ...
             'upcrossing_example', ...
             'wave_example', ...
diff --git a/mhkit/tests/runMatlabOnlyTests.m b/mhkit/tests/runMatlabOnlyTests.m
index 4b2b99df6..4b9d55229 100644
--- a/mhkit/tests/runMatlabOnlyTests.m
+++ b/mhkit/tests/runMatlabOnlyTests.m
@@ -24,6 +24,9 @@
         'Dolfyn_Test_VAP.m', ...
         'Mooring_TestMooring.m', ...
         'Power_TestDevice.m', ...
+        'River_TestDevice.m', ...
+        'River_TestIO.m', ...
+        'River_TestResource.m', ...
         'SignalProcessing_TestHilbert.m', ...
         'Tidal_TestIO.m', ...
         'upcrossing_Test.m', ...

From 87a7601caf68627965e47a89cb256bf9f6079465 Mon Sep 17 00:00:00 2001
From: "Simms, Andrew" <andrew.simms@nlr.gov>
Date: Wed, 30 Sep 2026 16:37:17 -0600
Subject: [PATCH 10/10] River: Add retries and timeouts to request_usgs_data

---
 mhkit/river/IO/request_usgs_data.m | 82 ++++++++++++++++++++++++------
 1 file changed, 67 insertions(+), 15 deletions(-)

diff --git a/mhkit/river/IO/request_usgs_data.m b/mhkit/river/IO/request_usgs_data.m
index fd4568d0d..a9af7895e 100644
--- a/mhkit/river/IO/request_usgs_data.m
+++ b/mhkit/river/IO/request_usgs_data.m
@@ -37,6 +37,10 @@
 %     write_json : str or None (optional)
 %         Name of json file to write data
 %         to call: request_usgs_data(station,parameter,start_date,end_date,"write_json",write_json)
+%
+%     timeout : double (optional)
+%         Timeout in seconds for each HTTP request. Default = 30
+%         to call: request_usgs_data(station,parameter,start_date,end_date,"timeout",timeout)
 %         
 % Returns
 % -------
@@ -58,10 +62,9 @@
     end_date string
     options.data_type = 'Daily';
     options.write_json string = "";
+    options.timeout (1,1) {mustBeNumeric, mustBePositive} = 30;
 end
 
-MAX_RETRIES = 5;                         % number of query retries if error
-
 % Formulate query
 if strcmp(options.data_type, 'Daily')
     data_url = "https://waterservices.usgs.gov/nwis/dv";
@@ -84,19 +87,8 @@
 disp("Data request URL: " + data_url + api_query)
 
 % Submit query and get data
-for i = 0:MAX_RETRIES
-    try
-        response = webread(data_url + api_query);
-        break;
-    catch ME
-        if i == MAX_RETRIES
-            disp(['MATLAB:request_usgs_data: ', ME.identifier]);
-            rethrow(ME)
-        else
-            pause(1);   % pause(seconds) and retry query
-        end
-    end
-end
+web_options = weboptions('Timeout', options.timeout);
+response = webread_with_retry(data_url + api_query, web_options);
 
 % Parse data label and units
 v = response.value.timeSeries.variable.variableDescription;
@@ -131,3 +123,63 @@
     fid = fopen(options.write_json, 'w');
     fprintf(fid, jsonencode(response));
 end
+
+end
+
+
+function response = webread_with_retry(url, web_options)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% Run a USGS web request, retrying transient failures with exponential
+% backoff and jitter
+%
+% Server errors (HTTP 5xx), rate limiting (HTTP 429), request timeouts
+% (HTTP 408), and connection failures are transient, so the request is
+% retried after 1, 2, 4, 8, 16, and 32 seconds, each scaled by a random
+% factor between 0.5 and 1.5 so repeated clients do not retry in
+% lockstep. Other client errors (HTTP 4xx, e.g. an invalid station or
+% parameter) are not transient and are raised immediately.
+%
+% Parameters
+% ------------
+%     url : string
+%         USGS request URL
+%     web_options : weboptions
+%         Options passed to webread
+%
+% Returns
+% ---------
+%     response : struct
+%         Decoded JSON response
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+MAX_RETRIES = 6;
+BASE_DELAY = 1;   % seconds, doubled after each failed attempt
+
+for attempt = 0:MAX_RETRIES
+    try
+        response = webread(url, web_options);
+        return
+    catch ME
+        status = regexp(ME.identifier, 'HTTP(\d+)StatusCodeError', 'tokens', 'once');
+        if isempty(status)
+            is_transient = true;   % timeout or connection failure
+        else
+            status = str2double(status{1});
+            is_transient = status >= 500 || status == 429 || status == 408;
+        end
+
+        if ~is_transient || attempt == MAX_RETRIES
+            disp(['MATLAB:request_usgs_data: ', ME.identifier]);
+            rethrow(ME)
+        end
+
+        delay = BASE_DELAY * 2^attempt * (0.5 + rand());
+        fprintf('USGS request failed (%s), retry %d of %d in %.1f s\n', ...
+            ME.identifier, attempt + 1, MAX_RETRIES, delay);
+        pause(delay);
+    end
+end
+end