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: This function is called by:

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