Home > wafo > spec > wallop.m

# wallop

## PURPOSE

Calculates (and plots) a Wallop spectral density.

## SYNOPSIS

S1=wallop(w1,sdata,plotflag)

## DESCRIPTION

 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

## CROSS-REFERENCE INFORMATION

This function calls:
 createspec Spectrum structure constructor w2k Translates from frequency to wave number wspecplot Plot a spectral density gamma Gamma function. linspace Linearly spaced vector. num2str Convert number to string. (Fast version)
This function is called by:

## SOURCE CODE

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 %
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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003