function [BS_we,M0,M1,M2,M3,M4]=BS(wm,we,H_sig)
%BS calculates a Bretschneider wave spectrum based on:
%modal frequency= wm
%significant wave height= H

w=we;

BS_we= (1.25/4).*((wm^4)./(w.^5)).*(H_sig^2).*exp(-1.25.*(wm./w).^4);


M0 = trapz(w,BS_we);
M1 = trapz(w,BS_we.*w);
M2 = trapz(w,BS_we.*w.^2);
M3 = trapz(w,BS_we.*w.^3);
M4 = trapz(w,BS_we.*w.^4);