Home > wafo > wavemodels > jhvpdf.m

# jhvpdf

## PURPOSE

Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum.

## SYNOPSIS

[f,Hrms,Vrms,fA,fB] = jhvpdf(Hd,Vcf,Hm0,Tp,gam,normalizedInput,condon)

## DESCRIPTION

``` 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

## CROSS-REFERENCE INFORMATION

This function calls:
 createpdf PDF class constructor getjonswappeakedness Peakedness factor Gamma given Hm0 and Tp for JONSWAP jhwparfun Wave height, Hd, distribution parameters for Jonswap spectrum. jonswap Calculates (and plots) a JONSWAP spectral density range Calculates the difference between the maximum and minimum values. smooth Calculates a smoothing spline. spec2char Evaluates spectral characteristics and their covariance wtweibpdf Truncated Weibull probability density function wweibcdf Weibull cumulative distribution function wweibpdf Weibull probability density function comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. interp2 2-D interpolation (table lookup). linspace Linearly spaced vector. meshgrid X and Y arrays for 3-D plots. polyval Evaluate polynomial. unique Set unique. warning Display warning message; disable or enable warning messages.
This function is called by:
 jhvcdf Joint (Vcf,Hd) CDF for linear waves with JONSWAP spectrum. jhvpdf2 Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum.

## SOURCE CODE

```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 %
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)
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```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003