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
Beta probability density function | |
Logarithm of beta function. | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Infinity. | |
Not-a-Number. |
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
Comments or corrections to the WAFO group