THSNLCDF Joint (Scf,Hd) CDF for nonlinear waves with Torsethaugen spectra. CALL: f = thsnlcdf(Hd,Scf,Hm0,Tp,tail) f = CDF evaluated at (Scf,Hd) Hd = zero down crossing wave height (vector) Scf = crest front steepness (vector) Hm0 = significant wave height [m] Tp = Spectral peak period [s] tail = 1 if upper tail is calculated 0 if lower tail is calulated (default) THSNLCDF approximates the joint CDF of (Scf, Hd), i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order nonlinear waves with a Torsethaugen spectral density. The empirical parameters of the model is fitted by least squares to simulated (Scf,Hd) data for 600 classes of Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves were simulated for each class of Hm0 and Tp. THSNLCDF is restricted to the following range for Hm0 and Tp: 0.5 < Hm0 [m] < 12, 3.5 < Tp [s] < 20, and Hm0 < (Tp-2)*12/11. Example: Hm0 = 6;Tp = 8; sc = 0.25; hc = 3; lowerTail = 0; upperTail = ~lowerTail thsnlcdf(hc,sc,Hm0,Tp) % Prob(Hd<hc,Scf<sc) thsnlcdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>hc,Scf>sc) % 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,10); Tp = linspace(4.5,19.5,16); [T,H] = meshgrid(Tp,Hs); p = thsnlcdf(hc,sc,H,T,upperTail); v = 10.^(-6:-1); contourf(Tp,Hs,log10(p),log10(v)) xlabel('Tp') ylabel('Hs') fcolorbar(log10(v)) See also thsnlpdf
Numerically evaluates a integral using a Gauss quadrature. | |
Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
2-D interpolation (table lookup). | |
X and Y arrays for 3-D plots. |
001 function f = thsnlcdf(Hd,Scf,Hm0,Tp,tail) 002 %THSNLCDF Joint (Scf,Hd) CDF for nonlinear waves with Torsethaugen spectra. 003 % 004 % CALL: f = thsnlcdf(Hd,Scf,Hm0,Tp,tail) 005 % 006 % f = CDF evaluated at (Scf,Hd) 007 % Hd = zero down crossing wave height (vector) 008 % Scf = crest front steepness (vector) 009 % Hm0 = significant wave height [m] 010 % Tp = Spectral peak period [s] 011 % tail = 1 if upper tail is calculated 012 % 0 if lower tail is calulated (default) 013 % 014 % THSNLCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 015 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order 016 % nonlinear waves with a Torsethaugen spectral density. 017 % The empirical parameters of the 018 % model is fitted by least squares to simulated (Scf,Hd) data for 600 019 % classes of Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves 020 % were simulated for each class of Hm0 and Tp. 021 % THSNLCDF is restricted to the following range for Hm0 and Tp: 022 % 0.5 < Hm0 [m] < 12, 3.5 < Tp [s] < 20, and Hm0 < (Tp-2)*12/11. 023 % 024 % Example: 025 % Hm0 = 6;Tp = 8; 026 % sc = 0.25; 027 % hc = 3; 028 % lowerTail = 0; 029 % upperTail = ~lowerTail 030 % thsnlcdf(hc,sc,Hm0,Tp) % Prob(Hd<hc,Scf<sc) 031 % thsnlcdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>hc,Scf>sc) 032 % 033 % % Conditional probability of steep and high waves given seastates 034 % % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz) 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 = thsnlcdf(hc,sc,H,T,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 thsnlpdf 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 % revised pab 09.08.2003 056 % changed input 057 % validated 20.11.2002 058 % By pab 20.12.2000 059 060 061 error(nargchk(3,5,nargin)) 062 063 if (nargin < 5|isempty(tail)),tail = 0; end 064 if (nargin < 4|isempty(Tp)),Tp = 8; end 065 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end 066 067 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 068 if multipleSeaStates 069 [errorcode, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp); 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 077 global THSNLPAR 078 if isempty(THSNLPAR) 079 %THSNLPAR = load('thsnlpar.mat'); 080 THSNLPAR = load('thsnlpar06-Aug-2004.mat') 081 end 082 083 Tpp = THSNLPAR.Tp; 084 Hm00 = THSNLPAR.Hm0; 085 Tm020 = THSNLPAR.Tm02; 086 % Interpolation method 087 method = '*cubic';% Faster interpolation 088 089 [Tp1,Hs1] = meshgrid(Tpp,Hm00); 090 if 1, 091 Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 092 else 093 094 Tm02 = Tp; 095 for ix= 1:100 096 dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method)); 097 Tp = Tp+dTp; 098 if all(abs(dTp)<0.01) 099 %dTp 100 ix 101 break 102 end 103 end 104 end 105 106 107 % w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 108 %w = linspace(0,10,2*1024+1).'; 109 % St = torsethaugen(w,[Hm0,Tp]); 110 % ch = spec2char(St,{'Tm02','eps2'}); 111 % Tm02 = ch(1); 112 % eps2 = ch(2); 113 Hrms = Hm0/sqrt(2); 114 Erms = 1.25*Hm0./(Tm02.^2); % Erms 115 116 s = Scf./Erms; 117 hMax = 5; 118 h = min(Hd./Hrms,hMax); 119 120 eps2 = 1e-6; 121 122 f = zeros(size(Hd)); 123 124 % Only compute within valid range 125 %k0 = find((2<=Tp) & (Tp<=21) & (Hm0<=(Tp-2)*12/11) & (Hm0<=12)); 126 upLimit = 6.5; 127 loLimit = 2.5; 128 k0 = find((2<=Tp) & (Tp<=21) & (loLimit*sqrt(Hm0)<Tp) & (Hm0<=12)); 129 if any(k0) 130 if multipleSeaStates 131 h = h(k0); 132 s = s(k0); 133 Hm0 = Hm0(k0); 134 Tp = Tp(k0); 135 else 136 k0 = 1:prod(size(Hd)); 137 end 138 139 hlim = h; 140 141 normalizedInput = 1; 142 lowerTail = 0; 143 144 if 0 145 % This is a trick to get the html documentation correct. 146 k = thsnlpdf(1,1,2,3); 147 end 148 149 if tail==lowerTail 150 k = find(h>2.5); 151 hlim(k) = 2.5; 152 f(k0) = gaussq('thsnlpdf',0,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,5)... 153 + gaussq('thsnlpdf',hlim,h,eps2/2,[],s,Hm0,Tp,normalizedInput,5); 154 else % upper tail 155 k = find(h<2.5); 156 hlim(k) = 2.5; 157 f(k0) = gaussq('thsnlpdf',h,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,7)... 158 + gaussq('thsnlpdf',hlim,hMax,eps2/2,[],s,Hm0,Tp,normalizedInput,7); 159 end 160 end 161 return 162 163
Comments or corrections to the WAFO group