TAY81PDF Tayfun (1981) PDF of breaking limited wave heights /inf f = x/a^2 * | u* J_o(u/b)^(b^2)*J_o(x/a* u) du *(1 - 4/pi*acos(b*a/x)*H(x/a-b) ) /0 where J_o and H(.) is the zero order Bessel function and the heavyside function CALL: [f tol]= tay81pdf(H,A,B,reltol) f = pdf x = quantiles 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(x,tay81pdf(x,a,b),x,wraylpdf(x,Hm0/2),'r') See also tay81cdf, braylpdf
Numerically evaluates a integral using a Gauss quadrature. | |
internal function to tay81pdf and tay81cdf. | |
Rayleigh probability density function | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Create figure window. | |
Not-a-Number. |
Tayfun (1981) CDF of breaking limited wave heights |
001 function [y, tol1] = tay81pdf(x,a,b,tol) 002 %TAY81PDF Tayfun (1981) PDF of breaking limited wave heights 003 % /inf 004 % f = x/a^2 * | u* J_o(u/b)^(b^2)*J_o(x/a* u) du *(1 - 4/pi*acos(b*a/x)*H(x/a-b) ) 005 % /0 006 % where J_o and H(.) is the zero order Bessel function and the heavyside function 007 % 008 % CALL: [f tol]= tay81pdf(H,A,B,reltol) 009 % f = pdf 010 % x = quantiles 0 <= x <= sqrt(2)*B 011 % A = scale parameter 012 % B = maximum limiting value (B->inf => pdf->Rayleigh ) 013 % reltol = relative tolerance (default 1e-3) 014 % tol = absolute tolerance abs(int-intold) 015 % 016 % Tayfun set the scale parameter and breaking limiting value to 017 % A = Hrms = 2*sqrt(2*m0)*sqrt(1-0.734*eps2^4) 018 % B = pi*tanh(k0*h)/(7*k0*2*sqrt(m0)) 019 % where 020 % k0 = the mean apparent wave number 021 % mi = i'th spectral moment (i=0 => variance of the process) 022 % eps2 = sqrt(m0*m2/m1^2-1) spectral bandwidth 023 % h = water depth 024 % 025 % 026 % The size of F is the common size of X, A and B. A scalar input 027 % functions as a constant matrix of the same size as the other input. 028 % 029 % Example: 030 % Hm0 = 7; Tp = 10; h = 150; 031 % S = jonswap([],[Hm0,Tp]); 032 % eps2 = spec2char(S,'eps2'); 033 % Tm02 = spec2char(S,'Tm02'); 034 % a = Hm0/sqrt(2)*sqrt(1-0.734*eps2^4); 035 % k0 = w2k(2*pi/Tm02,[],h); 036 % b = pi*tanh(k0*h)/(7*k0*Hm0/2); 037 % x = linspace(0,Hm0)'; 038 % plot(x,tay81pdf(x,a,b),x,wraylpdf(x,Hm0/2),'r') 039 % 040 % See also tay81cdf, braylpdf 041 042 % Reference: 043 % M. Aziz Tayfun (1981) "Breaking limited waveheights" 044 % Journal of the Waterway, port and coastal and ocean division 045 % vol 107 No.2 pp. 59-69 046 047 %tested on: 048 %history: 049 % Per A. Brodtkorb 01.02.99 050 051 % TODO % not tested 052 if (nargin <4)| isempty(tol), 053 tol=1e-3; %relative accuracy of the estimates % sqrt(eps);% 054 end 055 if nargin < 3, 056 error('Requires at least three input arguments.'); 057 end 058 059 [errorcode x a b] = comnsize(x,a,b.^2); 060 061 if errorcode > 0 062 error('Requires non-scalar arguments to match in size.'); 063 end 064 [N M]=size(x); 065 x=x(:); a=a(:);b=b(:); 066 067 % Initialize Y to zero. 068 y=zeros(size(x)); 069 tol1=y; 070 % Return NaN if A or B is not positive. 071 k1 = find(a <= 0|b<= 0); 072 if any(k1) 073 tmp = NaN; 074 y(k1) = tmp(ones(size(k1))); 075 end 076 077 078 count_limit=20; 079 gn=2^3;% # of of base points to start the integration with 080 trace=logical(0);%used for debugging 081 082 k=find(a > 0 & x >0 &x<a.*sqrt(2*b) & b<inf&b> 0); 083 if any(k), 084 yk=y(k);xk = x(k); ak = a(k); bk = b(k); 085 tmpy=yk; tmpt=tmpy; %size(xk) 086 uk_inf = traylfunzeros(xk,ak,bk,count_limit*10);% 087 uk_inf = [zeros(size(k)) uk_inf(:,3:4:end)]; 088 089 if 0 090 uk_inf= uk_inf(:,3:4:end); 091 k1=find( uk_inf(:,1)<=7); 092 if any(k1), 093 tmpy=uk_inf(:,1); 094 uk_inf(k1,:)=uk_inf(k1,:)+7-tmpy(k1,ones(1,size(uk_inf,2))); 095 end 096 %uk_inf=.5*ak.*sqrt(2.*bk); 097 %uk_inf(1:10,1:4), size(uk_inf),size(k) 098 [yk, tol1(k)]= gaussq('tay81fun',0,uk_inf(k), tol/2,trace,gn,xk./ak,bk); 099 end 100 if 0 101 % This is a trick to get the html documentation correct. 102 k = tay81fun(1,1,2); 103 end 104 105 106 % Break out of the iteration loop for three reasons: 107 % 1) the last update is very small (compared to int) 108 % 2) the last update is very small (compared to tol) 109 % 3) There are more than 20 iterations. This should NEVER happen. 110 k1=find(k);ix=1; 111 while(any(k1) & ix <= count_limit), 112 if trace,figure(ix),end % for debugging 113 [tmpy(k1), tmpt(k1)]= gaussq('tay81fun',uk_inf(k1,ix),uk_inf(k1,ix+1),... 114 tol/2/count_limit,trace,gn,xk(k1)./ak(k1),bk(k1)); 115 116 yk(k1)=yk(k1)+tmpy(k1); 117 tol1(k(k1))=tol1(k(k1))+tmpt(k1); 118 %find integrals which did not converge 119 k1=find(tol1(k) < abs(tmpy) & tol/count_limit <2*abs(tmpy)); 120 121 %k=find(tol1 > abs(tol*int)| tol1 > abs(tol)); 122 %indices to integrals which did not converge 123 ix = ix + 1; 124 end 125 %yk=yk.*xk./(ak.^2); 126 127 k2=find(yk<0 ); 128 if any(k2) 129 yk(k2)= zeros(size(k2)); 130 end 131 132 k3=find(xk>ak.*sqrt(bk) ); 133 if any(k3) 134 yk(k3)= yk(k3).*(1-4/pi.*acos(ak(k3).*sqrt(bk(k3))./xk(k3))); 135 end 136 y(k)=yk.*xk./ak.^2; 137 end 138 139 k4=find(a > 0 & x >= 0 & b==inf); 140 if any(k4), % special case b==inf -> rayleigh distribution 141 y(k4) = wraylpdf(x(k4),a(k4)/sqrt(2)); 142 end 143 y=reshape(y,N,M); 144 145 146 function z=traylfunzeros(x,a,b,n), 147 % internal function which calculates the approximate zerocrossings of 148 % the integrand: u* J_o(u/sqrt(b))^b*J_o(x/a* u) 149 nx=1:n; 150 z1=a./x; 151 z0=sqrt(b); 152 k=find(mod(b,2)==0); % remove 153 z0(k)=pi/2*(n+1)*max(z1);%size(x) 154 z0=z0(:,ones(1,n)).* pi.*(nx(ones(size(x)),:)-0.25); 155 z1=z1(:,ones(1,n)).* pi.*(nx(ones(size(x)),:)-0.25); 156 z=sort([z0,z1],2); 157 return 158
Comments or corrections to the WAFO group