Home > wafo > wavemodels > b04pdf.m

# b04pdf

## PURPOSE

Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data.

## SYNOPSIS

f = b04pdf(Hd,Scf,Hm0,Tm02, condon,normalizedInput)

## DESCRIPTION

``` B04PDF Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data.

CALL:  f = b04pdf(Hd,Scf,Hm0,Tm02)

f  = density
Hd  = zero down crossing wave height
Scf = crest front steepness
Hm0 = significant wave height.
Tm02 = average zero down crossing period.

B04PDF returns the joint PDF of (Scf, Hd) given Hm0 and Tm02,
i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height,
given the seastate. The root mean square values of Hd and Scf
(Hrms,Erms) are related to the significant waveheight and the
average zero down crossing period by:
Hrms = Hm0/sqrt(2);
Erms = 5/4*Hm0/(Tm02^2);
This distribution  is fitted to laboratory storm waves.
The size of f is the common size of  the input arguments

Example:
Hs = 7;Tz=10;
h  = linspace(0,3*Hs)';
s  = linspace(0,4*Hs/Tz^2)';
[S,H] = meshgrid(s,h);
contour(s,h,b04pdf(H,S,Hs,Tz))

## CROSS-REFERENCE INFORMATION

This function calls:
 wgamcdf Gamma cumulative distribution function wgampdf Gamma probability density function wraylpdf Rayleigh probability density function comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. polyval Evaluate polynomial.
This function is called by:
 b04cdf Brodtkorb (2004) joint (Scf,Hd) CDF of laboratory data. b04pdf2 Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data.

## SOURCE CODE

```001 function f = b04pdf(Hd,Scf,Hm0,Tm02, condon,normalizedInput)
002 %B04PDF Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data.
003 %
004 %    CALL:  f = b04pdf(Hd,Scf,Hm0,Tm02)
005 %
006 %       f  = density
007 %      Hd  = zero down crossing wave height
008 %      Scf = crest front steepness
009 %      Hm0 = significant wave height.
010 %     Tm02 = average zero down crossing period.
011 %
012 %    B04PDF returns the joint PDF of (Scf, Hd) given Hm0 and Tm02,
013 %    i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height,
014 %    given the seastate. The root mean square values of Hd and Scf
015 %    (Hrms,Erms) are related to the significant waveheight and the
016 %    average zero down crossing period by:
017 %                Hrms = Hm0/sqrt(2);
018 %                Erms = 5/4*Hm0/(Tm02^2);
019 %    This distribution  is fitted to laboratory storm waves.
020 %    The size of f is the common size of  the input arguments
021 %
022 %    Example:
023 %     Hs = 7;Tz=10;
024 %     h  = linspace(0,3*Hs)';
025 %     s  = linspace(0,4*Hs/Tz^2)';
026 %     [S,H] = meshgrid(s,h);
027 %     contour(s,h,b04pdf(H,S,Hs,Tz))
028 %
030
031 % Reference
032 % P. A. Brodtkorb (2004),
033 % The probability of Occurrence of dangerous Wave Situations at Sea.
034 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
035 % Trondheim, Norway.
036
037 % By pab 15 July 2004
038
039
040
041
042 error(nargchk(3,6,nargin))
043
044 if nargin < 6|isempty(normalizedInput),  normalizedInput=0; end
045 if nargin < 5|isempty(condon),  condon=0; end % regular pdf is returned
046 if nargin < 4|isempty(Tm02),  Tm02=8;end
047 if nargin < 3|isempty(Hm0),  Hm0=6;end
048
049 [errorcode, Scf,Hd,Tm02,Hm0] = comnsize(Scf,Hd,Tm02,Hm0);
050 if errorcode > 0
051     error('Requires non-scalar arguments to match in size.');
052 end
053 if (normalizedInput>0)
054   Hrms = 1;
055   Erms = 1;
056 else
057   %Hm00 = 6.7;
058   %Tm020 = 8.3
059   Hrms = Hm0/sqrt(2);
060   Erms = 5/4*Hm0./(Tm02.^2);
061   %Erms = (0.0202 + 0.826*Hm0./(Tm02.^2))/2;
062 end
063 s = Scf./Erms;
064 h = Hd./Hrms;
065
066 % Conditinal gamma dist.  parameters for steepness
067 % Error 0.02 for cAx and 0.005 for cBx
068 cA1 = [   3.04621335253845,...
069       -3.63962229567100,...
070         2.92552351922864];
071 cA2 = [  -0.00166846704576,...
072      0.02393968243269,...
073      -0.14693767011981,...
074      0.54351431859392,...
075      -0.72674041059717,...
076      1.00000000000000];
077
078 cB1 = [  -4.17697804696004,...
079      2.70403599166940,...
080      -1.96751680710343];
081 cB2 = [  -0.00261584642057,...
082      0.03498298832599,...
083      -0.08180403982217,...
084      2.05476574371033,...
085      -0.77207518455457,...
086      1.00000000000000];
087
088 A1 = polyval(cA1,h)./polyval(cA2,h);
089 B1 = exp(polyval(cB1,h)./polyval(cB2,h));
090 h0 = 4.49999999999989;
091 k = find(h>h0);
092 if any(k)
093   % Linear extrapolation
094   slope1 = 5.37234232285048;
095   slope2 = 0.00587463478855;
096   A1(k) =  (h(k)-h0)*slope1 + 23.14319996232766;
097   B1(k) = (h(k)-h0)*slope2 + 0.16439526054778;
098 end
099 % Rayleigh parameter
100 B0 = 0.70079568211392;
101
102
103 f   = zeros(size(h));
104
105 switch condon,
106  case 0, % regular pdf is returned
107   f = wraylpdf(h,B0).*wgampdf(s,A1,B1);
108  case 1, %pdf conditioned on x1 ie. p(x2|x1)
109   f = wgampdf(s,A1,B1);
110  case 3, % secret option  used by mk87stat: returns x2*p(x2|x1)
111   f = s.*wgampdf(s,A1,B1);
112  case 4, % secret option  used by mk87stat: returns x2.^2*p(x2|x1)
113   f = s.^2.*wgampdf(s,A1,B1);
114  case 5, % p(h)*P(V|h) is returned special case used by bmr00cdf
115   f = wraylpdf(h,B0).*wgamcdf(s,A1,B1);
116  case 6, % P(V|h) is returned special case used by bmr00cdf
117   f = wgamcdf(s,A1,B1);
118  case 7, % p(h)*(1-P(V|h)) is returned special case used by bmr00cdf
119   f = wraylpdf(h,B0).*(1-wgamcdf(s,A1,B1));
120  otherwise error('unknown option')
121 end
122 if condon~=6
123   f = f./Hrms./Erms;
124 end
125 return
126
127  return
128
129  % Code used for finding the rational polynomials
130 % \$\$\$ NFAC=8;
131 % \$\$\$ BIG=1e30;
132 % \$\$\$ MAXIT=5;
133 % \$\$\$
134 % \$\$\$ dev=BIG;
135 % \$\$\$
136 % \$\$\$
137 % \$\$\$
138 % \$\$\$ m = 7;
139 % \$\$\$ k = 7;
140 % \$\$\$ ncof=m+k+1;
141 % \$\$\$ npt=NFAC*ncof;
142 % \$\$\$ % Number of points where function is avaluated, i.e. fineness of mesh
143 % \$\$\$
144 % \$\$\$ a = 0; b = 2.5
145 % \$\$\$
146 % \$\$\$ x = zeros(npt,1);
147 % \$\$\$ ix1=1:floor(npt/2-1);
148 % \$\$\$ x(ix1)=a+(b-a).*sin(pi/2*(ix1-1)./(npt-1)).^2;
149 % \$\$\$ ix2=floor(npt/2):npt;
150 % \$\$\$ x(ix2)=b-(b-a).*sin(pi/2*(npt-ix2)./(npt-1)).^2;
151 % \$\$\$ [Av , Bv, Cv]=dist2dsmfun(sphat,x);
152 % \$\$\$
153 % \$\$\$
154 % \$\$\$ [cA1, cA2] = ratlsq([x,Av],m,k,a,b);
155 % \$\$\$
156 % \$\$\$ [cB1, cB2] = ratlsq([x,log(Bv)],m,k,a,b);
157 % \$\$\$ %[cA1, cA2] = ratlsq([x,Av],a,b,m,k);
158 % \$\$\$ subplot(2,1,1)
159 % \$\$\$  plot(x,polyval(cA1,x)./polyval(cA2,x),'g',x,Av,'r')
160 % \$\$\$  subplot(2,1,2)
161 % \$\$\$  plot(x,exp(polyval(cB1,x)./polyval(cB2,x)),'g',x,Bv,'r')
162 % \$\$\$
163 % \$\$\$
164 % \$\$\$
165 % \$\$\$ m = 5;
166 % \$\$\$ k = m+1;
167 % \$\$\$ ncof=m+k+1;
168 % \$\$\$ npt=NFAC*ncof;
169 % \$\$\$ % Number of points where function is avaluated, i.e. fineness of mesh
170 % \$\$\$
171 % \$\$\$ a = 0; b = 2
172 % \$\$\$
173 % \$\$\$ x = zeros(npt,1);
174 % \$\$\$ ix1=1:floor(npt/2-1);
175 % \$\$\$ x(ix1)=a+(b-a).*sin(pi/2*(ix1-1)./(npt-1)).^2;
176 % \$\$\$ ix2=floor(npt/2):npt;
177 % \$\$\$ x(ix2)=b-(b-a).*sin(pi/2*(npt-ix2)./(npt-1)).^2;
178 % \$\$\$ [Av , Bv, Cv]=dist2dsmfun(sphat,x);
179 % \$\$\$
180 % \$\$\$
181 % \$\$\$ [cA1, cA2] = ratlsq([x,Av],m,k,a,b);
182 % \$\$\$
183 % \$\$\$ [cB1, cB2] = ratlsq([x,log(Bv)],m,k,a,b);
184 % \$\$\$ %[cA1, cA2] = ratlsq([x,Av],a,b,m,k);
185 % \$\$\$ subplot(2,1,1)
186 % \$\$\$  plot(x,polyval(cA1,x)./polyval(cA2,x),'g',x,Av,'r')
187 % \$\$\$  subplot(2,1,2)
188 % \$\$\$  plot(x,exp(polyval(cB1,x)./polyval(cB2,x)),'g',x,Bv,'r')
189 % \$\$\$
190 % \$\$\$  N = 12
191 % \$\$\$  x1 = chebroot(N).'*(b-a)/2+(b+a)/2 ;
192 % \$\$\$  [Av , Bv, Cv]=dist2dsmfun(sphat,x1);
193 % \$\$\$  cA1 =chebfit([x1 Av],n,a,b);
194 % \$\$\$  cB1 =chebfit([x1 Bv],n,a,b);
195 % \$\$\$
196 % \$\$\$  dist2dparamplot(phat,sphat)
197 % \$\$\$  subplot(2,1,1), hold on
198 % \$\$\$  plot(x,chebval(x,cA1,a,b),'g')
199 % \$\$\$  subplot(2,1,2), hold on
200 % \$\$\$  plot(x,chebval(x,cB1,a,b),'g')
201 % \$\$\$
202 % \$\$\$
203```

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