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)) 
  
  See also  tay90cdf, spec2char, gaussq

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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

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