THSCDF Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra. CALL: f = thscdf(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) THSCDF approximates the joint CDF of (Scf, Hd), i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process 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. THSCDF 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 thscdf(hc,sc,Hm0,Tp) % Prob(Hd<Hc,Scf<Ec) thscdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>Hc,Scf>Ec) % 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 = thscdf(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 thspdf
Numerically evaluates a integral using a Gauss quadrature. | |
Joint (Scf,Hd) PDF for linear 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 = thscdf(Hd,Scf,Hm0,Tp,tail) 002 %THSCDF Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra. 003 % 004 % CALL: f = thscdf(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 % THSCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 015 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process 016 % with a Torsethaugen spectral density. The empirical parameters of the 017 % model is fitted by least squares to simulated (Scf,Hd) data for 600 018 % classes of Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves 019 % were simulated for each class of Hm0 and Tp. 020 % THSCDF is restricted to the following range for Hm0 and Tp: 021 % 0.5 < Hm0 [m] < 12, 3.5 < Tp [s] < 20, and Hm0 < (Tp-2)*12/11. 022 % 023 % Example: 024 % Hm0 = 6;Tp = 8; 025 % sc = 0.25; 026 % hc = 3; 027 % lowerTail = 0; 028 % upperTail = ~lowerTail 029 % thscdf(hc,sc,Hm0,Tp) % Prob(Hd<Hc,Scf<Ec) 030 % thscdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>Hc,Scf>Ec) 031 % 032 % % Conditional probability of steep and high waves given seastates 033 % % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz) 034 % upperTail = 1; 035 % Hs = linspace(2.5,11.5,10); 036 % Tp = linspace(4.5,19.5,16); 037 % [T,H] = meshgrid(Tp,Hs); 038 % p = thscdf(hc,sc,H,T,upperTail); 039 % v = 10.^(-6:-1); 040 % contourf(Tp,Hs,log10(p),log10(v)) 041 % xlabel('Tp') 042 % ylabel('Hs') 043 % fcolorbar(log10(v)) 044 % 045 % See also thspdf 046 047 % Reference 048 % P. A. Brodtkorb (2004), 049 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 050 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 051 % Trondheim, Norway. 052 053 % History 054 % revised pab 09.08.2003 055 % changed input 056 % validated 20.11.2002 057 % By pab 20.12.2000 058 059 060 error(nargchk(3,5,nargin)) 061 062 if (nargin < 5|isempty(tail)),tail = 0; end 063 if (nargin < 4|isempty(Tp)),Tp = 8; end 064 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end 065 066 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 067 if multipleSeaStates 068 [errorcode, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp); 069 else 070 [errorcode, Scf,Hd] = comnsize(Scf,Hd); 071 end 072 if errorcode > 0 073 error('Requires non-scalar arguments to match in size.'); 074 end 075 076 global THSPAR 077 if isempty(THSPAR) 078 %THSPAR = load('thspar.mat'); 079 THSPAR = load('thspar19-Jul-2004.mat'); 080 end 081 082 Tpp = THSPAR.Tp; 083 Hm00 = THSPAR.Hm0; 084 Tm020 = THSPAR.Tm02; 085 % Interpolation method 086 method = '*cubic';% Faster interpolation 087 088 [Tp1,Hs1] = meshgrid(Tpp,Hm00); 089 if 1, 090 Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 091 else 092 Tm02 = Tp; 093 for ix= 1:100 094 dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method)); 095 Tp = Tp+dTp; 096 if all(abs(dTp)<0.01) 097 %dTp 098 ix 099 break 100 end 101 end 102 end 103 % w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 104 %w = linspace(0,10,2*1024+1).'; 105 % St = torsethaugen(w,[Hm0,Tp]); 106 % ch = spec2char(St,{'Tm02','eps2'}); 107 % Tm02 = ch(1); 108 % eps2 = ch(2); 109 Hrms = Hm0/sqrt(2); 110 Erms = 1.25*Hm0./(Tm02.^2); % Erms 111 112 s = Scf./Erms; 113 hMax = 5; 114 h = min(Hd./Hrms,hMax); 115 116 eps2 = 1e-6; 117 118 f = zeros(size(Hd)); 119 120 % Only compute within valid range 121 %k0 = find((2<=Tp) & (Tp<=21) & (Hm0<=(Tp-2)*12/11) & (Hm0<=12)); 122 upLimit = 6.5; 123 loLimit = 2.5; 124 k0 = find((2<=Tp) & (Tp<=21) & (loLimit*sqrt(Hm0)<Tp) & (Hm0<=12)); 125 if any(k0) 126 if multipleSeaStates 127 h = h(k0); 128 s = s(k0); 129 Hm0 = Hm0(k0); 130 Tp = Tp(k0); 131 else 132 k0 = 1:prod(size(Hd)); 133 end 134 135 hlim = h; 136 137 normalizedInput = 1; 138 lowerTail = 0; 139 if 0 140 % This is a trick to get the html documentation correct. 141 k = thspdf(1,1,2,3); 142 end 143 144 if tail==lowerTail 145 k = find(h>2.5); 146 hlim(k) = 2.5; 147 f(k0) = gaussq('thspdf',0,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,5)... 148 + gaussq('thspdf',hlim,h,eps2/2,[],s,Hm0,Tp,normalizedInput,5); 149 else % upper tail 150 k = find(h<2.5); 151 hlim(k) = 2.5; 152 f(k0) = gaussq('thspdf',h,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,7)... 153 + gaussq('thspdf',hlim,hMax,eps2/2,[],s,Hm0,Tp,normalizedInput,7); 154 end 155 end 156 return 157 158
Comments or corrections to the WAFO group