From 1e405d4a02b912fa4d32fbce5e130385be45a234 Mon Sep 17 00:00:00 2001
From: MShabara <mshabara@nrel.gov>
Date: Tue, 29 Jul 2025 19:58:56 -0600
Subject: [PATCH 1/9] modifies the optimalGainCalc for the passive case

---
 Controls/Passive (P)/optimalGainCalc.m | 62 ++++++++++++++++----------
 1 file changed, 38 insertions(+), 24 deletions(-)

diff --git a/Controls/Passive (P)/optimalGainCalc.m b/Controls/Passive (P)/optimalGainCalc.m
index f51f09a5..230e6fe8 100644
--- a/Controls/Passive (P)/optimalGainCalc.m	
+++ b/Controls/Passive (P)/optimalGainCalc.m	
@@ -4,9 +4,10 @@
 
 close all; clear all; clc;
 
+dof = 3;            % Caluclate for heave motion
 % Inputs (from wecSimInputFile)
 simu = simulationClass();
-body(1) = bodyClass('../hydroData/sphere.h5');
+body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
 waves.height = 2.5;
 waves.period = 9.6664; % One of periods from BEM
 
@@ -19,20 +20,30 @@
 T = waves.period;
 omega = (1/T)*(2*pi);
 
+% Extend the freq vector to include the wave frequency
+hydro.simulation_parameters.w_extended = sort([hydro.simulation_parameters.w omega]);
+omegaIndex = find(hydro.simulation_parameters.w_extended == omega, 1, 'first');
+
 % Define excitation force based on wave conditions
