B04CDF Brodtkorb (2004) joint (Scf,Hd) CDF of laboratory data. CALL: F = b04cdf(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 = b04cdf(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,18.5,17); Tz = linspace(4.5,19.5,16); [T,H] = meshgrid(Tz,Hs); p = b04cdf(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 of laboratory data. | |
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 %B04CDF Brodtkorb (2004) joint (Scf,Hd) CDF of laboratory data. 003 % 004 % CALL: F = b04cdf(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 = b04cdf(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,18.5,17); 035 % Tz = linspace(4.5,19.5,16); 036 % [T,H] = meshgrid(Tz,Hs); 037 % p = b04cdf(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 %tested on matlab 5.2 053 %history: 054 % by Per A. brodtkorb July 2004 055 056 057 error(nargchk(3,5,nargin)) 058 059 if (nargin < 5|isempty(tail)),tail = 0; end 060 if (nargin < 4|isempty(Tz)),Tz = 8; end 061 if (nargin < 3|isempty(Hs)), Hs = 6; end 062 063 multipleSeaStates = any(prod(size(Hs))>1|prod(size(Tz))>1); 064 if multipleSeaStates 065 [errorcode, Scf,Hd,Tz,Hs] = comnsize(Scf,Hd,Tz,Hs); 066 else 067 [errorcode, Scf,Hd] = comnsize(Scf,Hd); 068 end 069 if errorcode > 0 070 error('Requires non-scalar arguments to match in size.'); 071 end 072 073 Hrms = Hs/sqrt(2); 074 Erms = 5/4*Hs./(Tz.^2); 075 %Erms = (0.0202 + 0.826*Hs./(Tz.^2)); 076 077 s = Scf./Erms; 078 hMax = 20; 079 h = min(Hd./Hrms,hMax); 080 081 eps2 = 1e-6; 082 083 084 p = zeros(size(Hd)); 085 %k0 = find((Hs<=(Tz-4)*13/6+4)); 086 upLimit = 6.5/1.4; 087 loLimit = 2.5/1.26;; 088 k0 = find((loLimit*sqrt(Hs)<Tz)); 089 if any(k0) 090 if multipleSeaStates 091 h = h(k0); 092 s = s(k0); 093 Hs = Hs(k0); 094 Tz = Tz(k0); 095 else 096 k0 = 1:prod(size(Hd)); 097 end 098 hlim = h; 099 100 101 h0 = 2.00528163239112; 102 normalizedInput = 1; 103 lowerTail = 0; 104 105 106 if 0 107 % This is a trick to get the html documentation correct. 108 k = b04pdf(1,1,2,3); 109 end 110 111 if (tail == lowerTail) 112 k = find(h>h0); 113 hlim(k) = h0; 114 p(k0) = gaussq('b04pdf',0,hlim,eps2/2,[],s,Hs,Tz,5,normalizedInput)... 115 + gaussq('b04pdf',hlim,h,eps2/2,[],s,Hs,Tz,5,normalizedInput); 116 else 117 k = find(h<h0); 118 hlim(k) = h0; 119 p(k0) = gaussq('b04pdf',h,hlim,eps2/2,[],s,Hs,Tz,7,normalizedInput)... 120 + gaussq('b04pdf',hlim,hMax,eps2/2,[],s,Hs,Tz,7,normalizedInput); 121 end 122 end 123 return 124 125
Comments or corrections to the WAFO group