THPDF Marginal wave height, Hd, PDF for Torsethaugen spectra. CALL: f = thpdf(h,Hm0,Tp,dim) f = pdf evaluated at h. h = vectors of evaluation points. Hm0 = significant wave height [m]. Tp = Spectral peak period [s]. dim = 'time' : Hd distribution in time (default) 'space' : Hd distribution in space THPDF approximates the marginal distribution of Hd, i.e., zero-downcrossing wave height, for a Gaussian process with a Torsethaugen spectral density. The empirical parameters of the model is fitted by least squares to simulated Hd data for 600 classes of Hm0 and Tp. Between 50000 and 150000 zero-downcrossing waves were simulated for each class of Hm0 and Tp. THPDF 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; h = linspace(0,4*Hm0/sqrt(2))'; f = thpdf(h,Hm0,Tp); dt = 0.4; w = linspace(0,2*pi/dt,256)'; S = torsethaugen(w,[Hm0 Tp]); xs = spec2sdat(S,20000,dt); rate=8; method=1; [S,H] = dat2steep(xs,rate,method); fk = kdebin(H,'epan',[],[],.5,128); plot(h,f), hold on, pdfplot(fk,'r'), hold off See also thvpdf
Wave height, Hd, distribution parameters for Torsethaugen spectra. | |
Wave height, Hd, distribution parameters for Torsethaugen spectra. | |
Generalized Gamma probability density function | |
Truncated Weibull probability density function | |
Display message and abort function. | |
Compare first N characters of strings ignoring case. |
001 function f = thpdf(h,Hm0,Tp,dim) 002 %THPDF Marginal wave height, Hd, PDF for Torsethaugen spectra. 003 % 004 % CALL: f = thpdf(h,Hm0,Tp,dim) 005 % 006 % f = pdf evaluated at h. 007 % h = vectors of evaluation points. 008 % Hm0 = significant wave height [m]. 009 % Tp = Spectral peak period [s]. 010 % dim = 'time' : Hd distribution in time (default) 011 % 'space' : Hd distribution in space 012 % 013 % THPDF approximates the marginal distribution of Hd, i.e., 014 % zero-downcrossing wave height, for a Gaussian process with a Torsethaugen 015 % spectral density. The empirical parameters of the model is fitted by 016 % least squares to simulated Hd data for 600 classes of Hm0 and 017 % Tp. Between 50000 and 150000 zero-downcrossing waves were simulated for 018 % each class of Hm0 and Tp. 019 % THPDF 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 % 022 % Example: 023 % Hm0 = 6;Tp = 8; 024 % h = linspace(0,4*Hm0/sqrt(2))'; 025 % f = thpdf(h,Hm0,Tp); 026 % dt = 0.4; w = linspace(0,2*pi/dt,256)'; 027 % S = torsethaugen(w,[Hm0 Tp]); 028 % xs = spec2sdat(S,20000,dt); rate=8; method=1; 029 % [S,H] = dat2steep(xs,rate,method); 030 % fk = kdebin(H,'epan',[],[],.5,128); 031 % plot(h,f), hold on, pdfplot(fk,'r'), hold off 032 % 033 % See also thvpdf 034 035 036 % Reference 037 % P. A. Brodtkorb (2004), 038 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 039 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 040 % Trondheim, Norway. 041 042 % History 043 % revised pab jan2004 044 % By pab 20.12.2000 045 046 047 048 error(nargchk(3,4,nargin)) 049 050 if nargin<4|isempty(dim), 051 dim = 'time';% dim='time'->wtweibpdf, dim='space'->wggampdf 052 end 053 054 if Hm0>12| Hm0>(Tp-2)*12/11 055 disp('Warning: Hm0 is outside the valid range') 056 disp('The validity of the Hd distribution is questionable') 057 end 058 if Tp>20|Tp<3 059 disp('Warning: Tp is outside the valid range') 060 disp('The validity of the Hd distribution is questionable') 061 end 062 Hrms = Hm0/sqrt(2); 063 064 [a b c] = thwparfun(Hm0,Tp,dim); 065 f = wtweibpdf(h/Hrms,a,b,c)/Hrms; 066 067 return 068 % old code kept just in case 069 if strncmpi(dim,'t',1) 070 [a b c] = thwparfun(Hm0,Tp,dim); 071 f = wtweibpdf(h/Hrms,a,b,c)/Hrms; 072 else 073 [a b c] = thgparfun(Hm0,Tp,dim); 074 f = wggampdf(h/Hrms,a,b,c)/Hrms; 075 end 076 077 return 078 079 % old code: saved just in case 080 % Weibull distribution parameters as a function of e2 and h2 081 % Hrms = Hm0/sqrt(2); 082 % if nargin<1 |isempty(h), h=linspace(0,4*Hrms).'; end 083 % if nargin>4, 084 % h = h*Hrms; 085 % if pdef==1 086 % f = hggam(h,Hm0,Tp,eps2,1); 087 % else 088 % f = htweib(h,Hm0,Tp,eps2,1); 089 % end 090 % elseif pdef==1 091 % f = hggam(h,Hm0,Tp,eps2); 092 % else 093 % f = htweib(h,Hm0,Tp,eps2); 094 % end 095 096 % f=f.f; 097 % return 098 099 % function f = hggam(h,Hm0,Tp,eps2,norm) 100 % pardef =2; 101 % if pardef == 1, 102 % if nargin<4|isempty(eps2), 103 % w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 104 % %w = linspace(0,1000/Tp,8*1024+1).'; % adaptiv spacing 105 % eps2 = spec2char(torsethaugen(w,[Hm0,Tp]),{'eps2'}); 106 % end 107 % if eps2<0.4 108 % disp('Warning: eps2 is less than 0.4. The pdf returned is questionable') 109 % elseif eps2>1.3 110 % disp('Warning: eps2 is larger than 1.3. The pdf returned is questionable') 111 % end 112 % end 113 114 % switch pardef 115 % case 1,% Wggampdf distribution parameters as a function of eps2 116 % % Best fit by smoothing spline forcing a=1, b=1 and c=2 for eps2=0. 117 % % Then approximate the spline with a rational polynomial 118 119 % da1 = [0.3579 0.9601 -1.8152 1.0009]; 120 % da2 = [1.5625 -0.1537 -1.1390 1.0000]; 121 % db1 = [-1.8555 3.5375 -3.4056 1.0026]; 122 % db2 = [-2.0855 4.1624 -3.6399 1.0000]; 123 % dc1 = [6.9572 -6.7932 2.1760 1.9974]; 124 % dc2 = [3.0979 -3.0781 0.5369 1.0000]; 125 % A0 = polyval(da1,eps2)./polyval(da2,eps2); 126 % B0 = polyval(db1,eps2)./polyval(db2,eps2); 127 % C0 = polyval(dc1,eps2)./polyval(dc2,eps2); 128 % case 2, % wave height distribution in time 129 % global THGPAR 130 % if isempty(THGPAR) 131 % THGPAR = load('thgpar.mat'); 132 % end 133 % % Generalized Gamma distribution parameters as a function of Tp, Hm0 134 % A00 = THGPAR.A00s; 135 % B00 = THGPAR.B00s; 136 % C00 = THGPAR.C00s; 137 138 % Tpp = THGPAR.Tp; 139 % Hm00 = THGPAR.Hm0; 140 % [E1, H1] = meshgrid(Tpp,Hm00); 141 % method = '*cubic'; 142 % A0 = interp2(E1,H1,A00,Tp,Hm0,method); 143 % B0 = interp2(E1,H1,B00,Tp,Hm0,method); 144 % C0 = interp2(E1,H1,C00,Tp,Hm0,method); 145 % case 3,% Waveheight distribution in space 146 % global THSSPAR 147 % if isempty(THSSPAR) 148 % THSSPAR = load('thsspar.mat'); 149 % end 150 % % Generalized Gamma distribution parameters as a function of Tp, Hm0 151 % A00 = THSSPAR.A00s; 152 % B00 = THSSPAR.B00s; 153 % C00 = THSSPAR.C00s; 154 155 % Tpp = THSSPAR.Tp; 156 % Hm00 = THSSPAR.Hm0; 157 % [E1, H1] = meshgrid(Tpp,Hm00); 158 % method = '*cubic'; 159 % A0 = interp2(E1,H1,A00,Tp,Hm0,method); 160 % B0 = interp2(E1,H1,B00,Tp,Hm0,method); 161 % C0 = interp2(E1,H1,C00,Tp,Hm0,method); 162 % end 163 % Hrms = Hm0/sqrt(2); 164 % f.f = wggampdf(h/Hrms,A0,B0,C0)/Hrms; 165 % return 166 167 % function f = htweib(h,Hm0,Tp,eps2,norm) 168 % pardef =7; 169 % if pardef < 7, 170 % if nargin<4|isempty(eps2), 171 % w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 172 % %w = linspace(0,1000/Tp,8*1024+1).'; % adaptiv spacing 173 % S = torsethaugen(w,[Hm0,Tp]); 174 % eps2 = spec2char(S,{'eps2'}); 175 % end 176 % if eps2<0.4 177 % disp('Warning: eps2 is less than 0.4. The pdf returned is questionable') 178 % elseif eps2>1.3 179 % disp('Warning: eps2 is larger than 1.3. The pdf returned is questionable') 180 % end 181 % end 182 % switch pardef, 183 % case 1,% Wtweibpdf distribution parameters as a function of eps2 184 % % Best fit by smoothing spline forcing a=1, b=2 and c=0 for eps2=0. 185 % brks = [0 .1 .2 .4 ,.6, .8, 1, 1.1 1.2]'; 186 % coefa = [0 0 0.02260415153596 0.99807186986167; ... 187 % 2.19065400617385 0 0.02260415153596 1.00033228501527; ... 188 % 4.34015195156053 0.65719620185215 0.03709199393185 1.00478335417504; ... 189 % -1.59533089716870 3.26128737278847 0.80983543882910 1.07321081664798; ... 190 % -6.81273221810880 2.30408883448726 1.92291068028425 1.35286675214799; ... 191 % -3.69498826658975 -1.78355049637802 2.09369407217829 1.77511058383946; ... 192 % 13.33514485443956 -4.00054345633187 0.94471547491809 2.09294747228728; ... 193 % 0 0 0.54466112928490 2.16074873007021]; 194 195 % coefb = [ 0 0 0.32503235228616 1.99054481866418; ... 196 % 3.28321899128157 0 0.32503235228616 2.02304805389280; ... 197 % 5.67672309005450 0.98496569738447 0.37964649056830 2.05883450811270; ... 198 % -5.29907238080822 4.39099955141717 1.43842344537222 2.21957621884217; .... 199 % -5.89663569823287 1.21155612293224 2.55893458024211 2.64050831092684; ... 200 % -6.21824739906323 -2.32642529600749 2.43691697455115 3.15358438630669; ... 201 % 20.19124578481806 -6.05737373544542 0.77134599291113 3.49816479018411; ... 202 % 0 0 0.16560861936659 3.53491689790559]; 203 % coefc =[ 0 0 0.04818579357214 -0.00817761487085; ... 204 % 2.94432030165157 0 0.04818579357214 -0.00335903551363; ... 205 % 4.77660844045250 0.88329609049547 0.09917317190900 0.00440386414523; ... 206 % -1.24578770271258 3.74926115476697 1.01096301945323 0.09778320967047; ... 207 % -7.70868155645400 3.00178853313943 2.36117295703451 0.43997995813009; ... 208 % -3.98346578867600 -1.62342040073298 2.70373824808144 0.97061663841094; ... 209 % 13.37833291312857 -4.01349987393858 1.57933014446810 1.41455974568850; ... 210 % 0 0 1.17798015707424 1.54573609430905]; 211 % pa = mkpp(brks,coefa); 212 % pb = mkpp(brks,coefb); 213 % pc = mkpp(brks,coefc); 214 % A0 = ppval(pa,eps2); 215 % B0 = ppval(pb,eps2); 216 % C0 = ppval(pc,eps2); 217 % case 2,% Wtweibpdf distribution parameters as a function of eps2 218 % % Best fit by smoothing spline forcing a=1, b=2 and c=0 for eps2=0. 219 % % Then approximate the spline with a rational polynomial. 220 221 % da1 = [1.8262 -2.1141 1.0034]; 222 % da2 = [-0.1371 0.1626 1.3189 -2.0016 1.0000]; 223 % db1 = [2.4360 -1.6331 1.9932]; 224 % db2 = [-1.3266 4.3755 -4.6795 2.4763 -1.0450 1.0000]; 225 % dc1 = [0.3162 -0.1647 0.0803 -0.0104]; 226 % dc2 = [-1.4325 6.3773 -11.1967 10.1490 -4.7401 1.0000]; 227 % A0 = polyval(da1,eps2)./polyval(da2,eps2); 228 % B0 = polyval(db1,eps2)./polyval(db2,eps2); 229 % C0 = polyval(dc1,eps2)./polyval(dc2,eps2); 230 % case 3, % Wtweibpdf distribution parameters as a function of eps2 231 % % LS fit to data 232 % % best fit to torsethaugen for eps2 = 0.4 to 1.2; 233 % % and forcing through A0=1 B0=2 C0=0 for eps2==0) 234 % A0 = sqrt(1+2.0094*eps2.^1.8492); 235 % B0 = 1+sqrt(1+3.4117*eps2.^1.4972); 236 % C0 = 0.9994*eps2.^1.6728; 237 % case 4,% Wtweibpdf distribution parameters as a function of eps2 238 % % best fit to torsethaugen eps2 = 0.4 to 1.2; 239 % A0 = 0.7315+eps2.^0.9916; 240 % B0 = 2.0628+eps2.^1.1927; 241 % C0 = 0.9994*eps2.^1.6728; 242 % case 5,% Naess (1985) 243 % R = spec2cov(S); 244 % % A0 = sqrt((1-min(R.R)/R.R(1))/2);% Naess (1985) 245 % A0 = sqrt((1-min(R.R)/R.R(1))/2)+0.03;% Modified approach broadbanded time 246 % % A0 = sqrt((1-min(R.R)/R.R(1))/2)+0.1;% Modified approach broadbanded space 247 248 % B0 = 2; 249 % C0 = 0; 250 % case 6,% Naess (1985) wave height in space 251 % R = spec2cov(spec2spec(specinterp(S,0.55),'k1d')); 252 % %A0 = sqrt((1-min(R.R)/R.R(1))/2); % Naess (1985) 253 % A0 = sqrt((1-min(R.R)/R.R(1))/2)+0.03; % modified approach 254 % B0 = 2; 255 % C0 = 0; 256 % case 7,% wave height distribution in time 257 % global THWPAR 258 % if isempty(THWPAR) 259 % THWPAR = load('thwpar.mat'); 260 % end 261 % % Truncated Weibull distribution parameters as a function of Tp, Hm0 262 % A00 = THWPAR.A00s; 263 % B00 = THWPAR.B00s; 264 % C00 = THWPAR.C00s; 265 266 % Tpp = THWPAR.Tp; 267 % Hm00 = THWPAR.Hm0; 268 % [E1, H1] = meshgrid(Tpp,Hm00); 269 % method = '*cubic'; 270 % A0 = interp2(E1,H1,A00,Tp,Hm0,method); 271 % B0 = interp2(E1,H1,B00,Tp,Hm0,method); 272 % C0 = interp2(E1,H1,C00,Tp,Hm0,method); 273 274 % end 275 276 % Hrms = Hm0/sqrt(2); 277 % f = wtweibpdf(h/Hrms,A0,B0,C0)/Hrms; 278 279 % return 280
Comments or corrections to the WAFO group