-ampSpect = zeros(length(hydro.simulation_parameters.w),1);
-[~,closestIndOmega] = min(abs(omega-hydro.simulation_parameters.w));
-ampSpect(closestIndOmega) = A;
-FeRao = squeeze(hydro.hydro_coeffs.excitation.re(3,:))'*simu.rho*simu.gravity + squeeze(hydro.hydro_coeffs.excitation.im(3,:))'*simu.rho*simu.gravity*1j;
-Fexc = ampSpect.*FeRao;
+ampSpect = zeros(length(hydro.simulation_parameters.w_extended),1);
+ampSpect(omegaIndex) = A;
+Fe_re = squeeze(hydro.hydro_coeffs.excitation.re(dof, 1, :)) * simu.rho * simu.gravity;
+Fe_im = squeeze(hydro.hydro_coeffs.excitation.im(dof, 1, :)) * simu.rho * simu.gravity;
+
+Fe_interp = interp1(hydro.simulation_parameters.w, Fe_re + 1j * Fe_im, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+Fexc = ampSpect.*Fe_interp;
 
 % Define the intrinsic mechanical impedance for the device
-mass = simu.rho*hydro.properties.volume;
-addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(3,3,:))*simu.rho;
-radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(3,3,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
-hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(3,3)*simu.rho*simu.gravity;
-Gi = -((hydro.simulation_parameters.w)'.^2.*(mass+addedMass)) + 1j*hydro.simulation_parameters.w'.*radiationDamping + hydrostaticStiffness;
-Zi = Gi./(1j*hydro.simulation_parameters.w');
+mass = simu.rho * hydro.properties.volume;
+addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(dof, dof, :)) * simu.rho;
+addedMass = interp1(hydro.simulation_parameters.w, addedMass, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+
+radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(dof,dof,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
+radiationDamping = interp1(hydro.simulation_parameters.w, radiationDamping, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+
+hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(dof, dof) * simu.rho * simu.gravity;
+Gi = -((hydro.simulation_parameters.w_extended)'.^2 .* (mass+addedMass)) + 1j * hydro.simulation_parameters.w_extended'.*radiationDamping + hydrostaticStiffness;
+Zi = Gi./(1j*hydro.simulation_parameters.w_extended');
 
 % Calculate magnitude and phase for bode plot
 Mag = 20*log10(abs(Zi));
@@ -40,36 +51,39 @@
 
 % Determine natural frequency based on the phase of the impedance
 [~,closestIndNat] = min(abs(imag(Zi)));
-natFreq = hydro.simulation_parameters.w(closestIndNat);
+natFreq = hydro.simulation_parameters.w_extended(closestIndNat);
 T0 = (2*pi)/natFreq;
 
 % Create bode plot for impedance
 figure()
 subplot(2,1,1)
-semilogx((hydro.simulation_parameters.w)/(2*pi),Mag)
-xlabel('freq (hz)')
-ylabel('mag (dB)')
+semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Mag)
+xlabel('freq (hz)','interpreter','latex')
+ylabel('mag (dB)','interpreter','latex')
 grid on
 xline(natFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','southwest')
+legend('','Natural Frequency','Wave Frequency','Location','southwest','interpreter','latex')
 
 subplot(2,1,2)
-semilogx((hydro.simulation_parameters.w)/(2*pi),Phase)
-xlabel('freq (hz)')
-ylabel('phase (deg)')
+semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
+xlabel('freq (hz)','interpreter','latex')
+ylabel('phase (deg)','interpreter','latex')
 grid on
 xline(natFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','northwest')
+legend('','Natural Frequency','Wave Frequency','Location','northwest','interpreter','latex')
 
 % Calculate the maximum potential power
-P_max = -sum(abs(Fexc).^2./(8*real(Zi)))
+P_max = -sum(abs(Fexc).^2./(8*real(Zi)));
+fprintf('Maximum potential power P_max = %f\n', P_max);
 
 % Optimal proportional gain for passive control:
-KpOpt = sqrt(radiationDamping(closestIndOmega)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(closestIndOmega)))^2)
+KpOpt = sqrt(radiationDamping(omegaIndex)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(omegaIndex)))^2)
 Ki = 0;
+fprintf('Optimal proportional gain for passive control KpOpt = %f\n', KpOpt);
 
 % Calculate expected power with optimal passive control
 Zpto = KpOpt + Ki/(1j*omega);
-P = -sum(0.5*real(Zpto).*((abs(Fexc)).^2./(abs(Zpto+Zi)).^2))
+P = -sum(0.5*real(Zpto).*((abs(Fexc)).^2./(abs(Zpto+Zi)).^2));
+fprintf('Expected power with optimal passive control P = %f\n', P);

From 82ad59ba6f65146b931a6e799385f99d559a6359 Mon Sep 17 00:00:00 2001
From: MShabara <mshabara@nrel.gov>
Date: Tue, 29 Jul 2025 20:05:52 -0600
Subject: [PATCH 2/9] modifies the optimalGainCalc for the reactive case

---
 Controls/Reactive (PI)/optimalGainCalc.m | 69 +++++++++++++++---------
 1 file changed, 43 insertions(+), 26 deletions(-)

diff --git a/Controls/Reactive (PI)/optimalGainCalc.m b/Controls/Reactive (PI)/optimalGainCalc.m
index 28bad950..d24909c4 100644
--- a/Controls/Reactive (PI)/optimalGainCalc.m	
+++ b/Controls/Reactive (PI)/optimalGainCalc.m	
@@ -1,10 +1,13 @@
 % This script identifies the dynamics of the float in the respective wave 
 % conditions and determines the optimal proportional gain value for a 
-% passive controller (for regular waves)
+% reactive controller (for regular waves)
+close all; clear all; clc;
+
+dof = 3;            % Caluclate for heave motion
 
 % Inputs (from wecSimInputFile)
 simu = simulationClass();
-body(1) = bodyClass('../hydroData/sphere.h5');
+body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
 waves.height = 2.5;
 waves.period = 9.6664;
 
@@ -17,20 +20,30 @@
 T = waves.period;
 omega = (1/T)*(2*pi);
 
+% Extend the freq vector to include the wave frequency
+hydro.simulation_parameters.w_extended = sort([hydro.simulation_parameters.w omega]);
+omegaIndex = find(hydro.simulation_parameters.w_extended == omega, 1, 'first');
+
 % Define excitation force based on wave conditions
-ampSpect = zeros(length(hydro.simulation_parameters.w),1);
-[~,closestIndOmega] = min(abs(omega-hydro.simulation_parameters.w));
-ampSpect(closestIndOmega) = A;
-FeRao = squeeze(hydro.hydro_coeffs.excitation.re(3,:))'*simu.rho*simu.gravity + squeeze(hydro.hydro_coeffs.excitation.im(3,:))'*simu.rho*simu.gravity*1j;
-Fexc = ampSpect.*FeRao;
+ampSpect = zeros(length(hydro.simulation_parameters.w_extended),1);
+ampSpect(omegaIndex) = A;
+Fe_re = squeeze(hydro.hydro_coeffs.excitation.re(dof, 1, :)) * simu.rho * simu.gravity;
+Fe_im = squeeze(hydro.hydro_coeffs.excitation.im(dof, 1, :)) * simu.rho * simu.gravity;
+
+Fe_interp = interp1(hydro.simulation_parameters.w, Fe_re + 1j * Fe_im, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+Fexc = ampSpect.*Fe_interp;
 
 % Define the intrinsic mechanical impedance for the device
-mass = simu.rho*hydro.properties.volume;
-addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(3,3,:))*simu.rho;
-radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(3,3,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
-hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(3,3)*simu.rho*simu.gravity;
-Gi = -((hydro.simulation_parameters.w)'.^2.*(mass+addedMass)) + 1j*hydro.simulation_parameters.w'.*radiationDamping + hydrostaticStiffness;
-Zi = Gi./(1j*hydro.simulation_parameters.w');
+mass = simu.rho * hydro.properties.volume;
+addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(dof, dof, :)) * simu.rho;
+addedMass = interp1(hydro.simulation_parameters.w, addedMass, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+
+radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(dof,dof,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
+radiationDamping = interp1(hydro.simulation_parameters.w, radiationDamping, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+
+hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(dof, dof) * simu.rho * simu.gravity;
+Gi = -((hydro.simulation_parameters.w_extended)'.^2 .* (mass+addedMass)) + 1j * hydro.simulation_parameters.w_extended'.*radiationDamping + hydrostaticStiffness;
+Zi = Gi./(1j*hydro.simulation_parameters.w_extended');
 
 % Calculate magnitude and phase for bode plot
 Mag = 20*log10(abs(Zi));
@@ -38,35 +51,39 @@
 
 % Determine natural frequency based on the phase of the impedance
 [~,closestIndNat] = min(abs(imag(Zi)));
-natFreq = hydro.simulation_parameters.w(closestIndNat);
+natFreq = hydro.simulation_parameters.w_extended(closestIndNat);
 T0 = (2*pi)/natFreq;
 
 % Create bode plot for impedance
 figure()
 subplot(2,1,1)
-semilogx((hydro.simulation_parameters.w)/(2*pi),Mag)
-xlabel('freq (hz)')
-ylabel('mag (dB)')
+semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Mag)
+xlabel('freq (hz)','interpreter','latex')
+ylabel('mag (dB)','interpreter','latex')
 grid on
 xline(natFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','southwest')
+legend('','Natural Frequency','Wave Frequency','Location','southwest','interpreter','latex')
 
 subplot(2,1,2)
-semilogx((hydro.simulation_parameters.w)/(2*pi),Phase)
-xlabel('freq (hz)')
-ylabel('phase (deg)')
+semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
+xlabel('freq (hz)','interpreter','latex')
+ylabel('phase (deg)','interpreter','latex')
 grid on
 xline(natFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','northwest')
+legend('','Natural Frequency','Wave Frequency','Location','northwest','interpreter','latex')
 
 % Calculate the maximum potential power
-P_max = -sum(abs(Fexc).^2./(8*real(Zi)))
+P_max = -sum(abs(Fexc).^2./(8*real(Zi)));
+fprintf('Maximum potential power P_max = %f\n', P_max);
 
 % Calculate optimal Kp and Ki gains for reactive control
-KpOpt = radiationDamping(closestIndOmega)
-KiOpt = -(-omega^2*(mass + addedMass(closestIndOmega)) + hydrostaticStiffness)
+KpOpt = radiationDamping(omegaIndex);
+KiOpt = -(-omega^2*(mass + addedMass(omegaIndex)) + hydrostaticStiffness);
+fprintf('Optimal gains for reactive control Kp = %f and Ki = %f\n', KpOpt, KiOpt);
 
+% Calculate expected power with optimal reactive control
 Zpto = KpOpt + KiOpt/(1j*omega);
-P = -sum(0.5*real(Zpto).*((abs(Fexc)).^2./(abs(Zpto+Zi)).^2))
+P = -sum(0.5*real(Zpto).*((abs(Fexc)).^2./(abs(Zpto+Zi)).^2));
+fprintf('Expected power with optimal passive control P = %f\n', P);

From 8274e105f086e0a1fc3439f86766596ecf026695 Mon Sep 17 00:00:00 2001
From: Mohamed Shabara <84589678+MShabara@users.noreply.github.com>
Date: Tue, 29 Jul 2025 20:35:26 -0600
Subject: [PATCH 3/9] Update optimalGainCalc.m

Adds semicolon to one of the lines
---
 Controls/Passive (P)/optimalGainCalc.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Controls/Passive (P)/optimalGainCalc.m b/Controls/Passive (P)/optimalGainCalc.m
index 230e6fe8..3767869e 100644
--- a/Controls/Passive (P)/optimalGainCalc.m	
+++ b/Controls/Passive (P)/optimalGainCalc.m	
@@ -79,7 +79,7 @@
 fprintf('Maximum potential power P_max = %f\n', P_max);
 
 % Optimal proportional gain for passive control:
-KpOpt = sqrt(radiationDamping(omegaIndex)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(omegaIndex)))^2)
+KpOpt = sqrt(radiationDamping(omegaIndex)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(omegaIndex)))^2);
 Ki = 0;
 fprintf('Optimal proportional gain for passive control KpOpt = %f\n', KpOpt);
 

From 9d1ba8d067be2dca0115e6aae004224bad0d9fe2 Mon Sep 17 00:00:00 2001
From: MShabara <mshabara@nrel.gov>
Date: Fri, 8 Aug 2025 21:22:00 -0600
Subject: [PATCH 4/9] modified the latching control gain calc

---
 Controls/Latching/optimalTimeCalc.m | 52 ++++++++++++++++++-----------
 1 file changed, 32 insertions(+), 20 deletions(-)

diff --git a/Controls/Latching/optimalTimeCalc.m b/Controls/Latching/optimalTimeCalc.m
index 21fd3b8a..e8c42553 100644
--- a/Controls/Latching/optimalTimeCalc.m
+++ b/Controls/Latching/optimalTimeCalc.m
@@ -6,33 +6,47 @@
 
 % Inputs (from wecSimInputFile)
 simu = simulationClass();
-body(1) = bodyClass('../hydroData/sphere.h5');
+body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
 waves.height = 2.5;
 waves.period = 9.6664;
 
+
 % Load hydrodynamic data for float from BEM
 hydro = readBEMIOH5(body.h5File{1}, 1, body.meanDrift);
 
+dof = 3;            % Caluclate for heave motion
+mass = simu.rho*hydro.properties.volume;
+
 % Define wave conditions
 H = waves.height;
 A = H/2;
 T = waves.period;
 omega = (1/T)*(2*pi);
 
+% Extend the freq vector to include the wave frequency
+hydro.simulation_parameters.w_extended = sort([hydro.simulation_parameters.w omega]);
+omegaIndex = find(hydro.simulation_parameters.w_extended == omega, 1, 'first');
+
 % Define excitation force based on wave conditions
-ampSpect = zeros(length(hydro.simulation_parameters.w),1);
-[~,closestIndOmega] = min(abs(omega-hydro.simulation_parameters.w));
-ampSpect(closestIndOmega) = A;
-FeRao = squeeze(hydro.hydro_coeffs.excitation.re(3,:))'*simu.rho*simu.gravity + squeeze(hydro.hydro_coeffs.excitation.im(3,:))'*simu.rho*simu.gravity*1j;
-Fexc = ampSpect.*FeRao;
+ampSpect = zeros(length(hydro.simulation_parameters.w_extended),1);
+[~,closestIndOmega] = min(abs(omega-hydro.simulation_parameters.w_extended));
+ampSpect(omegaIndex) = A;
+Fe_re = squeeze(hydro.hydro_coeffs.excitation.re(dof, 1, :)) * simu.rho * simu.gravity;
+Fe_im = squeeze(hydro.hydro_coeffs.excitation.im(dof, 1, :)) * simu.rho * simu.gravity;
+
+Fe_interp = interp1(hydro.simulation_parameters.w, Fe_re + 1j * Fe_im, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+Fexc = ampSpect.*Fe_interp;
 
 % Define the intrinsic mechanical impedance for the device
-mass = simu.rho*hydro.properties.volume;
-addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(3,3,:))*simu.rho;
-radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(3,3,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
-hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(3,3)*simu.rho*simu.gravity;
-Gi = -((hydro.simulation_parameters.w)'.^2.*(mass+addedMass)) + 1j*hydro.simulation_parameters.w'.*radiationDamping + hydrostaticStiffness;
-Zi = Gi./(1j*hydro.simulation_parameters.w');
+addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(dof, dof, :)) * simu.rho;
+addedMass = interp1(hydro.simulation_parameters.w, addedMass, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+
+radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(dof,dof,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
+radiationDamping = interp1(hydro.simulation_parameters.w, radiationDamping, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+
+hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(dof, dof) * simu.rho * simu.gravity;
+Gi = -((hydro.simulation_parameters.w_extended)'.^2 .* (mass+addedMass)) + 1j * hydro.simulation_parameters.w_extended'.*radiationDamping + hydrostaticStiffness;
+Zi = Gi./(1j*hydro.simulation_parameters.w_extended');
 
 % Calculate magnitude and phase for bode plot
 Mag = 20*log10(abs(Zi));
@@ -40,13 +54,13 @@
 
 % Determine natural frequency based on the phase of the impedance
 [~,closestIndNat] = min(abs(imag(Zi)));
-natFreq = hydro.simulation_parameters.w(closestIndNat);
+natFreq = hydro.simulation_parameters.w_extended(closestIndNat);
 natT = (2*pi)/natFreq;
 
 % Create bode plot for impedance
 figure()
 subplot(2,1,1)
-semilogx((hydro.simulation_parameters.w)/(2*pi),Mag)
+semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Mag)
 xlabel('freq (hz)')
 ylabel('mag (dB)')
 grid on
@@ -55,7 +69,7 @@
 legend('','Natural Frequency','Wave Frequency','Location','southwest')
 
 subplot(2,1,2)
-semilogx((hydro.simulation_parameters.w)/(2*pi),Phase)
+semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
 xlabel('freq (hz)')
 ylabel('phase (deg)')
 grid on
@@ -65,10 +79,8 @@
 
 % Determine optimal latching time
 optLatchTime = 0.5*(T - natT)
-KpOpt = radiationDamping(closestIndOmega)
-force = 80*(mass+addedMass(closestIndOmega))
+KpOpt = radiationDamping(omegaIndex)
+force = 80*(mass+addedMass(omegaIndex))
 
 % Calculate the maximum potential power
-P_max = -sum(abs(Fexc).^2./(8*real(Zi)))
-
-
+P_max = -sum(abs(Fexc).^2./(8*real(Zi)))
\ No newline at end of file

From 51c7372b80e72a22bbc72db44ad87c40365c68e2 Mon Sep 17 00:00:00 2001
From: MShabara <mshabara@nrel.gov>
Date: Fri, 8 Aug 2025 21:22:57 -0600
Subject: [PATCH 5/9] corrects the function description

---
 Controls/Latching/optimalTimeCalc.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Controls/Latching/optimalTimeCalc.m b/Controls/Latching/optimalTimeCalc.m
index e8c42553..f1221a3f 100644
--- a/Controls/Latching/optimalTimeCalc.m
+++ b/Controls/Latching/optimalTimeCalc.m
@@ -1,6 +1,6 @@
 % This script identifies the dynamics of the float in the respective wave 
-% conditions and determines the optimal proportional gain value for a 
-% passive controller (for regular waves)
+% conditions and determines the optimal proportional gain and latching time
+% value for a regular wave
 
 close all; clear all; clc;
 

From f38de3259b8d54c01a0ee0d9b5b5b525abfd59f6 Mon Sep 17 00:00:00 2001
From: MShabara <mshabara@nrel.gov>
Date: Fri, 8 Aug 2025 21:30:32 -0600
Subject: [PATCH 6/9] modified the declutching control gain calc

---
 Controls/Declutching/optimalTimeCalc.m | 45 ++++++++++++++++----------
 1 file changed, 28 insertions(+), 17 deletions(-)

diff --git a/Controls/Declutching/optimalTimeCalc.m b/Controls/Declutching/optimalTimeCalc.m
index 66baa166..ccb803e8 100644
--- a/Controls/Declutching/optimalTimeCalc.m
+++ b/Controls/Declutching/optimalTimeCalc.m
@@ -6,9 +6,10 @@
 
 % Inputs (from wecSimInputFile)
 simu = simulationClass();
-body(1) = bodyClass('../hydroData/sphere.h5');
-waves.height = 1;
-waves.period = 3.5;
+body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
+waves.height = 2.5;                          % Wave Height [m]
+waves.period = 9.6664;                       % Wave Period [s]
+dof = 3;        % changing this value reguires changing the device mass
 
 % Load hydrodynamic data for float from BEM
 hydro = readBEMIOH5(body.h5File{1}, 1, body.meanDrift);
@@ -19,20 +20,30 @@
 T = waves.period;
 omega = (1/T)*(2*pi);
 
+% Extend the freq vector to include the wave frequency
+hydro.simulation_parameters.w_extended = sort([hydro.simulation_parameters.w omega]);
+omegaIndex = find(hydro.simulation_parameters.w_extended == omega, 1, 'first');
+
 % Define excitation force based on wave conditions
-ampSpect = zeros(length(hydro.simulation_parameters.w),1);
-[~,closestIndOmega] = min(abs(omega-hydro.simulation_parameters.w));
-ampSpect(closestIndOmega) = A;
-FeRao = squeeze(hydro.hydro_coeffs.excitation.re(3,:))'*simu.rho*simu.gravity + squeeze(hydro.hydro_coeffs.excitation.im(3,:))'*simu.rho*simu.gravity*1j;
-Fexc = ampSpect.*FeRao;
+ampSpect = zeros(length(hydro.simulation_parameters.w_extended),1);
+ampSpect(omegaIndex) = A;
+Fe_re = squeeze(hydro.hydro_coeffs.excitation.re(dof, 1, :)) * simu.rho * simu.gravity;
+Fe_im = squeeze(hydro.hydro_coeffs.excitation.im(dof, 1, :)) * simu.rho * simu.gravity;
+
+Fe_interp = interp1(hydro.simulation_parameters.w, Fe_re + 1j * Fe_im, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+Fexc = ampSpect.*Fe_interp;
 
 % Define the intrinsic mechanical impedance for the device
 mass = simu.rho*hydro.properties.volume;
-addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(3,3,:))*simu.rho;
-radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(3,3,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
-hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(3,3)*simu.rho*simu.gravity;
-Gi = -((hydro.simulation_parameters.w)'.^2.*(mass+addedMass)) + 1j*hydro.simulation_parameters.w'.*radiationDamping + hydrostaticStiffness;
-Zi = Gi./(1j*hydro.simulation_parameters.w');
+addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(dof, dof, :)) * simu.rho;
+addedMass = interp1(hydro.simulation_parameters.w, addedMass, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+
+radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(dof,dof,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
+radiationDamping = interp1(hydro.simulation_parameters.w, radiationDamping, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+
+hydrostaticStiffness = hydro.hydro_coeffs.linear_restoring_stiffness(dof, dof) * simu.rho * simu.gravity;
+Gi = -((hydro.simulation_parameters.w_extended)'.^2 .* (mass+addedMass)) + 1j * hydro.simulation_parameters.w_extended'.*radiationDamping + hydrostaticStiffness;
+Zi = Gi./(1j*hydro.simulation_parameters.w_extended');
 
 % Calculate magnitude and phase for bode plot
 Mag = 20*log10(abs(Zi));
@@ -40,13 +51,13 @@
 
 % Determine natural frequency based on the phase of the impedance
 [~,closestIndNat] = min(abs(imag(Zi)));
-natFreq = hydro.simulation_parameters.w(closestIndNat);
+natFreq = hydro.simulation_parameters.w_extended(closestIndNat);
 natT = (2*pi)/natFreq;
 
 % Create bode plot for impedance
 figure()
 subplot(2,1,1)
-semilogx((hydro.simulation_parameters.w)/(2*pi),Mag)
+semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Mag)
 xlabel('freq (hz)')
 ylabel('mag (dB)')
 grid on
@@ -55,7 +66,7 @@
 legend('','Natural Frequency','Wave Frequency','Location','southwest')
 
 subplot(2,1,2)
-semilogx((hydro.simulation_parameters.w)/(2*pi),Phase)
+semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
 xlabel('freq (hz)')
 ylabel('phase (deg)')
 grid on
@@ -65,7 +76,7 @@
 
 % Determine optimal latching time
 optDeclutchTime = 0.5*(natT - T)
-KpOpt = sqrt(radiationDamping(closestIndOmega)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(closestIndOmega)))^2)
+KpOpt = sqrt(radiationDamping(omegaIndex)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(omegaIndex)))^2)
 
 % Calculate the maximum potential power
 P_max = -sum(abs(Fexc).^2./(8*real(Zi)))

From da8d4ad8cbe625df594517a820ea5b516b751fa4 Mon Sep 17 00:00:00 2001
From: jtgrasb <87095491+jtgrasb@users.noreply.github.com>
Date: Wed, 10 Sep 2025 16:56:52 -0600
Subject: [PATCH 7/9] update frequency calculation

---
 Controls/Declutching/optimalTimeCalc.m   | 20 +++++++++-----------
 Controls/Latching/optimalTimeCalc.m      | 24 ++++++++++--------------
 Controls/Passive (P)/optimalGainCalc.m   | 18 ++++++++----------
 Controls/Reactive (PI)/optimalGainCalc.m | 16 +++++++---------
 4 files changed, 34 insertions(+), 44 deletions(-)

diff --git a/Controls/Declutching/optimalTimeCalc.m b/Controls/Declutching/optimalTimeCalc.m
index ccb803e8..775a4ad7 100644
--- a/Controls/Declutching/optimalTimeCalc.m
+++ b/Controls/Declutching/optimalTimeCalc.m
@@ -1,15 +1,14 @@
 % This script identifies the dynamics of the float in the respective wave 
 % conditions and determines the optimal proportional gain value for a 
 % passive controller (for regular waves)
-
 close all; clear all; clc;
 
+dof = 3;            % Caluclate for heave motion
 % Inputs (from wecSimInputFile)
 simu = simulationClass();
 body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
 waves.height = 2.5;                          % Wave Height [m]
 waves.period = 9.6664;                       % Wave Period [s]
-dof = 3;        % changing this value reguires changing the device mass
 
 % Load hydrodynamic data for float from BEM
 hydro = readBEMIOH5(body.h5File{1}, 1, body.meanDrift);
@@ -49,10 +48,9 @@
 Mag = 20*log10(abs(Zi));
 Phase = (angle(Zi))*(180/pi);
 
-% Determine natural frequency based on the phase of the impedance
-[~,closestIndNat] = min(abs(imag(Zi)));
-natFreq = hydro.simulation_parameters.w_extended(closestIndNat);
-natT = (2*pi)/natFreq;
+% Determine resonant frequency based on the phase of the impedance
+resonantFreq = interp1(Phase, hydro.simulation_parameters.w_extended, 0, 'spline','extrap');
+resonantPeriod = (2*pi)/resonantFreq;
 
 % Create bode plot for impedance
 figure()
@@ -61,21 +59,21 @@
 xlabel('freq (hz)')
 ylabel('mag (dB)')
 grid on
-xline(natFreq/(2*pi))
+xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','southwest')
+legend('','resonant Frequency','Wave Frequency','Location','southwest')
 
 subplot(2,1,2)
 semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
 xlabel('freq (hz)')
 ylabel('phase (deg)')
 grid on
-xline(natFreq/(2*pi))
+xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','northwest')
+legend('','resonant Frequency','Wave Frequency','Location','northwest')
 
 % Determine optimal latching time
-optDeclutchTime = 0.5*(natT - T)
+optDeclutchTime = 0.5*(resonantPeriod - T)
 KpOpt = sqrt(radiationDamping(omegaIndex)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(omegaIndex)))^2)
 
 % Calculate the maximum potential power
diff --git a/Controls/Latching/optimalTimeCalc.m b/Controls/Latching/optimalTimeCalc.m
index f1221a3f..b6fbc11d 100644
--- a/Controls/Latching/optimalTimeCalc.m
+++ b/Controls/Latching/optimalTimeCalc.m
@@ -1,22 +1,18 @@
 % This script identifies the dynamics of the float in the respective wave 
 % conditions and determines the optimal proportional gain and latching time
 % value for a regular wave
-
 close all; clear all; clc;
 
+dof = 3;            % Caluclate for heave motion
 % Inputs (from wecSimInputFile)
 simu = simulationClass();
 body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
 waves.height = 2.5;
 waves.period = 9.6664;
 
-
 % Load hydrodynamic data for float from BEM
 hydro = readBEMIOH5(body.h5File{1}, 1, body.meanDrift);
 
-dof = 3;            % Caluclate for heave motion
-mass = simu.rho*hydro.properties.volume;
-
 % Define wave conditions
 H = waves.height;
 A = H/2;
@@ -38,6 +34,7 @@
 Fexc = ampSpect.*Fe_interp;
 
 % Define the intrinsic mechanical impedance for the device
+mass = simu.rho*hydro.properties.volume;
 addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(dof, dof, :)) * simu.rho;
 addedMass = interp1(hydro.simulation_parameters.w, addedMass, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
 
@@ -52,10 +49,9 @@
 Mag = 20*log10(abs(Zi));
 Phase = (angle(Zi))*(180/pi);
 
-% Determine natural frequency based on the phase of the impedance
-[~,closestIndNat] = min(abs(imag(Zi)));
-natFreq = hydro.simulation_parameters.w_extended(closestIndNat);
-natT = (2*pi)/natFreq;
+% Determine resonant frequency based on the phase of the impedance
+resonantFreq = interp1(Phase, hydro.simulation_parameters.w_extended, 0, 'spline','extrap');
+resonantPeriod = (2*pi)/resonantFreq;
 
 % Create bode plot for impedance
 figure()
@@ -64,21 +60,21 @@
 xlabel('freq (hz)')
 ylabel('mag (dB)')
 grid on
-xline(natFreq/(2*pi))
+xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','southwest')
+legend('','resonant Frequency','Wave Frequency','Location','southwest')
 
 subplot(2,1,2)
 semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
 xlabel('freq (hz)')
 ylabel('phase (deg)')
 grid on
-xline(natFreq/(2*pi))
+xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','northwest')
+legend('','resonant Frequency','Wave Frequency','Location','northwest')
 
 % Determine optimal latching time
-optLatchTime = 0.5*(T - natT)
+optLatchTime = 0.5*(T - resonantPeriod)
 KpOpt = radiationDamping(omegaIndex)
 force = 80*(mass+addedMass(omegaIndex))
 
diff --git a/Controls/Passive (P)/optimalGainCalc.m b/Controls/Passive (P)/optimalGainCalc.m
index 3767869e..b98a4c21 100644
--- a/Controls/Passive (P)/optimalGainCalc.m	
+++ b/Controls/Passive (P)/optimalGainCalc.m	
@@ -1,7 +1,6 @@
 % This script identifies the dynamics of the float in the respective wave 
 % conditions and determines the optimal proportional gain value for a 
 % passive controller (for regular waves)
-
 close all; clear all; clc;
 
 dof = 3;            % Caluclate for heave motion
@@ -9,7 +8,7 @@
 simu = simulationClass();
 body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
 waves.height = 2.5;
-waves.period = 9.6664; % One of periods from BEM
+waves.period = 4.5; % One of periods from BEM
 
 % Load hydrodynamic data for float from BEM
 hydro = readBEMIOH5(body.h5File{1}, 1, body.meanDrift);
@@ -49,10 +48,9 @@
 Mag = 20*log10(abs(Zi));
 Phase = (angle(Zi))*(180/pi);
 
-% Determine natural frequency based on the phase of the impedance
-[~,closestIndNat] = min(abs(imag(Zi)));
-natFreq = hydro.simulation_parameters.w_extended(closestIndNat);
-T0 = (2*pi)/natFreq;
+% Determine resonant frequency based on the phase of the impedance
+resonantFreq = interp1(Phase, hydro.simulation_parameters.w_extended, 0, 'spline','extrap');
+resonantPeriod = (2*pi)/resonantFreq;
 
 % Create bode plot for impedance
 figure()
@@ -61,18 +59,18 @@
 xlabel('freq (hz)','interpreter','latex')
 ylabel('mag (dB)','interpreter','latex')
 grid on
-xline(natFreq/(2*pi))
+xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','southwest','interpreter','latex')
+legend('','Resonant Frequency','Wave Frequency','Location','southwest','interpreter','latex')
 
 subplot(2,1,2)
 semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
 xlabel('freq (hz)','interpreter','latex')
 ylabel('phase (deg)','interpreter','latex')
 grid on
-xline(natFreq/(2*pi))
+xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','northwest','interpreter','latex')
+legend('','Resonant Frequency','Wave Frequency','Location','northwest','interpreter','latex')
 
 % Calculate the maximum potential power
 P_max = -sum(abs(Fexc).^2./(8*real(Zi)));
diff --git a/Controls/Reactive (PI)/optimalGainCalc.m b/Controls/Reactive (PI)/optimalGainCalc.m
index d24909c4..6b56bd6f 100644
--- a/Controls/Reactive (PI)/optimalGainCalc.m	
+++ b/Controls/Reactive (PI)/optimalGainCalc.m	
@@ -4,7 +4,6 @@
 close all; clear all; clc;
 
 dof = 3;            % Caluclate for heave motion
-
 % Inputs (from wecSimInputFile)
 simu = simulationClass();
 body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
@@ -49,10 +48,9 @@
 Mag = 20*log10(abs(Zi));
 Phase = (angle(Zi))*(180/pi);
 
-% Determine natural frequency based on the phase of the impedance
-[~,closestIndNat] = min(abs(imag(Zi)));
-natFreq = hydro.simulation_parameters.w_extended(closestIndNat);
-T0 = (2*pi)/natFreq;
+% Determine resonant frequency based on the phase of the impedance
+resonantFreq = interp1(Phase, hydro.simulation_parameters.w_extended, 0, 'spline','extrap');
+resonantPeriod = (2*pi)/resonantFreq;
 
 % Create bode plot for impedance
 figure()
@@ -61,18 +59,18 @@
 xlabel('freq (hz)','interpreter','latex')
 ylabel('mag (dB)','interpreter','latex')
 grid on
-xline(natFreq/(2*pi))
+xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','southwest','interpreter','latex')
+legend('','resonant Frequency','Wave Frequency','Location','southwest','interpreter','latex')
 
 subplot(2,1,2)
 semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
 xlabel('freq (hz)','interpreter','latex')
 ylabel('phase (deg)','interpreter','latex')
 grid on
-xline(natFreq/(2*pi))
+xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','Natural Frequency','Wave Frequency','Location','northwest','interpreter','latex')
+legend('','resonant Frequency','Wave Frequency','Location','northwest','interpreter','latex')
 
 % Calculate the maximum potential power
 P_max = -sum(abs(Fexc).^2./(8*real(Zi)));

From 9592114d282b96bef8d85cef366e19eb8aa9164a Mon Sep 17 00:00:00 2001
From: jtgrasb <87095491+jtgrasb@users.noreply.github.com>
Date: Wed, 10 Sep 2025 16:58:31 -0600
Subject: [PATCH 8/9] revert period

---
 Controls/Passive (P)/optimalGainCalc.m | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Controls/Passive (P)/optimalGainCalc.m b/Controls/Passive (P)/optimalGainCalc.m
index b98a4c21..be96a613 100644
--- a/Controls/Passive (P)/optimalGainCalc.m	
+++ b/Controls/Passive (P)/optimalGainCalc.m	
@@ -8,7 +8,7 @@
 simu = simulationClass();
 body(1) = bodyClass('../../_Common_Input_Files/Sphere/hydroData/sphere.h5');
 waves.height = 2.5;
-waves.period = 4.5; % One of periods from BEM
+waves.period = 9.6664; % One of periods from BEM
 
 % Load hydrodynamic data for float from BEM
 hydro = readBEMIOH5(body.h5File{1}, 1, body.meanDrift);

From 2c173873881d6cfe49ec600a6686122ca7ffeb7a Mon Sep 17 00:00:00 2001
From: jtgrasb <87095491+jtgrasb@users.noreply.github.com>
Date: Wed, 17 Sep 2025 08:40:32 -0600
Subject: [PATCH 9/9] final cleanup

---
 Controls/Declutching/optimalTimeCalc.m           | 4 ++--
 Controls/Latching/optimalTimeCalc.m              | 4 ++--
 Controls/Passive (P)/optimalGainCalc.m           | 3 ++-
 Controls/Passive (P)/userDefinedFunctionsMCR.m   | 2 +-
 Controls/Reactive (PI)/optimalGainCalc.m         | 4 ++--
 Controls/Reactive (PI)/userDefinedFunctionsMCR.m | 2 +-
 6 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/Controls/Declutching/optimalTimeCalc.m b/Controls/Declutching/optimalTimeCalc.m
index 775a4ad7..f3285da0 100644
--- a/Controls/Declutching/optimalTimeCalc.m
+++ b/Controls/Declutching/optimalTimeCalc.m
@@ -61,7 +61,7 @@
 grid on
 xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','resonant Frequency','Wave Frequency','Location','southwest')
+legend('','Resonant Frequency','Wave Frequency','Location','southwest')
 
 subplot(2,1,2)
 semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
@@ -70,7 +70,7 @@
 grid on
 xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','resonant Frequency','Wave Frequency','Location','northwest')
+legend('','Resonant Frequency','Wave Frequency','Location','northwest')
 
 % Determine optimal latching time
 optDeclutchTime = 0.5*(resonantPeriod - T)
diff --git a/Controls/Latching/optimalTimeCalc.m b/Controls/Latching/optimalTimeCalc.m
index b6fbc11d..cb611636 100644
--- a/Controls/Latching/optimalTimeCalc.m
+++ b/Controls/Latching/optimalTimeCalc.m
@@ -62,7 +62,7 @@
 grid on
 xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','resonant Frequency','Wave Frequency','Location','southwest')
+legend('','Resonant Frequency','Wave Frequency','Location','southwest')
 
 subplot(2,1,2)
 semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
@@ -71,7 +71,7 @@
 grid on
 xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','resonant Frequency','Wave Frequency','Location','northwest')
+legend('','Resonant Frequency','Wave Frequency','Location','northwest')
 
 % Determine optimal latching time
 optLatchTime = 0.5*(T - resonantPeriod)
diff --git a/Controls/Passive (P)/optimalGainCalc.m b/Controls/Passive (P)/optimalGainCalc.m
index be96a613..48c6f8b5 100644
--- a/Controls/Passive (P)/optimalGainCalc.m	
+++ b/Controls/Passive (P)/optimalGainCalc.m	
@@ -36,6 +36,7 @@
 mass = simu.rho * hydro.properties.volume;
 addedMass = squeeze(hydro.hydro_coeffs.added_mass.all(dof, dof, :)) * simu.rho;
 addedMass = interp1(hydro.simulation_parameters.w, addedMass, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
+addedMass = squeeze(hydro.hydro_coeffs.added_mass.inf_freq(dof, dof, :)) * simu.rho;
 
 radiationDamping = squeeze(hydro.hydro_coeffs.radiation_damping.all(dof,dof,:)).*squeeze(hydro.simulation_parameters.w')*simu.rho;
 radiationDamping = interp1(hydro.simulation_parameters.w, radiationDamping, hydro.simulation_parameters.w_extended, 'spline', 'extrap')';
@@ -77,7 +78,7 @@
 fprintf('Maximum potential power P_max = %f\n', P_max);
 
 % Optimal proportional gain for passive control:
-KpOpt = sqrt(radiationDamping(omegaIndex)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass(omegaIndex)))^2);
+KpOpt = sqrt(radiationDamping(omegaIndex)^2 + ((hydrostaticStiffness/omega) - omega*(mass + addedMass))^2);
 Ki = 0;
 fprintf('Optimal proportional gain for passive control KpOpt = %f\n', KpOpt);
 
diff --git a/Controls/Passive (P)/userDefinedFunctionsMCR.m b/Controls/Passive (P)/userDefinedFunctionsMCR.m
index 2eb7f4a8..3cb8481e 100644
--- a/Controls/Passive (P)/userDefinedFunctionsMCR.m	
+++ b/Controls/Passive (P)/userDefinedFunctionsMCR.m	
@@ -17,7 +17,7 @@
 mcr.maxPower(imcr) = max(controllersOutput.power(startInd:endInd,3));
 mcr.maxForce(imcr) = max(controllersOutput.force(startInd:endInd,3));
 
-if imcr == 9
+if imcr == length(mcr.cases)
     figure()
     plot(mcr.cases,mcr.meanPower)
     title('Mean Power vs. Proportional Gain')
diff --git a/Controls/Reactive (PI)/optimalGainCalc.m b/Controls/Reactive (PI)/optimalGainCalc.m
index 6b56bd6f..18bfab4c 100644
--- a/Controls/Reactive (PI)/optimalGainCalc.m	
+++ b/Controls/Reactive (PI)/optimalGainCalc.m	
@@ -61,7 +61,7 @@
 grid on
 xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','resonant Frequency','Wave Frequency','Location','southwest','interpreter','latex')
+legend('','Resonant Frequency','Wave Frequency','Location','southwest','interpreter','latex')
 
 subplot(2,1,2)
 semilogx((hydro.simulation_parameters.w_extended)/(2*pi),Phase)
@@ -70,7 +70,7 @@
 grid on
 xline(resonantFreq/(2*pi))
 xline(1/T,'--')
-legend('','resonant Frequency','Wave Frequency','Location','northwest','interpreter','latex')
+legend('','Resonant Frequency','Wave Frequency','Location','northwest','interpreter','latex')
 
 % Calculate the maximum potential power
 P_max = -sum(abs(Fexc).^2./(8*real(Zi)));
diff --git a/Controls/Reactive (PI)/userDefinedFunctionsMCR.m b/Controls/Reactive (PI)/userDefinedFunctionsMCR.m
index 92c62f69..ba79650f 100644
--- a/Controls/Reactive (PI)/userDefinedFunctionsMCR.m	
+++ b/Controls/Reactive (PI)/userDefinedFunctionsMCR.m	
@@ -17,7 +17,7 @@
 mcr.maxPower(imcr) = max(controllersOutput.power(startInd:endInd,3));
 mcr.maxForce(imcr) = max(controllersOutput.force(startInd:endInd,3));
 
-if imcr == 81
+if imcr == length(mcr.cases)
 
     % Kp and Ki gains
     kps = unique(mcr.cases(:,1));