Home > wafo > wavemodels > jhvcdf.m

# jhvcdf

## PURPOSE

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

## SYNOPSIS

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

## DESCRIPTION

``` JHVCDF Joint (Vcf,Hd) CDF for linear waves with JONSWAP spectrum.

CALL: f = jhvcdf(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)

JHVCDF 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. About 100000 zero-downcrossing waves were
simulated for each class of GAMMA.
JHVCDF 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
jhvcdf(hc,vc,Hm0,Tp,gam)           % Prob(Hd<Hc,Vcf<Vc)
jhvcdf(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);
p = jhvcdf(hc,vc,H,T,gam,upperTail);
v = 10.^(-6:-1);
contourf(Tp,Hs,log10(p),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 jhvpdf 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 = jhvcdf(Hd,Vcf,Hm0,Tp,gam,tail)
002 %JHVCDF Joint (Vcf,Hd) CDF for linear waves with JONSWAP spectrum.
003 %
004 %  CALL: f = jhvcdf(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 % JHVCDF 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. About 100000 zero-downcrossing waves were
020 % simulated for each class of GAMMA.
021 % JHVCDF 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 % jhvcdf(hc,vc,Hm0,Tp,gam)           % Prob(Hd<Hc,Vcf<Vc)
031 % jhvcdf(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 %  p = jhvcdf(hc,vc,H,T,gam,upperTail);
040 %  v = 10.^(-6:-1);
041 %  contourf(Tp,Hs,log10(p),log10(v))
042 %  xlabel('Tp')
043 %  ylabel('Hs')
044 %  fcolorbar(log10(v))
045 %
047
048
049 % Reference
050 % P. A. Brodtkorb (2004),
051 % The Probability of Occurrence of Dangerous Wave Situations at Sea.
052 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
053 % Trondheim, Norway.
054
055 % History
056 % revised pab 09.08.2003
057 % changed input
058 % validated 20.11.2002
059 % By pab 20.12.2000
060
061
062 error(nargchk(2,6,nargin))
063
064 if (nargin < 6|isempty(tail)),tail = 0; end
065 if (nargin < 4|isempty(Tp)),Tp = 8; end
066 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end
067 if (nargin < 5|isempty(gam)),
068    gam = getjonswappeakedness(Hm0,Tp);
069 end
070
071 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1);
072 if multipleSeaStates
073   [errorcode, Vcf,Hd,Hm0,Tp,gam] = comnsize(Vcf,Hd,Hm0,Tp,gam);
074 else
075   [errorcode, Vcf,Hd] = comnsize(Vcf,Hd);
076 end
077 if errorcode > 0
078   error('Requires non-scalar arguments to match in size.');
079 end
080 displayWarning = 0;
081 if displayWarning
082   if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0)))
083     disp('Warning: Hm0,Tp is outside the JONSWAP range')
084     disp('The validity of the parameters returned are questionable')
085   end
086 end
087 %dev = 2e-5;
088 c1 =[ 0.16183666835624   1.53691936441548   1.55852759524555];
089 c2 =[ 0.15659478203944   1.15736959972513   1];
090 Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam));
091
092 Hrms = Hm0/sqrt(2);
093 Vrms = 2*Hm0./Tm02; % Erms
094
095 v = Vcf./Vrms;
096
097 hMax = 6;
098 eps2 = 1e-6;
099 normalizedInput = 1;
100
101
102 h = min(Hd./Hrms,hMax);
103 f = zeros(size(Hd));
104 % Only compute
105 k0 = find((Tp<5*sqrt(Hm0)) & (3.6*sqrt(Hm0)<Tp));
106 if any(k0)
107   if multipleSeaStates
108     h = h(k0);
109     v = v(k0);
110     Hm0 = Hm0(k0);
111     Tp = Tp(k0);
112     gam = gam(k0);
113   else
114     k0 = 1:prod(size(Hd));
115   end
116   utprb = gaussq('jhvpdf',hMax,2*hMax,eps2/2,[],mean(v(:)),mean(Hm0(:)),mean(Tp(:)),mean(gam(:)),normalizedInput,7);
117   if eps2<utprb
118     warning('Check the accuracy of integration!')
119   end
120
121   if 0
122     % This is a trick to get the html documentation correct.
123     k = jhvpdf(1,1,2,3);
124   end
125
126   hlim  = h;
127
128
129   lowerTail = 0;
130   if tail==lowerTail,
131     k       = find(h>2*v);
132     hlim(k) = 2*v(k);
133     f(k0) = gaussq('jhvpdf',0,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5)...
134         + gaussq('jhvpdf',hlim,h,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5);
135   else % upper tail
136     k       = find(h<2*v);
137     hlim(k) = 2*v(k);
138     f(k0) = gaussq('jhvpdf',h,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7)...
139         + gaussq('jhvpdf',hlim,hMax,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7);
140   end
141 end
142 return
143
144```

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