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)) See also tay90cdf, spec2char, gaussq
Numerically evaluates a integral using a Gauss quadrature. | |
is an internal integrand function to tay90pdf | |
Rayleigh probability density function | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Not-a-Number. |
Tayfun (1990) CDF of large wave heights |
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 % 037 % See also tay90cdf, spec2char, gaussq 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
Comments or corrections to the WAFO group