THSPDF2 Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. CALL: f = thspdf2(Hd,Scf,Hm0,Tp) f = pdf struct evaluated at meshgrid(Scf,Hd) Hd = zero down crossing wave height (vector) Scf = crest front steepness (vector) Hm0 = significant wave height [m] Tp = Spectral peak period [s] THSPDF2 approximates the joint distribution of (Scf, Hd), i.e., crest 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. THSPDF 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)); s = linspace(0,6*1.25*Hm0/Tp^2); f = thspdf2(h,s,Hm0,Tp); w = linspace(0,40,5*1024+1).'; S = torsethaugen(w,[Hm0 Tp]); dt = 0.3; x = spec2sdat(S,80000,.2); rate = 8; [si,hi] = dat2steep(x,rate,2); fk = kdebin([si,hi],'epan',[],[],.5,128); fk.title = f.title; fk.labx = f.labx; plot(si,hi,'.'), hold on pdfplot(f),pdfplot(fk,'r'),hold off See also thsspdf
PDF class constructor | |
Calculates quantile levels which encloses P% of PDF | |
Calculates the difference between the maximum and minimum values. | |
Calculates a smoothing spline. | |
Evaluates spectral characteristics and their covariance | |
Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. | |
Wave height, Hd, distribution parameters for Torsethaugen spectra. | |
Calculates a double peaked (swell + wind) spectrum | |
Gamma probability density function | |
Truncated Weibull probability density function | |
Display message and abort function. | |
2-D interpolation (table lookup). | |
3-D interpolation (table lookup). | |
Linearly spaced vector. | |
X and Y arrays for 3-D plots. | |
Make piecewise polynomial. | |
Convert number to string. (Fast version) | |
Evaluate piecewise polynomial. |
001 function [f,varargout] = thspdf2(Hd,Scf,Hm0,Tp,normalizedInput) 002 %THSPDF2 Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. 003 % 004 % CALL: f = thspdf2(Hd,Scf,Hm0,Tp) 005 % 006 % f = pdf struct evaluated at meshgrid(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 % 012 % THSPDF2 approximates the joint distribution of (Scf, Hd), i.e., crest 013 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process with a 014 % Torsethaugen spectral density. The empirical parameters of the model is 015 % fitted by least squares to simulated (Scf,Hd) data for 600 classes of 016 % Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves were 017 % simulated for each class of Hm0 and Tp. 018 % THSPDF is restricted to the following range for Hm0 and Tp: 019 % 0.5 < Hm0 [m] < 12, 3.5 < Tp [s] < 20, and Hm0 < (Tp-2)*12/11. 020 % 021 % Example: 022 % Hm0 = 6;Tp = 8; 023 % h = linspace(0,4*Hm0/sqrt(2)); 024 % s = linspace(0,6*1.25*Hm0/Tp^2); 025 % f = thspdf2(h,s,Hm0,Tp); 026 % w = linspace(0,40,5*1024+1).'; 027 % S = torsethaugen(w,[Hm0 Tp]); 028 % dt = 0.3; 029 % x = spec2sdat(S,80000,.2); rate = 8; 030 % [si,hi] = dat2steep(x,rate,2); 031 % fk = kdebin([si,hi],'epan',[],[],.5,128); 032 % fk.title = f.title; fk.labx = f.labx; 033 % plot(si,hi,'.'), hold on 034 % pdfplot(f),pdfplot(fk,'r'),hold off 035 % 036 % See also thsspdf 037 038 039 % Reference 040 % P. A. Brodtkorb (2004), 041 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 042 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 043 % Trondheim, Norway. 044 045 % History 046 % revised pab 09.08.2003 047 % changed input 048 % validated 20.11.2002 049 % By pab 20.12.2000 050 051 error(nargchk(4,5,nargin)) 052 053 if (nargin < 5|isempty(normalizedInput)), normalizedInput = 0;end 054 if (nargin < 4|isempty(Tp)), Tp = 8;end 055 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end 056 057 058 [V,H] = meshgrid(Scf,Hd); 059 060 f = createpdf(2); 061 [f.f,Hrms,Vrms,varargout{1:nargout-1}] = thspdf(H,V,Hm0,Tp,normalizedInput); 062 063 f.x = {Scf(:),Hd(:)}; 064 065 if (normalizedInput) 066 f.labx={'Scf', 'Hd'}; 067 f.norm = 1; 068 else 069 f.norm=0; 070 f.labx={'Scf', 'Hd [m]'}; 071 end 072 f.title = 'Joint distribution of (Hd,Scf) in time'; 073 f.note = ['Torsethaugen Hm0=' num2str(Hm0) ' Tp = ' num2str(Tp)]; 074 [f.cl,f.pl] = qlevels(f.f); 075 076 return 077 % old call 078 if Hm0>11| Hm0>(Tp-2)*12/11 079 disp('Warning: Hm0 is outside the valid range') 080 disp('The validity of the Joint (Hd,Scf) distribution is questionable') 081 end 082 if Tp>20|Tp<3 083 disp('Warning: Tp is outside the valid range') 084 disp('The validity of the Joint (Hd,Scf) distribution is questionable') 085 end 086 087 global THSPAR 088 if isempty(THSPAR) 089 THSPAR = load('thspar.mat'); 090 end 091 % Gamma distribution parameters as a function of Tp Hm0 and h2 092 A11 = THSPAR.A11s; 093 B11 = THSPAR.B11s; 094 095 % Waveheight distribution in time 096 097 if 0, 098 % Truncated Weibull distribution parameters as a function of Tp, Hm0 099 global THWPAR 100 if isempty(THWPAR) 101 THWPAR = load('thwpar.mat'); 102 end 103 A00 = THWPAR.A00s; 104 B00 = THWPAR.B00s; 105 C00 = THWPAR.C00s; 106 else 107 108 % Truncated Weibull distribution parameters as a function of Tp, Hm0 109 A00 = THSPAR.A00s; 110 B00 = THSPAR.B00s; 111 C00 = THSPAR.C00s; 112 end 113 114 Tpp = THSPAR.Tp; 115 Hm00 = THSPAR.Hm0; 116 h2 = THSPAR.h2; 117 118 119 120 121 w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 122 %w = linspace(0,10,2*1024+1).'; 123 S = torsethaugen(w,[Hm0,Tp]); 124 ch = spec2char(S,{'Tm02','eps2'}); 125 Tm02 = ch(1); 126 eps2 = ch(2); 127 128 Hrms = Hm0/sqrt(2); 129 Vrms = 1.25*Hm0/(Tm02^2); % Erms 130 131 if nargin<1 |isempty(v), v=linspace(0,4*Vrms); end 132 if nargin<2 |isempty(h), h=linspace(0,4*Hrms); end 133 134 if nargin>4, 135 h = h*Hrms; 136 v = v*Vrms; 137 end 138 139 %Fh = thpdf(h(:)/Hrms,Hm0,Tp,eps2,1); 140 141 [A0, B0, C0] = thwparfun(Hm0,Tp,'time'); 142 143 method = '*cubic';% Faster interpolation 144 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2); 145 A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tp(ones(size(h2))),... 146 Hm0(ones(size(h2))) ,h2,method),1,h/Hrms,1)); 147 B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tp(ones(size(h2))),... 148 Hm0(ones(size(h2))),h2,method),1,h/Hrms,1)); 149 150 [V1 H1] = meshgrid(v/Vrms,h/Hrms); 151 [V1 A1] = meshgrid(v/Vrms,A1); 152 [V1 B1] = meshgrid(v/Vrms,B1); 153 154 f = createpdf(2); 155 f.title = 'Joint distribution of (Hd,Scf) in time'; 156 157 if nargin<5 158 159 f.f = wtweibpdf(H1,A0,B0,C0).*wgampdf(V1,A1,B1)/Vrms/Hrms; 160 % f.f = repmat(Fh/Hrms,[1 length(v)]).*wgampdf(V1,A1,B1)/Vrms; 161 f.x = {v,h}; 162 f.norm=0; 163 f.labx={'Scf', 'Hd [m]'}; 164 else 165 f.f = wtweibpdf(H1,A0,B0,C0).*wgampdf(V1,A1,B1); 166 %f.f = repmat(Fh,[1 length(v)]).*wgampdf(V1,A1,B1); 167 f.x = {v/Vrms,h/Hrms}; 168 f.labx={'Scf', 'Hd'}; 169 f.norm = 1; 170 end 171 f.f(find(isnan(f.f)|isinf(f.f) ))=0; 172 173 f.note = ['Torsethaugen Hm0=' num2str(Hm0) ' Tp = ' num2str(Tp)]; 174 [f.cl,f.pl] = qlevels(f.f); 175 if nargout>1, 176 fA = createpdf(2); 177 fA.x = {Tpp,Hm00}; 178 fA.labx = {'Tp', 'Hm0'}; 179 fA(3) = fA(1); 180 fA(2) = fA(1); 181 182 fA(1).f = A00; 183 fA(2).f = B00; 184 fA(3).f = C00; 185 186 fA(1).title = 'wtweibpdf parameter A'; 187 fA(2).title = 'wtweibpdf parameter B'; 188 fA(3).title = 'wtweibpdf parameter C'; 189 190 txt1 = ['The Wtweibpdf distribution Parameter ']; 191 txt2=[' for wave heigth in time as a function of Tp and Hm0 for' ... 192 'the Torsethaugen spectrum']; 193 fA(1).note =[txt1 'A' txt2]; 194 fA(2).note =[txt1 'B' txt2]; 195 fA(3).note =[txt1 'C' txt2]; 196 197 tmp = [A00(:) B00(:) C00(:)]; 198 ra = range(tmp); 199 st = round(min(tmp)*100)/100; 200 en = max(tmp); 201 for ix = 1:3, 202 fA(ix).cl = st(ix):ra(ix)/20:en(ix); 203 end 204 end 205 if nargout>2, 206 fB = createpdf(3); 207 fB.x = {Tpp,Hm00,h2}; 208 fB.labx = {'Tp','Hm0', 'h'}; 209 fB(2) = fB(1); 210 211 fB(1).f = A11; 212 fB(2).f = B11; 213 214 txt11 = 'The conditonal Wgampdf distribution Parameter '; 215 txt22 = [' for Scf given h=Hd/Hrms in time as function of Tp' ... 216 ' and Hm0 for the Torsethaugen spectrum']; 217 fB(1).title = 'wgampdf parameter A'; 218 fB(2).title = 'wgampdf parameter B'; 219 fB(1).note = [txt11,'A',txt22]; 220 fB(2).note = [txt11,'B',txt22]; 221 222 %fB(2).note = ['The conditonal Wggampdf distribution Parameter B(h)/Hrms', ...% ' for crest front steepness as a function of Tp,Hm0 and',... 223 % ' h=Hd/Hrms for the Torsethaugen spectrum']; 224 tmp= [A11(:) B11(:)]; 225 ra = range(tmp); 226 st = round(min(tmp)*100)/100; 227 en = max(tmp); 228 for ix=1:2 229 fB(ix).cl = st(ix):ra(ix)/20:en(ix); 230 end 231 end 232 return 233 234 235 236 237 238 239 %Old calls 240 %method ='spline'; 241 switch 1 242 case 3,% Best fit by smoothing spline 243 brks = [0 .1 .2 .4 ,.6, .8, 1, 1.1 1.2]'; 244 coefa = [0 0 0.02260415153596 0.99807186986167; ... 245 2.19065400617385 0 0.02260415153596 1.00033228501527; ... 246 4.34015195156053 0.65719620185215 0.03709199393185 1.00478335417504; ... 247 -1.59533089716870 3.26128737278847 0.80983543882910 1.07321081664798; ... 248 -6.81273221810880 2.30408883448726 1.92291068028425 1.35286675214799; ... 249 -3.69498826658975 -1.78355049637802 2.09369407217829 1.77511058383946; ... 250 13.33514485443956 -4.00054345633187 0.94471547491809 2.09294747228728; ... 251 0 0 0.54466112928490 2.16074873007021]; 252 253 coefb = [ 0 0 0.32503235228616 1.99054481866418; ... 254 3.28321899128157 0 0.32503235228616 2.02304805389280; ... 255 5.67672309005450 0.98496569738447 0.37964649056830 2.05883450811270; ... 256 -5.29907238080822 4.39099955141717 1.43842344537222 2.21957621884217; .... 257 -5.89663569823287 1.21155612293224 2.55893458024211 2.64050831092684; ... 258 -6.21824739906323 -2.32642529600749 2.43691697455115 3.15358438630669; ... 259 20.19124578481806 -6.05737373544542 0.77134599291113 3.49816479018411; ... 260 0 0 0.16560861936659 3.53491689790559]; 261 coefc =[ 0 0 0.04818579357214 -0.00817761487085; ... 262 2.94432030165157 0 0.04818579357214 -0.00335903551363; ... 263 4.77660844045250 0.88329609049547 0.09917317190900 0.00440386414523; ... 264 -1.24578770271258 3.74926115476697 1.01096301945323 0.09778320967047; ... 265 -7.70868155645400 3.00178853313943 2.36117295703451 0.43997995813009; ... 266 -3.98346578867600 -1.62342040073298 2.70373824808144 0.97061663841094; ... 267 13.37833291312857 -4.01349987393858 1.57933014446810 1.41455974568850; ... 268 0 0 1.17798015707424 1.54573609430905]; 269 pa = mkpp(brks,coefa); 270 pb = mkpp(brks,coefb); 271 pc = mkpp(brks,coefc); 272 A0 = ppval(pa,eps2); 273 B0 = ppval(pb,eps2); 274 C0 = ppval(pc,eps2); 275 case 1, 276 277 [E1, H1] = meshgrid(Tpp,Hm00); 278 A0 = interp2(E1,H1,A00,Tp,Hm0,method); 279 B0 = interp2(E1,H1,B00,Tp,Hm0,method); 280 C0 = interp2(E1,H1,C00,Tp,Hm0,method); 281 end
Comments or corrections to the WAFO group