Home > wafo > wavemodels > tay81pdf.m

tay81pdf

PURPOSE ^

Tayfun (1981) PDF of breaking limited wave heights

SYNOPSIS ^

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

DESCRIPTION ^

 TAY81PDF  Tayfun (1981) PDF of  breaking limited wave heights 
               /inf 
   f = x/a^2 * | u* J_o(u/b)^(b^2)*J_o(x/a* u) du *(1 - 4/pi*acos(b*a/x)*H(x/a-b) ) 
               /0 
  where J_o and H(.) is the zero order Bessel function and the heavyside function 
  
   CALL:  [f tol]= tay81pdf(H,A,B,reltol)  
             f = pdf   
             x = quantiles  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(x,tay81pdf(x,a,b),x,wraylpdf(x,Hm0/2),'r') 
  
  See also  tay81cdf, braylpdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function [y, tol1] = tay81pdf(x,a,b,tol) 
002 %TAY81PDF  Tayfun (1981) PDF of  breaking limited wave heights 
003 %              /inf 
004 %  f = x/a^2 * | u* J_o(u/b)^(b^2)*J_o(x/a* u) du *(1 - 4/pi*acos(b*a/x)*H(x/a-b) ) 
005 %              /0 
006 % where J_o and H(.) is the zero order Bessel function and the heavyside function 
007 % 
008 %  CALL:  [f tol]= tay81pdf(H,A,B,reltol)  
009 %            f = pdf   
010 %            x = quantiles  0 <= x <= sqrt(2)*B 
011 %            A = scale parameter 
012 %            B = maximum limiting value (B->inf => pdf->Rayleigh ) 
013 %       reltol = relative tolerance (default 1e-3) 
014 %          tol = absolute tolerance abs(int-intold) 
015 % 
016 % Tayfun set the scale parameter and  breaking limiting value to 
017 %     A = Hrms = 2*sqrt(2*m0)*sqrt(1-0.734*eps2^4) 
018 %     B = pi*tanh(k0*h)/(7*k0*2*sqrt(m0)) 
019 % where 
020 %    k0 = the mean apparent wave number 
021 %    mi = i'th spectral moment (i=0 => variance of the process) 
022 %  eps2 = sqrt(m0*m2/m1^2-1) spectral bandwidth 
023 %     h = water depth 
024 % 
025 % 
026 %  The size of F is the common size of X, A and B.  A scalar input    
027 %  functions as a constant matrix of the same size as the other input. 
028 % 
029 % Example: 
030 %  Hm0  = 7; Tp = 10; h = 150; 
031 %  S    = jonswap([],[Hm0,Tp]); 
032 %  eps2 = spec2char(S,'eps2'); 
033 %  Tm02 = spec2char(S,'Tm02'); 
034 %  a = Hm0/sqrt(2)*sqrt(1-0.734*eps2^4); 
035 %  k0   = w2k(2*pi/Tm02,[],h); 
036 %  b = pi*tanh(k0*h)/(7*k0*Hm0/2); 
037 %  x = linspace(0,Hm0)'; 
038 %  plot(x,tay81pdf(x,a,b),x,wraylpdf(x,Hm0/2),'r') 
039 % 
040 % See also  tay81cdf, braylpdf 
041  
042 %   Reference: 
043 %  M. Aziz Tayfun (1981) "Breaking limited waveheights" 
044 % Journal of the Waterway, port and coastal and ocean division  
045 %      vol 107 No.2 pp. 59-69 
046  
047 %tested on:  
048 %history: 
049 %  Per A. Brodtkorb 01.02.99 
050  
051 % TODO % not tested 
052 if (nargin <4)| isempty(tol),  
053   tol=1e-3; %relative accuracy of the estimates % sqrt(eps);% 
054 end 
055 if nargin <  3,  
056     error('Requires at least three input arguments.');  
057 end 
058  
059 [errorcode x a b] = comnsize(x,a,b.^2); 
060  
061 if errorcode > 0 
062     error('Requires non-scalar arguments to match in size.'); 
063 end 
064 [N M]=size(x); 
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 or B is not positive. 
071 k1 = find(a <= 0|b<= 0); 
072 if any(k1)  
073     tmp   = NaN; 
074     y(k1) = tmp(ones(size(k1))); 
075 end 
076  
077  
078 count_limit=20; 
079 gn=2^3;% # of of base points to start the integration with 
080 trace=logical(0);%used for debugging 
081  
082 k=find(a > 0 & x >0 &x<a.*sqrt(2*b) & b<inf&b> 0); 
083 if any(k), 
084   yk=y(k);xk = x(k); ak = a(k); bk = b(k); 
085   tmpy=yk;    tmpt=tmpy; %size(xk) 
086   uk_inf = traylfunzeros(xk,ak,bk,count_limit*10);% 
087   uk_inf = [zeros(size(k)) uk_inf(:,3:4:end)]; 
088                
089   if 0   
090     uk_inf= uk_inf(:,3:4:end); 
091     k1=find( uk_inf(:,1)<=7); 
092     if any(k1), 
093       tmpy=uk_inf(:,1); 
094       uk_inf(k1,:)=uk_inf(k1,:)+7-tmpy(k1,ones(1,size(uk_inf,2))); 
095     end 
096     %uk_inf=.5*ak.*sqrt(2.*bk); 
097     %uk_inf(1:10,1:4), size(uk_inf),size(k) 
098     [yk, tol1(k)]= gaussq('tay81fun',0,uk_inf(k),  tol/2,trace,gn,xk./ak,bk); 
099   end 
100   if 0 
101     % This is a trick to get the html documentation correct. 
102     k = tay81fun(1,1,2); 
103   end 
104    
105    
106   % Break out of the iteration loop for three reasons: 
107   %  1) the last update is very small (compared to int) 
108   %  2) the last update is very small (compared to tol) 
109   %  3) There are more than 20 iterations. This should NEVER happen.  
110   k1=find(k);ix=1; 
111   while(any(k1) & ix <= count_limit),  
112     if trace,figure(ix),end % for debugging  
113     [tmpy(k1), tmpt(k1)]= gaussq('tay81fun',uk_inf(k1,ix),uk_inf(k1,ix+1),... 
114                  tol/2/count_limit,trace,gn,xk(k1)./ak(k1),bk(k1)); 
115      
116     yk(k1)=yk(k1)+tmpy(k1);        
117     tol1(k(k1))=tol1(k(k1))+tmpt(k1); 
118      %find integrals which did not converge   
119     k1=find(tol1(k) < abs(tmpy) & tol/count_limit <2*abs(tmpy));  
120     
121     %k=find(tol1 > abs(tol*int)| tol1 > abs(tol)); 
122     %indices to integrals which did not converge 
123     ix = ix + 1; 
124   end 
125   %yk=yk.*xk./(ak.^2); 
126    
127   k2=find(yk<0 ); 
128   if any(k2) 
129     yk(k2)= zeros(size(k2)); 
130   end 
131  
132   k3=find(xk>ak.*sqrt(bk) ); 
133   if any(k3) 
134     yk(k3)= yk(k3).*(1-4/pi.*acos(ak(k3).*sqrt(bk(k3))./xk(k3)));  
135   end 
136   y(k)=yk.*xk./ak.^2; 
137 end 
138  
139 k4=find(a > 0 & x >= 0 & b==inf); 
140 if any(k4), % special case b==inf -> rayleigh distribution 
141  y(k4) = wraylpdf(x(k4),a(k4)/sqrt(2)); 
142 end  
143 y=reshape(y,N,M); 
144  
145  
146 function z=traylfunzeros(x,a,b,n), 
147 % internal function which calculates the approximate zerocrossings of  
148 %  the integrand:  u* J_o(u/sqrt(b))^b*J_o(x/a* u) 
149 nx=1:n; 
150 z1=a./x;  
151 z0=sqrt(b); 
152 k=find(mod(b,2)==0); % remove 
153 z0(k)=pi/2*(n+1)*max(z1);%size(x) 
154 z0=z0(:,ones(1,n)).* pi.*(nx(ones(size(x)),:)-0.25); 
155 z1=z1(:,ones(1,n)).* pi.*(nx(ones(size(x)),:)-0.25); 
156 z=sort([z0,z1],2); 
157 return 
158

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