SPEC2SKEW Estimates the moments of 2'nd order non-linear waves CALL: [skew, kurt, mean, sigma] = spec2skew(S,h,method); skew, kurt, mean, sigma = skewness, kurtosis, mean and standard deviation, respectively, of 2'nd order waves to the leading order. (skew=kurt=0 for a Gaussian process) S = spectral density structure h = water depth (default S.h) method = 'approximate' method due to Marthinsen & Winterstein (default) 'eigenvalue' method due to Kac and Siegert Skewness = kurtosis-3 = 0 for a Gaussian process. The mean, sigma, skewness and kurtosis are determined as follows: method == 'approximate': due to Marthinsen and Winterstein mean = 2 * int Hd(w1,w1)*S(w1) dw1 sigma = sqrt(int S(w1) dw1) skew = 6 * int int [Hs(w1,w2)+Hd(w1,w2)]*S(w1)*S(w2) dw1*dw2/m0^(3/2) kurt = (4*skew/3)^2 where Hs = sum frequency effects and Hd = difference frequency effects method == 'eigenvalue' mean = sum(E); sigma = sqrt(sum(C^2)+2*sum(E^2)); skew = sum((6*C^2+8*E^2).*E)/sigma^3; kurt = 3+48*sum((C^2+E^2).*E^2)/sigma^4; where h1 = sqrt(S*dw/2); C = (ctranspose(V)*[h1;h1]); and E and V is the eigenvalues and eigenvectors, respectively, of the 2'order transfer matrix. S is the spectrum and dw is the frequency spacing of S. Example: Simulate a Transformed Gaussian process: Hm0=7;Tp=11; S = jonswap([],[Hm0 Tp]); [sk, ku, me]=spec2skew(S); g=hermitetr([],[Hm0/4 sk ku me]); g2=[g(:,1), g(:,2)*Hm0/4]; ys = spec2sdat(S,15000); % Simulated in the Gaussian world xs = gaus2dat(ys,g2); % Transformed to the real world See also hermitetr, ochitr, lc2tr, dat2tr
returns the constant acceleration of gravity | |
Transforms between different types of spectra | |
Toggle Transform between angular frequency and frequency spectrum | |
Translates from frequency to wave number | |
' Complex conjugate transpose. | |
Difference and approximate derivative. | |
Find a few eigenvalues and eigenvectors of a matrix using ARPACK. | |
Display message and abort function. | |
Convert string to lowercase. | |
X and Y arrays for 3-D plots. | |
Trapezoidal numerical integration. | |
Display warning message; disable or enable warning messages. |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
% CHAPTER4 contains the commands used in Chapter 4 of the tutorial | |
Calculate transformation, g, proposed by Winterstein |
001 function [skew, kurt, ma, sa, Hs ,Hd]=spec2skew(S,h,method) 002 %SPEC2SKEW Estimates the moments of 2'nd order non-linear waves 003 % 004 % CALL: [skew, kurt, mean, sigma] = spec2skew(S,h,method); 005 % 006 % skew, kurt, 007 % mean, sigma = skewness, kurtosis, mean and standard deviation, 008 % respectively, of 2'nd order waves to the leading 009 % order. (skew=kurt=0 for a Gaussian process) 010 % S = spectral density structure 011 % h = water depth (default S.h) 012 % method = 'approximate' method due to Marthinsen & Winterstein (default) 013 % 'eigenvalue' method due to Kac and Siegert 014 % 015 % Skewness = kurtosis-3 = 0 for a Gaussian process. 016 % The mean, sigma, skewness and kurtosis are determined as follows: 017 % 018 % method == 'approximate': due to Marthinsen and Winterstein 019 % mean = 2 * int Hd(w1,w1)*S(w1) dw1 020 % sigma = sqrt(int S(w1) dw1) 021 % skew = 6 * int int [Hs(w1,w2)+Hd(w1,w2)]*S(w1)*S(w2) dw1*dw2/m0^(3/2) 022 % kurt = (4*skew/3)^2 023 % 024 % where Hs = sum frequency effects and Hd = difference frequency effects 025 % 026 % method == 'eigenvalue' 027 % 028 % mean = sum(E); 029 % sigma = sqrt(sum(C^2)+2*sum(E^2)); 030 % skew = sum((6*C^2+8*E^2).*E)/sigma^3; 031 % kurt = 3+48*sum((C^2+E^2).*E^2)/sigma^4; 032 % 033 % where 034 % h1 = sqrt(S*dw/2); 035 % C = (ctranspose(V)*[h1;h1]); 036 % and E and V is the eigenvalues and eigenvectors, respectively, of the 2'order 037 % transfer matrix. S is the spectrum and dw is the frequency spacing of S. 038 % 039 % Example: Simulate a Transformed Gaussian process: 040 % Hm0=7;Tp=11; 041 % S = jonswap([],[Hm0 Tp]); [sk, ku, me]=spec2skew(S); 042 % g=hermitetr([],[Hm0/4 sk ku me]); g2=[g(:,1), g(:,2)*Hm0/4]; 043 % ys = spec2sdat(S,15000); % Simulated in the Gaussian world 044 % xs = gaus2dat(ys,g2); % Transformed to the real world 045 % 046 % See also hermitetr, ochitr, lc2tr, dat2tr 047 048 % References: 049 % Langley, RS (1987) 050 % 'A statistical analysis of nonlinear random waves' 051 % Ocean Engineering, Vol 14, No 5, pp 389-407 052 % 053 % Marthinsen, T. and Winterstein, S.R (1992) 054 % 'On the skewness of random surface waves' 055 % In proceedings of the 2nd ISOPE Conference, San Francisco, 14-19 june. 056 % 057 % Winterstein, S.R, Ude, T.C. and Kleiven, G. (1994) 058 % "Springing and slow drift responses: 059 % predicted extremes and fatigue vs. simulation" 060 % In Proc. 7th International behaviour of Offshore structures, (BOSS) 061 % Vol. 3, pp.1-15 062 063 % tested on: matlab 5.2 064 % History 065 % revised pab 22.07.2002 066 % -fixed a bug in the calculaton of the mean for the 2'nd order process. 067 % - added method 068 % -updated the help header accordingly 069 % revised pab 12.11.2001 070 % - added Langley's version of the quadratic transfer functions. 071 % revised pab 09.01.2001 072 % -simplified calculation of quadratic transfer functions when h=inf 073 % revised pab 09.01.2001 074 % - changed kurtosis so that kurtosis correspond to what wkurtosis 075 % measures, i.e., kurtosis-3=0 for a Gaussian process 076 % by pab 01.03.2000 077 078 079 error(nargchk(1,3,nargin)) 080 081 % default options 082 opts.disp = 0; 083 084 if nargin<2|isempty(h), h = S.h; end 085 if nargin<3|isempty(method), method = 'approximate'; end 086 087 S = spec2spec(S,'freq'); 088 S = ttspec(S,'w'); 089 090 091 S1 = S.S(:); 092 m0 = trapz(S.w,S1); 093 Nw = length(S.w); 094 095 [Hs, Hd,Hdii] = qtf(S.w(:),h); 096 097 %return 098 %skew=6/sqrt(m0)^3*simpson(S.w,simpson(S.w,(Hs+Hd).*S1(:,ones(1,Nw))).*S1.'); 099 100 Hspd = trapz(S.w,trapz(S.w,(Hs+Hd).*S1(:,ones(1,Nw))).*S1.'); 101 switch lower(method(1)) 102 case 'a', %approx : Marthinsen, T. and Winterstein, S.R (1992) method 103 if nargout>2 104 ma = 2*trapz(S.w,Hdii.*S1); 105 end 106 sa = sqrt(m0); 107 skew = 6/sa^3*Hspd; 108 kurt = (4*skew/3).^2+3; 109 110 case 'q', % quasi method 111 dw = diff(S.w(1:2)); 112 tmp1 =sqrt(S1(:,ones(1,Nw)).*(S1(:,ones(1,Nw)).'))*dw; 113 Hd = Hd.*tmp1; 114 Hs = Hs.*tmp1; 115 k = 6; 116 stop = 0; 117 while (~stop) 118 E = eigs([Hd,Hs;Hs,Hd],[],k); 119 stop = (length(find(abs(E)<1e-4))>0 | k>1200); 120 k = 2*k; 121 end 122 123 124 m02=2*sum(E.^2); % variance of 2'nd order contribution 125 126 %Hstd = 16*trapz(S.w,(Hdii.*S1).^2); 127 %Hstd = trapz(S.w,trapz(S.w,((Hs+Hd)+ 2*Hs.*Hd).*S1(:,ones(1,Nw))).*S1.'); 128 ma = 2*trapz(S.w,Hdii.*S1); 129 %m02 = Hstd-ma^2% variance of second order part 130 sa = sqrt(m0+m02); 131 skew = 6/sa^3*Hspd; 132 kurt = (4*skew/3).^2+3; 133 case 'e', % Kac and Siegert eigenvalue analysis 134 dw = diff(S.w(1:2)); 135 tmp1 =sqrt(S1(:,ones(1,Nw)).*(S1(:,ones(1,Nw)).'))*dw; 136 Hd = Hd.*tmp1; 137 Hs = Hs.*tmp1; 138 k = 6; 139 stop = 0; 140 E2 = 1; 141 142 while (~stop) 143 [V,D] = eigs([Hd,Hs;Hs,Hd],[],k); 144 E = diag(D); 145 stop = (length(find(abs(E)<1e-4))>0 | k>=min(2*Nw,1200)); 146 k = min(2*k,2*Nw); 147 end 148 149 150 h1 = sqrt(S1*dw/2); 151 C = (ctranspose(V)*[h1;h1]) 152 153 E2 = E.^2; 154 C2 = C.^2; 155 156 ma = sum(E); % mean 157 sa = sqrt(sum(C2)+2*sum(E2)); % standard deviation 158 skew = sum((6*C2+8*E2).*E)/sa^3; % skewness 159 kurt = 3+48*sum((C2+E2).*E2)/sa^4; % kurtosis 160 otherwise error('Method is not available') 161 end 162 163 164 165 return 166 167 function [Hs, Hd,Hdii]=qtf(w,h) 168 % QTF Quadratic Transfer Function 169 % 170 % CALL: [Hs, Hd, Hdii]=qtf(w,h) 171 % 172 % Hs = sum frequency effects 173 % Hd = difference frequency effects 174 % Hdii = diagonal of Hd 175 % w = angular frequencies 176 % h = water depth 177 178 179 Nw = length(w); 180 g = gravity; 181 kw = w2k(w,0,h,g); 182 [k1, k2] = meshgrid(kw); 183 184 185 186 Hd = zeros(size(k1)); 187 if h==inf,% go here for faster calculations 188 Hs = 0.25*(abs(k1)+abs(k2)); 189 Hd = -0.25*abs(abs(k1)-abs(k2)); 190 Hdii = zeros(size(w)); 191 return 192 end 193 194 [w1, w2]=meshgrid(w); 195 196 warning off %off % suppress warnings on division by zero 197 198 w12 = (w1.*w2); 199 w1p2 = (w1+w2); 200 w1m2 = (w1-w2); 201 k12 = (k1.*k2); 202 k1p2 = (k1+k2); 203 k1m2 = abs(k1-k2); 204 if 0, % Langley 205 p1 = (-2*w1p2.*(k12*g^2-w12.^2)+... 206 w1.*(w2.^4-g^2*k2.^2)+w2.*(w1.^4-g^2*k1.^2))./(4.*w12); 207 p2= w1p2.^2.*cosh((k1p2).*h)-g*(k1p2).*sinh((k1p2).*h); 208 209 Hs = -p1./p2.*w1p2.*cosh((k1p2).*h)/g-... 210 (k12*g^2-w12.^2)./(4*g*w12)+(w1.^2+w2.^2)/(4*g); 211 212 p3 = (-2*w1m2.*(k12*g^2+w12.^2)-... 213 w1.*(w2.^4-g^2*k2.^2)+w2.*(w1.^4-g^2*k1.^2))./(4.*w12); 214 p4= w1m2.^2.*cosh(k1m2.*h)-g*(k1m2).*sinh((k1m2).*h); 215 216 217 Hd = -p3./p4.*(w1m2).*cosh((k1m2).*h)/g-... 218 (k12*g^2+w12.^2)./(4*g*w12)+(w1.^2+w2.^2)/(4*g); 219 220 else % Marthinsen & Winterstein 221 tmp1 = 0.5*g*k12./w12; 222 tmp2 = 0.25/g*(w1.^2+w2.^2+w12); 223 Hs = (tmp1-tmp2+0.25*g*(w1.*k2.^2+w2.*k1.^2)./(w12.*(w1p2)))..... 224 ./(1-g*(k1p2)./(w1p2).^2.*tanh((k1p2).*h))+tmp2-0.5*tmp1; % OK 225 226 tmp2 = 0.25/g*(w1.^2+w2.^2-w12); %OK 227 Hd = (tmp1-tmp2-0.25*g*(w1.*k2.^2-w2.*k1.^2)./(w12.*(w1m2)))..... 228 ./(1-g*(k1m2)./(w1m2).^2.*tanh((k1m2).*h))+tmp2-0.5*tmp1; % OK 229 end 230 231 tmp1 = 0.5*g*kw./(w.*sqrt(g*h)); 232 tmp2 = 0.25*w.^2/g; 233 234 235 Cg = 0.5*g*(tanh(kw*h) +kw*h.*(1- tanh(kw*h).^2))./w; %Wave group velocity 236 Hdii = 0.5*(0.5*g*(kw./w).^2-0.5*w.^2/g+g*kw./(w.*Cg)).... 237 ./(1-g*h./Cg.^2)-0.5*kw./sinh(2*kw*h); % OK 238 Hd(1:Nw+1:end) = Hdii; 239 240 %k = find(w1==w2); 241 %Hd(k) = Hdii; 242 243 % The NaN's occur due to division by zero. => Set the isnans to zero 244 Hdii(isnan(Hdii))=0; 245 Hd(isnan(Hd))=0; 246 Hs(isnan(Hs))=0; 247 248 warning on 249 return 250 251 252
Comments or corrections to the WAFO group