diff --git a/examples/shapeStudy/Performance.m b/examples/shapeStudy/Performance.m
new file mode 100644
index 0000000..fb71a18
--- /dev/null
+++ b/examples/shapeStudy/Performance.m
@@ -0,0 +1,297 @@
+classdef Performance < handle
+    % Data type for storage of the controlled WEC performance.
+    %
+    % This data type defines a set of parameters that are common to
+    % analysis of controlled WEC performance.
+    %
+    % The following parameters must be provided within the input struct
+    % (any additional parameters given will also be stored within the
+    % created object:
+    %
+    %     * powPerFreq
+    %
+    % The following additional parameters are available following
+    % instantiation of the object:
+    %
+    %     * pow
+    %
+    % Once created, the parameters in the object are read only, but the
+    % object can be converted to a struct, and then modified.
+    %
+    % Arguments:
+    %    input (struct):
+    %        A struct (not array) whose fields represent the parameters
+    %        to be stored.
+    %
+    % Attributes:
+    %     powPerFreq (array of float):
+    %         The exported power of the WEC per wave frequency
+    %     pow (float):
+    %         The total exported power of the WEC.
+    %
+    % Methods:
+    %    struct(): convert to struct
+    %
+    % Note:
+    %    To create an array of Performance objects see the
+    %    :mat:func:`+WecOptTool.types` function.
+    %
+    % --
+    %
+    %  Performance Properties:
+    %     powPerFreq - The exported power of the WEC per wave frequency
+    %     pow - The total exported power of the WEC
+    %
+    %  Performance Methods:
+    %    struct - convert to struct
+    %
+    % See also WecOptTool.types
+    %
+    % --
+    
+    % Copyright 2020 National Technology & Engineering Solutions of Sandia,
+    % LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the
+    % U.S. Government retains certain rights in this software.
+    %
+    % This file is part of WecOptTool.
+    %
+    %     WecOptTool is free software: you can redistribute it and/or
+    %     modify it under the terms of the GNU General Public License as
+    %     published by the Free Software Foundation, either version 3 of
+    %     the License, or (at your option) any later version.
+    %
+    %     WecOptTool is distributed in the hope that it will be useful,
+    %     but WITHOUT ANY WARRANTY; without even the implied warranty of
+    %     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    %     GNU General Public License for more details.
+    %
+    %     You should have received a copy of the GNU General Public
+    %     License along with WecOptTool.  If not, see
+    %     <https://www.gnu.org/licenses/>.
+    
+    properties
+        w (:,:) double {mustBeFinite,mustBeReal,mustBePositive}
+        ph (:,:) double {mustBeFinite,mustBeReal}
+        eta (:,:) double {mustBeFinite}
+        F0 (:,:) double {mustBeFinite}
+        u (:,:) double {mustBeFinite}
+        pos (:,:) double {mustBeFinite}
+        Zpto (:,:) double {}
+        Fpto (:,:) double {mustBeFinite}
+        pow (:,:) double {mustBeFinite}
+        name (1,:) char = 'tmp'
+        date (1,1) double {mustBeFinite,mustBePositive} = now
+        mass (:,:) double {mustBeFinite}
+    end
+    
+    methods
+        
+        function plotTime(obj ,t)
+            
+            if nargin < 2
+                trep = obj(1).getRepeatPer();
+                t = 0:0.05:trep;
+            end
+            
+            fig = figure('Name','Performance.plotTime');
+            fig.Position = fig.Position.*[1 1 1 1.5];
+            
+            % fields for plotting
+            fns = {'eta','F0','pos','u','Fpto','pow'};
+            
+            for ii = 1:length(fns)
+                ax(ii) = subplot(length(fns), 1, ii);
+                hold on
+                grid on
+            end
+            
+            for jj = 1:length(obj)
+                
+                for ii = 1:length(fns)
+                    timeRes.(fns{ii}) = getTimeRes(obj(jj),fns{ii}, t);
+                    plot(ax(ii),t,timeRes.(fns{ii}))
+                    ylabel(ax(ii),fns{ii})
+                end
+                
+                for ii = 1:length(ax) - 1
+                    set(ax(ii),'XTickLabel',[])
+                end
+                linkaxes(ax,'x')
+                xlabel(ax(end),'Time [s]')
+            end
+            xlim([t(1), t(end)])
+            
+            if length(obj) > 1
+                legend(ax(1),{obj.name})
+            end 
+            
+        end
+        
+        function plotFreq(obj,fig)
+            
+            if nargin < 2 || isempty(fig)
+                fig = figure;
+            end
+            set(fig,'Name','Performance.plotFreq');
+            
+            fns = {'F0','u','Fpto'};
+            mrks = {'o','.','+','s'};
+            
+            n = length(obj);
+            for jj = 1:n
+                for ii = 1:length(fns)
+                    
+                    fv = obj(jj).(fns{ii})(:,1); % use the first column if this is PS
+                    
+                    % mag plot
+                    ax(jj,1) = subplot(2,n,sub2ind([n,2],jj,1));
+                    title(obj(jj).name,'interpreter','none')
+                    hold on
+                    grid on
+                    
+                    stem(ax(jj,1),obj(jj).w, mag2db(abs(fv))...
+                        ,mrks{ii},...
+                        'DisplayName',fns{ii},...
+                        'MarkerSize',8,...
+                        'Color','b')
+                    
+                    % phase plot
+                    ax(jj,2) = subplot(2,n,sub2ind([n,2],jj,2));
+                    hold on
+                    grid on
+                    
+                    stem(ax(jj,2),obj(jj).w, angle(fv)...
+                        ,mrks{ii},...
+                        'DisplayName',fns{ii},...
+                        'MarkerSize',8,...
+                        'Color','b')
+                    
+                    ylim(ax(jj,2),[-pi,pi])
+                end
+                xlabel(ax(jj,2),'Frequency [rad/s]')
+            end
+            ylabel(ax(1,1),'Magnitude [dB]')
+            ylabel(ax(1,2),'Angle [rad]')
+            legend(ax(n,1))
+            linkaxes(ax,'x')
+            linkaxes(ax(:,1),'y')
+            
+        end
+        
+        function T = summary(obj)
+            
+            if length(obj) > 1
+                for ii = 1:length(obj)
+                    Tr(ii,:) = summary(obj(ii));
+                end
+                
+                % augment names if they are the same
+                if any(strcmp(obj(1).name, {obj(2:end).name}))
+                    for ii = 1:length(obj)
+                        rnames{ii} = [obj(ii).name, '_', num2str(ii)];
+                    end
+                else
+                    rnames = {obj.name};
+                end
+                
+                Tr.Properties.RowNames = rnames;
+                mT = Tr;
+                
+                if nargout
+                    T = mT;
+                else
+                    disp(mT)
+                end
+                
+                return
+                
+            else
+                rnames = {obj.name};
+            end
+            
+            trep = obj.getRepeatPer();
+            t = linspace(0,trep,1e3);
+
+            for jj = 1:size(obj.ph,2) % for each phase in PS cases
+
+                tmp.pow_avg(jj) = sum(real(obj.pow(:,jj)));
+
+                pow_t = getTimeRes(obj, 'pow', t, jj);
+                tmp.pow_max(jj) = max(abs(pow_t));
+
+                try
+                    tmp.pow_thd(jj) = thd(pow_t);
+                catch ME
+                    warning(ME.message)
+                    tmp.pow_thd(jj) = NaN;
+                end
+
+                pos_t = getTimeRes(obj, 'pos', t, jj);
+                tmp.pos_max(jj) = max(abs(pos_t));
+
+                vel_t = getTimeRes(obj, 'u', t, jj);
+                tmp.vel_max(jj) = max(abs(vel_t));
+
+                Fpto_t = getTimeRes(obj, 'Fpto', t, jj);
+                tmp.Fpto_max(jj) = max(abs(Fpto_t));
+            end
+
+            fn = fieldnames(tmp);
+            for kk = 1:length(fn)
+                out.(fn{kk}) = mean(tmp.(fn{kk}), 2);
+            end
+                
+            rnames = reshape(rnames,[],1);
+            
+            mT = table(out.pow_avg(:),out.pow_max(:),out.pow_thd(:),...
+                out.pos_max(:),out.vel_max(:),out.Fpto_max(:),...
+                'VariableNames',...
+                {'AvgPow','|MaxPow|','PowTHD_dBc','MaxPos','MaxVel','MaxPTO'},...
+                'RowNames',rnames);
+            
+            if nargout
+                T = mT;
+            else
+                disp(mT)
+            end
+            
+        end
+        
+    end
+    
+    methods (Access=protected)
+        
+        function [tRep] = getRepeatPer(obj)
+            tRep = 2*pi/(obj.w(2) - obj.w(1));
+        end
+        
+        function [timeRes] = getTimeRes(obj, fn, t_vec, ph_idx)
+            if nargin < 4
+                ph_idx = 1;
+            end
+            
+            if strcmp(fn,'pow')
+                vel = obj.getTimeRes('u',t_vec);
+                f = obj.getTimeRes('Fpto',t_vec);
+                timeRes = vel .* f;
+            else
+                timeRes = zeros(size(t_vec));
+                fv = obj.(fn)(:,ph_idx); % use the first column if this is PS
+                for ii = 1:length(obj.w) % for each freq. TODO - use IFFT
+                    timeRes = timeRes ...
+                        + real(fv(ii) * exp(1i * obj.w(ii) * t_vec));
+                end
+            end
+        end
+        
+%         function checkSizes(varargin) % TODO
+%             n = length(varargin);
+%             for ii = 1:n
+%                 if ~isequal(varargin(varargin{ii}),size(varargin{1}))
+%                     error('Frequency vectors must have same size')
+%                 end
+%             end
+%         end
+        
+    end
+end
\ No newline at end of file
diff --git a/examples/shapeStudy/ShapeStudy.m b/examples/shapeStudy/ShapeStudy.m
new file mode 100644
index 0000000..cb79dcf
--- /dev/null
+++ b/examples/shapeStudy/ShapeStudy.m
@@ -0,0 +1,475 @@
+% ShapeStudy
+%
+% This file is part of WecOptTool.
+%
+%     WecOptTool is free software: you can redistribute it and/or modify it
+%     under the terms of the GNU General Public License as published by the
+%     Free Software Foundation, either version 3 of the License, or (at
+%     your option) any later version.
+%
+%     WecOptTool is distributed in the hope that it will be useful, but
+%     WITHOUT ANY WARRANTY; without even the implied warranty of
+%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+%     General Public License for more details.
+%
+%     You should have received a copy of the GNU General Public License
+%     along with WecOptTool.  If not, see <https://www.gnu.org/licenses/>.
+constant='Spect';
+%constant='Force';
+%constant='Power';
+
+%% define sea state of interest
+dw = 0.3142;
+nf = 50;
+%nf = 30;
+w = dw * (1:nf)';
+S = jonswap(w,[0.25, 3.5, 1]);
+SS = WecOptTool.SeaState(S);
+
+% Amp = 0.125/2;
+% wp = w(6);
+% fp = wp/(2*pi);
+% Tp = 1/fp; 
+% SS = WecOptTool.SeaState.regularWave(w,[Amp,Tp]);
+
+%% set up device types
+
+controlType{1} = 'CC';
+controlType{2} = 'P';
+controlType{3} = 'PS';
+zmax = 0.6;
+fmax = 1e10;
+
+folder = WecOptTool.AutoFolder();
+
+%% Create Devices & Plot Geometries
+
+%AR = [ 0.1, 0.125, 0.15,  0.2, 0.3, 0.4, 0.5, 0.7, 1.0, 2, 3, 6, 10];
+AR = [0.3, 0.4, 0.5, 0.7, 1.0, 2, 3, 6, 10];
+%AR = [ 5, 10, 50, 100, 200];
+%AR = [0.1:0.5:10];
+volume=0.875;
+heights = (volume./(pi * AR.^2)).^(1/3);
+radii = AR.*heights;
+
+fig = figure('name','Shape Study');
+%fig.Position = fig.Position .*[1,1,1.5,0.75];
+hold on
+grid on
+set(gca, 'fontsize',16);
+for ii = 1:length(AR)
+    radius = radii(ii);
+    height = heights(ii);
+    xCoords =  [0, radius, radius, 0];
+    yCoords = [0.2, 0.2, -height, -height];
+    ax = gca;
+    p(ii) = plot(ax, xCoords, yCoords, 'bo-','DisplayName',num2str(ii));    
+end
+axis equal
+% xlabel('Radius [m]')
+% ylabel('Height [m]')
+% ylim([-inf 0])
+% x = [0.7 0.5];
+% y = [0.2 0.2];
+% annotation('textarrow',x,y,'String','Low AR')
+% x = [0.85 0.85];
+% y = [0.4 0.7];
+% annotation('textarrow',x,y,'String','High AR')
+% saveas(fig, 'plots/aspectRatios.pdf')
+
+%% Running hydrodynamics
+
+N = length(AR);
+% Determine the value of the constant force
+iAR = N; 
+remaingARIndicies = setdiff(find(AR), iAR);
+
+[deviceHydro, mesh] = designDevice('parametric', folder.path, ...
+                                   radii(iAR), heights(iAR), w);
+% Tune mass to desired natural frequency
+wdes = 0.625*2*pi;  
+dynModel = getDynamicsModel(deviceHydro, ...
+                            SS,          ...
+                            'linear',    ...
+                            wdes); 
+%mean(dynModel.F0)
+% Constant Force? 
+if all(constant == 'Force')
+    Fconstant = ones(length(w),1) * mean(dynModel.F0);
+    dynModel.F0 = Fconstant;
+elseif all(constant == 'Power')
+    Fconstant = sqrt(8*real([dynModel.Zi]));
+    dynModel.F0 = Fconstant(:,1);
+end
+if ~all(constant == 'Spect')
+    dynModel.eta_fd = dynModel.F0 ./ dynModel.Hex;
+end
+perform = simulateDevice(dynModel(1), controlType{1},     ...
+                         'interpMethod','nearest',    ...
+                         'Zmax',zmax, 'Fmax',fmax);       
+
+deviceHydro = repmat(deviceHydro, N, 1 );
+meshes = repmat(mesh, N, 1 );
+dynModel = repmat(dynModel, N, 1 );
+perform  = repmat(perform, N, 1 );
+
+for i = remaingARIndicies
+    radius = radii(i);
+    height = heights(i);
+    [deviceHydro(i), meshes(i)] = designDevice('parametric', ...
+                                      folder.path, radius, height, w);    
+
+    %% Tune Device Mass to a Natural Frequency
+    dynModel(i) = getDynamicsModel(deviceHydro(i), ...
+                                   SS,       ...
+                                   'linear', ...
+                                   wdes);      
+                           
+    if all(constant == 'Force')
+        Fconstant = ones(length(w),1) * mean(dynModel(iAR).F0);%Global Const
+        %Fconstant = ones(length(w),1) * mean(dynModel(i).F0); % AR const
+        dynModel(i).F0 = Fconstant;
+    elseif all(constant == 'Power')
+        Fconstant = sqrt(8*real([dynModel(i).Zi]));
+        dynModel(i).F0 = Fconstant(:,1);        
+    end
+    if ~all(constant == 'Spect')
+        dynModel(i).eta_fd = dynModel(i).F0 ./ dynModel(i).Hex;
+    end
+    perform(i) = simulateDevice(dynModel(i), controlType{1}, ...
+                                'interpMethod','nearest',    ...
+                                'Zmax',zmax, 'Fmax',fmax);
+end
+
+%% Plot the parametric meshes
+plotMesh=false;
+if plotMesh == true
+%WecOptTool.plot.plotMesh(mesh)
+for ii = 1:length(AR)
+    WecOptTool.plot.plotMesh(meshes(ii))
+end    
+end
+ 
+
+%% Plot mass ratio
+fig = figure('name','Mass Ratio');
+hold on
+grid on
+set(gca, 'fontsize',16);
+semilogx(AR,[dynModel.mass] ./( [deviceHydro.Vo] .* [deviceHydro.rho] ), ...
+        '--xk' )
+set(gca,'XScale','log')
+%ylim([0 1])
+xlabel('Aspect Ratio [-]')
+ylabel('Mass ratio, m^\prime/m [-]')
+saveas(fig, 'plots/massRatios.pdf')
+
+
+%% Plot hydro
+%close all
+clear ax
+
+[xp,yp] = meshgrid(w,AR);
+
+flag3D=false;
+%flag3D=true;
+
+% Figure 1: Added Mass
+nFig = 1;
+fh(nFig) = figure('name','Radiation added mass (A)');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+A = [dynModel.A]-[dynModel.Ainf];
+if flag3D == 1
+    surf(xp',yp',(A));
+    rotate3d on
+else
+    contourf(xp',yp',(A));
+    c(nFig) = colorbar;
+end
+
+
+% Figure 2: Radiation Damping
+nFig = nFig + 1;
+fh(nFig) = figure('name','Radiation wave damping (B)');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+if flag3D == 1
+    surf(xp',yp',[dynModel.B])
+    rotate3d on
+else
+    contourf(xp',yp',[dynModel.B])
+    c(nFig) = colorbar;
+end
+
+
+% Figure 3: Position
+nFig = nFig + 1;
+fh(nFig) = figure('name','|Position|');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+
+if flag3D == 1
+    surf(xp',yp', abs([perform.pos]) )
+    rotate3d on
+else
+    contourf(xp',yp', abs([perform.pos]))% ./max(max(abs([perform.pos]))))
+    c(nFig) = colorbar;
+    set(gca,'ColorScale','log')
+end
+
+
+% Figure 4: Power
+nFig = nFig + 1;
+fh(nFig) = figure('name','|\Re(Power)|');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+% surf(xp',yp',abs(Ex)./max(max(abs(Ex))))
+Popt = abs([dynModel.Hex].*SS.S).^2./(8*real([dynModel.Zi]));
+Popt(~isfinite(Popt)) = 0;
+%surf(xp',yp',Popt./max(max(Popt)))
+
+if flag3D == 1
+    surf(xp',yp',abs(real([perform.pow])) )
+    rotate3d on
+else
+    contourf(xp',yp',abs(real([perform.pow])) )
+    c(nFig) = colorbar;
+    set(gca,'ColorScale','log')
+end
+
+
+% Figure 5: Zpto
+nFig = nFig + 1;
+fh(nFig) = figure('name','|Zpto|');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+%surf(xp',yp',abs([perform.Zpto])./max(max(abs([perform.Zpto]))))
+if flag3D == 1
+    surf(xp',yp',abs([perform.Zpto]))
+    rotate3d on
+else
+    contourf(xp',yp',abs([perform.Zpto]))
+    c(nFig) = colorbar;
+    set(gca,'ColorScale','log')
+end
+
+
+% Figure 6: Real Zpto
+nFig = nFig + 1;
+fh(nFig) = figure('name','Real(Zpto)');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+if flag3D == 1
+    surf(xp',yp',real([perform.Zpto]) )
+    rotate3d on
+else
+    contourf(xp',yp',real([perform.Zpto]) )
+    c(nFig) = colorbar;
+end
+
+
+% Figure 7: Imag Zpto
+nFig = nFig + 1;
+fh(nFig) = figure('name','Im(Zpto)');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+
+if flag3D == 1
+    surf(xp',yp',imag([perform.Zpto]) )
+    surf(xp',yp',zeros(size(xp')),'FaceAlpha',0.25,'EdgeColor','none',...
+     'FaceColor','blue')
+    rotate3d on
+else
+    contourf(xp',yp',imag([perform.Zpto]) )
+    c(nFig) = colorbar;
+end
+
+% Figure 8: Velocity
+nFig = nFig + 1;
+fh(nFig) = figure('name','|Velocity|');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+
+if flag3D == 1
+    surf(xp',yp',abs([perform.u]) )
+    rotate3d on
+else
+    contourf(xp',yp',abs([perform.u]) )
+    c(nFig) = colorbar;
+    set(gca,'ColorScale','log')
+end
+
+
+% Figure 9: eta_fd
+nFig = nFig + 1;
+fh(nFig) = figure('name','|\eta_{fd}|');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid on
+eta = [dynModel.eta_fd];
+if flag3D == 1
+    surf(xp(:,:)',yp(:,:)',abs(eta(:,:) ) )% ./ max(max((abs([eta])))))
+    set(ax(nFig),'zscale','log')
+    rotate3d on
+else
+    contourf(xp(:,:)',yp(:,:)',abs(eta(:,:) ) )% ./ max(max((abs([eta])))))
+    c(nFig) = colorbar;    
+    set(gca,'ColorScale','log')
+end
+
+
+% Figure 10: Fex
+nFig = nFig + 1;
+fh(nFig) = figure('name','|F_{Ex}|');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+
+if flag3D == 1
+    surf(xp',yp',abs([perform.F0]))
+    set(ax(nFig),'zscale','log')
+    rotate3d on
+else
+    contourf(xp(:,:)',yp(:,:)',abs([perform.F0]) )
+    c(nFig) = colorbar;    
+    set(gca,'ColorScale','log')
+end
+
+
+% Figure 11: Zi
+nFig = nFig + 1;
+fh(nFig) = figure('name','|Z_{i}|');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+
+if flag3D == 1
+    surf(xp',yp',abs([dynModel.Zi]))
+    set(ax(nFig),'zscale','log')
+    rotate3d on
+else
+    contourf(xp(:,:)',yp(:,:)',abs([dynModel.Zi]) )
+    c(nFig) = colorbar;    
+    set(gca,'ColorScale','log')
+end
+
+
+% Figure 12: Hex
+nFig = nFig + 1;
+fh(nFig) = figure('name','H_{ex}');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+
+if flag3D == 1
+    surf(xp',yp',abs([dynModel.Hex]))
+    set(ax(nFig),'zscale','log')
+    rotate3d on
+else
+    contourf(xp(:,:)',yp(:,:)',abs([dynModel.Hex]) )
+    c(nFig) = colorbar;    
+    set(gca,'ColorScale','log')
+end
+
+% Figure 13: Reactive Power
+nFig = nFig + 1;
+fh(nFig) = figure('name','|\Im(Power)|');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+
+if flag3D == 1
+    surf(xp',yp',imag([perform.pow])) 
+    rotate3d on
+else
+    contourf(xp',yp',abs(imag([perform.pow])))
+    c(nFig) = colorbar;
+    set(gca,'ColorScale','log')
+end
+
+
+if flag3D == 1
+    Link = linkprop(ax, ...
+           {'CameraUpVector', 'CameraPosition', 'CameraTarget'});
+    setappdata(gcf, 'StoreTheLink', Link);
+    view([20, 12])
+end
+
+
+for ii = 1:nFig
+    xlabel(ax(ii),'Freq. [rad/s]')
+    ylabel(ax(ii),'Aspect ratio [ ]')
+    if flag3D==1
+        zlabel(ax(ii),fh(ii).Name)
+    else
+        c(ii).Label.String = fh(ii).Name;
+    end
+    %fname = strrep(strrep(strrep(strrep(fh(ii).Name,'|',''), '\',''), '{',''),'}','');
+    fname = strrep(fh(ii).Name,'|','');
+    path = "plots/%s/%s.pdf";
+    
+    str = sprintf(path,constant,fname);
+    saveas(fh(ii), str)
+end
+
+
+%% Plot omega, Zi for each AR (2D)
+nFig = nFig + 1;
+fh(nFig) = figure('name','Z_{i}');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+hold on
+grid 
+for ii = 1:length(AR)
+    plot(xp(ii,:)',abs([dynModel(ii).Zi]), 'DisplayName',string(AR(ii)))    
+end   
+xlabel('Freq. [rad/s]')
+ylabel('Z_i [\Omega]')
+legend()
+
+
+%% Plot omega, Zi for each AR (2D)
+nFig = nFig + 1;
+fh(nFig) = figure('name','Z_{i}');
+ax(nFig) = gca;
+set(gca, 'fontsize',16);
+set(gca, 'yscale','log');
+hold on
+grid 
+for ii = 1:length(AR)
+    plot(w',abs(imag([dynModel(ii).Zi])), 'DisplayName',string(AR(ii)))    
+end   
+xlabel('Freq. [rad/s]')
+ylabel('')
+legend()
+
+%% Bode 
+hold on
+
+for ii = 1:length(AR)
+    bode(frd(dynModel(ii).Zi,w)),  
+end   
+xlabel('Freq. [rad/s]')
+xlim([0.3,11])
+grid 
+legendCell = cellstr(num2str(AR'));
+legend(legendCell)
+
diff --git a/examples/shapeStudy/designDevice.m b/examples/shapeStudy/designDevice.m
new file mode 100644
index 0000000..ed80e68
--- /dev/null
+++ b/examples/shapeStudy/designDevice.m
@@ -0,0 +1,98 @@
+function [hydro, meshes] = designDevice(type, varargin) 
+    % WaveBot   WEC based on the Sandia "WaveBot" device.
+    %
+    % The WaveBot is a model-scale wave energy converter (WEC) tested in
+    % the Navy's Manuevering and Sea Keeping (MASK) basin. Reports and
+    % papers about the WaveBot are available at advweccntrls.sandia.gov.
+    
+    switch type
+        
+        case 'existing'
+            hydro = WecOptTool.geometry.existingNEMOH(varargin{:});
+        case 'scalar'
+            hydro = getHydroScalar(varargin{:});
+        case 'parametric'
+            [hydro, meshes] = getHydroParametric(varargin{:});
+        otherwise
+            error('WecOptTool:UnknownGeometryType',...
+                'Invalid geometry type')
+    end
+    
+end
+
+
+function [hydro, meshes] = getHydroParametric(folder, r1, d1, w)
+               
+    if w(1) == 0
+        w = w(2:end);
+    end
+    
+    z0 = 0;
+    r = [0,    r1,          r1,  r1,   0];
+    z = [z0, z0, (-d1-z0)/2, -d1, -d1];
+
+    % Mesh
+    ntheta = 20;
+    nfobj = 200;
+    zG = 0;
+    
+    meshes = WecOptTool.mesh("AxiMesh",    ...
+                             folder,       ...
+                             r,            ...
+                             z,            ...
+                             ntheta,       ...
+                             nfobj,        ...
+                             zG,           ...
+                             1);
+
+    hydro = WecOptTool.solver("NEMOH", folder, meshes, w);
+           
+end
+
+% function hydro = getHydroScalar(folder, lambda, w)
+%                    
+%     if w(1) == 0
+%     error('WecOptTool:UnknownGeometryType',...
+%                 'Invalid frequency vector')     % TODO - more checks
+%     end
+%     
+%     r = lambda * [0, 0.88, 0.88, 0.35, 0];
+%     z = lambda * [0.2, 0.2, -0.16, -0.53, -0.53];
+% 
+%     % Mesh
+%     ntheta = 20;
+%     nfobj = 200;
+%     zG = 0;
+%     
+%     meshes = WecOptTool.mesh("AxiMesh",    ...
+%                              folder,       ...
+%                              r,            ...
+%                              z,            ...
+%                              ntheta,       ...
+%                              nfobj,        ...
+%                              zG,           ...
+%                              1);
+%     
+%     hydro = WecOptTool.solver("NEMOH", folder, meshes, w);
+%            
+% end
+
+
+% Copyright 2020 National Technology & Engineering Solutions of Sandia, 
+% LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the 
+% U.S. Government retains certain rights in this software.
+%
+% This file is part of WecOptTool.
+% 
+%     WecOptTool is free software: you can redistribute it and/or modify
+%     it under the terms of the GNU General Public License as published by
+%     the Free Software Foundation, either version 3 of the License, or
+%     (at your option) any later version.
+% 
+%     WecOptTool is distributed in the hope that it will be useful,
+%     but WITHOUT ANY WARRANTY; without even the implied warranty of
+%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%     GNU General Public License for more details.
+% 
+%     You should have received a copy of the GNU General Public License
+%     along with WecOptTool.  If not, see <https://www.gnu.org/licenses/>.
diff --git a/examples/shapeStudy/getDynamicsModel.m b/examples/shapeStudy/getDynamicsModel.m
new file mode 100644
index 0000000..e095640
--- /dev/null
+++ b/examples/shapeStudy/getDynamicsModel.m
@@ -0,0 +1,89 @@
+function dynModel = getDynamicsModel(hydro, SS, interpMethod, wdes)
+    
+    % Restoring
+    K = hydro.C(3,3) * hydro.g * hydro.rho;
+
+    function result = interp_mass(hydro, dof1, dof2, w)
+        result = interp1(hydro.w,                           ...
+                         squeeze(hydro.A(dof1, dof2, :)),   ...
+                         w,                                 ...
+                         interpMethod,                          ...
+                         0);
+    end
+
+    function result = interp_rad(hydro, dof1, dof2, w)
+        result = interp1(hydro.w,                           ...
+                         squeeze(hydro.B(dof1, dof2, :)),   ...
+                         w,                                 ...
+                         interpMethod,                          ...
+                         0);
+    end
+
+    function result = interp_ex(hydro, dof, w)
+
+        h = squeeze(hydro.ex(dof, 1, :));
+        result = interp1(hydro.w, h ,w, interpMethod, 0);
+
+    end
+    
+    w = hydro.w(:);
+    dw = w(2) - w(1);
+    
+    % Calculate wave amplitude
+    waveAmpSS = SS.getAmplitudeSpectrum;
+    waveAmp = interp1(SS.w, waveAmpSS, w, interpMethod, 'extrap');
+
+    % Row vector of random phases
+    ph = rand(size(waveAmp))*2*pi;
+
+    % Wave height in frequency domain
+    eta_fd = waveAmp .* exp(1i * ph);
+    eta_fd = eta_fd(:);
+
+    % radiation damping FRF
+    B = interp_rad(hydro, 3, 3, w) * hydro.rho .* w;
+%   B(:,ii) = squeeze(hydro.B(3,3,:)).*w'*rho;
+
+    % added mass FRF
+    A = interp_mass(hydro, 3, 3, w) * hydro.rho;
+    Ainf = hydro.Ainf(3,3)*hydro.rho;
+
+    % friction
+    Bf = max(B) * 0.1;      % TODO - make this adjustable 
+
+    % Tune device mass to desired natural frequency
+    m = hydro.Vo * hydro.rho;
+    fun = @(m) tune_wdes(wdes,m,K,w',A);
+    mass = fminsearch(fun,m);
+    
+    % intrinsic impedance
+    Zi = B + Bf + 1i * (w .* (mass + A) - K ./ w);
+
+    % Excitation Forces
+    Hex = interp_ex(hydro, 3, w) * hydro.g * hydro.rho;
+    F0 = Hex .* eta_fd;   
+    
+    dynModel.mass = mass;
+    dynModel.K = K;
+    dynModel.w = w;
+    dynModel.eta_fd = eta_fd;
+    dynModel.dw = dw;
+    dynModel.wave_amp = waveAmp;
+    dynModel.ph = ph;
+    dynModel.B = B;
+    dynModel.A = A;
+    dynModel.Ainf = Ainf;
+    dynModel.Bf = Bf;
+    dynModel.Zi = Zi;
+    dynModel.Hex = Hex;
+    dynModel.F0 = F0;
+    
+    
+end
+
+function [err] = tune_wdes(wdes,m,k,w,A)
+    fun1 = @(w1) w1.^2.*(m + interp1(w,A,w1)) - k;
+    options = optimset('Display','off');
+    w0 = fsolve(fun1,5, options);
+    err = (w0 - wdes).^2;
+end
\ No newline at end of file
diff --git a/examples/shapeStudy/simulateDevice.m b/examples/shapeStudy/simulateDevice.m
new file mode 100644
index 0000000..2361baa
--- /dev/null
+++ b/examples/shapeStudy/simulateDevice.m
@@ -0,0 +1,400 @@
+function performance = simulateDevice(dynModel, controlType,...
+                                      options)
+    % simulateDevice   WEC based on the Sandia "WaveBot" device.
+    %
+    % The WaveBot is a model-scale wave energy converter (WEC) tested in
+    % the Navy's Manuevering and Sea Keeping (MASK) basin. Reports and
+    % papers about the WaveBot are available at advweccntrls.sandia.gov.
+    %
+    % Arguments:
+    %  dynModel      
+    %  controlType  controller type:
+    %                   complex conjugate:      'CC'
+    %                   proportional damping:   'P'
+    %                   pseudo-spectral:        'PS'
+    %  name-value pairs
+    %  interpMethod (optional) method to use for linear interpolation
+    %  Zmax         (only valid for controlType == 'controlType') maximum
+    %               displacement
+    %  Fmax         (only valid for controlType == 'controlType') maximum
+    %               PTO force
+    %
+    % See also WecOptTool.SeaState, interp1
+    %hydro (1,1) WecOptTool.Hydrodynamics
+    arguments
+        dynModel 
+        controlType (1,1) string
+        options.Zmax (1,:) double  = Inf % TODO - can be assymetric, need to check throughout
+        options.Fmax (1,:) double = Inf
+        options.interpMethod (1,1) string = 'linear'
+    end
+    
+   
+    switch controlType
+        case 'CC'
+            performance = complexCongugateControl(dynModel);
+        case 'P'
+            performance = dampingControl(dynModel);
+        case 'PS'
+            performance = psControl(dynModel,options.Zmax, options.Fmax);            
+    end
+end
+        
+
+
+function myPerf = complexCongugateControl(dynModel,~)
+    
+    myPerf = Performance();
+            
+    myPerf.Zpto = conj(dynModel.Zi);
+    
+    % velocity
+    myPerf.u = dynModel.F0 ./ (myPerf.Zpto + dynModel.Zi);
+    
+    % position
+    myPerf.pos = myPerf.u ./ (1i * dynModel.w);
+    
+    % PTO force
+    myPerf.Fpto = -1 * myPerf.Zpto .* myPerf.u;
+    
+    % power
+    myPerf.pow = 0.5 * myPerf.Fpto .* conj(myPerf.u);
+    
+    myPerf.ph = dynModel.ph;
+    myPerf.w = dynModel.w;
+    myPerf.eta = dynModel.eta_fd;
+    myPerf.F0 = dynModel.F0;
+    myPerf.mass = dynModel.mass;
+
+end
+
+function myPerf = dampingControl(dynModel,~)
+    
+    myPerf = Performance();
+            
+    P_max = @(b) -0.5*b*sum(abs(dynModel.F0 ./ ...
+                                (dynModel.Zi + b)).^2);
+                            
+    % Solve for damping to produce most power (can do analytically for a
+    % single frequency, but must use numerical solution for spectrum). Note
+    % that fval is the sum of power absorbed (negative being "good") - the
+    % following should be true: -1 * fval = sum(pow), where pow is the
+    % frequency dependent array calculated below.
+    [B_opt, ~] = fminsearch(P_max, max(real(dynModel.Zi)));
+
+    % PTO impedance
+    myPerf.Zpto = complex(B_opt * ones(size(dynModel.Zi)),0);
+    
+    % velocity
+    myPerf.u = dynModel.F0 ./ (myPerf.Zpto + dynModel.Zi);
+    
+    % position
+    myPerf.pos = myPerf.u ./ (1i * dynModel.w);
+    
+    % PTO force
+    myPerf.Fpto = -1 * myPerf.Zpto .* myPerf.u;
+    
+    % power
+    myPerf.pow = 0.5 * myPerf.Fpto .* conj(myPerf.u);
+    
+    
+    myPerf.ph = dynModel.ph;
+    myPerf.w = dynModel.w;
+    myPerf.eta = dynModel.eta_fd;
+    myPerf.F0 = dynModel.F0;
+    myPerf.mass = dynModel.mass;
+
+end
+
+function myPerf = psControl(dynModel,delta_Zmax,delta_Fmax)
+%     motion = getPSCoefficients(motion, delta_Zmax, delta_Fmax);
+%     ps.wave_amp = waveAmp; % TODO
+%     
+%     % Use mutliple phase realizations for PS at the model
+%     % is nonlinear (note that we use the original phasing
+%     % from the other cases)
+%     n_ph = 5;
+%     ph_mat = [ph, rand(length(ps.w), n_ph-1)];
+%     
+%     n_freqs = length(motion.w);
+%     phasePowMat = zeros(n_ph, 1);
+%     powPerFreqMat = zeros(n_freqs, n_ph);
+%     
+%     for ind_ph = 1 : n_ph
+%         
+%         ph = ph_mat(:, ind_ph);
+% %         [powTot, fRes(ind_ph), tRes(ind_ph)] = getPSPhasePower(ps, ph);
+%         [pow, powPerFreq] = getPSPhasePower(motion, ph)
+%         phasePowMat(ind_ph) = powTot;
+%         powPerFreqMat(:, ind_ph) = fRes(ind_ph).pow;
+%         
+%     end
+%     
+%     ph = ph_mat(:,1);
+%     u = fRes(1).vel;
+%     pos = fRes(1).pos;
+%     Zpto = nan(size(motion.hydro.Zi)); % TODO
+%     Fpto = fRes(1).u;
+%     pow = powPerFreqMat(:,1);
+
+    arguments
+        dynModel
+        delta_Zmax (1,:) double {mustBeFinite,mustBeReal,mustBePositive}
+        delta_Fmax (1,:) double {mustBeFinite,mustBeReal,mustBePositive}
+    end
+        
+    % Fix random seed <- Do we want this???
+    rng(1);
+    
+    % Reformulate equations of motion
+    dynModel = getPSCoefficients(dynModel, delta_Zmax, delta_Fmax);
+    
+    % Add phase realizations
+    n_ph = 5;
+    ph_mat = [dynModel.ph, rand(length(dynModel.w), n_ph-1)];
+
+    for ind_ph = 1 : n_ph
+        
+        ph = ph_mat(:, ind_ph);
+        [phasePowMat(ind_ph), fRes(ind_ph), tRes(ind_ph)] = ...
+            getPSPhasePower(dynModel, ph);
+        
+        
+        pos(:, ind_ph) = fRes(ind_ph).pos;
+        u(:, ind_ph) = fRes(ind_ph).vel;
+        Zpto(:, ind_ph) = fRes(ind_ph).Zpto;
+        Fpto(:, ind_ph) = fRes(ind_ph).u;
+        pow(:, ind_ph) = fRes(ind_ph).pow;
+        eta(:, ind_ph) = fRes(ind_ph).eta;
+        F0(:, ind_ph) = fRes(ind_ph).F0;
+        
+    end
+    
+    % assemble results
+    myPerf = Performance();
+    myPerf.w = dynModel.w;
+    myPerf.eta = eta;
+    myPerf.F0 = F0;
+    myPerf.ph = ph_mat;
+    myPerf.u = u;
+    myPerf.pos = pos;
+    myPerf.Zpto = Zpto;
+    myPerf.Fpto = Fpto;
+    myPerf.pow = pow;
+    myPerf.mass = dynModel.mass;
+    
+end
+
+function dynModel = getPSCoefficients(dynModel, delta_Zmax, delta_Fmax)
+    % getPSCoefficients   constructs the necessary coefficients and
+    % matrices used in the pseudospectral control optimization
+    % problem
+    %
+    % Note that these coefficients are not sea state dependent,
+    % thus it is beneficial to find them once only when doing a
+    % study involving multiple sea states.
+    %
+    % Bacelli 2014: Background Chapter 4.1, 4.2; RM3 in section 6.1
+    
+    % Number of frequency - half the number of Fourier coefficients
+    Nf = length(dynModel.w);
+    
+    % Collocation points uniformly distributed between 0 and T
+    % note that we have 2*Nf collocation points since we will have
+    % two Fourier coefficients for each frequency
+    Nc = (2*Nf) + 2;
+    
+    % Rebuild frequency vector to ensure monotonically increasing
+    % with w(1) = w0
+    w0 = dynModel.dw;                    % fundamental frequency
+    T = 2 * pi/w0;                  % '' period
+    
+    % Building cost function component
+    % we will form the cost function as transpose(x) * H x, where x
+    % is a vector of [vel, u]; we want the product above to result
+    % in power (u*vel)
+    H = [0,1;1,0];
+    H_mat = 0.5 * kron(H, eye(2*Nf));
+    
+    % Building matrices B33 and A33
+    Adiag33 = zeros(2*Nf-1,1);
+    Bdiag33 = zeros(2*Nf,1);
+    
+    Adiag33(1:2:end) = dynModel.w.* dynModel.A;
+    Bdiag33(1:2:end) = dynModel.B;
+    Bdiag33(2:2:end) = Bdiag33(1:2:end);
+    
+    Bmat = diag(Bdiag33);
+    Amat = diag(Adiag33,1);
+    Amat = Amat - Amat';
+    
+    G = Amat + Bmat;
+    
+    B = dynModel.Bf * eye(2*Nf);
+    C = blkdiag(dynModel.K * eye(2*Nf));
+    M = blkdiag(dynModel.mass * eye(2*Nf));
+    
+    % Building derivative matrix
+    d = [dynModel.w(:)'; zeros(1, length(dynModel.w))];
+    Dphi1 = diag(d(1:end-1), 1);
+    Dphi1 = (Dphi1 - Dphi1');
+    Dphi = blkdiag(Dphi1);
+    
+    % scaling factor to improve optimization performance
+    m_scale = dynModel.mass;
+    
+    % equality constraints for EOM
+    P =  (M*Dphi + B + G + (C / Dphi)) / m_scale;
+    Aeq = [P, -eye(2*Nf) ];
+    Aeq = [Aeq,            zeros(2*Nf,2);
+        zeros(1,4*Nf), dynModel.K / m_scale, -1];
+    
+    % Calculating collocation points for constraints
+    tkp = linspace(0, T, 4*(Nc));
+    tkp = tkp(1:end);
+    Wtkp = dynModel.w*tkp;
+    Phip1 = zeros(2*size(Wtkp,1),size(Wtkp,2));
+    Phip1(1:2:end,:) = cos(Wtkp);
+    Phip1(2:2:end,:) = sin(Wtkp);
+    
+    Phip = blkdiag(Phip1);
+    
+    A_ineq =  [kron([1 0], Phip1' / Dphi1), ones(4*Nc,1), zeros(4*Nc,1)];
+    A_ineq = [A_ineq; -A_ineq];
+    
+    % position constraints
+    if length(delta_Zmax)==1
+        B_ineq = [ones(size(A_ineq, 1),1) * delta_Zmax];
+    else
+        B_ineq = [ones(size(A_ineq, 1)/2,1) * max(delta_Zmax);
+            -ones(size(A_ineq, 1)/2,1) * min(delta_Zmax)];
+    end
+    
+    % force constraints
+    siz = size(A_ineq);
+    forc =  [kron([0 1], Phip'), zeros(4*Nc,1), ones(4*Nc,1)];
+    if length(delta_Fmax)==1
+        B_ineq = [B_ineq; ones(siz(1),1) * delta_Fmax/m_scale];
+    else
+        B_ineq = [B_ineq; ones(siz(1)/2,1) * max(delta_Fmax)/m_scale;
+            -ones(siz(1)/2,1) * min(delta_Fmax)/m_scale];
+    end
+    A_ineq = [A_ineq; forc; -forc];
+    
+    dynModel.Nf = Nf;
+    dynModel.T = T;
+    dynModel.H_mat = H_mat;
+    dynModel.tkp = tkp;
+    dynModel.Aeq = Aeq;
+    dynModel.A_ineq = A_ineq;
+    dynModel.B_ineq = B_ineq;
+    dynModel.Phip = Phip;
+    dynModel.Phip1 = Phip1;
+    dynModel.Dphi = Dphi;
+    dynModel.mass_scale = m_scale;
+end
+
+function [powTot, fRes, tRes] = getPSPhasePower(dynModel, ph)
+    % getPSPhasePower   calculates power using the pseudospectral
+    % method given a phase and a descrption of the body movement.
+    % Returns total phase power and power per frequency
+
+    eta_fd = dynModel.wave_amp .* exp(1i*ph);
+    E3 = dynModel.Hex .* eta_fd;
+    
+    fef3 = zeros(2*dynModel.Nf,1);
+    
+    fef3(1:2:end) =  real(E3);
+    fef3(2:2:end) = -imag(E3);
+    
+    Beq = [fef3; 0] / dynModel.mass_scale;
+    
+    % constrained optimization settings
+    qp_options = optimoptions('fmincon',  ...
+        'Algorithm', 'sqp',               ...
+        'Display', 'off',                 ...
+        'MaxIterations', 1e3,             ...
+        'MaxFunctionEvaluations', 1e5,    ...
+        'OptimalityTolerance', 1e-8,      ...
+        'StepTolerance', 1e-8);
+    
+    siz = size(dynModel.A_ineq);
+    X0 = zeros(siz(2),1);
+    [y, fval, exitflag, output] = fmincon(@pow_calc,...
+        X0,...
+        dynModel.A_ineq,...
+        dynModel.B_ineq,...
+        dynModel.Aeq,...         % Aeq and Beq are the hydrodynamic model
+        Beq,...
+        [], [], [],...
+        qp_options);
+    
+    %     if exitflag ~= 1      % for debugging
+    %         disp(exitflag)
+    %         disp(output)
+    %     end
+    
+    % y is a column vector containing [vel; u] of the
+    % pseudospectral coefficients
+    tmp = reshape(y(1:end-2),[],2);
+    x1hat = tmp(:,1);
+    uhat = tmp(:,2);
+    
+    % find the spectra
+    ps2spec = @(x) (x(1:2:end) - 1i * x(2:2:end));  % TODO - probably make this a global function
+    velFreq = ps2spec(x1hat);
+    posFreq = velFreq ./ (1i * dynModel.w);
+    uFreq = dynModel.mass_scale * ps2spec(uhat);
+    powFreq = 1/2 * uFreq .* conj(velFreq);
+    zFreq = uFreq ./ velFreq;
+    
+    % find time histories
+    spec2time = @(x) dynModel.Phip' * x;              % TODO - probably make this a global function
+    velT = spec2time(x1hat);
+    posT = y(end-1) + (dynModel.Phip' / dynModel.Dphi) * x1hat;
+    uT = dynModel.mass_scale * (y(end) + spec2time(uhat));
+    powT = 1 * velT .* uT;
+    
+    powTot = trapz(dynModel.tkp, powT) / (dynModel.tkp(end) - dynModel.tkp(1));
+    assert(WecOptTool.math.isClose(powTot, sum(real(powFreq)),...
+        'rtol', eps*1e2),...
+        sprintf('Mismatch in PS results\n\tpowTot: %.3e\n\tpowFreq: %.3e',...
+        powTot,sum(real(powFreq))))
+    
+    % assemble outputs
+    fRes.pos = posFreq;
+    fRes.vel = velFreq;
+    fRes.u = uFreq;
+    fRes.pow = powFreq;
+    fRes.Zpto = zFreq;
+    fRes.eta = eta_fd;
+    fRes.F0 = E3;
+    
+    tRes.pos = posT;
+    tRes.vel = velT;
+    tRes.u = uT;
+    tRes.pow = powT;
+    
+    function P = pow_calc(X)
+        P = X(1:end-2)' * dynModel.H_mat * X(1:end-2); % 1/2 factor dropped for simplicity
+    end
+end
+
+% Copyright 2020 National Technology & Engineering Solutions of Sandia, 
+% LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the 
+% U.S. Government retains certain rights in this software.
+%
+% This file is part of WecOptTool.
+% 
+%     WecOptTool is free software: you can redistribute it and/or modify
+%     it under the terms of the GNU General Public License as published by
+%     the Free Software Foundation, either version 3 of the License, or
+%     (at your option) any later version.
+% 
+%     WecOptTool is distributed in the hope that it will be useful,
+%     but WITHOUT ANY WARRANTY; without even the implied warranty of
+%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+%     GNU General Public License for more details.
+% 
+%     You should have received a copy of the GNU General Public License
+%     along with WecOptTool.  If not, see <https://www.gnu.org/licenses/>.
diff --git a/toolbox/+WecOptTool/SeaState.m b/toolbox/+WecOptTool/SeaState.m
index cc85979..52702ea 100644
--- a/toolbox/+WecOptTool/SeaState.m
+++ b/toolbox/+WecOptTool/SeaState.m
@@ -418,7 +418,7 @@ function checkSpectrum(S)
             %    >>> Hm0 = 5;
             %    >>> Tp = 8;
             %    >>> S = bretschneider([],[Hm0,Tp]);
-            %    >>> WecOptTool.utils.checkSpectrum(S)
+            %    >>> WecOptTool.SeaState.checkSpectrum(S)
             %
             
             arguments