JHVPDF Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. CALL: f = jhvpdf(Hd,Vcf,Hm0,Tp,gamma) f = pdf evaluated at (Vcf,Hd) Hd = zero down crossing wave height [m] Vcf = crest front velocity [m/s] Hm0 = significant wave height [m] Tp = Spectral peak period [s] Gamma = Peakedness parameter of the JONSWAP spectrum JHVPDF approximates the joint distribution of (Vcf, Hd), i.e., crest front velocity (Ac/Tcf) and wave height, for a Gaussian process with a JONSWAP spectral density. The empirical parameters of the model is fitted by least squares to simulated (Vcf,Hd) data for 13 classes of GAMMA between 1 and 7. About 100000 zero-downcrossing waves were simulated for each class. JHVPDF is restricted to the following range for GAMMA: 1 <= GAMMA <= 7 NOTE: The size of f is the common size of the input arguments. Example: Hm0 = 6;Tp = 9; gam=3.5 h = linspace(0,4*Hm0/sqrt(2))'; v = linspace(0,4*2*Hm0/Tp)'; [V,H] = meshgrid(v,h); f = jhvpdf(H,V,Hm0,Tp,gam); w = linspace(0,40,5*1024+1).'; S = jonswap(w,[Hm0, Tp, gam]); dt = .3; x = spec2sdat(S,80000,dt); rate = 4; [vi,hi] = dat2steep(x,rate,1); fk = kdebin([vi,hi],'epan',[],[],.5,128); plot(vi,hi,'.'), hold on contour(v,h,f,fk.cl), pdfplot(fk,'r'),hold off See also thvpdf
PDF class constructor | |
Peakedness factor Gamma given Hm0 and Tp for JONSWAP | |
Wave height, Hd, distribution parameters for Jonswap spectrum. | |
Calculates (and plots) a JONSWAP spectral density | |
Calculates the difference between the maximum and minimum values. | |
Calculates a smoothing spline. | |
Evaluates spectral characteristics and their covariance | |
Truncated Weibull probability density function | |
Weibull cumulative distribution function | |
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). | |
Linearly spaced vector. | |
X and Y arrays for 3-D plots. | |
Evaluate polynomial. | |
Set unique. | |
Display warning message; disable or enable warning messages. |
Joint (Vcf,Hd) CDF for linear waves with JONSWAP spectrum. | |
Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. |
001 function [f,Hrms,Vrms,fA,fB] = jhvpdf(Hd,Vcf,Hm0,Tp,gam,normalizedInput,condon) 002 %JHVPDF Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. 003 % 004 % CALL: f = jhvpdf(Hd,Vcf,Hm0,Tp,gamma) 005 % 006 % f = pdf evaluated at (Vcf,Hd) 007 % Hd = zero down crossing wave height [m] 008 % Vcf = crest front velocity [m/s] 009 % Hm0 = significant wave height [m] 010 % Tp = Spectral peak period [s] 011 % Gamma = Peakedness parameter of the JONSWAP spectrum 012 % 013 % JHVPDF approximates the joint distribution of (Vcf, Hd), i.e., crest 014 % front velocity (Ac/Tcf) and wave height, for a Gaussian process with a 015 % JONSWAP spectral density. The empirical parameters of the model is 016 % fitted by least squares to simulated (Vcf,Hd) data for 13 classes of 017 % GAMMA between 1 and 7. About 100000 zero-downcrossing waves 018 % were simulated for each class. 019 % JHVPDF is restricted to the following range for GAMMA: 020 % 1 <= GAMMA <= 7 021 % 022 % NOTE: The size of f is the common size of the input arguments. 023 % 024 % Example: 025 % Hm0 = 6;Tp = 9; gam=3.5 026 % h = linspace(0,4*Hm0/sqrt(2))'; 027 % v = linspace(0,4*2*Hm0/Tp)'; 028 % [V,H] = meshgrid(v,h); 029 % f = jhvpdf(H,V,Hm0,Tp,gam); 030 % w = linspace(0,40,5*1024+1).'; 031 % S = jonswap(w,[Hm0, Tp, gam]); 032 % dt = .3; 033 % x = spec2sdat(S,80000,dt); rate = 4; 034 % [vi,hi] = dat2steep(x,rate,1); 035 % fk = kdebin([vi,hi],'epan',[],[],.5,128); 036 % plot(vi,hi,'.'), hold on 037 % contour(v,h,f,fk.cl), 038 % pdfplot(fk,'r'),hold off 039 % 040 % See also thvpdf 041 042 % Reference 043 % P. A. Brodtkorb (2004), 044 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 045 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 046 % Trondheim, Norway. 047 048 % History 049 % revised pab 10 Jan 2004 050 % By pab 20.12.2000 051 052 error(nargchk(4,7,nargin)) 053 if (nargin < 7|isempty(condon)), condon = 0; end 054 if (nargin < 6|isempty(normalizedInput)), normalizedInput = 0;end 055 if (nargin < 5|isempty(gam)) 056 gam = getjonswappeakedness(Hm0,Tp); 057 end 058 059 multipleSeaStates = any(prod(size(Hm0))>1|... 060 prod(size(Tp)) >1|... 061 prod(size(gam))>1); 062 if multipleSeaStates 063 [errorcode, Vcf,Hd,Hm0,Tp,gam] = comnsize(Vcf,Hd,Hm0,Tp,gam); 064 else 065 [errorcode, Vcf,Hd] = comnsize(Vcf,Hd); 066 end 067 if errorcode > 0 068 error('Requires non-scalar arguments to match in size.'); 069 end 070 071 displayWarning = 0; 072 if displayWarning 073 if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0))) 074 disp('Warning: Hm0,Tp is outside the JONSWAP range') 075 disp('The validity of the parameters returned are questionable') 076 end 077 end 078 k = find(gam<1); 079 if any(k) 080 if displayWarning, 081 warning('Peakedness parameter less than 1. Must be larger than 1.') 082 end 083 gam(k)=1; 084 end 085 k1 = find(7<gam); 086 if any(k1) 087 if displayWarning 088 warning('Peakedness parameter larger than 7. The pdf returned is questionable') 089 end 090 gam(k1) = 7; 091 end 092 093 global JHVPAR 094 if isempty(JHVPAR) 095 JHVPAR = load('jhvpar.mat'); 096 end 097 % Weibull distribution parameters as a function of e2 and h2 098 A11 = JHVPAR.A11s; 099 B11 = JHVPAR.B11s; 100 e2 = JHVPAR.gam; % gamma 101 h2 = JHVPAR.h2; % Hd/Hrms 102 [E2 H2] = meshgrid(e2,h2); 103 104 if 1, 105 106 %Tm02 = Tp./(1.30301-0.01698*gam+0.12102./gam); 107 %dev = 2e-5; 108 c1 =[ 0.16183666835624 1.53691936441548 1.55852759524555]; 109 c2 =[ 0.15659478203944 1.15736959972513 1]; 110 Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam)); 111 112 %dev = 2e-4; 113 %EPS2cof = [0.00068263671017 -0.01802256231624 0.44176198490431]; 114 %eps2 = polyval(EPS2cof,gam); 115 else 116 w = linspace(0,100,16*1024+1).'; % jonswap original spacing 117 %Hm0 = 6; 118 %gam = linspace(1,7,32); 119 Tm02 = zeros(size(gam)); 120 eps2 = Tm02; 121 for ix=1:length(gam(:)) 122 ch = spec2char(jonswap(w,[Hm0(ix),Tp(ix) gam(ix)]),{'Tm02','eps2'}); 123 Tm02(ix) = ch(1); 124 eps2(ix) = ch(2); 125 end 126 end 127 128 129 if normalizedInput 130 Hrms = 1; 131 Vrms = 1; 132 else 133 Hrms = Hm0/sqrt(2); 134 Vrms = 2*Hm0./Tm02; % Erms 135 end 136 137 v = Vcf./Vrms; 138 h = Hd./Hrms; 139 cSize = size(h); % common size of input 140 141 142 143 method ='*cubic'; 144 Nh2 = length(h2); 145 if multipleSeaStates 146 h = h(:); 147 v = v(:); 148 Tp = Tp(:); 149 Hm0 = Hm0(:); 150 gam = gam(:); 151 % eps2 = eps2(:); 152 A1 = zeros(length(h),1); 153 B1 = A1; 154 [gamu,ix,jx] = unique(gam); 155 numSeaStates = length(gamu); 156 gami = zeros(Nh2,1); 157 for iz=1:numSeaStates 158 k = find(jx==iz); 159 gami(:) = gamu(iz); 160 A1(k) = smooth(h2,interp2(E2,H2,A11,gami,h2,method),1,h(k),1); 161 B1(k) = smooth(h2,interp2(E2,H2,B11,gami,h2,method),1,h(k),1); 162 end 163 else 164 A1 = smooth(h2,interp2(E2,H2,A11,gam(ones(size(h2))),h2,method),1,h,1); 165 B1 = smooth(h2,interp2(E2,H2,B11,gam(ones(size(h2))),h2,method),1,h,1); 166 end 167 168 % Waveheight distribution in time 169 % Truncated Weibull distribution parameters as a function of Tp, Hm0, gam 170 [A0, B0, C0] = jhwparfun(Hm0,Tp,gam,'time'); 171 172 switch condon, 173 case 0, % regular pdf is returned 174 f = wtweibpdf(h,A0,B0,C0).*wweibpdf(v,A1,B1); 175 case 1, %pdf conditioned on x1 ie. p(x2|x1) 176 f = wweibpdf(v,A1,B1); 177 case 3, % secret option used by XXstat: returns x2*p(x2|x1) 178 f = v.*wweibpdf(v,A1,B1); 179 case 4, % secret option used by XXstat: returns x2.^2*p(x2|x1) 180 f = v.^2.*wweibpdf(v,A1,B1); 181 case 5, % p(h)*P(V|h) is returned special case used by jhvcdf2 182 f = wtweibpdf(h,A0,B0,C0).*wweibcdf(v,A1,B1); 183 case 6, % P(V|h) is returned special case used by jhvcdf2 184 f = wweibcdf(v,A1,B1); 185 case 7,% p(h)*(1-P(V|h)) is returned special case used by jhvcdf2 186 f = wtweibpdf(h,A0,B0,C0).*(1-wweibcdf(v,A1,B1)); 187 otherwise error('unknown option') 188 end 189 190 if multipleSeaStates 191 f = reshape(f,cSize); 192 end 193 194 if condon~=6 195 f = f./Hrms./Vrms; 196 end 197 f(find(isnan(f)|isinf(f) ))=0; 198 if any(size(f)~=cSize) 199 disp('Wrong size') 200 end 201 202 203 if nargout>3, 204 fA = createpdf(2); 205 fA.f = A11; 206 fA.x = {e2,h2}; 207 fA.labx = {'Gamma', 'h'}; 208 fA.note = ['The conditonal Weibull distribution Parameter A(h,gamma)/Hrms' ... 209 'for Vcf as a function of h=Hd/Hrms and gamma for' ... 210 'the Jonswap spectrum']; 211 212 ra = range(A11(:)); 213 st = round(min(A11(:))*100)/100; 214 en = max(A11(:)); 215 fA.cl = st:ra/20:en; 216 end 217 if nargout>4, 218 fB = createpdf(2); 219 fB.f = B11; 220 fB.x = {e2,h2}; 221 fB.labx = {'Gamma', 'h'}; 222 fB.note = ['The conditonal Weibull distribution Parameter B(h,gamma)/Hrms' ... 223 'for Vcf as a function of h=Hd/Hrms and gamma for' ... 224 'the Jonswap spectrum']; 225 ra = range(B11(:)); 226 st = round(min(B11(:))*100)/100; 227 en = max(B11(:)); 228 fB.cl = st:ra/20:en; 229 end 230 return 231 232 233 234 235 236 237
Comments or corrections to the WAFO group