Home > wafo > wavemodels > bmr00pdf.m

# bmr00pdf

## PURPOSE

Brodtkorb et.al (2000) joint (Scf,Hd) PDF from North Sea.

## SYNOPSIS

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

## DESCRIPTION

``` BMR00PDF Brodtkorb et.al (2000) joint (Scf,Hd) PDF from North Sea.

CALL:  f = bmr00pdf(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.

BMR00PDF 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 is a revised distribution of MK87 and is fitted to storm waves
from 1995 obtained from the Draupner field in the North Sea.
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,bmr00pdf(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:
 bmr00cdf Brodtkorb et.al. (2000) joint (Scf,Hd) CDF from North Sea. bmr00pdf2 Brodtkorb et.al (2000) joint (Scf,Hd) PDF from North Sea.

## SOURCE CODE

```001 function f = bmr00pdf(Hd,Scf,Hm0,Tm02, condon,normalizedInput)
002 %BMR00PDF Brodtkorb et.al (2000) joint (Scf,Hd) PDF from North Sea.
003 %
004 %    CALL:  f = bmr00pdf(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 %    BMR00PDF 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 is a revised distribution of MK87 and is fitted to storm waves
020 %    from 1995 obtained from the Draupner field in the North Sea.
021 %      The size of f is the common size of  the input arguments
022 %
023 %    Example:
024 %     Hs = 7;Tz=10;
025 %     h  = linspace(0,3*Hs)';
026 %     s  = linspace(0,4*Hs/Tz^2)';
027 %     [S,H] = meshgrid(s,h);
028 %     contour(s,h,bmr00pdf(H,S,Hs,Tz))
029 %
031
032 % Reference
033 %  Brodtkorb, P.A. and Myrhaug, D. and Rue, H. (2000)
034 %  "Joint Distributions of Wave Height and Wave Steepness Parameters",
035 % In Proc. 27'th Int. Conf. on Coastal Eng., ICCE, Sydney, Australia },
036 %  vol. 1, pp. 545--558, Paper No. 162
037
038 % By pab 15 July 2004
039
040
041
042
043 error(nargchk(3,6,nargin))
044
045 if nargin < 6|isempty(normalizedInput),  normalizedInput=0; end
046 if nargin < 5|isempty(condon),  condon=0; end % regular pdf is returned
047 if nargin < 4|isempty(Tm02),  Tm02=8;end
048 if nargin < 3|isempty(Hm0),  Hm0=6;end
049
050 [errorcode, Scf,Hd,Tm02,Hm0] = comnsize(Scf,Hd,Tm02,Hm0);
051 if errorcode > 0
052     error('Requires non-scalar arguments to match in size.');
053 end
054 if (normalizedInput>0)
055   Hrms = 1;
056   Erms = 1;
057 else
058   %Hm00 = 6.7;
059   %Tm020 = 8.3
060   Hrms = Hm0/sqrt(2);
061   Erms = 5/4*Hm0./(Tm02.^2);
062   %Erms = (0.0202 + 0.826*Hm0./(Tm02.^2))/2;
063 end
064 s = Scf./Erms;
065 h = Hd./Hrms;
066
067 % Conditinal gamma dist.  parameters for steepness
068 if 1,
069   % Error 0.02 for cAx and 0.02 for cBx
070   cA1 = [ 4.19193631070257,...
071       -12.03484717258579,...
072       15.95695603911886,...
073       -7.31618437661838,...
074       3.03590705490789];
075   cA2 = [  -0.02838644340635,
076        0.29409707273437,
077        -0.79138075176131,
078        1.20024905431293,
079        -0.58192823600509,
080        1.00000000000000];
081   cB1 =[ -5.96970985003487,...
082   23.24982082461128,...
083  -37.69215199928053,...
084    5.18520245368619,...
085   -1.87015799594426];
086   cB2 = [ -0.47236526233572,
087       5.71130418016660,
088       -17.45884089250511,
089       21.60781934904390,
090       0.11061694237412,
091       1.00000000000000];
092 else
093   % Error 0.08 for cAx and 0.03 for cBx
094   cA1 = [3.39863744229230, 0.42853512533872,2.97401909675030];
095   cA2 = [0.26749624884949, -1.47019799044842,2.17197357932780,1];
096   cB1 = [-15.50443060097295, -0.21474910268201,  -1.87709428385351];
097   cB2 = [-0.49720581030185, 6.35119691685354, 3.88968350685379, 1];
098  end
099 A1 = polyval(cA1,h)./polyval(cA2,h);
100 B1 = exp(polyval(cB1,h)./polyval(cB2,h));
101 h0 = 3.79937912772197;
102 k = find(h>h0);
103 if any(k)
104   % Linear extrapolation
105  slope1 = 14.03628528937437;
106  slope2 = -0.37497306844631;
107   A1(k) =  (h(k)-h0)*slope1 + 36.15168757587687;
108   B1(k) = exp((h(k)-h0)*slope2 + log(0.05673842727364));
109   %3.37291820522257  30.10141068490231   0.06979731246760
110   %3.41168737999524  30.65143585680909   0.06861014108633
111 end
112 % Rayleigh parameter
113 B0 = 0.69233682309018;
114
115
116
117 f   = zeros(size(h));
118
119 switch condon,
120  case 0, % regular pdf is returned
121   f = wraylpdf(h,B0).*wgampdf(s,A1,B1);
122  case 1, %pdf conditioned on x1 ie. p(x2|x1)
123   f = wgampdf(s,A1,B1);
124  case 3, % secret option  used by mk87stat: returns x2*p(x2|x1)
125   f = s.*wgampdf(s,A1,B1);
126  case 4, % secret option  used by mk87stat: returns x2.^2*p(x2|x1)
127   f = s.^2.*wgampdf(s,A1,B1);
128  case 5, % p(h)*P(V|h) is returned special case used by bmr00cdf
129   f = wraylpdf(h,B0).*wgamcdf(s,A1,B1);
130  case 6, % P(V|h) is returned special case used by bmr00cdf
131   f = wgamcdf(s,A1,B1);
132  case 7, % p(h)*(1-P(V|h)) is returned special case used by bmr00cdf
133   f = wraylpdf(h,B0).*(1-wgamcdf(s,A1,B1));
134  otherwise error('unknown option')
135 end
136 if condon~=6
137   f = f./Hrms./Erms;
138 end
139 return
140
141  return
142
143  % Code used for finding the rational polynomials
144 % \$\$\$ NFAC=8;
145 % \$\$\$ BIG=1e30;
146 % \$\$\$ MAXIT=5;
147 % \$\$\$
148 % \$\$\$ dev=BIG;
149 % \$\$\$
150 % \$\$\$
151 % \$\$\$
152 % \$\$\$ m = 7;
153 % \$\$\$ k = 7;
154 % \$\$\$ ncof=m+k+1;
155 % \$\$\$ npt=NFAC*ncof;
156 % \$\$\$ % Number of points where function is avaluated, i.e. fineness of mesh
157 % \$\$\$
158 % \$\$\$ a = 0; b = 2.5
159 % \$\$\$
160 % \$\$\$ x = zeros(npt,1);
161 % \$\$\$ ix1=1:floor(npt/2-1);
162 % \$\$\$ x(ix1)=a+(b-a).*sin(pi/2*(ix1-1)./(npt-1)).^2;
163 % \$\$\$ ix2=floor(npt/2):npt;
164 % \$\$\$ x(ix2)=b-(b-a).*sin(pi/2*(npt-ix2)./(npt-1)).^2;
165 % \$\$\$ [Av , Bv, Cv]=dist2dsmfun(sphat,x);
166 % \$\$\$
167 % \$\$\$
168 % \$\$\$ [cA1, cA2] = ratlsq([x,Av],m,k,a,b);
169 % \$\$\$
170 % \$\$\$ [cB1, cB2] = ratlsq([x,log(Bv)],m,k,a,b);
171 % \$\$\$ %[cA1, cA2] = ratlsq([x,Av],a,b,m,k);
172 % \$\$\$ subplot(2,1,1)
173 % \$\$\$  plot(x,polyval(cA1,x)./polyval(cA2,x),'g',x,Av,'r')
174 % \$\$\$  subplot(2,1,2)
175 % \$\$\$  plot(x,exp(polyval(cB1,x)./polyval(cB2,x)),'g',x,Bv,'r')
176 % \$\$\$
177 % \$\$\$
178 % \$\$\$
179 % \$\$\$ m = 3;
180 % \$\$\$ k = m;
181 % \$\$\$ ncof=m+k+1;
182 % \$\$\$ npt=NFAC*ncof;
183 % \$\$\$ % Number of points where function is avaluated, i.e. fineness of mesh
184 % \$\$\$
185 % \$\$\$ a = 0; b = 2
186 % \$\$\$
187 % \$\$\$ x = zeros(npt,1);
188 % \$\$\$ ix1=1:floor(npt/2-1);
189 % \$\$\$ x(ix1)=a+(b-a).*sin(pi/2*(ix1-1)./(npt-1)).^2;
190 % \$\$\$ ix2=floor(npt/2):npt;
191 % \$\$\$ x(ix2)=b-(b-a).*sin(pi/2*(npt-ix2)./(npt-1)).^2;
192 % \$\$\$ [Av , Bv, Cv]=dist2dsmfun(sphat,x);
193 % \$\$\$
194 % \$\$\$
195 % \$\$\$ [cA1, cA2] = ratlsq([x,Av],m,k,a,b);
196 % \$\$\$
197 % \$\$\$ [cB1, cB2] = ratlsq([x,log(Bv)],m,k,a,b);
198 % \$\$\$ %[cA1, cA2] = ratlsq([x,Av],a,b,m,k);
199 % \$\$\$ subplot(2,1,1)
200 % \$\$\$  plot(x,polyval(cA1,x)./polyval(cA2,x),'g',x,Av,'r')
201 % \$\$\$  subplot(2,1,2)
202 % \$\$\$  plot(x,exp(polyval(cB1,x)./polyval(cB2,x)),'g',x,Bv,'r')
203 % \$\$\$
204 % \$\$\$  N = 12
205 % \$\$\$  x1 = chebroot(N).'*(b-a)/2+(b+a)/2 ;
206 % \$\$\$  [Av , Bv, Cv]=dist2dsmfun(sphat,x1);
207 % \$\$\$  cA1 =chebfit([x1 Av],n,a,b);
208 % \$\$\$  cB1 =chebfit([x1 Bv],n,a,b);
209 % \$\$\$
210 % \$\$\$  dist2dparamplot(phat,sphat)
211 % \$\$\$  subplot(2,1,1), hold on
212 % \$\$\$  plot(x,chebval(x,cA1,a,b),'g')
213 % \$\$\$  subplot(2,1,2), hold on
214 % \$\$\$  plot(x,chebval(x,cB1,a,b),'g')```

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