Home > wafo > wavemodels > thvcdf.m

# thvcdf

## PURPOSE

Joint (Vcf,Hd) CDF for linear waves with Torsethaugen spectra.

## SYNOPSIS

f = thvcdf(Hd,Vcf,Hm0,Tp,tail)

## DESCRIPTION

``` THVCDF Joint (Vcf,Hd) CDF for linear waves with Torsethaugen spectra.

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

f   = CDF evaluated at (Vcf,Hd)
Hd  = zero down crossing wave height
Vcf = crest front velocity
Hm0 = significant wave height [m]
Tp  = Spectral peak period    [s]
tail = 1 if upper tail is calculated
0 if lower tail is calulated (default)

THVCDF approximates the joint CDF of (Vcf, Hd), i.e., crest front
velocity (Ac/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 (Vcf,Hd) data for 600 classes of
Hm0 and Tp. Between 50000 and 150000 zero-downcrossing waves were
simulated for each class of Hm0 and Tp.
THVCDF 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;
vc = 3;
hc = 3;
lowerTail = 0;
upperTail = ~lowerTail
thvcdf(hc,vc,Hm0,Tp)           % Prob(Hd<Hc,Vcf<Vc)
thvcdf(hc,vc,Hm0,Tp,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 = thvcdf(hc,vc,H,T,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. thvpdf Joint (Vcf,Hd) PDF for linear waves with Torsethaugen spectra. 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). mean Average or mean value. meshgrid X and Y arrays for 3-D plots. warning Display warning message; disable or enable warning messages.
This function is called by:

## SOURCE CODE

```001 function f = thvcdf(Hd,Vcf,Hm0,Tp,tail)
002 %THVCDF Joint (Vcf,Hd) CDF for linear waves with Torsethaugen spectra.
003 %
004 %  CALL: f = thvcdf(Hd,Vcf,Hm0,Tp,tail)
005 %
006 %   f   = CDF evaluated at (Vcf,Hd)
007 %   Hd  = zero down crossing wave height
008 %   Vcf = crest front velocity
009 %   Hm0 = significant wave height [m]
010 %   Tp  = Spectral peak period    [s]
011 %  tail = 1 if upper tail is calculated
012 %         0 if lower tail is calulated (default)
013 %
014 % THVCDF approximates the joint CDF of (Vcf, Hd), i.e., crest front
015 % velocity (Ac/Tcf) and wave height, for a Gaussian process with a
016 % Torsethaugen spectral density. The empirical parameters of the model is
017 % fitted by least squares to simulated (Vcf,Hd) data for 600 classes of
018 % Hm0 and Tp. Between 50000 and 150000 zero-downcrossing waves were
019 % simulated for each class of Hm0 and Tp.
020 % THVCDF is restricted to the following range for Hm0 and Tp:
021 %  0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11.
022 %
023 % Example:
024 % Hm0 = 6;Tp = 8;
025 % vc = 3;
026 % hc = 3;
027 % lowerTail = 0;
028 % upperTail = ~lowerTail
029 % thvcdf(hc,vc,Hm0,Tp)           % Prob(Hd<Hc,Vcf<Vc)
030 % thvcdf(hc,vc,Hm0,Tp,upperTail) % Prob(Hd>Hc,Vcf>Vc)
031 %
032 %  % Conditional probability of steep and high waves given seastates
033 %  % i.e., Prob(Hd>hc,Vcf>vc|Hs,Tp)
034 %  upperTail = 1;
035 %  Hs = linspace(2.5,11.5,10);
036 %  Tp = linspace(4.5,19.5,16);
037 %  [T,H] = meshgrid(Tp,Hs);
038 %  p = thvcdf(hc,vc,H,T,upperTail);
039 %  v = 10.^(-6:-1);
040 %  contourf(Tp,Hs,log10(p),log10(v))
041 %  xlabel('Tp')
042 %  ylabel('Hs')
043 %  fcolorbar(log10(v))
044 %
046
047 % Reference
048 % P. A. Brodtkorb (2004),
049 % The Probability of Occurrence of Dangerous Wave Situations at Sea.
050 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
051 % Trondheim, Norway.
052
053 % History
054 % revised pab 09.08.2003
055 % changed input
056 % validated 20.11.2002
057 % By pab 20.12.2000
058
059
060 error(nargchk(3,5,nargin))
061
062 if (nargin < 5|isempty(tail)),tail = 0; end
063 if (nargin < 4|isempty(Tp)),Tp = 8; end
064 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end
065
066 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1);
067 if multipleSeaStates
068   [errorcode, Vcf,Hd,Hm0,Tp] = comnsize(Vcf,Hd,Hm0,Tp);
069 else
070   [errorcode, Vcf,Hd] = comnsize(Vcf,Hd);
071 end
072 if errorcode > 0
073   error('Requires non-scalar arguments to match in size.');
074 end
075
076 global THVPAR
077 if isempty(THVPAR)
079 end
080
081 Tpp  = THVPAR.Tp;
082 Hm00 = THVPAR.Hm0;
083 Tm020 = THVPAR.Tm02;
084 % Interpolation method
085 method = '*cubic';% Faster interpolation
086
087 [Tp1,Hs1] = meshgrid(Tpp,Hm00);
088 Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method);
089 %  w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing
090   %w    = linspace(0,10,2*1024+1).';
091 %  St = torsethaugen(w,[Hm0,Tp]);
092 %  ch   = spec2char(St,{'Tm02','eps2'});
093 %  Tm02 = ch(1);
094 %  eps2 = ch(2);
095 Hrms = Hm0/sqrt(2);
096 Vrms = 2*Hm0./Tm02; % Erms
097
098 v = Vcf./Vrms;
099 f = zeros(size(Hd));
100
101 % Only compute within valid range
102 k0 = find((2<=Tp) & (Tp<=21) & (Hm0<=(Tp-2)*12/11) & (Hm0<=12));
103 if any(k0)
104   hMax = 5;
105   eps2 = 1e-6;
106   h = min(Hd./Hrms,hMax);
107
108   if multipleSeaStates
109     h = h(k0);
110     v = v(k0);
111     Hm0 = Hm0(k0);
112     Tp = Tp(k0);
113   else
114     k0 = 1:prod(size(Hd));
115   end
116   if 0
117     % This is a trick to get the html documentation correct.
118     k = thvpdf(1,1,2,3);
119   end
120   normalizedInput = 1;
121   utprb = gaussq('thvpdf',hMax,2*hMax,eps2/2,[],mean(v(:)),mean(Hm0(:)),mean(Tp(:)),normalizedInput,7);
122   if eps2<utprb
123     warning('Check the accuracy of integration!')
124   end
125
126
127
128   hlim    = h;
129
130
131   lowerTail = 0;
132   if tail==lowerTail,
133     k       = find(h>2.5);
134     hlim(k) = 2.5;
135     k       = find(h>1.3*v);
136     hlim(k) = 1.3*v(k);
137     f(k0) = gaussq('thvpdf',0,hlim,eps2/2,[],v,Hm0,Tp,normalizedInput,5)...
138     + gaussq('thvpdf',hlim,h,eps2/2,[],v,Hm0,Tp,normalizedInput,5);
139   else % upper tail
140     k       = find(h<1.3*v);
141     hlim(k) = 1.3*v(k);
142     f(k0) = gaussq('thvpdf',h,hlim,eps2/2,[],v,Hm0,Tp,normalizedInput,7)...
143     + gaussq('thvpdf',hlim,hMax,eps2/2,[],v,Hm0,Tp,normalizedInput,7);
144   end
145 end
146 return
147
148```

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