B04JCDF Brodtkorb (2004) joint (Scf,Hd) CDF from Japan Sea. CALL: F = b04jcdf(Hd,Scf,Hs,Tz,tail) F = CDF Hd = zero down crossing wave height. Scf = crest front steepness. Hs = significant wave height. Tz = average zero down crossing period. tail = 1 if upper tail is calculated 0 if lower tail is calulated (default) B04CDF returns the joint CDF of (Scf, Hd) given Hs and Tz, 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 = Hs/sqrt(2); Erms = 5/4*Hs/(Tz^2); The size of F is the common size of the input arguments Examples: Hs = 5.5; Tz = 8.5; sc = 0.25; hc = 4; p = b04jcdf(hc,sc,Hs,Tz) % Prob(Hd<=hc,Scf<=sc|Hs,Tz) = 0.66 % Conditional probability of steep and high waves given seastates % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz) upperTail = 1; Hs = linspace(2.5,11.5,17); Tz = linspace(4.5,11,16); [T,H] = meshgrid(Tz,Hs); p = b04jcdf(hc,sc,H,T,upperTail); v = 10.^(-6:-1); contourf(Tz,Hs,log10(p),log10(v)) xlabel('Tz') ylabel('Hs') fcolorbar(log10(v)) See also mk87cdf, gaussq
Brodtkorb (2004) joint (Scf,Hd) PDF from Japan Sea. | |
Numerically evaluates a integral using a Gauss quadrature. | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. |
001 function [p, eps2, a, b] = b04cdf(Hd,Scf,Hs,Tz,tail) 002 %B04JCDF Brodtkorb (2004) joint (Scf,Hd) CDF from Japan Sea. 003 % 004 % CALL: F = b04jcdf(Hd,Scf,Hs,Tz,tail) 005 % 006 % F = CDF 007 % Hd = zero down crossing wave height. 008 % Scf = crest front steepness. 009 % Hs = significant wave height. 010 % Tz = average zero down crossing period. 011 % tail = 1 if upper tail is calculated 012 % 0 if lower tail is calulated (default) 013 % 014 % B04CDF returns the joint CDF of (Scf, Hd) given Hs and Tz, 015 % i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, given 016 % the seastate. The root mean square values of Hd and Scf (Hrms,Erms) are 017 % related to the significant waveheight and the average zero down 018 % crossing period by: 019 % Hrms = Hs/sqrt(2); 020 % Erms = 5/4*Hs/(Tz^2); 021 % 022 % The size of F is the common size of the input arguments 023 % 024 % Examples: 025 % Hs = 5.5; 026 % Tz = 8.5; 027 % sc = 0.25; 028 % hc = 4; 029 % p = b04jcdf(hc,sc,Hs,Tz) % Prob(Hd<=hc,Scf<=sc|Hs,Tz) = 0.66 030 % 031 % % Conditional probability of steep and high waves given seastates 032 % % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz) 033 % upperTail = 1; 034 % Hs = linspace(2.5,11.5,17); 035 % Tz = linspace(4.5,11,16); 036 % [T,H] = meshgrid(Tz,Hs); 037 % p = b04jcdf(hc,sc,H,T,upperTail); 038 % v = 10.^(-6:-1); 039 % contourf(Tz,Hs,log10(p),log10(v)) 040 % xlabel('Tz') 041 % ylabel('Hs') 042 % fcolorbar(log10(v)) 043 % 044 % See also mk87cdf, gaussq 045 046 % Reference 047 % P. A. Brodtkorb (2004), 048 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 049 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 050 % Trondheim, Norway. 051 052 053 %tested on matlab 5.2 054 %history: 055 % by Per A. brodtkorb July 2004 056 057 058 error(nargchk(3,5,nargin)) 059 060 if (nargin < 5|isempty(tail)),tail = 0; end 061 if (nargin < 4|isempty(Tz)),Tz = 8; end 062 if (nargin < 3|isempty(Hs)), Hs = 6; end 063 064 multipleSeaStates = any(prod(size(Hs))>1|prod(size(Tz))>1); 065 if multipleSeaStates 066 [errorcode, Scf,Hd,Tz,Hs] = comnsize(Scf,Hd,Tz,Hs); 067 else 068 [errorcode, Scf,Hd] = comnsize(Scf,Hd); 069 end 070 if errorcode > 0 071 error('Requires non-scalar arguments to match in size.'); 072 end 073 074 Hrms = Hs/sqrt(2); 075 Erms = 5/4*Hs./(Tz.^2); 076 %Erms = (0.0202 + 0.826*Hs./(Tz.^2)); 077 078 s = Scf./Erms; 079 hMax = 20; 080 h = min(Hd./Hrms,hMax); 081 082 eps2 = 1e-6; 083 084 085 p = zeros(size(Hd)); 086 %k0 = find((Hs<=(Tz-4)*13/6+4)); 087 upLimit = 6.5/1.4; 088 loLimit = 2.5/1.26;; 089 k0 = find((loLimit*sqrt(Hs)<Tz)); 090 if any(k0) 091 if multipleSeaStates 092 h = h(k0); 093 s = s(k0); 094 Hs = Hs(k0); 095 Tz = Tz(k0); 096 else 097 k0 = 1:prod(size(Hd)); 098 end 099 hlim = h; 100 101 102 h0 = 2.50704530736889; 103 normalizedInput = 1; 104 lowerTail = 0; 105 106 if 0 107 % This is a trick to get the html documentation correct. 108 k = b04jpdf(1,1,2,3); 109 end 110 111 if (tail == lowerTail) 112 k = find(h>h0); 113 hlim(k) = h0; 114 p(k0) = gaussq('b04jpdf',0,hlim,eps2/2,[],s,Hs,Tz,5,normalizedInput)... 115 + gaussq('b04jpdf',hlim,h,eps2/2,[],s,Hs,Tz,5,normalizedInput); 116 else 117 k = find(h<h0); 118 hlim(k) = h0; 119 p(k0) = gaussq('b04jpdf',h,hlim,eps2/2,[],s,Hs,Tz,7,normalizedInput)... 120 + gaussq('b04jpdf',hlim,hMax,eps2/2,[],s,Hs,Tz,7,normalizedInput); 121 end 122 end 123 return 124 125
Comments or corrections to the WAFO group