Home > wafo > trgauss > specq2lc.m

# specq2lc

## 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,
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
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:
 createpdf PDF class constructor gravity returns the constant acceleration of gravity jonswap Calculates (and plots) a JONSWAP spectral density spec2mom Calculates spectral moments from spectrum w2k Translates from frequency to wave number wnormcdf Normal cumulative distribution function wnormpdf Normal probability density function det Determinant. eig Eigenvalues and eigenvectors. error Display message and abort function. inv Matrix inverse. meshgrid X and Y arrays for 3-D plots. trapz Trapezoidal numerical integration.
This function is called by:

## 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,
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