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
Spectrum structure constructor | |
Peakedness factor Gamma given Hm0 and Tp for JONSWAP | |
returns the constant acceleration of gravity | |
Numerical integration with the Simpson method | |
Plot a spectral density | |
Linearly spaced vector. | |
Convert number to string. (Fast version) |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
% CHAPTER4 contains the commands used in Chapter 4 of the tutorial | |
Loads a precreated spectrum of chosen type | |
gives parameters in non-central CHI-TWO process for directional Stokes waves. | |
Script to computer exercises 3 | |
Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. | |
Make a directional spectrum | |
Saddlepoint approximation of the crossing intensity for the quadratic sea. | |
Quick test of the routines in module 'cycles' | |
Calculates a JONSWAP spectral density for finite water depth | |
Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac | |
Joint distribution (pdf) of crest wavelength, Lc, and crest amplitude, Ac for extremal waves |
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
Comments or corrections to the WAFO group