%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.017 Homework 1. Bretschnieder Spectrum.
% FSH MIT Mechanical Engineering
clear all;
figure(1);clf;hold off; hold on; % note: hold is on so we can
% figures
% show the data from the table in the problem
SSvec = [2 3 4 5 6] ; % sea states
wmvec = 2*pi*ones(size(SSvec)) ./ [6.3 7.5 8.8 9.7 12.4] ;
% modal frequencies
Hsigvec = [0.3 0.9 1.9 3.3 5.0] ; % significant wave heights
% vector of frequencies for the spectrum calculation
wvec = [.1:.01:2];
% step through the different seastates
for i = 1:length(wmvec),
    wm = wmvec(i) ;
    Hsig = Hsigvec(i) ;
    SS = SSvec(i) ;
    for j = 1:length(wvec),
        w = wvec(j) ;
        S(j) = 5/16 * wm^4 / w^5 * Hsig^2 * exp(-5 * wm^4 / 4 / w^4) ;
    end;
    % check that we got the right formula!
    disp(sprintf(...
        "Square Root of Integral of Area of S: %g; 4*Hsig: %g", ...
    sqrt(sum(S)*mean(diff(wvec))), 1/4*Hsig));
    % add on the curves and some labels
    plot(wvec,S,"LineWidth",2);
    if SS > 3,
        text(wm-.18,max(S),sprintf("SS %g", SS));
    end;
end;
% finish off the plot
grid;
title("Sea Wave Spectra for Sea States 2-6");
    xlabel("frequency \omega, rad/s");
    ylabel("S(\omega)");