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

## CROSS-REFERENCE INFORMATION

This function calls:
 wbetapdf Beta probability density function betaln Logarithm of beta function. comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. inf Infinity. nan Not-a-Number.
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 %
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
046 % revised pab 14.10.1999
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