DIRSP2CHITWO gives parameters in non-central CHI-TWO process for directional Stokes waves. CALL: [gamma,beta,S12,S22]= dirsp2chitwo(s,w,L0,L2,th,h,eps,dwth); s,w,th = spectral density s(w,th), where w is a vector with angular frequencies th wave directions (frequences and angles are equally spaced) Note that usually the spectrum s is truncated to exclude low and high frequencies. The degree of truncation is measured by using the spectral moments L0, L2 for untruncated spectrum. L0,L2 = Two first spectral moments of the linear spectrum. These can be smaller than the corresponding spectral moments of s. h = water depth (default: 5000 [m]). eps = all eigenvalues which have absolutvalue below eps*variance of the CHI2 process are replaced by zero. dwth = spacing in w and th vectors dtwh=dw*dth gamma,beta,S12,S22 parameters in CHI-TWO model.
Computes values of the quadratic transfer function E, for quadratic sea | |
returns the constant acceleration of gravity | |
Calculates (and plots) a JONSWAP spectral density | |
Eigenvalues and eigenvectors. |
001 function [gam,bet,S12,S22]= dirsp2chitwo(s,w,L0,L2,th,h,eps,dthdw) 002 % DIRSP2CHITWO gives parameters in non-central CHI-TWO process for directional Stokes waves. 003 % 004 % CALL: [gamma,beta,S12,S22]= dirsp2chitwo(s,w,L0,L2,th,h,eps,dwth); 005 % 006 % s,w,th = spectral density s(w,th), where w is a vector with angular frequencies 007 % th wave directions (frequences and angles are equally spaced) Note that 008 % usually the spectrum s is truncated to exclude low and high frequencies. 009 % The degree of truncation is measured by using the spectral 010 % moments L0, L2 for untruncated spectrum. 011 % L0,L2 = Two first spectral moments of the linear spectrum. These can 012 % be smaller than the corresponding spectral moments of s. 013 % h = water depth (default: 5000 [m]). 014 % eps = all eigenvalues which have absolutvalue below eps*variance of the 015 % CHI2 process are replaced by zero. 016 % dwth = spacing in w and th vectors dtwh=dw*dth 017 % 018 % gamma,beta,S12,S22 parameters in CHI-TWO model. 019 % 020 021 022 % References: U. Machado, I. Rychlik (2002) "Wave statistics in nonlinear sea" to 023 % appear in Extremes. 024 % Butler, R., Machado, U. Rychlik, I. (2002) "Distribution of wave crests in non- 025 % linear random sea - application of saddlepoint methods" by , presented at ISOPE 2003. 026 % Calls: ewwdir 027 % By I.R 24.10.04 028 % 029 %------------------------------------------------------------------------------------ 030 031 032 033 if nargin<2 034 Spec=jonswap; 035 s=Spec.S; 036 w=Spec.w; 037 th=zeros(size(w)); 038 h=Spec.h; 039 nb=0; 040 dthdw=w(3)-w(2); 041 end 042 if nargin<4 043 dthdw=w(3)-w(2); 044 L0=dthdw*sum(s); 045 L2=dthdw*sum(w.^2.*s); 046 end 047 if nargin<5 048 th=zeros(size(w)); 049 end 050 if nargin<6 051 h=5000; 052 end 053 if (h>5000) 054 h=5000; 055 end 056 if nargin<7 057 eps=0.00001; 058 end 059 060 if nargin<8 061 dth=th(3)-th(2); 062 dthdw=(w(3)-w(2))*dth; 063 end 064 if (dthdw<0.0000000001) 065 dthdw=(w(3)-w(2)); 066 end 067 068 g=gravity; 069 omega=w; 070 spec=sqrt(s); 071 spec=spec*spec'; 072 073 % if narrow-band 074 % [i,j]=max(Spec.S); 075 % w=Spec.w(j)*ones(size(Spec.w)); 076 077 078 %=== Computation of the transfer functions Q, R and S 079 080 W=diag(-w); 081 Q=(ewwdir(w,th,-w,th,h)+ewwdir(w,th,w,th,h)).*spec*dthdw; 082 R=(ewwdir(w,th,-w,th,h)-ewwdir(w,th,w,th,h)).*spec*dthdw; 083 084 085 %=== 086 Q=(Q+Q'); 087 R=(R+R'); 088 S=Q*W-W*R; 089 sigma=sqrt(s*dthdw); 090 N=length(W); 091 092 093 %=== 094 [P1c,Delta]=eig(Q); 095 P1=P1c'; 096 [P2c,Gama]=eig(R); 097 P2=P2c'; 098 099 100 101 gam=[diag(Delta,0)' diag(Gama,0)']/2; 102 variance=2*(sum(gam.^2)); 103 104 105 % computations of beta 106 107 sigma=sqrt(s*dthdw); 108 bet=[(P1*sigma)' zeros(1,N)]; 109 110 SS=eye(2*N); 111 Z=zeros(N); 112 SS12=[Z -W;W Z]; 113 SS22=[W.^2 Z; Z W.^2]; 114 PP=[P1 Z;Z P2]; 115 S12=PP*SS12*PP'; 116 S22=PP*SS22*PP'; 117 118 [bet*SS*bet' bet*S22*bet']; 119 120 % 121 % In this part of program we are removing the quadratic processes with 122 % negligable gamma_i coefficients.% 123 124 n=2*N; 125 if (eps>0) 126 [gammasort indexgamma]=sort(abs(gam)); 127 128 gammasort=gam(indexgamma); 129 betasort=bet(indexgamma); 130 S12sort=S12(indexgamma,indexgamma); 131 S22sort=S22(indexgamma,indexgamma); 132 betasort*S22sort*betasort'; 133 varapprox=2*(cumsum(gammasort.^2))/variance; 134 mm=sum(varapprox<eps); 135 if (mm>0) 136 beta1=betasort(1:mm); 137 bet=[sqrt(sum(beta1.^2)) betasort(mm+1:end)]; 138 gam=[0 gammasort(mm+1:end)]; 139 S12=zeros(n-mm+1,n-mm+1); 140 S22=zeros(n-mm+1,n-mm+1); 141 S22(1,1)=beta1*S22sort(1:mm,1:mm)*beta1'/bet(1)^2; 142 S22(1,2:end)=beta1*S22sort(1:mm,mm+1:end)/bet(1); 143 S22(2:end,2:end)=S22sort(mm+1:end,mm+1:end); 144 S12(2:end,2:end)=S12sort(mm+1:end,mm+1:end); 145 S12(1,2:end)=beta1*S12sort(1:mm,mm+1:end)/bet(1); 146 S12(2:end,1)=-S12(1,2:end)'; 147 S22(2:end,1)=S22(1,2:end)'; 148 end 149 end 150 dL0=L0-sum(s)*dthdw; 151 dL2=L2-(s'*w.^2)*dthdw; 152 bet(1)=sqrt(bet(1).^2+dL0); 153 if(dL2>0.001) 154 S22(1,1)=S22(1,1)+dL2/dL0; 155 end
Comments or corrections to the WAFO group