WALLOP Calculates (and plots) a Wallop spectral density. CALL: S = wallop(w,data,plotflag); S = wallop(wc,data,plotflag); S = a struct containing the spectral density, see datastructures. w = angular frequency (default linspace(0,3,257)) wc = angular cutoff frequency (default 33/Tp) data = [Hm0 Tp M] Hm0 = significant wave height (default 7 (m)) Tp = peak period (default 11 (sec)) M = shape factor, i.e. slope for the high frequency part (default depending on Hm0 and Tp, see below) plotflag = 0, do not plot the spectrum (default). 1, plot the spectrum. The WALLOP spectrum parameterization used is S(w)=Bw*Hm0^2/wp*(wp./w).^M.*exp(-M/4*(wp./w).^4); where Bw = normalization factor M = abs((log(2*pi^2)+2*log(Hm0/4)-2*log(Lp))/log(2)); Lp = wave length corresponding to the peak frequency, wp. If M=5 it becomes the same as the JONSWAP spectrum with peak enhancement factor gamma=1 or the Pierson-Moskowitz spectrum. Example: S = wallop(1.1,[6.5 10]), wspecplot(S) See also jonswap, torsethaugen, simpson
Spectrum structure constructor | |
Translates from frequency to wave number | |
Plot a spectral density | |
Gamma function. | |
Linearly spaced vector. | |
Convert number to string. (Fast version) |
001 function S1=wallop(w1,sdata,plotflag) 002 %WALLOP Calculates (and plots) a Wallop spectral density. 003 % 004 % CALL: S = wallop(w,data,plotflag); 005 % S = wallop(wc,data,plotflag); 006 % 007 % S = a struct containing the spectral density, see datastructures. 008 % w = angular frequency (default linspace(0,3,257)) 009 % wc = angular cutoff frequency (default 33/Tp) 010 % data = [Hm0 Tp M] 011 % Hm0 = significant wave height (default 7 (m)) 012 % Tp = peak period (default 11 (sec)) 013 % M = shape factor, i.e. slope for the high frequency 014 % part (default depending on Hm0 and Tp, see below) 015 % plotflag = 0, do not plot the spectrum (default). 016 % 1, plot the spectrum. 017 % 018 % The WALLOP spectrum parameterization used is 019 % 020 % S(w)=Bw*Hm0^2/wp*(wp./w).^M.*exp(-M/4*(wp./w).^4); 021 % where 022 % Bw = normalization factor 023 % M = abs((log(2*pi^2)+2*log(Hm0/4)-2*log(Lp))/log(2)); 024 % Lp = wave length corresponding to the peak frequency, wp. 025 % 026 % If M=5 it becomes the same as the JONSWAP spectrum with 027 % peak enhancement factor gamma=1 or the Pierson-Moskowitz spectrum. 028 % 029 % Example: 030 % S = wallop(1.1,[6.5 10]), wspecplot(S) 031 % 032 % See also jonswap, torsethaugen, simpson 033 034 % References: 035 % Huang, N.E., Long, S.R., Tung, C.C, Yuen, Y. and Bilven, L.F. (1981) 036 % "A unified two parameter wave spectral model for a generous sea state" 037 % J. Fluid Mechanics, Vol.112, pp 203-224 038 039 % Tested on: matlab 6.0, 5.3 040 % History: 041 % revised jr 03.04.2001 042 % - added wc to input 043 % - updated information 044 % revised pab 18.02.2000 045 % - normalization so that int S(w) dw = m0 046 % revised pab 24.01.2000 047 % - updated note to 'Wallop Hm0='.... 048 % by pab 01.12.99 049 050 monitor=0; 051 052 if nargin<3|isempty(plotflag) 053 plotflag=0; 054 end 055 056 % Old call 057 %if nargin<1|isempty(w1) 058 % w=linspace(0,3,257).'; 059 %end 060 061 M=[]; 062 if nargin<2|isempty(sdata) 063 sdata=[7 11]; 064 else 065 switch length(sdata) 066 case 1, sdata=[sdata 11]; 067 case 3, M=sdata(3); sdata=sdata(1:2); 068 end 069 end % 070 071 if nargin<1|isempty(w1), wc = 33/sdata(2) 072 elseif length(w1)==1, wc = w1; 073 else w = w1 ; end 074 nw = 257; 075 if isempty(w), w = linspace(0,wc,nw).'; end 076 077 n=length(w); 078 079 S1=createspec; 080 S1.S=zeros(n,1); 081 S1.w=w; 082 S1.norm=0; % The spectrum is not normalized 083 084 Hm0=sdata(1); 085 Tp=sdata(2); 086 S1.note=['Wallop, Hm0 = ' num2str(Hm0) ', Tp = ' num2str(Tp)]; 087 wp=2*pi/Tp; 088 089 if monitor 090 disp(['Hm0, Tp = ' num2str([Hm0 Tp])]) 091 end 092 093 if isempty(M), 094 kp=w2k(wp,0,inf); % wavenumber at peak frequency 095 Lp=2*pi/kp; % wave length at the peak frequency 096 M=abs((log(2*pi^2)+2*log(Hm0/4)-2*log(Lp))/log(2)); 097 end 098 099 %Bw=0.06238*M^((M-1)/4)/(4^((M-5)/4)*gamma((M-1)/4))*(1+0.7458*(M+2)^(-1.057)); 100 101 Bw = M^((M-1)/4)/(4^((M-5)/4)*gamma((M-1)/4))/16; 102 103 % for w>0 % avoid division by zero 104 k=find(w>0); 105 S1.S(k)=Bw*Hm0^2/wp*(wp./w(k)).^M.*exp(-M/4*(wp./w(k)).^4); 106 107 if plotflag 108 wspecplot(S1,plotflag) 109 end 110 111
Comments or corrections to the WAFO group