Home > wafo > wavemodels > b04jpdf.m

# b04jpdf

## PURPOSE

Brodtkorb (2004) joint (Scf,Hd) PDF from Japan Sea.

## SYNOPSIS

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

## DESCRIPTION

``` B04JPDF Brodtkorb (2004) joint (Scf,Hd) PDF from Japan Sea.

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

B04JPDF 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 storm waves from the Japan 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,b04jpdf(H,S,Hs,Tz))

See also mk87pdf```

## 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:
 b04jcdf Brodtkorb (2004) joint (Scf,Hd) CDF from Japan Sea. b04jpdf2 Brodtkorb (2004) joint (Scf,Hd) PDF from Japan Sea.

## SOURCE CODE

```001 function f = b04jpdf(Hd,Scf,Hm0,Tm02, condon,normalizedInput)
002 %B04JPDF Brodtkorb (2004) joint (Scf,Hd) PDF from Japan Sea.
003 %
004 %    CALL:  f = b04jpdf(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 %    B04JPDF 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 storm waves from the Japan Sea.
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,b04jpdf(H,S,Hs,Tz))
028 %
029 % See also mk87pdf
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 if 1,
068   % Error 0.07 for cAx and 0.01 for cBx
069   cA1 = [   2.42767540920461,...
070         -5.12101444148408,...
071         4.49060259316430];
072
073   cA2 =[   0.00572707656955,...
074        -0.10679331978077,...
075        0.81081632991506,...
076        -3.19389926947809,...
077        6.81410759646213,...
078        -7.14860241955959,...
079        2.40890619464040,...
080        1.00000000000000];
081   cB1 =[  -0.70310486398236,...
082       1.74679333441188,...
083       -2.47501962948775];
084
085   cB2 =[   0.00536662963627,...
086        -0.10086105895100,...
087        0.77434723641863,...
088        -3.10647077544176,...
089        6.88647211568827,...
090        -7.71847113553527,...
091        3.19755172037110,...
092        1.00000000000000];
093
094 end
095 A1 = polyval(cA1,h)./polyval(cA2,h);
096 B1 = exp(polyval(cB1,h)./polyval(cB2,h));
097 h0 = 4.49999999999989;
098 k = find(h>h0);
099 if any(k)
100   % Linear extrapolation
101  slope1 = 1.53878793498107;
102
103  slope2 =  0.01549012000025;
104
105   A1(k) =  (h(k)-h0)*slope1 + 17.80538369601295;
106   B1(k) = (h(k)-h0)*slope2  + 0.15016765927218;
107 end
108 % Rayleigh parameter
109 B0 = 0.67768680795379;
110
111
112 f   = zeros(size(h));
113
114 switch condon,
115  case 0, % regular pdf is returned
116   f = wraylpdf(h,B0).*wgampdf(s,A1,B1);
117  case 1, %pdf conditioned on x1 ie. p(x2|x1)
118   f = wgampdf(s,A1,B1);
119  case 3, % secret option  used by mk87stat: returns x2*p(x2|x1)
120   f = s.*wgampdf(s,A1,B1);
121  case 4, % secret option  used by mk87stat: returns x2.^2*p(x2|x1)
122   f = s.^2.*wgampdf(s,A1,B1);
123  case 5, % p(h)*P(V|h) is returned special case used by bmr00cdf
124   f = wraylpdf(h,B0).*wgamcdf(s,A1,B1);
125  case 6, % P(V|h) is returned special case used by bmr00cdf
126   f = wgamcdf(s,A1,B1);
127  case 7, % p(h)*(1-P(V|h)) is returned special case used by bmr00cdf
128   f = wraylpdf(h,B0).*(1-wgamcdf(s,A1,B1));
129  otherwise error('unknown option')
130 end
131 if condon~=6
132   f = f./Hrms./Erms;
133 end
134 return
135
136  return
137
138  % Code used for finding the rational polynomials
139 % \$\$\$ NFAC=8;
140 % \$\$\$ BIG=1e30;
141 % \$\$\$ MAXIT=5;
142 % \$\$\$
143 % \$\$\$ dev=BIG;
144 % \$\$\$
145 % \$\$\$
146 % \$\$\$
147 % \$\$\$ m = 7;
148 % \$\$\$ k = 7;
149 % \$\$\$ ncof=m+k+1;
150 % \$\$\$ npt=NFAC*ncof;
151 % \$\$\$ % Number of points where function is avaluated, i.e. fineness of mesh
152 % \$\$\$
153 % \$\$\$ a = 0; b = 2.5
154 % \$\$\$
155 % \$\$\$ x = zeros(npt,1);
156 % \$\$\$ ix1=1:floor(npt/2-1);
157 % \$\$\$ x(ix1)=a+(b-a).*sin(pi/2*(ix1-1)./(npt-1)).^2;
158 % \$\$\$ ix2=floor(npt/2):npt;
159 % \$\$\$ x(ix2)=b-(b-a).*sin(pi/2*(npt-ix2)./(npt-1)).^2;
160 % \$\$\$ [Av , Bv, Cv]=dist2dsmfun(sphat,x);
161 % \$\$\$
162 % \$\$\$
163 % \$\$\$ [cA1, cA2] = ratlsq([x,Av],m,k,a,b);
164 % \$\$\$
165 % \$\$\$ [cB1, cB2] = ratlsq([x,log(Bv)],m,k,a,b);
166 % \$\$\$ %[cA1, cA2] = ratlsq([x,Av],a,b,m,k);
167 % \$\$\$ subplot(2,1,1)
168 % \$\$\$  plot(x,polyval(cA1,x)./polyval(cA2,x),'g',x,Av,'r')
169 % \$\$\$  subplot(2,1,2)
170 % \$\$\$  plot(x,exp(polyval(cB1,x)./polyval(cB2,x)),'g',x,Bv,'r')
171 % \$\$\$
172 % \$\$\$
173 % \$\$\$
174 % \$\$\$ m = 5;
175 % \$\$\$ k = m+1;
176 % \$\$\$ ncof=m+k+1;
177 % \$\$\$ npt=NFAC*ncof;
178 % \$\$\$ % Number of points where function is avaluated, i.e. fineness of mesh
179 % \$\$\$
180 % \$\$\$ a = 0; b = 2
181 % \$\$\$
182 % \$\$\$ x = zeros(npt,1);
183 % \$\$\$ ix1=1:floor(npt/2-1);
184 % \$\$\$ x(ix1)=a+(b-a).*sin(pi/2*(ix1-1)./(npt-1)).^2;
185 % \$\$\$ ix2=floor(npt/2):npt;
186 % \$\$\$ x(ix2)=b-(b-a).*sin(pi/2*(npt-ix2)./(npt-1)).^2;
187 % \$\$\$ [Av , Bv, Cv]=dist2dsmfun(sphat,x);
188 % \$\$\$
189 % \$\$\$
190 % \$\$\$ [cA1, cA2] = ratlsq([x,Av],m,k,a,b);
191 % \$\$\$
192 % \$\$\$ [cB1, cB2] = ratlsq([x,log(Bv)],m,k,a,b);
193 % \$\$\$ %[cA1, cA2] = ratlsq([x,Av],a,b,m,k);
194 % \$\$\$ subplot(2,1,1)
195 % \$\$\$  plot(x,polyval(cA1,x)./polyval(cA2,x),'g',x,Av,'r')
196 % \$\$\$  subplot(2,1,2)
197 % \$\$\$  plot(x,exp(polyval(cB1,x)./polyval(cB2,x)),'g',x,Bv,'r')
198 % \$\$\$
199 % \$\$\$  N = 12
200 % \$\$\$  x1 = chebroot(N).'*(b-a)/2+(b+a)/2 ;
201 % \$\$\$  [Av , Bv, Cv]=dist2dsmfun(sphat,x1);
202 % \$\$\$  cA1 =chebfit([x1 Av],n,a,b);
203 % \$\$\$  cB1 =chebfit([x1 Bv],n,a,b);
204 % \$\$\$
205 % \$\$\$  dist2dparamplot(phat,sphat)
206 % \$\$\$  subplot(2,1,1), hold on
207 % \$\$\$  plot(x,chebval(x,cA1,a,b),'g')
208 % \$\$\$  subplot(2,1,2), hold on
209 % \$\$\$  plot(x,chebval(x,cB1,a,b),'g')
210 % \$\$\$
211 % \$\$\$
212```

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