Home > wafo > wavemodels > tay81cdf.m

tay81cdf

PURPOSE ^

Tayfun (1981) CDF of breaking limited wave heights

SYNOPSIS ^

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

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function [y, tol1] = tay81cdf(x,a,b,tol) 
002 %TAY81CDF Tayfun (1981) CDF of  breaking limited wave heights 
003 %                  /inf  
004 %   pdf = x/a^2 * | u* J_o(u/b)^(b^2)*J_o(xu) du *(1 - 4/pi*acos(b/x)*H(1-b/x) ) 
005 %                 /0 
006 % where J_o and H(.) is the zero order Bessel function and the heavyside function 
007 % 
008 %  CALL:  [F tol]= tay81cdf(x,A,B,reltol)  
009 % 
010 %            F = CDF 
011 %            x = wave heights  0 <= x <= sqrt(2)*B 
012 %            A = scale parameter 
013 %            B = maximum limiting value (B->inf => pdf->Rayleigh ) 
014 %       reltol = relative tolerance (default 1e-3) 
015 %           tol = absolute tolerance abs(int-intold) 
016 % 
017 % Tayfun set the scale parameter and  breaking limiting value to 
018 %     A = Hrms = 2*sqrt(2*m0)*sqrt(1-0.734*eps2^4) 
019 %     B = pi*tanh(k0*h)/(7*k0*2*sqrt(m0)) 
020 % where 
021 %    k0 = the mean apparent wave number 
022 %    mi = i'th spectral moment (i=0 => variance of the process) 
023 %  eps2 = sqrt(m0*m2/m1^2-1) spectral bandwidth 
024 %     h = water depth 
025 % 
026 % 
027 %The size of F is the common size of X, A and B.  A scalar input    
028 %    functions as a constant matrix of the same size as the other input. 
029 % 
030 % Example: 
031 %  Hm0  = 7; Tp = 10; h = 150; 
032 %  S    = jonswap([],[Hm0,Tp]); 
033 %  eps2 = spec2char(S,'eps2'); 
034 %  Tm02 = spec2char(S,'Tm02'); 
035 %  a = Hm0/sqrt(2)*sqrt(1-0.734*eps2^4); 
036 %  k0   = w2k(2*pi/Tm02,[],h); 
037 %  b = pi*tanh(k0*h)/(7*k0*Hm0/2); 
038 %  x = linspace(0,Hm0)'; 
039 %  plot(1:5,tay81cdf(1:5,a,b),x,wraylcdf(x,Hm0/2),'r') 
040 % 
041 % See also  tay81pdf, braylcdf 
042  
043 %   Reference: 
044 %  M. Aziz Tayfun (1981) "Breaking limited waveheights" 
045 % Journal of the Waterway, port and coastal and ocean division  
046 %      vol 107 No.2 pp. 59-69 
047  
048 %tested on: 
049 %history 
050 % 
051 %  Per A. Brodtkorb 01.02.99 
052  
053 % TODO % not tested 
054    
055 if (nargin <4)| isempty(tol),  
056   tol=1e-3;%sqrt(eps);%relative accuracy of the estimates 
057 end 
058 if nargin <  3,  
059     error('Requires at least three input arguments.');  
060 end 
061  
062 [errorcode x a b] = comnsize(x,a,b.^2); 
063  
064 if errorcode > 0 
065     error('Requires non-scalar arguments to match in size.'); 
066 end 
067  
068 % Initialize Y to zero. 
069 y=zeros(size(x)); 
070 tol1=y; 
071 % Return NaN if A or B is not positive. 
072 k1 = find(a <= 0|b<= 0); 
073 if any(k1)  
074     tmp   = NaN; 
075     y(k1) = tmp(ones(size(k1))); 
076 end 
077  
078 trace=logical(0);%[];%used for debugging 
079  
080 k=find(a > 0 & x > 0 &x<=a.*sqrt(2*b) & b<inf&b> 0); 
081 if any(k), 
082   yk = y(k); xk = x(k); ak = a(k); bk = b(k); 
083   yk = yk(:); xk = xk(:); ak = ak(:); bk = bk(:); 
084   if 0,all((ak==ak(1)).*(bk==bk(1))), 
085     [xk ind]=sort(xk); 
086     gn=2; 
087     [y(k), tol(k)]= gaussq('tay81pdf',[0;xk(1:end-1)] ,xk,tol,trace,gn,ak,bk,tol); 
088     %[y(k), tol(k)]= asimpson('tay81pdf',0,xk,tol,trace,ak,bk); 
089     y(k)=cumsum(y(k)); 
090     y(k)=y(k(ind)); %make sure the ordering is right 
091   else 
092     gn=2; 
093     [y(k), tol(k)]= gaussq('tay81pdf',0,xk,tol,trace,gn,ak,bk,tol); 
094     %[y(k), tol(k)]= asimpson('traylpdf',0,xk,tol,trace,ak,bk); 
095   end 
096 end 
097 if 0 
098   % This is a trick to get the html documentation correct. 
099   k = tay81pdf(1,1); 
100 end 
101  
102 k2=find((a > 0 & y>1) | (x>=a.*sqrt(2*b) )& b<inf & b> 0); 
103 if any(k2) 
104      y(k2)= ones(size(k2)); 
105 end 
106  
107 k4=find(a > 0 & x >= 0 & b==inf); 
108 if any(k4), % special case b==inf -> rayleigh distribution 
109  y(k4) = wraylcdf(x(k4),a(k4)/sqrt(2)); 
110 end  
111  
112 function z=traylfunzeros(x,a,b,n), 
113 % internal function which calculates the approximate zerocrossings of  
114 %  the integrand:  x/a* J_o(u/sqrt(b))^b*J_1(x/a* u) 
115 nx=1:n; 
116 z1=a./x;  
117 z0=sqrt(b); 
118 %k=find(mod(b,2)==0); % remove 
119 %z0(k)=pi/2*(n+1)*max(z1); 
120 z0=z0(:,ones(1,n)).* pi.*(nx(ones(size(x)),:)+0.25); 
121 z1=z1(:,ones(1,n)).* pi.*(nx(ones(size(x)),:)-0.25); 
122 z=sort([z0,z1],2); 
123 return 
124  
125  
126

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