TAY81CDF Tayfun (1981) CDF of breaking limited wave heights /inf pdf = x/a^2 * | u* J_o(u/b)^(b^2)*J_o(xu) du *(1 - 4/pi*acos(b/x)*H(1-b/x) ) /0 where J_o and H(.) is the zero order Bessel function and the heavyside function CALL: [F tol]= tay81cdf(x,A,B,reltol) F = CDF x = wave heights 0 <= x <= sqrt(2)*B A = scale parameter B = maximum limiting value (B->inf => pdf->Rayleigh ) reltol = relative tolerance (default 1e-3) tol = absolute tolerance abs(int-intold) Tayfun set the scale parameter and breaking limiting value to A = Hrms = 2*sqrt(2*m0)*sqrt(1-0.734*eps2^4) B = pi*tanh(k0*h)/(7*k0*2*sqrt(m0)) where k0 = the mean apparent wave number mi = i'th spectral moment (i=0 => variance of the process) eps2 = sqrt(m0*m2/m1^2-1) spectral bandwidth h = water depth The size of F is the common size of X, A and B. A scalar input functions as a constant matrix of the same size as the other input. Example: Hm0 = 7; Tp = 10; h = 150; S = jonswap([],[Hm0,Tp]); eps2 = spec2char(S,'eps2'); Tm02 = spec2char(S,'Tm02'); a = Hm0/sqrt(2)*sqrt(1-0.734*eps2^4); k0 = w2k(2*pi/Tm02,[],h); b = pi*tanh(k0*h)/(7*k0*Hm0/2); x = linspace(0,Hm0)'; plot(1:5,tay81cdf(1:5,a,b),x,wraylcdf(x,Hm0/2),'r') See also tay81pdf, braylcdf
Numerically evaluates a integral using a Gauss quadrature. | |
Tayfun (1981) PDF of breaking limited wave heights | |
Rayleigh cumulative distribution function | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Not-a-Number. |
001 function [y, tol1] = tay81cdf(x,a,b,tol) 002 %TAY81CDF Tayfun (1981) CDF of breaking limited wave heights 003 % /inf 004 % pdf = x/a^2 * | u* J_o(u/b)^(b^2)*J_o(xu) du *(1 - 4/pi*acos(b/x)*H(1-b/x) ) 005 % /0 006 % where J_o and H(.) is the zero order Bessel function and the heavyside function 007 % 008 % CALL: [F tol]= tay81cdf(x,A,B,reltol) 009 % 010 % F = CDF 011 % x = wave heights 0 <= x <= sqrt(2)*B 012 % A = scale parameter 013 % B = maximum limiting value (B->inf => pdf->Rayleigh ) 014 % reltol = relative tolerance (default 1e-3) 015 % tol = absolute tolerance abs(int-intold) 016 % 017 % Tayfun set the scale parameter and breaking limiting value to 018 % A = Hrms = 2*sqrt(2*m0)*sqrt(1-0.734*eps2^4) 019 % B = pi*tanh(k0*h)/(7*k0*2*sqrt(m0)) 020 % where 021 % k0 = the mean apparent wave number 022 % mi = i'th spectral moment (i=0 => variance of the process) 023 % eps2 = sqrt(m0*m2/m1^2-1) spectral bandwidth 024 % h = water depth 025 % 026 % 027 %The size of F is the common size of X, A and B. A scalar input 028 % functions as a constant matrix of the same size as the other input. 029 % 030 % Example: 031 % Hm0 = 7; Tp = 10; h = 150; 032 % S = jonswap([],[Hm0,Tp]); 033 % eps2 = spec2char(S,'eps2'); 034 % Tm02 = spec2char(S,'Tm02'); 035 % a = Hm0/sqrt(2)*sqrt(1-0.734*eps2^4); 036 % k0 = w2k(2*pi/Tm02,[],h); 037 % b = pi*tanh(k0*h)/(7*k0*Hm0/2); 038 % x = linspace(0,Hm0)'; 039 % plot(1:5,tay81cdf(1:5,a,b),x,wraylcdf(x,Hm0/2),'r') 040 % 041 % See also tay81pdf, braylcdf 042 043 % Reference: 044 % M. Aziz Tayfun (1981) "Breaking limited waveheights" 045 % Journal of the Waterway, port and coastal and ocean division 046 % vol 107 No.2 pp. 59-69 047 048 %tested on: 049 %history 050 % 051 % Per A. Brodtkorb 01.02.99 052 053 % TODO % not tested 054 055 if (nargin <4)| isempty(tol), 056 tol=1e-3;%sqrt(eps);%relative accuracy of the estimates 057 end 058 if nargin < 3, 059 error('Requires at least three input arguments.'); 060 end 061 062 [errorcode x a b] = comnsize(x,a,b.^2); 063 064 if errorcode > 0 065 error('Requires non-scalar arguments to match in size.'); 066 end 067 068 % Initialize Y to zero. 069 y=zeros(size(x)); 070 tol1=y; 071 % Return NaN if A or B is not positive. 072 k1 = find(a <= 0|b<= 0); 073 if any(k1) 074 tmp = NaN; 075 y(k1) = tmp(ones(size(k1))); 076 end 077 078 trace=logical(0);%[];%used for debugging 079 080 k=find(a > 0 & x > 0 &x<=a.*sqrt(2*b) & b<inf&b> 0); 081 if any(k), 082 yk = y(k); xk = x(k); ak = a(k); bk = b(k); 083 yk = yk(:); xk = xk(:); ak = ak(:); bk = bk(:); 084 if 0,all((ak==ak(1)).*(bk==bk(1))), 085 [xk ind]=sort(xk); 086 gn=2; 087 [y(k), tol(k)]= gaussq('tay81pdf',[0;xk(1:end-1)] ,xk,tol,trace,gn,ak,bk,tol); 088 %[y(k), tol(k)]= asimpson('tay81pdf',0,xk,tol,trace,ak,bk); 089 y(k)=cumsum(y(k)); 090 y(k)=y(k(ind)); %make sure the ordering is right 091 else 092 gn=2; 093 [y(k), tol(k)]= gaussq('tay81pdf',0,xk,tol,trace,gn,ak,bk,tol); 094 %[y(k), tol(k)]= asimpson('traylpdf',0,xk,tol,trace,ak,bk); 095 end 096 end 097 if 0 098 % This is a trick to get the html documentation correct. 099 k = tay81pdf(1,1); 100 end 101 102 k2=find((a > 0 & y>1) | (x>=a.*sqrt(2*b) )& b<inf & b> 0); 103 if any(k2) 104 y(k2)= ones(size(k2)); 105 end 106 107 k4=find(a > 0 & x >= 0 & b==inf); 108 if any(k4), % special case b==inf -> rayleigh distribution 109 y(k4) = wraylcdf(x(k4),a(k4)/sqrt(2)); 110 end 111 112 function z=traylfunzeros(x,a,b,n), 113 % internal function which calculates the approximate zerocrossings of 114 % the integrand: x/a* J_o(u/sqrt(b))^b*J_1(x/a* u) 115 nx=1:n; 116 z1=a./x; 117 z0=sqrt(b); 118 %k=find(mod(b,2)==0); % remove 119 %z0(k)=pi/2*(n+1)*max(z1); 120 z0=z0(:,ones(1,n)).* pi.*(nx(ones(size(x)),:)+0.25); 121 z1=z1(:,ones(1,n)).* pi.*(nx(ones(size(x)),:)-0.25); 122 z=sort([z0,z1],2); 123 return 124 125 126
Comments or corrections to the WAFO group