Home > wafo > wavemodels > tay90pdf.m

# tay90pdf

## PURPOSE

Tayfun (1990) PDF of large wave heights

## SYNOPSIS

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

## DESCRIPTION

``` TAY90PDF Tayfun (1990) PDF of large wave heights
/1
f = .5*H^3/A^2 *  |  z/sqrt(1-z)* exp(-(H/A)^2*(2-z))*Io(  (H/A)^2*z*b/2) dz
/0

where Io(.) is the modified bessel funcition of zero order

CALL:  [f tol]= tay90pdf(H,A,B,reltol)

f = pdf
H = wave height
A = scale parameter (=sqrt(m0))
B = correlation between trough and crest amplitude squared
i.e. Cov(At^2, Ac^2) approx -R(T/2)  (B->1 => pdf->Rayleigh )
reltol = relative tolerance default=1e-3
tol = absolute tolerance abs(int-intold)

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.

Tayfun calculates the correlation from the Spectral density as:

B = sqrt( [int S(w)*cos(w*T) dw]^2+[int S(w)*sin(w*T) dw]^2 )/m0

where T = Tm02, mean zero crossing period.
Note that this is the same as the groupiness factor, Ka, found in
spec2char

Example:
S = jonswap;
a = sqrt(spec2mom(S,1));
b = spec2char(S,'Ka');
h = linspace(0,7*a)';
plot(h,tay90pdf(h,a,b))

## CROSS-REFERENCE INFORMATION

This function calls:
 gaussq Numerically evaluates a integral using a Gauss quadrature. tay90fun is an internal integrand function to tay90pdf wraylpdf Rayleigh probability density 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:
 tay90cdf Tayfun (1990) CDF of large wave heights

## SOURCE CODE

```001 function [y, tol1] = tay90pdf(x,a,b,tol)
002 %TAY90PDF Tayfun (1990) PDF of large wave heights
003 %                      /1
004 %   f = .5*H^3/A^2 *  |  z/sqrt(1-z)* exp(-(H/A)^2*(2-z))*Io(  (H/A)^2*z*b/2) dz
005 %                     /0
006 %
007 %   where Io(.) is the modified bessel funcition of zero order
008 %
009 %  CALL:  [f tol]= tay90pdf(H,A,B,reltol)
010 %
011 %            f = pdf
012 %            H = wave height
013 %            A = scale parameter (=sqrt(m0))
014 %            B = correlation between trough and crest amplitude squared
015 %                i.e. Cov(At^2, Ac^2) approx -R(T/2)  (B->1 => pdf->Rayleigh )
016 %       reltol = relative tolerance default=1e-3
017 %          tol = absolute tolerance abs(int-intold)
018 %
019 % The size of f is the common size of X, A and B.  A scalar input
020 % functions as a constant matrix of the same size as the other input.
021 %
022 %  Tayfun calculates the correlation from the Spectral density as:
023 %
024 %  B = sqrt( [int S(w)*cos(w*T) dw]^2+[int S(w)*sin(w*T) dw]^2 )/m0
025 %
026 %  where T = Tm02, mean zero crossing period.
027 %  Note that this is the same as the groupiness factor, Ka, found in
028 %  spec2char
029 %
030 % Example:
031 %  S = jonswap;
032 %  a = sqrt(spec2mom(S,1));
033 %  b = spec2char(S,'Ka');
034 %  h = linspace(0,7*a)';
035 %  plot(h,tay90pdf(h,a,b))
036 %
038
039
040
041 %   Reference:
042 % [1]  M. Aziz Tayfun (1990) "Distribution of large waveheights"
043 % Journal of the Waterway, port and coastal and ocean division
044 %      vol 116 No.6 pp. 686-707
045 % [2]  M. Aziz Tayfun (1981) "Distribution of crest-to-trough waveheights"
046 % Journal of the Waterway, port and coastal and ocean division
047 %      vol 107 No.3 pp. 149-158
048
049 %tested on:
050 %history:
051 % revised pab 26.07.2001
052 % - added forgotten ) in error(nargchk(3,4,nargin))
053 % revised IR, JR, PAB
054 %  raylpdf -> wraylpdf
055 % revised pab 07.03.2000
056 % - updated error in help header
057 % revised pab 28.02.2000
058 %  - fixed some bugs
059 %  Per A. Brodtkorb 21.02.99
060
061
062 error(nargchk(3,4,nargin))
063 if (nargin <4)| isempty(tol),
064   tol=1e-3; %relative accuracy of the estimates % sqrt(eps);%
065 end
066
067 [errorcode x a b] = comnsize(x,a,b);
068
069 if errorcode > 0
070   error('h, a and b must be of common size or scalar.');
071 end
072 xsz = size(x);
073 x=x(:); a=a(:);b=b(:);
074 %a
075 %b
076 % Initialize Y to zero.
077 y=zeros(size(x));
078 tol1=y;
079 % Return NaN if A is not positive.
080 k1 = find(a <= 0| abs(b)>1);
081 if any(k1)
082     tmp   = NaN;
083     y(k1) = tmp(ones(size(k1)));
084 end
085
086 if 0
087   % This is a trick to get the html documentation correct.
088   k = tay90fun(1,1);
089 end
090
091 gn=2;% # of of base points to start the integration with
092 trace=0;%used for debugging
093 %size(x)
094 %size(a)
095 %size(y)
096 k=find(a > 0 & x >0 & abs(b)<1);
097 %size(k)
098 if any(k),
099   yk=y(k);xk = x(k); ak = a(k);
100
101  [yk, tol1(k)]= gaussq('tay90fun',0,ones(size(yk)),  [tol 4],[trace gn],xk./ak,b(k));
102   yk=yk./ak;
103
104   k2=find(yk<0 );
105    if any(k2)
106      disp('Some less than zero')
107      yk(k2)= zeros(size(k2));
108   end
109  y(k)=yk;
110 end
111
112 k4=find(a > 0 & x >= 0 & abs(b)==1);
113 if any(k4), % special case b==1 -> rayleigh distribution
114  y(k4) = wraylpdf(x(k4),2*a(k4));
115 end
116 y=reshape(y,xsz);
117
118
119
120
121
122
123```

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