B04PDF Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data. CALL: f = b04pdf(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. B04PDF 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 laboratory storm waves. 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,b04pdf(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 (2004) joint (Scf,Hd) CDF of laboratory data. | |
Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data. |
001 function f = b04pdf(Hd,Scf,Hm0,Tm02, condon,normalizedInput) 002 %B04PDF Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data. 003 % 004 % CALL: f = b04pdf(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 % B04PDF 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 laboratory storm waves. 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,b04pdf(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 % Error 0.02 for cAx and 0.005 for cBx 068 cA1 = [ 3.04621335253845,... 069 -3.63962229567100,... 070 2.92552351922864]; 071 cA2 = [ -0.00166846704576,... 072 0.02393968243269,... 073 -0.14693767011981,... 074 0.54351431859392,... 075 -0.72674041059717,... 076 1.00000000000000]; 077 078 cB1 = [ -4.17697804696004,... 079 2.70403599166940,... 080 -1.96751680710343]; 081 cB2 = [ -0.00261584642057,... 082 0.03498298832599,... 083 -0.08180403982217,... 084 2.05476574371033,... 085 -0.77207518455457,... 086 1.00000000000000]; 087 088 A1 = polyval(cA1,h)./polyval(cA2,h); 089 B1 = exp(polyval(cB1,h)./polyval(cB2,h)); 090 h0 = 4.49999999999989; 091 k = find(h>h0); 092 if any(k) 093 % Linear extrapolation 094 slope1 = 5.37234232285048; 095 slope2 = 0.00587463478855; 096 A1(k) = (h(k)-h0)*slope1 + 23.14319996232766; 097 B1(k) = (h(k)-h0)*slope2 + 0.16439526054778; 098 end 099 % Rayleigh parameter 100 B0 = 0.70079568211392; 101 102 103 f = zeros(size(h)); 104 105 switch condon, 106 case 0, % regular pdf is returned 107 f = wraylpdf(h,B0).*wgampdf(s,A1,B1); 108 case 1, %pdf conditioned on x1 ie. p(x2|x1) 109 f = wgampdf(s,A1,B1); 110 case 3, % secret option used by mk87stat: returns x2*p(x2|x1) 111 f = s.*wgampdf(s,A1,B1); 112 case 4, % secret option used by mk87stat: returns x2.^2*p(x2|x1) 113 f = s.^2.*wgampdf(s,A1,B1); 114 case 5, % p(h)*P(V|h) is returned special case used by bmr00cdf 115 f = wraylpdf(h,B0).*wgamcdf(s,A1,B1); 116 case 6, % P(V|h) is returned special case used by bmr00cdf 117 f = wgamcdf(s,A1,B1); 118 case 7, % p(h)*(1-P(V|h)) is returned special case used by bmr00cdf 119 f = wraylpdf(h,B0).*(1-wgamcdf(s,A1,B1)); 120 otherwise error('unknown option') 121 end 122 if condon~=6 123 f = f./Hrms./Erms; 124 end 125 return 126 127 return 128 129 % Code used for finding the rational polynomials 130 % $$$ NFAC=8; 131 % $$$ BIG=1e30; 132 % $$$ MAXIT=5; 133 % $$$ 134 % $$$ dev=BIG; 135 % $$$ 136 % $$$ 137 % $$$ 138 % $$$ m = 7; 139 % $$$ k = 7; 140 % $$$ ncof=m+k+1; 141 % $$$ npt=NFAC*ncof; 142 % $$$ % Number of points where function is avaluated, i.e. fineness of mesh 143 % $$$ 144 % $$$ a = 0; b = 2.5 145 % $$$ 146 % $$$ x = zeros(npt,1); 147 % $$$ ix1=1:floor(npt/2-1); 148 % $$$ x(ix1)=a+(b-a).*sin(pi/2*(ix1-1)./(npt-1)).^2; 149 % $$$ ix2=floor(npt/2):npt; 150 % $$$ x(ix2)=b-(b-a).*sin(pi/2*(npt-ix2)./(npt-1)).^2; 151 % $$$ [Av , Bv, Cv]=dist2dsmfun(sphat,x); 152 % $$$ 153 % $$$ 154 % $$$ [cA1, cA2] = ratlsq([x,Av],m,k,a,b); 155 % $$$ 156 % $$$ [cB1, cB2] = ratlsq([x,log(Bv)],m,k,a,b); 157 % $$$ %[cA1, cA2] = ratlsq([x,Av],a,b,m,k); 158 % $$$ subplot(2,1,1) 159 % $$$ plot(x,polyval(cA1,x)./polyval(cA2,x),'g',x,Av,'r') 160 % $$$ subplot(2,1,2) 161 % $$$ plot(x,exp(polyval(cB1,x)./polyval(cB2,x)),'g',x,Bv,'r') 162 % $$$ 163 % $$$ 164 % $$$ 165 % $$$ m = 5; 166 % $$$ k = m+1; 167 % $$$ ncof=m+k+1; 168 % $$$ npt=NFAC*ncof; 169 % $$$ % Number of points where function is avaluated, i.e. fineness of mesh 170 % $$$ 171 % $$$ a = 0; b = 2 172 % $$$ 173 % $$$ x = zeros(npt,1); 174 % $$$ ix1=1:floor(npt/2-1); 175 % $$$ x(ix1)=a+(b-a).*sin(pi/2*(ix1-1)./(npt-1)).^2; 176 % $$$ ix2=floor(npt/2):npt; 177 % $$$ x(ix2)=b-(b-a).*sin(pi/2*(npt-ix2)./(npt-1)).^2; 178 % $$$ [Av , Bv, Cv]=dist2dsmfun(sphat,x); 179 % $$$ 180 % $$$ 181 % $$$ [cA1, cA2] = ratlsq([x,Av],m,k,a,b); 182 % $$$ 183 % $$$ [cB1, cB2] = ratlsq([x,log(Bv)],m,k,a,b); 184 % $$$ %[cA1, cA2] = ratlsq([x,Av],a,b,m,k); 185 % $$$ subplot(2,1,1) 186 % $$$ plot(x,polyval(cA1,x)./polyval(cA2,x),'g',x,Av,'r') 187 % $$$ subplot(2,1,2) 188 % $$$ plot(x,exp(polyval(cB1,x)./polyval(cB2,x)),'g',x,Bv,'r') 189 % $$$ 190 % $$$ N = 12 191 % $$$ x1 = chebroot(N).'*(b-a)/2+(b+a)/2 ; 192 % $$$ [Av , Bv, Cv]=dist2dsmfun(sphat,x1); 193 % $$$ cA1 =chebfit([x1 Av],n,a,b); 194 % $$$ cB1 =chebfit([x1 Bv],n,a,b); 195 % $$$ 196 % $$$ dist2dparamplot(phat,sphat) 197 % $$$ subplot(2,1,1), hold on 198 % $$$ plot(x,chebval(x,cA1,a,b),'g') 199 % $$$ subplot(2,1,2), hold on 200 % $$$ plot(x,chebval(x,cB1,a,b),'g') 201 % $$$ 202 % $$$ 203
Comments or corrections to the WAFO group