THSSCDF Joint (Scf,Hd) CDF for linear waves in space with Torsethaugen spectra. CALL: f = thsscdf(Hd,Scf,Hm0,Tp,tail) f = CDF evaluated at (Scf,Hd) Hd = zero down crossing wave height Scf = crest front steepness Hm0 = significant wave height [m] Tp = Spectral peak period [s] tail = 1 if upper tail is calculated 0 if lower tail is calulated (default) THSSCDF approximates the joint CDF of (Scf, Hd), i.e., crest front steepness (Ac/Lcf) and wave height in space, 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 100000 and 1000000 zero-downcrossing waves were simulated for each class of Hm0 and Tp. THSSCDF 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; Ec = 0.25; Hc = 3; lowerTail = 0; upperTail = ~lowerTail thsscdf(Hc,Ec,Hm0,Tp) % Prob(Hd<Hc,Scf<Ec) thsscdf(Hc,Ec,Hm0,Tp,upperTail) % Prob(Hd>Hc,Scf>Ec) See also thscdf
Numerically evaluates a integral using a Gauss quadrature. | |
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 = thsscdf(Hd,Scf,Hm0,Tp,tail) 002 %THSSCDF Joint (Scf,Hd) CDF for linear waves in space with Torsethaugen spectra. 003 % 004 % CALL: f = thsscdf(Hd,Scf,Hm0,Tp,tail) 005 % 006 % f = CDF evaluated at (Scf,Hd) 007 % Hd = zero down crossing wave height 008 % Scf = crest front steepness 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 % THSSCDF approximates the joint CDF of (Scf, Hd), i.e., crest 015 % front steepness (Ac/Lcf) and wave height in space, for a Gaussian 016 % process with a Torsethaugen spectral density. The empirical parameters 017 % of the model is fitted by least squares to simulated (Scf,Hd) data for 018 % 600 classes of Hm0 and Tp. Between 100000 and 1000000 zero-downcrossing 019 % waves were simulated for each class of Hm0 and Tp. 020 % THSSCDF 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 % Ec = 0.25; 026 % Hc = 3; 027 % lowerTail = 0; 028 % upperTail = ~lowerTail 029 % thsscdf(Hc,Ec,Hm0,Tp) % Prob(Hd<Hc,Scf<Ec) 030 % thsscdf(Hc,Ec,Hm0,Tp,upperTail) % Prob(Hd>Hc,Scf>Ec) 031 % 032 % See also thscdf 033 034 % Reference 035 % P. A. Brodtkorb (2004), 036 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 037 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 038 % Trondheim, Norway. 039 040 % History 041 % revised pab 14.03.2004 042 % revised pab 09.08.2003 043 % changed input 044 % validated 20.11.2002 045 % By pab 20.12.2000 046 047 048 error(nargchk(3,5,nargin)) 049 050 if (nargin < 5|isempty(tail)),tail = 0; end 051 if (nargin < 4|isempty(Tp)),Tp = 8; end 052 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end 053 054 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 055 if multipleSeaStates 056 [errorcode, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp); 057 else 058 [errorcode, Scf,Hd] = comnsize(Scf,Hd); 059 end 060 if errorcode > 0 061 error('Requires non-scalar arguments to match in size.'); 062 end 063 useWeibull = 1; 064 if useWeibull 065 global THSSPARW 066 if isempty(THSSPARW) 067 THSSPARW= load('thsspar27-Jul-2004.mat'); 068 end 069 Tpp = THSSPARW.Tp; 070 Hm00 = THSSPARW.Hm0; 071 Tm020 = THSSPARW.Tm02; 072 else 073 global THSSPARG 074 if isempty(THSSPARG) 075 THSSPARG = load('thsspar.mat'); 076 end 077 078 Tpp = THSSPARG.Tp; 079 Hm00 = THSSPARG.Hm0; 080 Tm020 = THSSPARG.Tm02; 081 end 082 % Interpolation method 083 method = '*cubic';% Faster interpolation 084 085 [Tp1,Hs1] = meshgrid(Tpp,Hm00); 086 Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 087 % w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 088 %w = linspace(0,10,2*1024+1).'; 089 % St = torsethaugen(w,[Hm0,Tp]); 090 % ch = spec2char(St,{'Tm02','eps2'}); 091 % Tm02 = ch(1); 092 % eps2 = ch(2); 093 Hrms = Hm0/sqrt(2); 094 Erms = 2*Hm0./(Tm02); % Srms 095 096 s = Scf./Erms; 097 hMax = 10; 098 h = min(Hd./Hrms,hMax); 099 100 eps2 = 1e-6; 101 102 f = zeros(size(Hd)); 103 104 % Only compute within valid range 105 %k0 = find((2<=Tp) & (Tp<=21) & (Hm0<=(Tp-2)*12/11) & (Hm0<=12)); 106 upLimit = 6.5; 107 loLimit = 2.5; 108 k0 = find((2<=Tp) & (Tp<=21) & (loLimit*sqrt(Hm0)<Tp) & (Hm0<=12)); 109 if any(k0) 110 if multipleSeaStates 111 h = h(k0); 112 s = s(k0); 113 Hm0 = Hm0(k0); 114 Tp = Tp(k0); 115 else 116 k0 = 1:prod(size(Hd)); 117 end 118 hlim = h; 119 120 normalizedInput = 1; 121 lowerTail = 0; 122 if tail==lowerTail 123 k = find(h>2.5); 124 hlim(k) = 2.5; 125 f(k0) = gaussq('thsspdf',0,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,5)... 126 + gaussq('thsspdf',hlim,h,eps2/2,[],s,Hm0,Tp,normalizedInput,5); 127 else % upper tail 128 k = find(h<2.5); 129 hlim(k) = 2.5; 130 f(k0) = gaussq('thsspdf',h,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,7)... 131 + gaussq('thsspdf',hlim,hMax,eps2/2,[],s,Hm0,Tp,normalizedInput,7); 132 end 133 end 134 return 135 136
Comments or corrections to the WAFO group