THSNLPDF Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. CALL: f = thsnlpdf(Hd,Scf,Hm0,Tp) f = pdf Hd = zero down crossing wave height Scf = crest front steepness Hm0 = significant wave height Tp = Spectral peak period THSNLPDF approximates the joint PDF of (Scf, Hd), i.e., crest steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for for 2nd order non-linear 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. THSNLPDF 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. The size of f is the common size of the input arguments. Example: Hm0 = 6;Tp = 8; h = linspace(0,4*Hm0/sqrt(2)); s = linspace(0,6*1.25*Hm0/Tp^2,101); [S,H] = meshgrid(s,h); f = thsnlpdf(H,S,Hm0,Tp); contourf(s,h,f) See also thspdf2, thsspdf
Calculates a smoothing spline. | |
Wave height, Hd, distribution parameters for Torsethaugen spectra. | |
Gamma cumulative distribution function | |
Gamma probability density function | |
Truncated Weibull probability density function | |
Check if all input arguments are either scalar or of common size. | |
Display message and abort function. | |
2-D interpolation (table lookup). | |
3-D interpolation (table lookup). | |
X and Y arrays for 3-D plots. | |
Set unique. |
Joint (Scf,Hd) CDF for nonlinear waves with Torsethaugen spectra. | |
Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. |
001 function [f,Hrms,Vrms] = thsnlpdf(Hd,Scf,Hm0,Tp,normalizedInput,condon) 002 %THSNLPDF Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. 003 % 004 % CALL: f = thsnlpdf(Hd,Scf,Hm0,Tp) 005 % 006 % f = pdf 007 % Hd = zero down crossing wave height 008 % Scf = crest front steepness 009 % Hm0 = significant wave height 010 % Tp = Spectral peak period 011 % 012 % THSNLPDF approximates the joint PDF of (Scf, Hd), i.e., crest 013 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for for 2nd order 014 % non-linear waves with a Torsethaugen spectral density. The empirical 015 % parameters of the model is 016 % fitted by least squares to simulated (Scf,Hd) data for 600 classes of 017 % Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves were 018 % simulated for each class of Hm0 and Tp. 019 % THSNLPDF is restricted to the following range for Hm0 and Tp: 020 % 0.5 < Hm0 [m] < 12, 3.5 < Tp [s] < 20, and Hm0 < (Tp-2)*12/11. 021 % The size of f is the common size of the input arguments. 022 % 023 % Example: 024 % Hm0 = 6;Tp = 8; 025 % h = linspace(0,4*Hm0/sqrt(2)); 026 % s = linspace(0,6*1.25*Hm0/Tp^2,101); 027 % [S,H] = meshgrid(s,h); 028 % f = thsnlpdf(H,S,Hm0,Tp); 029 % contourf(s,h,f) 030 % 031 % See also thspdf2, thsspdf 032 033 % Reference 034 % P. A. Brodtkorb (2004), 035 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 036 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 037 % Trondheim, Norway. 038 039 040 % History 041 % revised pab 09.08.2003 042 % changed input and help header 043 % validated 20.11.2002 044 % By pab 20.12.2000 045 046 error(nargchk(3,6,nargin)) 047 048 049 if (nargin < 6|isempty(condon)), condon = 0; end 050 if (nargin < 5|isempty(normalizedInput)), normalizedInput = 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 displayWarning = 0; 064 if displayWarning, 065 if any(Hm0>11| Hm0>(Tp-2)*12/11) 066 disp('Warning: Hm0 is outside the valid range') 067 disp('The validity of the Joint (Hd,Scf) distribution is questionable') 068 end 069 if any(Tp>20|Tp<3) 070 disp('Warning: Tp is outside the valid range') 071 disp('The validity of the Joint (Hd,Scf) distribution is questionable') 072 end 073 end 074 075 global THSNLPAR 076 if isempty(THSNLPAR) 077 %THSNLPAR = load('thsnlpar.mat'); 078 THSNLPAR = load('thsnlpar06-Aug-2004.mat') 079 end 080 081 Tpp = THSNLPAR.Tp; 082 Hm00 = THSNLPAR.Hm0; 083 Tm020 = THSNLPAR.Tm02; 084 h2 = THSNLPAR.h2; 085 % Interpolation method 086 method = '*cubic';% Faster interpolation 087 088 if normalizedInput, 089 Hrms = 1; 090 Vrms = 1; 091 else 092 [Tp1,Hs1] = meshgrid(Tpp,Hm00); 093 094 if 1, 095 Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 096 %Tp/Tm02 097 else 098 099 Tm02 = Tp; 100 for ix= 1:100 101 dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method)); 102 Tp = Tp+dTp; 103 if all(abs(dTp)<0.01) 104 dTp 105 ix 106 break 107 end 108 end 109 %Tp, Tp/Tm02 110 end 111 112 % w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 113 %w = linspace(0,10,2*1024+1).'; 114 % St = torsethaugen(w,[Hm0,Tp]); 115 % ch = spec2char(St,{'Tm02','eps2'}); 116 % Tm02 = ch(1); 117 % eps2 = ch(2); 118 Hrms = Hm0/sqrt(2); 119 Vrms = 1.25*Hm0./(Tm02.^2); % Erms 120 end 121 122 h = Hd./Hrms; 123 v = Scf./Vrms; 124 cSize = size(h); % common size 125 126 % Gamma distribution parameters as a function of Tp Hm0 and h2 127 A11 = THSNLPAR.A11s; 128 B11 = THSNLPAR.B11s; 129 130 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2); 131 Nh2 = length(h2); 132 if multipleSeaStates 133 h = h(:); 134 v = v(:); 135 Tp = Tp(:); 136 Hm0 = Hm0(:); 137 A1 = zeros(length(h),1); 138 B1 = A1; 139 [TpHm0,ix,jx] = unique([Tp,Hm0],'rows'); 140 numSeaStates = length(ix); 141 Tpi = zeros(Nh2,1); 142 Hm0i = zeros(Nh2,1); 143 for iz=1:numSeaStates 144 k = find(jx==iz); 145 Tpi(:) = TpHm0(iz,1); 146 Hm0i(:) = TpHm0(iz,2); 147 A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),... 148 1,h(k),1)); 149 B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),... 150 1,h(k),1)); 151 end 152 else 153 Tpi = repmat(Tp,[Nh2,1]); 154 Hm0i = repmat(Hm0,[Nh2,1]); 155 A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1)); 156 B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1)); 157 end 158 % Waveheight distribution in time 159 % Truncated Weibull distribution parameters as a function of Tp, Hm0 160 [A0, B0, C0] = thwparfun(Hm0,Tp,'time'); 161 switch condon, 162 case 0, % regular pdf is returned 163 f = wtweibpdf(h,A0,B0,C0).*wgampdf(v,A1,B1); 164 case 1, %pdf conditioned on x1 ie. p(x2|x1) 165 f = wgampdf(v,A1,B1); 166 case 3, % secret option used by XXstat: returns x2*p(x2|x1) 167 f = v.*wgampdf(v,A1,B1); 168 case 4, % secret option used by XXstat: returns x2.^2*p(x2|x1) 169 f = v.^2.*wgampdf(v,A1,B1); 170 case 5, % p(h)*P(V|h) is returned special case used by thscdf 171 f = wtweibpdf(h,A0,B0,C0).*wgamcdf(v,A1,B1); 172 case 6, % P(V|h) is returned special case used by thscdf 173 f = wgamcdf(v,A1,B1); 174 case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf 175 f = wtweibpdf(h,A0,B0,C0).*(1-wgamcdf(v,A1,B1)); 176 otherwise error('unknown option') 177 end 178 if multipleSeaStates 179 f = reshape(f,cSize); 180 end 181 182 if condon~=6 183 f = f./Hrms./Vrms; 184 end 185 f(find(isnan(f)|isinf(f) ))=0; 186 if any(size(f)~=cSize) 187 disp('Wrong size') 188 end 189 190 return 191 192 193 194 195
Comments or corrections to the WAFO group