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)) 
  
  See also mk87pdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 % 
030 % See also mk87pdf 
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