Home > wafo > wavemodels > jhvnlcdf.m

# jhvnlcdf

## PURPOSE

Joint (Vcf,Hd) CDF for non-linear waves with JONSWAP spectrum.

## SYNOPSIS

f = jhvnlcdf(Hd,Vcf,Hm0,Tp,gam,tail)

## DESCRIPTION

``` JHVNLCDF Joint (Vcf,Hd) CDF for non-linear waves with JONSWAP spectrum.

CALL: f = jhvnlcdf(Hd,Vcf,Hm0,Tp,Gamma,tail)

f   = CDF 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
tail = 1 if upper tail is calculated
0 if lower tail is calulated (default)

JHVNLCDF approximates the joint CDF 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. Between 50000 and 150000 zero-downcrossing waves were
simulated for each class of GAMMA.
JHVNLCDF is restricted to the following range for GAMMA:
1 <= GAMMA <= 7

Example:
Hm0 = 6;Tp = 8; gam = 3.5;
vc = 3;
hc = 3;
lowerTail = 0;
upperTail = ~lowerTail
jhvnlcdf(hc,vc,Hm0,Tp,gam)           % Prob(Hd<Hc,Vcf<Vc)
jhvnlcdf(hc,vc,Hm0,Tp,gam,upperTail) % Prob(Hd>Hc,Vcf>Vc)

% Conditional probability of steep and high waves given seastates
% i.e., Prob(Hd>hc,Vcf>vc|Hs,Tp)
upperTail = 1;
Hs = linspace(2.5,11.5,10);
Tp = linspace(4.5,19.5,16);
[T,H] = meshgrid(Tp,Hs);
pl = jhvcdf(hc,vc,H,T,gam,upperTail);
p = jhvnlcdf(hc,vc,H,T,gam,upperTail);
v = 10.^(-6:-1);
contourf(Tp,Hs,log10(p),log10(v))
xlabel('Tp'), ylabel('Hs'),  fcolorbar(log10(v))
figure(2)
contourf(Tp,Hs,log10(pl),log10(v))
xlabel('Tp'), ylabel('Hs'),  fcolorbar(log10(v))

## CROSS-REFERENCE INFORMATION

This function calls:
 gaussq Numerically evaluates a integral using a Gauss quadrature. getjonswappeakedness Peakedness factor Gamma given Hm0 and Tp for JONSWAP jhvnlpdf Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. mean Average or mean value. polyval Evaluate polynomial. warning Display warning message; disable or enable warning messages.
This function is called by:

## SOURCE CODE

```001 function f = jhvnlcdf(Hd,Vcf,Hm0,Tp,gam,tail)
002 %JHVNLCDF Joint (Vcf,Hd) CDF for non-linear waves with JONSWAP spectrum.
003 %
004 %  CALL: f = jhvnlcdf(Hd,Vcf,Hm0,Tp,Gamma,tail)
005 %
006 %   f   = CDF 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 %  tail = 1 if upper tail is calculated
013 %         0 if lower tail is calulated (default)
014 %
015 % JHVNLCDF approximates the joint CDF of (Vcf, Hd), i.e., crest front
016 % velocity (Ac/Tcf) and wave height, for a Gaussian process with a
017 % JONSWAP spectral density. The empirical parameters of the model is
018 % fitted by least squares to simulated (Vcf,Hd) data for 13 classes of
019 % GAMMA between 1 and 7. Between 50000 and 150000 zero-downcrossing waves were
020 % simulated for each class of GAMMA.
021 % JHVNLCDF is restricted to the following range for GAMMA:
022 %  1 <= GAMMA <= 7
023 %
024 % Example:
025 % Hm0 = 6;Tp = 8; gam = 3.5;
026 % vc = 3;
027 % hc = 3;
028 % lowerTail = 0;
029 % upperTail = ~lowerTail
030 % jhvnlcdf(hc,vc,Hm0,Tp,gam)           % Prob(Hd<Hc,Vcf<Vc)
031 % jhvnlcdf(hc,vc,Hm0,Tp,gam,upperTail) % Prob(Hd>Hc,Vcf>Vc)
032 %
033 %  % Conditional probability of steep and high waves given seastates
034 %  % i.e., Prob(Hd>hc,Vcf>vc|Hs,Tp)
035 %  upperTail = 1;
036 %  Hs = linspace(2.5,11.5,10);
037 %  Tp = linspace(4.5,19.5,16);
038 %  [T,H] = meshgrid(Tp,Hs);
039 %  pl = jhvcdf(hc,vc,H,T,gam,upperTail);
040 %  p = jhvnlcdf(hc,vc,H,T,gam,upperTail);
041 %  v = 10.^(-6:-1);
042 %  contourf(Tp,Hs,log10(p),log10(v))
043 %  xlabel('Tp'), ylabel('Hs'),  fcolorbar(log10(v))
044 %  figure(2)
045 %  contourf(Tp,Hs,log10(pl),log10(v))
046 %  xlabel('Tp'), ylabel('Hs'),  fcolorbar(log10(v))
047 %
048 %
050
051 % Reference
052 % P. A. Brodtkorb (2004),
053 % The Probability of Occurrence of Dangerous Wave Situations at Sea.
054 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
055 % Trondheim, Norway.
056
057
058 % History
059 % revised pab 09.08.2003
060 % changed input
061 % validated 20.11.2002
062 % By pab 20.12.2000
063
064
065 error(nargchk(2,6,nargin))
066
067 if (nargin < 6|isempty(tail)),tail = 0; end
068 if (nargin < 4|isempty(Tp)),Tp = 8; end
069 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end
070 if (nargin < 5|isempty(gam)),
071    gam = getjonswappeakedness(Hm0,Tp);
072 end
073
074 multipleSeaStates = any(prod(size(Hm0))>1|...
075             prod(size(Tp))>1|...
076             prod(size(gam))>1);
077 if multipleSeaStates
078   [errorcode, Vcf,Hd,Hm0,Tp,gam] = comnsize(Vcf,Hd,Hm0,Tp,gam);
079 else
080   [errorcode, Vcf,Hd] = comnsize(Vcf,Hd);
081 end
082 if errorcode > 0
083   error('Requires non-scalar arguments to match in size.');
084 end
085 displayWarning = 0;
086 if displayWarning
087   if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0)))
088     disp('Warning: Hm0,Tp is outside the JONSWAP range')
089     disp('The validity of the parameters returned are questionable')
090   end
091 end
092 %dev = 2e-5;
093 c1 =[ 0.16183666835624   1.53691936441548   1.55852759524555];
094 c2 =[ 0.15659478203944   1.15736959972513   1];
095 Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam));
096
097 Hrms = Hm0/sqrt(2);
098 Vrms = 2*Hm0./Tm02; % Erms
099
100 v = Vcf./Vrms;
101
102 hMax = 5;
103 eps2 = 1e-6;
104 normalizedInput = 1;
105
106
107 h = min(Hd./Hrms,hMax);
108
109 f = zeros(size(Hd));
110 % Only compute
111 k0 = find((Tp<5*sqrt(Hm0)) & (3.6*sqrt(Hm0)<Tp));
112 if any(k0)
113   if multipleSeaStates
114     h = h(k0);
115     v = v(k0);
116     Hm0 = Hm0(k0);
117     Tp = Tp(k0);
118     gam = gam(k0);
119   else
120     k0 = 1:prod(size(Hd));
121   end
122
123
124   utprb = gaussq('jhvnlpdf',hMax,2*hMax,eps2/2,[],mean(v(:)),mean(Hm0(:)),mean(Tp(:)),mean(gam(:)),normalizedInput,7);
125   if eps2<utprb
126     warning('Check the accuracy of integration!')
127   end
128
129   if 0
130     % This is a trick to get the html documentation correct.
131     k = jhvnlpdf(1,1,2,3);
132   end
133   hlim  = h;
134
135   lowerTail = 0;
136   if tail==lowerTail,
137     k       = find(h>1.3*v);
138     hlim(k) = 1.3*v(k);
139
140     f(k0) = gaussq('jhvnlpdf',0,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5)...
141         + gaussq('jhvnlpdf',hlim,h,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5);
142   else % upper tail
143     k       = find(h<1.3*v);
144     hlim(k) = 1.3*v(k);
145     f(k0) = gaussq('jhvnlpdf',h,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7)...
146       + gaussq('jhvnlpdf',hlim,hMax,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7);
147   end
148 end
149 return
150
151```

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