Home > wafo > wavemodels > braylpdf.m

braylpdf

PURPOSE ^

Beta Rayleigh PDF of wave heigth

SYNOPSIS ^

y = braylpdf(x,a,b,c)

DESCRIPTION ^

 BRAYLPDF Beta Rayleigh PDF of wave heigth 
                       
    f = 2*(a+b-1)!/((a-1)! * (b-1)!)*h^(2a-1)*(1-(h/c)^2)^(b-1)/c^(2a) 
                         
   CALL:   f = braylpdf(h,a,b,c);  
  
        f = pdf 
        h = waveheigth (0 <= h <= c) 
        a = abs(k1*(k2-k1)/(k1^2-k2))  
        b = abs((1-k1)*(k2-k1)/(k1^2-k2))  
        c = Hb, breaking wave height approximated by water depth, d. 
  where 
       k1 = E(H^2)/Hb^2 
       k2 = E(H^4)/Hb^4 
   E(H^2) = .5*exp(0.00272*(d/g*Tp^2)^(-0.834))*Hm0^2 
   E(H^2) = .5*exp(0.00046*(d/g*Tp^2)^(-1.208))*Hm0^2 
      Hm0 = significant waveheight 
      Tp  = modal period of wave spectrum 
  
     The size of F is the common size of H, A, B and C.  A scalar input    
     functions as a constant matrix of the same size as the other input. 
  
  Example: % Compare with Rayleigh distribution 
   Hm0 = 7;Tp = 11;d = 50; g = gravity; 
   k1  = .5*exp(0.00272*(d/g*Tp^2)^(-0.834))*Hm0^2/d^2; 
   k2  = .5*exp(0.00046*(d/g*Tp^2)^(-1.208))*Hm0^2/d^4; 
   a   = abs(k1*(k2-k1)/(k1^2-k2));  
   b   = abs((1-k1)*(k2-k1)/(k1^2-k2)); 
   h   = linspace(0,2*Hm0)'; 
   plot(h,braylpdf(h,a,b,d),'r',h,wraylpdf(h,Hm0/2)) 
  
  See also  wbetapdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function y = braylpdf(x,a,b,c) 
002 %BRAYLPDF Beta Rayleigh PDF of wave heigth 
003 %                      
004 %   f = 2*(a+b-1)!/((a-1)! * (b-1)!)*h^(2a-1)*(1-(h/c)^2)^(b-1)/c^(2a) 
005 %                        
006 %  CALL:   f = braylpdf(h,a,b,c);  
007 % 
008 %       f = pdf 
009 %       h = waveheigth (0 <= h <= c) 
010 %       a = abs(k1*(k2-k1)/(k1^2-k2))  
011 %       b = abs((1-k1)*(k2-k1)/(k1^2-k2))  
012 %       c = Hb, breaking wave height approximated by water depth, d. 
013 % where 
014 %      k1 = E(H^2)/Hb^2 
015 %      k2 = E(H^4)/Hb^4 
016 %  E(H^2) = .5*exp(0.00272*(d/g*Tp^2)^(-0.834))*Hm0^2 
017 %  E(H^2) = .5*exp(0.00046*(d/g*Tp^2)^(-1.208))*Hm0^2 
018 %     Hm0 = significant waveheight 
019 %     Tp  = modal period of wave spectrum 
020 % 
021 %    The size of F is the common size of H, A, B and C.  A scalar input    
022 %    functions as a constant matrix of the same size as the other input. 
023 % 
024 % Example: % Compare with Rayleigh distribution 
025 %  Hm0 = 7;Tp = 11;d = 50; g = gravity; 
026 %  k1  = .5*exp(0.00272*(d/g*Tp^2)^(-0.834))*Hm0^2/d^2; 
027 %  k2  = .5*exp(0.00046*(d/g*Tp^2)^(-1.208))*Hm0^2/d^4; 
028 %  a   = abs(k1*(k2-k1)/(k1^2-k2));  
029 %  b   = abs((1-k1)*(k2-k1)/(k1^2-k2)); 
030 %  h   = linspace(0,2*Hm0)'; 
031 %  plot(h,braylpdf(h,a,b,d),'r',h,wraylpdf(h,Hm0/2)) 
032 % 
033 % See also  wbetapdf 
034  
035  
036 %  
037 %   Reference: 
038 %       Michel K. Ochi (1998), 
039 %      "OCEAN WAVES, The stochastic approach", 
040 %       OCEAN TECHNOLOGY series 6, Cambridge, pp 279. (pd of peaks to trough)  
041  
042 %tested on: matlab 5.2 
043 %History: 
044 % revised pab 31.03.2001 
045 % added example 
046 % revised pab 14.10.1999 
047 % updated help header 
048 % by  Per A. Brodtkorb 21.02.99 
049  
050 error(nargchk(4,4,nargin)) 
051  
052  
053 [errorcode, x, a ,b, c] = comnsize(x,a,b,c); 
054  
055 if errorcode > 0 
056     error('h, a, b and c must be of common size or scalar.'); 
057 end 
058  
059  
060 % Initialize Y to zero. 
061 y=zeros(size(x)); 
062  
063  
064 k=find(a > 0 & x >=0 & b>0 & c>0 & x<=c & ~((x == 0 & a < .5) | (x == c & b < 1))); 
065 if any(k), 
066   xk = x(k); ak = a(k); bk = b(k);ck=c(k); 
067  switch 1, % choose between different implementations 
068   case 1,  
069          y(k)=2*(xk./ck).^(2*ak-1)./ck.*(1-(xk./ck).^2).^(bk-1).*exp(-betaln(ak,bk)); 
070  case 2, 
071         tmp(k) = (2*ak - 1).*log(xk./ck)-log(ck) + (bk - 1).*log((1 - (xk./ck).^2)) - betaln(ak,bk); 
072         y(k) = 2*exp(tmp(k)); 
073  case 3, 
074       y(k)=2*wbetapdf((xk./ck).^2,ak,bk).*xk./ck.^2; 
075  end 
076 end 
077  
078 % Return Inf for x = 0 and a < 1 or x = 1 and b < 1. 
079 % Required for non-IEEE machines. 
080 k2 = find((x == 0 & a < .5) | (x == c & b < 1)); 
081 if any(k2) 
082     tmp = Inf; 
083     y(k2) = tmp(ones(size(k2)));  
084 end 
085  
086 % Return NaN if A,B or C  is not positive. 
087 k1 = find(a <= 0| b<=0|c<=0); 
088 if any(k1)  
089     tmp   = NaN; 
090     y(k1) = tmp(ones(size(k1))); 
091 end 
092

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