Home > wafo > spec > jonswap.m

# jonswap

## PURPOSE

Calculates (and plots) a JONSWAP spectral density

## SYNOPSIS

S1 = jonswap(w1,sdata,plotflag)

## DESCRIPTION

``` JONSWAP Calculates (and plots) a JONSWAP spectral density

CALL:  S = jonswap(w,sdata,plotflag);
S = jonswap(wc,sdata,plotflag);

S     = a struct containing the spectral density. See datastructures
w     = angular frequency        (default linspace(0,wc,257))
wc    = angular cutoff frequency (default 33/Tp)
sdata = [Hm0 Tp gamma sa sb A], where
Hm0   = significant wave height (default 7 (m))
Tp    = peak period (default 11 (sec))
gamma = peakedness factor determines the concentraton
of the spectrum on the peak frequency,  1 <= gamma <= 7.
(default depending on Hm0, Tp, see below)
sa,sb = spectral width parameters (default 0.07 0.09)
A     = alpha, normalization factor, (default -1)
A<0 : A calculated by integration so that int S dw =Hm0^2/16
A==0 : A = 5.061*Hm0^2/Tp^4*(1-0.287*log(gamma))
A>0  : A = A
plotflag = 0, do not plot the spectrum (default).
1, plot the spectrum.

For zero values, NaN's or parameters not specified in DATA the
default values are used.

S(w) = A*g^2/w^5*exp(-5/4(wp/w)^4)*j^exp(-.5*((w/wp-1)/s)^2)
where
s   = sa w<=wp
sb w>wp (wp = angular peak frequency)
j   = gamma,  (j=1, => Bretschneider spectrum)

This spectrum is assumed to be especially suitable for the North Sea,
and does not represent a fully developed sea. It is a reasonable model for
wind generated sea when 3.6*sqrt(Hm0) < Tp < 5*sqrt(Hm0)
A standard value for gamma is 3.3. However, a more correct approach is
to relate gamma to Hm0:
D = 0.036-0.0056*Tp/sqrt(Hm0);
gamma = exp(3.484*(1-0.1975*D*Tp^4/(Hm0^2)));
This parameterization is based on qualitative considerations of deep water
wave data from the North Sea, see Torsethaugen et. al. (1984)
Here gamma is limited to 1..7.

The relation between the peak period and mean zero-upcrossing period
may be approximated by
Tz = Tp/(1.30301-0.01698*gamma+0.12102/gamma)

Example:  % Bretschneider spectrum Hm0=7, Tp=11
S = jonswap([],[0 0 1])

## CROSS-REFERENCE INFORMATION

This function calls:
 createspec Spectrum structure constructor getjonswappeakedness Peakedness factor Gamma given Hm0 and Tp for JONSWAP gravity returns the constant acceleration of gravity simpson Numerical integration with the Simpson method wspecplot Plot a spectral density linspace Linearly spaced vector. num2str Convert number to string. (Fast version)
This function is called by:
 Chapter2 % CHAPTER2 Modelling random loads and stochastic waves Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics Chapter4 % CHAPTER4 contains the commands used in Chapter 4 of the tutorial demospec Loads a precreated spectrum of chosen type dirsp2chitwo gives parameters in non-central CHI-TWO process for directional Stokes waves. itmkurs_lab3 Script to computer exercises 3 jhvpdf Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. mkdspec Make a directional spectrum specq2lc Saddlepoint approximation of the crossing intensity for the quadratic sea. test_cycles Quick test of the routines in module 'cycles' tmaspec Calculates a JONSWAP spectral density for finite water depth wafofig7 Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac wafofig8 Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac for extremal waves

## SOURCE CODE

```001 function S1 = jonswap(w1,sdata,plotflag)
002 %JONSWAP Calculates (and plots) a JONSWAP spectral density
003 %
004 % CALL:  S = jonswap(w,sdata,plotflag);
005 %        S = jonswap(wc,sdata,plotflag);
006 %
007 %   S     = a struct containing the spectral density. See datastructures
008 %   w     = angular frequency        (default linspace(0,wc,257))
009 %   wc    = angular cutoff frequency (default 33/Tp)
010 %   sdata = [Hm0 Tp gamma sa sb A], where
011 %          Hm0   = significant wave height (default 7 (m))
012 %          Tp    = peak period (default 11 (sec))
013 %          gamma = peakedness factor determines the concentraton
014 %                  of the spectrum on the peak frequency,  1 <= gamma <= 7.
015 %                  (default depending on Hm0, Tp, see below)
016 %          sa,sb = spectral width parameters (default 0.07 0.09)
017 %          A     = alpha, normalization factor, (default -1)
018 %                  A<0 : A calculated by integration so that int S dw =Hm0^2/16
019 %                  A==0 : A = 5.061*Hm0^2/Tp^4*(1-0.287*log(gamma))
020 %                  A>0  : A = A
021 %       plotflag = 0, do not plot the spectrum (default).
022 %                  1, plot the spectrum.
023 %
024 %  For zero values, NaN's or parameters not specified in DATA the
025 %  default values are used.
026 %
027 %
028 %         S(w) = A*g^2/w^5*exp(-5/4(wp/w)^4)*j^exp(-.5*((w/wp-1)/s)^2)
029 %    where
030 %         s   = sa w<=wp
031 %               sb w>wp (wp = angular peak frequency)
032 %         j   = gamma,  (j=1, => Bretschneider spectrum)
033 %
034 %  This spectrum is assumed to be especially suitable for the North Sea,
035 %  and does not represent a fully developed sea. It is a reasonable model for
036 %  wind generated sea when 3.6*sqrt(Hm0) < Tp < 5*sqrt(Hm0)
037 %  A standard value for gamma is 3.3. However, a more correct approach is
038 %  to relate gamma to Hm0:
039 %        D = 0.036-0.0056*Tp/sqrt(Hm0);
040 %        gamma = exp(3.484*(1-0.1975*D*Tp^4/(Hm0^2)));
041 %  This parameterization is based on qualitative considerations of deep water
042 %  wave data from the North Sea, see Torsethaugen et. al. (1984)
043 %  Here gamma is limited to 1..7.
044 %
045 %  The relation between the peak period and mean zero-upcrossing period
046 %  may be approximated by
047 %         Tz = Tp/(1.30301-0.01698*gamma+0.12102/gamma)
048 %
049 % Example:  % Bretschneider spectrum Hm0=7, Tp=11
050 %      S = jonswap([],[0 0 1])
051 %
053
054 % References:
055 % Torsethaugen et al. (1984)
056 % Characteristica for extreme Sea States on the Norwegian continental shelf.
057 % Report No. STF60 A84123. Norwegian Hydrodyn. Lab., Trondheim
058 %
059 % Hasselman et al. (1973)
060 % Measurements of Wind-Wave Growth and Swell Decay during the Joint
061 % North Sea Project (JONSWAP).
062 % Ergansungsheft, Reihe A(8), Nr. 12, Deutschen Hydrografischen Zeitschrift.
063
064
065 % Tested on: matlab 6.0, 5.3
066 % History:
067 % revised pab June 2005
068 % -fixed a bug in help header: the jonswap range is now correct
069 % revised pab 11jan2004
070 % - replaced code with call to getjonswappeakedness
071 % revised jr 22.08.2001
072 % - correction in formula for S(w) in help: j^exp(-.5*((w-wp)/s*wp)^2)
073 %   (the first minus sign added)
074 % revised pab 01.04.2001
075 % - added wc to input
076 % revised jr 30.10 2000
077 %   - changed 'data' to 'sdata' in the function call
078 % revised pab 20.09.2000
079 %   - changed default w: made it dependent on Tp
080 % revised es 25.05.00
081 %   - revision of help text
082 % revised pab 16.02.2000
083 %  - fixed a bug for sa,sb and the automatic calculation of gamma.
084 %  - added sa, sb, A=alpha to data input
085 %  - restricted values of gamma to 1..7
086 % revised by pab 01.12.99
087 %  added gamma to data input
088 % revised by pab 11.08.99
089 % changed so that parameters are only dependent on the
090 % seastate parameters Hm0 and Tp.
091 % also checks if Hm0 and Tp are reasonable.
092
093
094 %NOTE: In order to calculate the short term statistics of the response,
095 %      it is extremely important that the resolution of the transfer
096 %      function is sufficiently good. In addition, the transfer function
097 %      must cover a sufficietn range of wave periods, especially in the
098 %      range where the wave spectrum contains most of its
099 %      energy. VIOLATION OF THIS MAY LEAD TO MEANINGLESS RESULTS FROM THE
100 %      CALCULATIONS OF SHORT TERM STATISTICS. The highest wave period
101 %      should therefore be at least 2.5 to 3 times the highest peak
102 %      period in the transfer function. The lowest period should be selected
103 %      so that the transfer function value is low. This low range is
104 %      especially important when studying velocities and accelerations.
105
106 monitor=0;
107
108 if nargin<3|isempty(plotflag),  plotflag=0;end
109
110 Hm0=7;Tp=11; gam=0; sa=0.07; sb=0.09; A=-1;% default values
111 data2=[Hm0 Tp gam sa sb A];
112 nd2=length(data2);
113 if (nargin>1) & ~isempty(sdata),
114   nd=length(sdata);
115   ind=find(~isnan(sdata(1:min(nd,nd2))));
116   if any(ind) % replace default values with those from input data
117     data2(ind)=sdata(ind);
118   end
119 end
120 if (nd2>0) & (data2(1)>0),
121   Hm0 = data2(1);
122 end
123 if (nd2>1) & (data2(2)>0),
124   Tp = data2(2);
125 end
126 if (nd2>2) & (data2(3)>=1) & (data2(3)<=7),
127   gam = data2(3);
128 end
129 if (nd2>3) & (data2(4)>0),
130   sa = data2(4);
131 end
132 if (nd2>4) & (data2(5)>0),
133   sb = data2(5);
134 end
135 if (nd2>5) ,
136   A = data2(6);
137 end
138
139 w = [];
140 if nargin<1|isempty(w1),
141   wc = 33/Tp;
142 elseif length(w1)==1,
143   wc = w1;
144 else
145   w = w1 ;
146 end
147 nw = 257;
148 if isempty(w),
149   w = linspace(0,wc,nw).';
150 end
151
152
153 n=length(w);
154 S1=createspec;
155 S1.S=zeros(n,1);
156 S1.w=w(:);
157 S1.norm=0; % The spectrum is not normalized
158 S1.note=['JONSWAP, Hm0 = ' num2str(Hm0)  ', Tp = ' num2str(Tp)];
159
160 M=4;
161 N=5;
162 wp=2*pi/Tp;
163
164
165 if gam<1
166   gam = getjonswappeakedness(Hm0,Tp);
167 end
168 S1.note=[S1.note ', gamma = ' num2str(gam)];
169 %end
170
171
172
173 if Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0)
174   disp('Warning: Hm0,Tp is outside the JONSWAP range')
175   disp('The validity of the spectral density is questionable')
176 end
177 if gam>7|gam<1
178   disp('Warning: gamma is outside the valid range')
179   disp('The validity of the spectral density is questionable')
180 end
181
182
183 % for w>wp
184 k=(w>wp);
185 S1.S(k)=1./(w(k).^N).*(gam.^(exp(-(w(k)/wp-1).^2 ...
186                 /(2*sb^2)))).*exp(-N/M*(wp./w(k)).^M);
187 % for 0<w<=wp
188 k=~k;
189 k(1)=k(1)*(w(1)>0); % avoid division by zero
190 S1.S(k)=1./(w(k).^N).*(gam.^(exp(-(w(k)/wp-1).^2 ...
191     /(2*sa^2)))).*exp(-N/M*(wp./w(k)).^M);
192
193
194 g=gravity; % acceleration of gravity
195
196 if A<0, % normalizing by integration
197   A=(Hm0/g)^2/16/simpson(w,S1.S);% make sure m0=Hm0^2/16=int S(w)dw
198 elseif A==0,% original normalization
199   % NOTE: that  Hm0^2/16 generally is not equal to intS(w)dw
200   %       with this definition of A if sa or sb are changed from the
201   %       default values
202   A=5.061*Hm0^2/Tp^4*(1-0.287*log(gam)); % approx D
203 end
204
205 S1.S=S1.S*A*g^2; %normalization
206
207 if monitor
208   D=max(0,0.036-0.0056*Tp/sqrt(Hm0)); % approx 5.061*Hm0^2/Tp^4*(1-0.287*log(gam));
209   disp(['sa, sb       = ' num2str([sa sb])])
210   disp(['alpha, gamma = ' num2str([A gam])])
211   disp(['Hm0, Tp      = ' num2str([Hm0 Tp])])
212   disp(['D            = ' num2str(D)])
213 end
214
215 if plotflag
216   wspecplot(S1,plotflag)
217 end
218```

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