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
Gamma cumulative distribution function | |
Gamma probability density function | |
Rayleigh probability density function | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Evaluate polynomial. |
Brodtkorb et.al. (2000) joint (Scf,Hd) CDF from North Sea. | |
Brodtkorb et.al (2000) joint (Scf,Hd) PDF from North Sea. |
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')
Comments or corrections to the WAFO group