JHSNLCDF Joint (Scf,Hd) CDF for non-linear waves with JONSWAP spectrum. CALL: f = jhsnlcdf(Hd,Scf,Hm0,Tp,Gamma,tail) f = CDF evaluated at (Scf,Hd) Hd = zero down crossing wave height [m] Scf = crest front steepness Hm0 = significant wave height [m] Tp = Spectral peak period [s] Gamma = Peakedness parameter of the JONSWAP spectrum tail = 1 if upper tail is calculated 0 if lower tail is calulated (default) JHSNLCDF approximates the joint CDF of (Scf, Hd), i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order Stokes waves with a JONSWAP spectral density. The empirical parameters of the model is fitted by least squares to simulated (Scf,Hd) data for 13 classes of GAMMA between 1 and 7. Between 47000 and 55000 zero-downcrossing waves were simulated for each class of GAMMA. JHSNLCDF is restricted to the following range for GAMMA: 1 <= GAMMA <= 7 Example: Hm0 = 6;Tp = 9; gam = 3.5; ec = 0.25; hc = 3; lowerTail = 0; upperTail = ~lowerTail jhsnlcdf(hc,ec,Hm0,Tp,gam) % Prob(Hd<Hc,Scf<ec) jhsnlcdf(hc,ec,Hm0,Tp,gam,upperTail) % Prob(Hd>Hc,Scf>ec) % Conditional probability of steep and high waves given seastates % i.e., Prob(Hd>hc,Scf>ec|Hs,Tp) upperTail = 1; Hs = linspace(2.5,11.5,10); Tp = linspace(4.5,19.5,16); [T,H] = meshgrid(Tp,Hs); p = jhsnlcdf(hc,ec,H,T,gam,upperTail); v = 10.^(-6:-1); contourf(Tp,Hs,log10(p),log10(v)) xlabel('Tp') ylabel('Hs') fcolorbar(log10(v)) See also jhscdf
Numerically evaluates a integral using a Gauss quadrature. | |
Peakedness factor Gamma given Hm0 and Tp for JONSWAP | |
Joint (Scf,Hd) PDF for nonlinear waves with a JONSWAP spectra. | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
Average or mean value. | |
Evaluate polynomial. | |
Display warning message; disable or enable warning messages. |
001 function f = jhsnlcdf(Hd,Scf,Hm0,Tp,gam,tail) 002 %JHSNLCDF Joint (Scf,Hd) CDF for non-linear waves with JONSWAP spectrum. 003 % 004 % CALL: f = jhsnlcdf(Hd,Scf,Hm0,Tp,Gamma,tail) 005 % 006 % f = CDF evaluated at (Scf,Hd) 007 % Hd = zero down crossing wave height [m] 008 % Scf = crest front steepness 009 % Hm0 = significant wave height [m] 010 % Tp = Spectral peak period [s] 011 % Gamma = Peakedness parameter of the JONSWAP spectrum 012 % tail = 1 if upper tail is calculated 013 % 0 if lower tail is calulated (default) 014 % 015 % JHSNLCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 016 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order Stokes 017 % waves with a JONSWAP spectral density. The empirical parameters of the 018 % model is fitted by least squares to simulated (Scf,Hd) data for 13 019 % classes of GAMMA between 1 and 7. Between 47000 and 55000 020 % zero-downcrossing waves were simulated for each class of GAMMA. 021 % JHSNLCDF is restricted to the following range for GAMMA: 022 % 1 <= GAMMA <= 7 023 % 024 % Example: 025 % Hm0 = 6;Tp = 9; gam = 3.5; 026 % ec = 0.25; 027 % hc = 3; 028 % lowerTail = 0; 029 % upperTail = ~lowerTail 030 % jhsnlcdf(hc,ec,Hm0,Tp,gam) % Prob(Hd<Hc,Scf<ec) 031 % jhsnlcdf(hc,ec,Hm0,Tp,gam,upperTail) % Prob(Hd>Hc,Scf>ec) 032 % 033 % % Conditional probability of steep and high waves given seastates 034 % % i.e., Prob(Hd>hc,Scf>ec|Hs,Tp) 035 % upperTail = 1; 036 % Hs = linspace(2.5,11.5,10); 037 % Tp = linspace(4.5,19.5,16); 038 % [T,H] = meshgrid(Tp,Hs); 039 % p = jhsnlcdf(hc,ec,H,T,gam,upperTail); 040 % v = 10.^(-6:-1); 041 % contourf(Tp,Hs,log10(p),log10(v)) 042 % xlabel('Tp') 043 % ylabel('Hs') 044 % fcolorbar(log10(v)) 045 % 046 % See also jhscdf 047 048 % Reference 049 % P. A. Brodtkorb (2004), 050 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 051 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 052 % Trondheim, Norway. 053 054 % History 055 % By pab 20.01.2003 056 057 058 error(nargchk(2,6,nargin)) 059 060 if (nargin < 6|isempty(tail)),tail = 0; end 061 if (nargin < 4|isempty(Tp)),Tp = 8; end 062 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end 063 if (nargin < 5|isempty(gam)), 064 gam = getjonswappeakedness(Hm0,Tp); 065 end 066 067 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 068 if multipleSeaStates 069 [errorcode, Scf,Hd,Hm0,Tp,gam] = comnsize(Scf,Hd,Hm0,Tp,gam); 070 else 071 [errorcode, Scf,Hd] = comnsize(Scf,Hd); 072 end 073 if errorcode > 0 074 error('Requires non-scalar arguments to match in size.'); 075 end 076 displayWarning = 0; 077 if displayWarning 078 if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0))) 079 disp('Warning: Hm0,Tp is outside the JONSWAP range') 080 disp('The validity of the parameters returned are questionable') 081 end 082 end 083 %dev = 2e-5; 084 c1 =[ 0.16183666835624 1.53691936441548 1.55852759524555]; 085 c2 =[ 0.15659478203944 1.15736959972513 1]; 086 if 0, 087 Tm02 = Tp 088 Tp = Tm02.*(polyval(c1,gam)./polyval(c2,gam)); 089 else 090 Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam)); 091 end 092 Hrms = Hm0/sqrt(2); 093 Vrms = 1.25*Hm0./(Tm02.^2); % Erms 094 095 v = Scf./Vrms; 096 097 hMax = 6; 098 eps2 = 1e-6; 099 normalizedInput = 1; 100 utprb = gaussq('jhsnlpdf',hMax,2*hMax,eps2/2,[],mean(v(:)),mean(Hm0(:)),mean(Tp(:)),mean(gam(:)),normalizedInput,7); 101 if eps2<utprb 102 warning('Check the accuracy of integration!') 103 end 104 105 h = min(Hd./Hrms,hMax); 106 107 f = zeros(size(Hd)); 108 % Only compute 109 if 0, % haver parametrization 110 loLimit = 3.6; 111 upLimit = 5; 112 else 113 loLimit = 2.5; 114 upLimit = 6.5; 115 end 116 117 k0 = find( (loLimit*sqrt(Hm0)<Tp) & (Tp<upLimit*sqrt(Hm0)) ); 118 if any(k0) 119 if multipleSeaStates 120 h = h(k0); 121 v = v(k0); 122 Hm0 = Hm0(k0); 123 Tp = Tp(k0); 124 gam = gam(k0); 125 else 126 k0 = 1:prod(size(Hd)); 127 end 128 else 129 return 130 end 131 132 if 0 133 % This is a trick to get the html documentation correct. 134 k = jhsnlpdf(1,1,2,3); 135 end 136 137 hlim = h; 138 139 lowerTail = 0; 140 if tail==lowerTail, 141 k = find(h>2.5); 142 hlim(k) = 2.5; 143 f(k0) = gaussq('jhsnlpdf',0,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5)... 144 + gaussq('jhsnlpdf',hlim,h,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5); 145 else % upper tail 146 k = find(h<2.5); 147 hlim(k) = 2.5; 148 f(k0) = gaussq('jhsnlpdf',h,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7)... 149 + gaussq('jhsnlpdf',hlim,hMax,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7); 150 end 151 return 152 153
Comments or corrections to the WAFO group