Loads: Be explicit about how column vs row vectors are handled/validated

Author: simmsa

mhkit/loads/extreme/mler_coefficients.m

@@ -34,16 +34,21 @@
     error('ERROR: response_desired must be an int or double')
 end
 
+% validate all inputs have same length
+N = length(RAO);
+if length(wave_spectrum.frequency) ~= N || length(wave_spectrum.spectrum) ~= N
+    error('MHKiT:loads:mler_coefficients: RAO, frequency, and spectrum must have same length');
+end
+
 % convert from Hz to rad/s
-freq = wave_spectrum.frequency * (2*pi);
-freq_hz = wave_spectrum.frequency;
-wave_spectrum = wave_spectrum.spectrum / (2*pi);
-dw = (2*pi - 0) / (length(freq)-1);
+freq_rad = wave_spectrum.frequency * (2*pi);
+wave_spectrum_rad = wave_spectrum.spectrum / (2*pi);
+dw = (2*pi - 0) / (N-1);
 
 % response spectrum
-R.spectrum = abs(RAO).^2 .* (2*wave_spectrum);
+R.spectrum = abs(RAO).^2 .* (2*wave_spectrum_rad);
 R.type = 'response';
-R.frequency = freq;
+R.frequency = freq_rad;
 
 % spectral moment calculations
 m0 = frequency_moment(R, 0);
@@ -52,16 +57,16 @@
 wBar = m1/m0;
 
 % calculate coefficient_a from Quon2016 Eqn.8
-coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum) .* ((m2 - freq.*m1) + wBar.*(freq.*m0 - m1)) ./ (m0*m2 - m1^2);
+coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum_rad) .* ((m2 - freq_rad.*m1) + wBar.*(freq_rad.*m0 - m1)) ./ (m0*m2 - m1^2);
 % phase delay should be positive number
-phase = unwrap(angle(RAO));
+phase = -unwrap(angle(RAO));
 % for negative values of Amp, add pi phase shift, flip sign
 phase(coeff_a_rn < 0) = phase(coeff_a_rn < 0) - pi;
 coeff_a_rn(coeff_a_rn < 0) = coeff_a_rn(coeff_a_rn < 0) * -1;
 
 % calculate conditioned spectrum [m^2-s/rad]
-S = wave_spectrum .* coeff_a_rn.^2 .* response_desired^2;
-S(isnan(S)) = 0; % replace nans with zero
+conditioned_spectrum = wave_spectrum_rad .* coeff_a_rn.^2 .* response_desired^2;
+conditioned_spectrum(isnan(conditioned_spectrum)) = 0; % replace nans with zero
 % if the response amplitude we ask for is negative, we will add
 % a pi phase shift to the phase information.  This is because
 % the sign of response_desired is lost in the squaring above.
@@ -73,8 +78,8 @@
 end
 
 % outputs
-mler.conditioned_spectrum = S;
+mler.conditioned_spectrum = conditioned_spectrum;
 mler.phase = phase;
-mler.frequency = freq_hz;
+mler.frequency = wave_spectrum.frequency;
 
 end