Full dir spectra (resolved conflicts)

docs/theory/theory.rst

@@ -512,7 +512,7 @@ formulation of :math:`D\left( \theta \right)` requires that
                                     \int_{0}^{\infty} S\left( f \right) df
 
 so that the total energy in the directional spectrum must be the same as the 
-total energy in the one-dimensional spectrum. 
+total energy in the one-dimensional spectrum. An example one-dimensional spectrum is.
 
 .. math::
 
@@ -532,7 +532,33 @@ significant wave height :math:`H_{m0}` (in m), as:
 
 .. math::
     H_{m0} = 4 \sqrt{m_{0}}~~
+
+Frequency-varying direction and spread
+""""""""""""""""""""""""""""""""""""""""""""""
+
+A more general description of a wave field can be expressed
+
+.. math::
+	S(\omega,\theta)=S(\omega)D(\omega,\theta)
 
+where :math:`D(\omega,\theta)` indicates that the directional components of a wave can also
+vary with frequency. The total energy in the directional spectrum must be the same at each frequency
+as the total energy in the one-dimensional spectrum. At each frequency,:math:`D(\theta)` is often parameterized by an
+analytical distribution about a mean value, called a spreading function. The default spreading function
+used by WEC-Sim is a Gaussian defined by a mean :math:`\theta` and a standard deviation `\psi`.
+
+.. math::
+	D(\overline{\theta},\psi) = \frac{1}{\psi \sqrt(2 \pi)}\exp(-\overline{\theta}^2/(2\psi^2)) 
+
+When treated computationally, :math:`D(\overline{\theta},\psi)` must be discretized into a finite number
+of bins spanning a non-infinite range: in order to ensure that the total energy is preserved, the resulting
+discretization of :math:`D(\overline{\theta},\psi)` must be normalized by the fraction of energy contained in the included bins. The wave elevation for frequency-varying direction and spread can then be expressed as
+
+.. math::
+	\eta(x,y,t)\ = \sum_i \sum_j\ \sqrt{2S(\omega_i)(P(\omega_i,\theta_j))(\Delta f)}(\cos(\omega_i\ t\ -\ k_i(x \cos(\theta_j) + y \sin(\theta_j))+ \phi_{ij})
+
+where :math:`(P(\omega_n,\theta_m))` is the appropriate normalization of :math:`D(\omega,\theta)` into :math:`j` directional bins. Note that the phase :math:`\phi_{ij}` must be specified in both frequency and direction. 
+
 Irregular Wave Power
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

docs/user/code_structure.rst

@@ -207,16 +207,17 @@ depending on which wave type is selected, as shown in the table below. A more
 detailed description of the available wave types is given in the following 
 sections. 
 
-=================== ==================================================================
-**Wave Type**       **Required Properties**                         
-``noWave``          ``waves.period``                                     
+========================== ==================================================================
+**Wave Type**              **Required Properties**                         
+``noWave``                 ``waves.period``                                     
 ``noWaveCIC``                                                      
-``regular``         ``waves.height``, ``waves.period``                        
-``regularCIC``      ``waves.height``, ``waves.period``                        
-``irregular``       ``waves.height``, ``waves.period``, ``waves.spectrumType``
-``spectrumImport``  ``waves.spectrumFile``                      
-``elevationImport`` ``waves.elevationFile``                           
-=================== ==================================================================
+``regular``                ``waves.height``, ``waves.period``                        
+``regularCIC``             ``waves.height``, ``waves.period``                        
+``irregular``              ``waves.height``, ``waves.period``, ``waves.spectrumType``
+``spectrumImport``         ``waves.spectrumFile`` 
+``spectrumImportFullDir``  ``waves.spectrumFile``, ``waves.freqDepDirection.nBins``                 
+``elevationImport``        ``waves.elevationFile``                           
+========================== ==================================================================
 
 Available wave class properties, default values, and functions can be found by 
 typing ``doc waveClass`` in the MATLAB command window, or by opening the 
@@ -329,11 +330,28 @@ in the input file::
     waves.spectrumFile ='<spectrumFile >.mat';
 
 .. Note::
-    When using the ``spectrumImport`` option, users must specify a sufficient 
+    When using the ``spectrumImport`` or ``spectrumImportFullDir`` option, users must specify a sufficient 
     number of wave frequencies (typically ~1000) to adequately describe the 
     wave spectra. These wave frequencies are the same that will be used to 
     define the wave forces on the WEC, for more information refer to the 
     :ref:`user-advanced-features-irregular-wave-binning` section.
+
+spectrumImportFullDir
+"""""""""""""""""""""
+
+The ``spectrumImportFullDir`` case is for irregular wave simulations where the imported wave spectrum
+has frequency-dependent directions and/or spread values (ex: from buoy data). The user-defined spectrum
+must be defined with the wave frequency (Hz) in the first column, the spectral energy density (m^2/Hz) in the second column,
+the mean direction (degrees) in the third column, and the spread (degress) in the fourth column. 
+
+If specified, wave phase must be a rectangular matrix of size [i,j], where :math:`i` is the number of wave frequencies and :math:`j`
+is the number of directional bins. To generate a random spectra, specifying a single phase seed value (e.g., waves.phaseSeed = 128) is sufficient
+to ensure that the generated wave spectra is repeatable. 
+
+.. Note::
+    The default spread function is a Gaussian discretized into ``waves.freqDepDirection.nBins`` defined at each frequency for a mean direction (third column) and standard deviation (fourth column) over a range defined by ``waves.freqDepDirection.spreadRange`` (default = 2, so that the Gaussian will be defined for +/-2 standard deviations). It is recommended that this range is at least 2, though it will be normalized in any case so that energy is not lost due to discretization. 
+
+	At this time, this wave spectra does **NOT** work with nonlinear hydrodynamics.
 
 elevationImport
 """""""""""""""

source/functions/initializeWecSim.m

@@ -448,9 +448,10 @@
     eval(['sv_regularWavesYaw_b' num2str(ii) '= Simulink.Variant(''typeNum>=10 && typeNum<20 && yaw_' num2str(ii) '==1'');'])
     eval(['sv_irregularWaves_b' num2str(ii) '= Simulink.Variant(''typeNum>=20 && typeNum<30 && yaw_' num2str(ii) '==0'');'])
     eval(['sv_irregularWavesYaw_b' num2str(ii) '= Simulink.Variant(''typeNum>=20 && typeNum<30 && yaw_' num2str(ii) '==1'');'])
+    eval(['sv_fullDirIrregularWaves_b' num2str(ii) '= Simulink.Variant(''typeNum>=35 && typeNum<40'');'])
 end; clear ii
 
-sv_udfWaves=Simulink.Variant('typeNum>=30');
+sv_udfWaves=Simulink.Variant('typeNum>=40');
 
 % Body2Body
 B2B = simu.b2b;
@@ -497,7 +498,7 @@
 end
 
 % Visualization Blocks
-if ~isempty(waves(1).marker.location) && typeNum < 30
+if ~isempty(waves(1).marker.location) && typeNum < 40
     visON = 1;
 else
     visON = 0;

source/functions/simulink/model/calcDispPhase.m

@@ -2,21 +2,18 @@
 
 % INPUTS:
 % disp: body displacement vector, x(1) and y(2) will be used
-% waveObj: the waveClass object.
-% dispLast: the previous displacement for which a phase correction was
-%   calculated.
-% phaseLast: the previous phase correction.
 % enable: boolean, simu.largeXYDisp. To calculate transformation vector.
-% dispThresh: simu.displacementThresh. If enabled, recalculation will occur
-%   if displacement exceeds this amount from the previous threshold. 
+% direction: the direction or directional bins of the wave
+% frequency: the frequencies for which the waves have been calculated
+% wavenumber: the wavenumber for which the waves have been calculated
 
 % OUTPUTS:
 % dispPhase: a transformation matrix for real(F_exc) and
 %     imag(F_exc) based on displacement. This is identity matrix if
 %     enable = 0.
 
 %% Initialize 
-%dispPhase = zeros(length(frequency),length(direction));
+dispPhase = zeros(length(frequency),length(direction));
 %dispNew = zeros(2,1);
 
 if enable == 1  

source/functions/simulink/model/irregExcFullDirF.m

@@ -0,0 +1,108 @@
+function [Fext, relYawLast, coeffsLastMD,coeffsLastRE,coeffsLastIM] = irregExcFullDirF(pYaw, nBins, A,w, nW, dofGRD,dirGRD,wGRD,qDofGRD,qWGRD,fEHRE,fEHIM, fEHMD, phaseRand,dw,time,spreadBins, dirBins, Disp, intThresh, prevYaw, prevCoeffMD, prevCoeffRE, prevCoeffIM)
+
+
+[M,N]=size(dirBins);
+A1=bsxfun(@plus,w*time,pi/2);
+A1=repmat(A1,1,nBins);
+A = repmat(A,1,nBins);
+%initialize outputs
+Fext = zeros(1,6);
+relYawLast=zeros(M,N);
+coeffsLastMD=zeros(N,M,6); %  dirBins, freq, dof
+coeffsLastRE=zeros(N,M,6);
+coeffsLastIM=zeros(N,M,6);
+
+relYaw = dirBins-(Disp(6)*180/pi); % relative yaw angle, size = dirBins = [length(w) nBins]
+
+
+% compare relYaw to available direction data Direction data must
+% span 360 deg, inclusive (i.e [-180 180])
+if max(relYaw,[],'all')>max(dirGRD,[],'all') % handle interpolation out of BEM range by wrapping other bound
+    relYaw=wrapTo180(relYaw);
+    [relYaw,I]=sort(relYaw,2); % grid must be in ascending order
+
+elseif min(relYaw,[],'all')<min(dirGRD,[],'all') 
+    relYaw=wrapTo180(relYaw);
+    [relYaw,I]=sort(relYaw,2); % grid must be in ascending order
+else
+    I=1:length(relYaw(1,:));
+end
+
+ relYawGRD = zeros([6 size(relYaw.')]);
+ for k=1:6
+     relYawGRD(k,:,:) = relYaw.';
+ end
+
+yawError = zeros(M,N);
+if pYaw ==1 % fewer cases here because you WILL have to interpolate for any new yaw position, no chance it aligns with precalc'd grid
+    yawDiff = abs(relYaw - prevYaw);
+    if max(yawDiff,[],'all')>intThresh% interpolate for nonlinear yaw
+
+        % performs 1D interpolation in wave direction
+        fExtMDint=interpn(dofGRD,dirGRD,wGRD,fEHMD,qDofGRD,relYawGRD,qWGRD);
+        fExtREint=interpn(dofGRD,dirGRD,wGRD,fEHRE,qDofGRD,relYawGRD,qWGRD);
+        fExtIMint=interpn(dofGRD,dirGRD,wGRD,fEHIM,qDofGRD,relYawGRD,qWGRD);
+
+        % permute for dimensional consistency with linear yaw
+        fExtMDint=permute(fExtMDint,[2 3 1]);
+        fExtREint=permute(fExtREint,[2 3 1]);
+        fExtIMint=permute(fExtIMint,[2 3 1]);
+
+        relYawLast=relYaw;
+        coeffsLastMD=fExtMDint;
+        coeffsLastRE=fExtREint;
+        coeffsLastIM=fExtIMint;
+
+    else  % significant yaw maybe present, but close to previous value
+        fExtMDint=prevCoeffMD;
+        fExtREint=prevCoeffRE;
+        fExtIMint=prevCoeffIM;
+
+        relYawLast=prevYaw;
+        coeffsLastMD=prevCoeffMD;
+        coeffsLastRE=prevCoeffRE;
+        coeffsLastIM=prevCoeffIM;
+
+    end
+else
+    if sum(sum(prevYaw)) ==0; %abs(relYaw-prevYaw)>intThresh % this should only execute on the first run through
+        fExtMDint=interpn(dofGRD,dirGRD,wGRD,fEHMD,qDofGRD,relYawGRD,qWGRD);
+        fExtREint=interpn(dofGRD,dirGRD,wGRD,fEHRE,qDofGRD,relYawGRD,qWGRD);
+        fExtIMint=interpn(dofGRD,dirGRD,wGRD,fEHIM,qDofGRD,relYawGRD,qWGRD);
+        fExtMDint=permute(fExtMDint,[2 3 1]);
+        fExtREint=permute(fExtREint,[2 3 1]);
+        fExtIMint=permute(fExtIMint,[2 3 1]);
+
+        coeffsLastMD=fExtMDint;
+        coeffsLastRE=fExtREint;
+        coeffsLastIM=fExtIMint;
+        relYawLast=relYaw;
+
+    else  % pass-through for all other cases
+        fExtMDint = prevCoeffMD;
+        fExtREint = prevCoeffRE;
+        fExtIMint = prevCoeffIM;
+        coeffsLastMD=fExtMDint;
+        coeffsLastRE=fExtREint;
+        coeffsLastIM=fExtIMint;
+        relYawLast=relYaw;
+    end
+
+
+end
+
+
+B1= sin(bsxfun(@plus,A1,phaseRand(:,:,1)));
+B11 = sin(bsxfun(@plus,w*time,phaseRand(:,:,1)));
+C0 = bsxfun(@times,A.*spreadBins,dw);
+C1 = sqrt(bsxfun(@times,A.*spreadBins.^2,dw));
+for k=1:6
+    D0 =bsxfun(@times,fExtMDint(:,:,k).',C0);
+    D1 =bsxfun(@times,fExtREint(:,:,k).',C1);
+    D11 = bsxfun(@times,fExtIMint(:,:,k).',C1);
+    E1 = D0+ bsxfun(@times,B1,D1);
+    E11 = bsxfun(@times,B11,D11);
+    Fext(k) = Fext(k) + sum(sum(bsxfun(@minus,E1,E11)));
+end
+
+end