Home > wafo > trgauss > specq2lc.m

specq2lc

PURPOSE ^

Saddlepoint approximation of the crossing intensity for the quadratic sea.

SYNOPSIS ^

[mu,fu,Fu,muG,mu1,mu2]= specq2lc(Spec,ds0,nb)

DESCRIPTION ^

 SPECQ2LC  Saddlepoint approximation of the crossing intensity for the quadratic sea.
 
   CALL:  [mu,fu,Fu]= specq2lc(S,ds);
   
      mu   = three column vector containing: levels in the first column, 
             the saddlepoint approximation of crossing intensity for quadratic sea  
             in the second column and for Gaussian sea in the third column.
             Note that the levels $u$ are not equally spaced. If the spacing 
             in the grid is not fine  enough, choose smaller value of parameter  ds. 
       fu  = the pdf structure: f.x contains levels u, f.f is a two collon matrix containing
             the Danniel's-saddlepoint approximation of the density (pdf) of the quadratic 
             sea in the first collon and pdf of the Gaussian sea in the second.  
       Fu  = the two collon matrix containing levels u and the Lugannani Rice's-saddlepoint 
             approximation of the cdf of the quadratic sea.  
      S    = spectral density structure, computed using WAFO. S.h contains the water depth.
             (default: Jonswap with depth 5000 [m]).
      ds   = parameter defining the levels: default value is ds=0.1.
 
    Example: S=Jonswap; S.h=40;
             [mu, fu, Fu]=specq2lc(S,0.05); semilogy(mu(:,1),mu(:,2:3))
             pdfplot(fu)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [mu,fu,Fu,muG,mu1,mu2]= specq2lc(Spec,ds0,nb)
