Home > wafo > wavemodels > tay81cdf.m

# tay81cdf

## PURPOSE

Tayfun (1981) CDF of breaking limited wave heights

## SYNOPSIS

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

## DESCRIPTION

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

## CROSS-REFERENCE INFORMATION

This function calls:
 gaussq Numerically evaluates a integral using a Gauss quadrature. tay81pdf Tayfun (1981) PDF of breaking limited wave heights wraylcdf Rayleigh cumulative distribution function comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. nan Not-a-Number.
This function is called by:

## SOURCE CODE

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

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