Home > wafo > wavemodels > tay81pdf.m

tay81pdf

PURPOSE

Tayfun (1981) PDF of breaking limited wave heights

SYNOPSIS

[y, tol1] = tay81pdf(x,a,b,tol)

DESCRIPTION

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')

CROSS-REFERENCE INFORMATION

This function calls:
 gaussq Numerically evaluates a integral using a Gauss quadrature. tay81fun internal function to tay81pdf and tay81cdf. wraylpdf Rayleigh probability density function comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. figure Create figure window. nan Not-a-Number.
This function is called by:
 tay81cdf Tayfun (1981) CDF of breaking limited wave heights

SOURCE CODE

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 %
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

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