Home > wafo > trgauss > dirsp2chitwo.m

dirsp2chitwo

PURPOSE ^

gives parameters in non-central CHI-TWO process for directional Stokes waves.

SYNOPSIS ^

[gam,bet,S12,S22]= dirsp2chitwo(s,w,L0,L2,th,h,eps,dthdw)

DESCRIPTION ^

  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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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