MK87CDF Myrhaug and Kjeldsen (1987) joint (Scf,Hd) CDF. CALL: F = mk87cdf(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) MK87CDF 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 = 0.715*Hs; Erms = 0.0202+0.826*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 = mk87cdf(hc,sc,Hs,Tz) % Prob(Hd<=hc,Scf<=sc|Hs,Tz) = 0.59 % 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 = mk87cdf(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 mk87pdf, gaussq
Numerically evaluates a integral using a Gauss quadrature. | |
Myrhaug and Kjeldsen (1987) joint (Scf,Hd) PDF. | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. |
001 function [p, eps2, a, b] = mk87cdf(Hd,Scf,Hs,Tz,tail) 002 %MK87CDF Myrhaug and Kjeldsen (1987) joint (Scf,Hd) CDF. 003 % 004 % CALL: F = mk87cdf(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 % MK87CDF 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 = 0.715*Hs; 020 % Erms = 0.0202+0.826*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 = mk87cdf(hc,sc,Hs,Tz) % Prob(Hd<=hc,Scf<=sc|Hs,Tz) = 0.59 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 = mk87cdf(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 mk87pdf, gaussq 045 046 % References: 047 % Myrhaug, D. and Kjelsen S.P. (1987) 048 % 'Prediction of occurences of steep and high waves in deep water'. 049 % Journal of waterway, Port, Coastal and Ocean Engineers, 050 % Vol. 113, pp 122--138 051 % 052 % Myrhaug & Dahle (1984) Parametric modelling of joint probability 053 % density distributions for steepness and asymmetry in deep water 054 % waves 055 % 056 057 %tested on matlab 5.2 058 %history: 059 % revised pab 09.08.2003 060 % Changed input + updated help header 061 % revised pab 04.11.2000 062 % - no dependency on stats toolbox anymore. 063 % by Per A. brodtkorb 1998 064 065 066 error(nargchk(3,5,nargin)) 067 068 if (nargin < 5|isempty(tail)),tail = 0; end 069 if (nargin < 4|isempty(Tz)),Tz = 8; end 070 if (nargin < 3|isempty(Hs)), Hs = 6; end 071 072 multipleSeaStates = any(prod(size(Hs))>1|prod(size(Tz))>1); 073 if multipleSeaStates 074 [errorcode, Scf,Hd,Tz,Hs] = comnsize(Scf,Hd,Tz,Hs); 075 else 076 [errorcode, Scf,Hd] = comnsize(Scf,Hd); 077 end 078 if errorcode > 0 079 error('Requires non-scalar arguments to match in size.'); 080 end 081 082 Hrms = 0.715*Hs; 083 Erms = 0.0202 + 0.826*Hs./(Tz.^2); 084 085 s = Scf./Erms; 086 hMax = 20; 087 h = min(Hd./Hrms,hMax); 088 089 eps2 = 1e-6; 090 091 092 p = zeros(size(Hd)); 093 %k0 = find((Hs<=(Tz-4)*13/6+4)); 094 upLimit = 6.5/1.4; 095 loLimit = 2.5/1.26;; 096 k0 = find((loLimit*sqrt(Hs)<Tz)); 097 if any(k0) 098 if multipleSeaStates 099 h = h(k0); 100 s = s(k0); 101 Hs = Hs(k0); 102 Tz = Tz(k0); 103 else 104 k0 = 1:prod(size(Hd)); 105 end 106 hlim = h; 107 108 normalizedInput = 1; 109 lowerTail = 0; 110 111 112 113 if 0 114 % This is a trick to get the html documentation correct. 115 k = mk87pdf(1,1,2,3); 116 end 117 118 if (tail == lowerTail) 119 k = find(h>2.5); 120 hlim(k) = 2.5; 121 p(k0) = gaussq('mk87pdf',0,hlim,eps2/2,[],s,Hs,Tz,5,normalizedInput)... 122 + gaussq('mk87pdf',hlim,h,eps2/2,[],s,Hs,Tz,5,normalizedInput); 123 else 124 k = find(h<2.5); 125 hlim(k) = 2.5; 126 p(k0) = gaussq('mk87pdf',h,hlim,eps2/2,[],s,Hs,Tz,7,normalizedInput)... 127 + gaussq('mk87pdf',hlim,hMax,eps2/2,[],s,Hs,Tz,7,normalizedInput); 128 end 129 end 130 return 131 132
Comments or corrections to the WAFO group