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])
 
  See also  pmspec, torsethaugen, simpson

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 %
052 % See also  pmspec, torsethaugen, simpson
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