TAY90CDF Tayfun (1990) CDF of large wave heights /H /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 dH /0 /0 where Io(.) is the modified bessel funcition of zero order CALL: [F, tol]= tay90cdf(H,A,B,reltol); F = CDF H = wave height A = scale parameter (=sqrt(m0)) B = correlation between At and Ac (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,tay90cdf(h,a,b)) See also tay90pdf, spec2char, gaussq
Numerically evaluates a integral using a Gauss quadrature. | |
Tayfun (1990) PDF of large wave heights | |
Rayleigh cumulative distribution function | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Not-a-Number. |
001 function [y, tol1] = tay90cdf(x,a,b,tol) 002 %TAY90CDF Tayfun (1990) CDF of large wave heights 003 % /H /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 dH 005 % /0 /0 006 % 007 % where Io(.) is the modified bessel funcition of zero order 008 % 009 % CALL: [F, tol]= tay90cdf(H,A,B,reltol); 010 % 011 % F = CDF 012 % H = wave height 013 % A = scale parameter (=sqrt(m0)) 014 % B = correlation between At and Ac (B->1 => pdf->Rayleigh ) 015 % reltol = relative tolerance (default 1e-3) 016 % tol = absolute tolerance abs(int-intold) 017 % 018 % The size of F is the common size of X, A and B. A scalar input 019 % functions as a constant matrix of the same size as the other input. 020 % 021 % Tayfun calculates the correlation from the Spectral density as: 022 % 023 % B = sqrt( [int S(w)*cos(w*T) dw]^2+[int S(w)*sin(w*T) dw]^2 )/m0 024 % 025 % where T = Tm02, mean zero crossing period. 026 % Note that this is the same as the groupiness factor, Ka, found in 027 % spec2char 028 % 029 % Example: 030 % S = jonswap; 031 % a = sqrt(spec2mom(S,1)); 032 % b = spec2char(S,'Ka'); 033 % h = linspace(0,7*a)'; 034 % plot(h,tay90cdf(h,a,b)) 035 % 036 % See also tay90pdf, spec2char, gaussq 037 038 039 % References: 040 % [1] M. Aziz Tayfun (1990) "Distribution of large waveheights" 041 % Journal of the Waterway, port and coastal and ocean division 042 % vol 116 No.6 pp. 686-707 043 % [2] M. Aziz Tayfun (1981) "Distribution of crest-to-trough waveheights" 044 % Journal of the Waterway, port and coastal and ocean division 045 % vol 107 No.3 pp. 149-158 046 047 %tested 048 % history: 049 % revised pab 01.12.2000 050 % -fixed a bug in the normalization 051 % Per A. Brodtkorb 21.02.99 052 053 error(nargchk(3,4,nargin)) 054 055 if (nargin <4)| isempty(tol), 056 tol=1e-3; %relative accuracy of the estimates % sqrt(eps);% 057 end 058 059 [errorcode, x, a, b] =comnsize(x,a,b); 060 061 if errorcode > 0 062 error('h, a and b must be of common size or scalar.'); 063 end 064 xsiz = size(x); % remember old size 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 is not positive. 071 k1 = find(a <= 0| abs(b)>1); 072 if any(k1) 073 tmp = NaN; 074 y(k1) = tmp(ones(size(k1))); 075 end 076 077 078 079 gn = 2;% # of of base points to start the integration with 080 trace = 0;%used for debugging 081 082 083 if 0 084 % This is a trick to get the html documentation correct. 085 k = tay90pdf(1,1,0.1); 086 end 087 088 k=find(a > 0 & x >0 & abs(b)<1); 089 if any(k), 090 yk=y(k);xk = x(k); ak = a(k); bk = b(k); 091 092 [yk, tol1(k)]= gaussq('tay90pdf',0,xk, tol,[trace gn],ak,bk,tol); 093 % yk=yk./ak; 094 095 k2=find(yk>1 ); 096 if any(k2) 097 yk(k2)= ones(size(k2)); 098 end 099 y(k)=yk; 100 end 101 102 k4=find(a > 0 & x >= 0 & abs(b)==1); 103 if any(k4), % special case b==1 -> rayleigh distribution 104 y(k4) = wraylcdf(x(k4),2*a(k4)); 105 end 106 y=reshape(y,xsiz); 107 108
Comments or corrections to the WAFO group