0002 %SPECQ2LC  Saddlepoint approximation of the crossing intensity for the quadratic sea.
0003 %
0004 %  CALL:  [mu,fu,Fu]= specq2lc(S,ds);
0005 %  
0006 %     mu   = three column vector containing: levels in the first column, 
0007 %            the saddlepoint approximation of crossing intensity for quadratic sea  
0008 %            in the second column and for Gaussian sea in the third column.
0009 %            Note that the levels $u$ are not equally spaced. If the spacing 
0010 %            in the grid is not fine  enough, choose smaller value of parameter  ds. 
0011 %      fu  = the pdf structure: f.x contains levels u, f.f is a two collon matrix containing
0012 %            the Danniel's-saddlepoint approximation of the density (pdf) of the quadratic 
0013 %            sea in the first collon and pdf of the Gaussian sea in the second.  
0014 %      Fu  = the two collon matrix containing levels u and the Lugannani Rice's-saddlepoint 
0015 %            approximation of the cdf of the quadratic sea.  
0016 %     S    = spectral density structure, computed using WAFO. S.h contains the water depth.
0017 %            (default: Jonswap with depth 5000 [m]).
0018 %     ds   = parameter defining the levels: default value is ds=0.1.
0019 %
0020 %   Example: S=Jonswap; S.h=40;
0021 %            [mu, fu, Fu]=specq2lc(S,0.05); semilogy(mu(:,1),mu(:,2:3))
0022 %            pdfplot(fu)
0023 %
0024 
0025 
0026 % References: U. Machado, I. Rychlik (2002) "Wave statistics ib nonlinear sea" to
0027 %             appear in Extremes.
0028 %             Butler, R., Machado, U. Rychlik, I. (2002)"Distribution of wave crests in non-
0029 %             linear random sea - application of saddlepoint methods" by , presented at ISOPE 2003. 
0030 % Calls: kwaves.m  which evaluates the cumulant generating function K(s,t) as defined in Eq. 13
0031 % By U.M. 04.11.02
0032 % Revised by I.R 22.11.02
0033 %
0034 %------------------------------------------------------------------------------------
0035 
0036 
0037 
0038 if nargin<1 
0039   S=jonswap;
0040 end
0041 
0042 if nargin<2 | ds0<=0
0043   ds0=0.1;
0044 end
0045 
0046 if nargin<3
0047   nb=0;
0048 end
0049 
0050 h=Spec.h;
0051    if (h>5000)
0052        h=5000;
0053    end
0054 
0055 g=gravity;
0056 omega=Spec.w;
0057 spec=sqrt(Spec.S);
0058 spec=spec*spec';
0059 dw=Spec.w(2,1)-Spec.w(1,1);
0060 w=Spec.w;
0061 
0062 if nb==1 % if narrow-band
0063  [i,j]=max(Spec.S);
0064  w=Spec.w(j)*ones(size(Spec.w));
0065 end 
0066 
0067 %=== Computation of the transfer functions Q, R and S
0068 
0069 W=diag(-Spec.w);
0070 Q=(eww(w,-w,h)+eww(w,w,h)).*spec*dw;
0071 R=(eww(w,-w,h)-eww(w,w,h)).*spec*dw;
0072       
0073 %===
0074 Q=(Q+Q')/2;
0075 R=(R+R')/2;
0076 S=Q*W-W*R;
0077 sigma=sqrt(Spec.S*dw);
0078 N=length(W);
0079 
0080 %===
0081 [P1c,Delta]=eig(Q);
0082 P1=P1c';
0083 % eigenvectors per line
0084 %check: mesh(Q-P1'*Delta*P1)
0085 
0086 [P2c,Gama]=eig(R);
0087 P2=P2c';
0088 %check: mesh(R-P2'*Delta*P2)
0089 
0090 %==============================================
0091 %              Delta P1 Gama P2
0092 %
0093 % ORDER from the small to the largest values
0094 %==============================================
0095 %
0096 % small letters are vectors:
0097 
0098 lambda=diag(real(Delta));
0099 [Vl, Il]=sort(abs(lambda));
0100 lambdaord(:,1)=lambda(Il(:,1),1);
0101 
0102 %Vl: absolute values of eigenvalues in ascending order
0103 %Il: positions where they were
0104 
0105 Deltaord=diag(lambdaord);
0106 
0107 for i=1:length(P1)
0108  P1ord(i,:)=P1(Il(i,1),:);
0109 end
0110 %====
0111 
0112 gamda=diag(real(Gama));
0113 [Vl, Ill]=sort(abs(gamda));
0114 gamdaord(:,1)=gamda(Ill(:,1),1);
0115 
0116 Gamaord=diag(gamdaord);
0117 
0118 for i=1:length(P2)
0119  P2ord(i,:)=P2(Ill(i,1),:);
0120 end
0121 
0122 %=== Computation of P now ordered
0123 P=[P1ord zeros(size(P1)) ; zeros(size(P2)) P2ord];
0124 
0125 %==============================================
0126 % REDUCE i.e. take away lines (replace mm lines by zeros)
0127 %==============================================
0128 
0129 variance=2*(sum(lambdaord.^2) + sum(gamdaord.^2));
0130 varapprox=2*(cumsum(lambdaord.^2) + cumsum(gamdaord.^2))/variance;
0131 mm=max(1,sum(varapprox<0.00001));
0132 
0133 lambdatilde=lambdaord;
0134 lambdatilde(1:mm,1)=zeros(mm,1);
0135 
0136 gamdatilde=gamdaord;
0137 gamdatilde(1:mm,1)=zeros(mm,1);
0138 
0139 Deltatilde=diag(lambdatilde);
0140 Gamatilde=diag(gamdatilde);
0141 
0142 % position of the first element that is not zero
0143 np=mm+1;
0144 
0145 
0146 % the same for both
0147 nn=mm;
0148 npp=np;
0149 
0150 
0151 % ====================================================================================
0152 
0153 % for the computation r
0154 
0155 sigma=sqrt(Spec.S*dw);
0156 N=length(W);
0157 
0158 % Here P1 is ordered and without zeros
0159 
0160 rr1=P1ord*sigma;
0161 rr2=P2ord*W*sigma;
0162 
0163 r=[rr1;rr2];
0164 
0165 x1=zeros(mm,1);
0166 x2=zeros(N-mm,1);
0167 x3=zeros(nn,1);
0168 x4=zeros(N-nn,1);
0169 
0170 x1=r(1:mm,1);
0171 x2=r(mm+1:N,1);
0172 x3=r(N+1:N+nn,1);
0173 x4=r(N+nn+1:2*N,1);
0174 
0175 sum1=sum(x1.^2)/2;
0176 sum2=sum(x3.^2)/2;
0177 
0178 Sn=Deltatilde*P1ord*W*P2ord'-P1ord*W*P2ord'*Gamatilde;
0179 % ====================================================================================
0180 %check: mesh(P1ord'*Sn*P2ord-S)
0181 
0182 A22=Deltatilde(mm+1:N,mm+1:N);
0183 A44=Gamatilde(nn+1:N,nn+1:N);
0184 A24=Sn(mm+1:N,mm+1:N);
0185 %A14=Sn(1:mm,mm+1:end);
0186 
0187 C=zeros(2*mm,2*N-2*mm);
0188 C(1:mm,N-mm+1:end)=Sn(1:mm,mm+1:end);
0189 C(mm+1:end,1:N-mm)=Sn(mm+1:end,1:mm)';
0190 
0191 CC=C'*C;
0192 
0193 
0194 %%%%%%%%%%%%%%%%%%%%%%% CHECK of  accuracy of K %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0195 if nb==99
0196   s1=-3;
0197   s2=2;
0198   t=[s1*sigma ; s2*W*sigma];
0199   A=[s1*Q s2*S;  s2*S' s1*R];
0200   IA=eye(size(A));
0201   InvA=inv(IA-A);
0202   IQ=inv(eye(size(Q))-s1*Q);
0203   %Note that we have  minus -S' giving the plus in the following formula.
0204   r=s2*W*sigma+s2*S'*IQ*sigma*s1;
0205   % Here I am checking Lemma 6.
0206   K0=0.5*s1^2*sigma'*IQ*sigma +0.5*r'*inv(eye(size(R))-s1*R-s2*S'*IQ*S*s2)*r-0.5*log(det(eye(size(Q))-s1*Q))-0.5*log(det(eye(size(R))-s1*R-s2*S'*IQ*S*s2));
0207     lamb = spec2mom(Spec);
0208     [i,j]=max(Spec.S);
0209     wp=Spec.w(j);
0210     Const=wp*wp/(2*g);
0211 [ K0 0.5*t'*InvA*t-0.5*log(det(IA-A)) kwaves(s1,s2,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4)   knb(s1,s2,lamb(1),lamb(2),lamb(3),Const)]
0212 end
0213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0214 
0215 
0216 
0217 % ====================================================================================
0218 
0219 
0220 K=[];
0221 DK=[];
0222 DDK=[];
0223 DK111=[];
0224 
0225 DK1122=[];
0226 DK122=[];
0227 
0228 DDKt=[];
0229 
0230 DK2222=[];
0231 
0232 S=[];
0233 
0234 s=0;
0235 j=0;
0236 flag=0;
0237 
0238 ds=0.01;
0239 dt=ds;
0240 
0241 % In these formulas dt must be equal to ds
0242 
0243 %==>positive side
0244 while flag~=1
0245  
0246   j=j+1; 
0247   
0248      
0249  [K0]=kwaves(s,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0250  [K1]=kwaves(s-ds,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0251  [K2]=kwaves(s+ds,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);   
0252  [K3]=kwaves(s-2*ds,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0253  [K4]=kwaves(s+2*ds,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0254  
0255  [K1t]=kwaves(s,-dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0256  [K2t]=kwaves(s,dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0257  
0258  [K11]=kwaves(s-ds,-dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0259  [K22]=kwaves(s+ds,dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0260  
0261  [K12t]=kwaves(s-ds,dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0262  [K21t]=kwaves(s+ds,-dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0263  
0264  [K2tt]=kwaves(s,2*dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0265  [K1tt]=kwaves(s,-2*dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0266 
0267  dK=(K2-K1)/(2*ds);
0268  ddK=(K1+K2-2*K0)/(ds*ds);
0269  
0270  dddK=(K4-2*K2+2*K1-K3)/(2*ds*ds*ds);
0271  d2sd2t=(K22-2*K2t+K12t-2*K2+4*K0-2*K1+K21t-2*K1t+K11)/(ds*ds*ds*ds); 
0272  dsd2t=(K22-2*K2+K21t-K12t+2*K1-K11)/(2*ds*ds*ds);
0273  
0274  ddKt=(K1t+K2t-2*K0)/(dt*dt);
0275  d4t=(K2tt-4*K2t+6*K0-4*K1t+K1tt)/(dt*dt*dt*dt);
0276  
0277  if (j~=1) & ((s*dK-K0)>25 | 0>(s*dK-K0) )
0278 %if (j~=1) & (abs(K0-s*dK)>20)
0279 %    if (j~=1) & (abs(K0-s*dK)>10) 
0280   flag=1;
0281 else
0282 
0283  S=[S; s]; 
0284  s=s+ds0/2; 
0285  
0286  K=[K; K0];  
0287  DK=[DK; dK];
0288  DDK=[DDK; ddK];
0289  
0290  DK111=[DK111; dddK];
0291  DDKt=[DDKt; ddKt];
0292  DK1122=[DK1122; d2sd2t]; 
0293  DK122=[DK122; dsd2t]; 
0294  DK2222=[DK2222; d4t];
0295  end    
0296 end     
0297 
0298 %=======================================================
0299 
0300 s=-ds0/2;
0301 j=0;
0302 flag=0;
0303 
0304 %==>negative side
0305 
0306 while flag~=1
0307  j=j+1;
0308  [K0]=kwaves(s,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0309  [K1]=kwaves(s-ds,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0310  [K2]=kwaves(s+ds,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0311  [K3]=kwaves(s-2*ds,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0312  [K4]=kwaves(s+2*ds,0,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0313   
0314  [K1t]=kwaves(s,-dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0315  [K2t]=kwaves(s,dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4); 
0316  
0317  [K11]=kwaves(s-ds,-dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0318  [K22]=kwaves(s+ds,dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0319 
0320  [K12t]=kwaves(s-ds,dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0321  [K21t]=kwaves(s+ds,-dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0322  
0323  [K2tt]=kwaves(s,2*dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0324  [K1tt]=kwaves(s,-2*dt,A22,A44,A24,C',CC,sum1,sum2,x1,x2,x3,x4);
0325  
0326  dK=(K2-K1)/(2*ds);
0327  ddK=(K1+K2-2*K0)/(ds*ds);
0328  
0329  dddK=(K4-2*K2+2*K1-K3)/(2*ds*ds*ds); 
0330  d2sd2t=(K22-2*K2t+K12t-2*K2+4*K0-2*K1+K21t-2*K1t+K11)/(ds*ds*ds*ds);
0331  dsd2t=(K22-2*K2+K21t-K12t+2*K1-K11)/(2*ds*ds*ds);
0332  
0333  ddKt=(K1t+K2t-2*K0)/(dt*dt);
0334  d4t=(K2tt-4*K2t+6*K0-4*K1t+K1tt)/(dt*dt*dt*dt);
0335  
0336  if (j~=1) & ((s*dK-K0)>25 | 0>(s*dK-K0) )
0337 %         if (j~=1) & (abs(K0-s*dK)>10)
0338     flag=1;
0339 else
0340  
0341  S=[s;S];
0342  s=s-ds0/2;
0343  K=[K0;K];
0344  DK=[dK;DK];
0345  DDK=[ddK;DDK];
0346  
0347  DK111=[dddK;DK111];
0348  DDKt=[ddKt; DDKt];
0349  DK1122=[d2sd2t; DK1122]; 
0350  DK122=[dsd2t; DK122]; 
0351  DK2222=[d4t; DK2222];
0352  end
0353 end     
0354 
0355 % ====================================================================================
0356 % ====================================================================================
0357 % marginal density
0358 
0359 f=1/sqrt(2*pi)*exp(K-S.*DK)./sqrt(DDK);
0360 %(S.*DK-K)
0361 wh=sign(S).*sqrt(2*(S.*DK-K));
0362 %f1=wnormpdf(wh)./sqrt(DDK);
0363 L0=spec2mom(Spec,2);
0364 f1=1/sqrt(2*pi)*exp(-DK.^2/2/L0(1))./sqrt(L0(1));
0365 
0366 [area]=trapz(DK,f);
0367 fu=[DK f/area]; 
0368 fu=createpdf;
0369   Htxt = ['Density of X(0) for guadratic sea'];
0370   xtxt = ['X(0) [m]'];
0371 fu.title=Htxt;
0372 fu.labx{1}=xtxt;
0373 fu.x{1}=DK;
0374 fu.f=[f/area f1];   
0375 
0376 %fu1=[DK f1/area];
0377 F1=wnormcdf(wh)-wnormpdf(wh).*((1./S)./sqrt(DDK)-1./wh);
0378 Fu=[DK F1];
0379 
0380 % ====================================================================================
0381 % 1-TERM SADDLEPONT APRROXIMATION (GAUSSIAN APPR)
0382 % ====================================================================================
0383 
0384 muG=f.*sqrt(DDKt)/sqrt(2*pi)/area;
0385 %figure(3)
0386 %plot(DK,muG,'k:o')
0387 
0388 % ====================================================================================
0389 % 2-TERMS SADDLEPONT APRROXIMATION (GAUSSIAN APPR)
0390 % ====================================================================================
0391 
0392 % equation 22
0393 eq22=-1/4*(DK1122.*DDK-DK111.*DK122)./((DDK.^2).*DDKt);
0394 
0395 %mu_2=muG.*(1+eq22);
0396 
0397 %figure(3)
0398 %hold on
0399 %plot(DK,mu_2,'k:*')
0400 
0401 % ====================================================================================
0402 % 3-TERMS SADDLEPONT APRROXIMATION (GAUSSIAN APPR)
0403 % ====================================================================================
0404 
0405 % equation 23
0406 eq23=(DK2222.*DDK-3*DK122.^2)./(DDK.*(DDKt.^2));
0407 
0408 mu1=muG.*(1+eq22);
0409 mu=muG.*(1+eq22-(1/24)*eq23);
0410 mu2=muG.*(1+eq22-(1/24)*eq23.*(1-eq22));
0411 %figure(3)
0412 %hold on
0413 %plot(DK,mu,'k:d')
0414 
0415 %mu1=DDKt;
0416 %mu2=(DK2222.*DDK-3*DK122.^2)./DDK;
0417 
0418 mu=[DK mu 1/2/pi*exp(-DK.^2/2/L0(1))*sqrt(L0(2)/L0(1))];
0419 muG=[DK muG];
0420 mu1=[DK mu1];
0421 mu2=[DK mu2];
0422 
0423 
0424 
0425 return % main
0426 
0427 % ====================================================================================
0428 
0429 function [out]=eww(omega,omegat,h);
0430 % EWW: Computes values of the quadratic transfer function E, given in Eq. 6. in
0431 %             R. Butler, U. Machado, I. Rychlik (2002). 
0432 %
0433 %  CALL:  [Eww]= eww(w,wt,h);
0434 %  
0435 %     Eww  = a matrix with the quadratic transfer function E(w,wt). 
0436 %     w,wt = two equally long vectors with angular frequencies.
0437 %     h    = water depth (default 5000 [m]).
0438 %
0439 %   Example: S=Jonswap; w=[-flipud(S.w) ;S.w]; Eww=eww(w,w); mesh(w,w,Eww)
0440 
0441 %  Bugs:   E(w,-w)=0. (This should be checked?)
0442 %
0443 %----------------------------------------------------------------------
0444 % References: R. Butler, U. Machado, I. Rychlik (2002) Distribution of waves crests in nonlinear 
0445 %             random sea - applications of saddlepoint methods. ISOPE 2003. 
0446 % Calls: kwaves.m
0447 % By U.M. 10.11.02
0448 % Revised by I.R 22.11.02
0449 
0450 g=gravity;
0451 eps0=0.000001;
0452 
0453 if (length(omega)~=length(omegat))
0454    error('error in input to eww')    
0455 end
0456 
0457 if nargin<3 | h<=0
0458   h=5000;
0459 end
0460 
0461 [w wt]=meshgrid(omega,omegat);
0462 wpl=w+wt;
0463 
0464 
0465 ind=find(abs(w.*wt.*wpl)<eps0);
0466 wpl(ind)=1;
0467 w(ind)=1;
0468 wt(ind)=1;
0469      
0470 kw=w2k(w,[],h,g);
0471 kwt=w2k(wt,[],h,g);
0472 
0473 Etop=g*(kw.*kwt)./(w.*wt)-(w.^2+wt.^2+w.*wt)/(2*g)+...
0474             (g/2)*(w.*kwt.^2+wt.*kw.^2)./(w.*wt.*wpl);
0475 
0476 Ebottom=1-g*(kw+kwt)./wpl.^2.*tanh((kw+kwt)*h);
0477 
0478 out=Etop./Ebottom-(g*kw.*kwt)./(2*w.*wt)+(w.^2+wt.^2+w.*wt)/(2*g);
0479 
0480 out(ind)=0;
0481 return %eww
0482 
0483 %-----------------------------------------------------------------------
0484 
0485 function [K]=kwaves(s1,s2,A22,A44,A24,CT,CC,sum1,sum2,x1,x2,x3,x4);
0486 % KWAVES: help function which computes the cumulant generating function 
0487 %         formula (36) in Machado and Rychlik (2002) 
0488 %         "Wave statistics in nonlinear random sea" to appear in Extremes.
0489 %_____________________________________________________________________________ 
0490 %   Note that s1,s2 must be constants (vectors are not allowed):
0491 
0492 
0493 % tested on matlab 6.1
0494 % By U Machado, last update 20.10.02
0495 % revised IR 19.11.02, changed name and added some checks and comments.
0496 
0497 
0498 
0499 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0500 % We tried to keep notation of the Appendix in Machado&Rychlik
0501 % CT means C',  CC=C'*C
0502 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0503 % 
0504 % r=r(s_1,s_2) and B=B(s_1,s_2)   are defined in Eq. 37.
0505 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0506 
0507 r=[s1*x2; s2*x4]+s2*CT*[s1*x1; s2*x3];
0508 
0509 %r=[s1*x2; s2*x4]-s2*CT*[s1*x1; s2*x3];
0510 
0511 %===
0512 
0513 Aux=s2^2*CC;
0514 
0515 I=eye(size(Aux));
0516 
0517 B=I-([s1*A22 s2*A24; s2*A24' s1*A44])-Aux;
0518 
0519 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0520 % K=K(s_1,s_2) This is formulas (36) in Machado&Rychlik          
0521 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0522 
0523 K=-1/2*log(det(B))+s1^2*sum1+s2^2*sum2+1/2*r'*inv(B)*r;
0524 
0525 %[B rn]
0526 %[K s1^2*sum1+s2^2*sum2 1/2*rn'*inv(B)*rn]
0527 return % kwaves
0528 
0529 
0530 
0531 
0532 
0533 
0534 
0535 
0536 
0537

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