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)
PDF class constructor | |
returns the constant acceleration of gravity | |
Calculates (and plots) a JONSWAP spectral density | |
Calculates spectral moments from spectrum | |
Translates from frequency to wave number | |
Normal cumulative distribution function | |
Normal probability density function | |
Determinant. | |
Eigenvalues and eigenvectors. | |
Display message and abort function. | |
Matrix inverse. | |
X and Y arrays for 3-D plots. | |
Trapezoidal numerical integration. |
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
Comments or corrections to the WAFO group