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