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
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 from Japan Sea. | |
Brodtkorb (2004) joint (Scf,Hd) PDF from Japan Sea. |
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
Comments or corrections to the WAFO group