Home > wafo > wavemodels > thscdf.m

# thscdf

## PURPOSE

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

## SYNOPSIS

f = thscdf(Hd,Scf,Hm0,Tp,tail)

## DESCRIPTION

``` THSCDF Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra.

CALL: f = thscdf(Hd,Scf,Hm0,Tp,tail)

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

THSCDF approximates the joint CDF of (Scf, Hd), i.e., crest front
steepness (2*pi*Ac/(g*Td*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 (Scf,Hd) data for 600
classes of Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves
were simulated for each class of Hm0 and Tp.
THSCDF 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;
sc = 0.25;
hc = 3;
lowerTail = 0;
upperTail = ~lowerTail
thscdf(hc,sc,Hm0,Tp)           % Prob(Hd<Hc,Scf<Ec)
thscdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>Hc,Scf>Ec)

% Conditional probability of steep and high waves given seastates
% i.e., Prob(Hd>hc,Scf>sc|Hs,Tz)
upperTail = 1;
Hs = linspace(2.5,11.5,10);
Tp = linspace(4.5,19.5,16);
[T,H] = meshgrid(Tp,Hs);
p = thscdf(hc,sc,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. thspdf Joint (Scf,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). meshgrid X and Y arrays for 3-D plots.
This function is called by:

## SOURCE CODE

```001 function f = thscdf(Hd,Scf,Hm0,Tp,tail)
002 %THSCDF Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra.
003 %
004 %  CALL: f = thscdf(Hd,Scf,Hm0,Tp,tail)
005 %
006 %   f   = CDF evaluated at (Scf,Hd)
007 %   Hd  = zero down crossing wave height (vector)
008 %   Scf = crest front steepness (vector)
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 % THSCDF approximates the joint CDF of (Scf, Hd), i.e., crest front
015 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process
016 % with a Torsethaugen spectral density. The empirical parameters of the
017 % model is fitted by least squares to simulated (Scf,Hd) data for 600
018 % classes of Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves
019 % were simulated for each class of Hm0 and Tp.
020 % THSCDF 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 % sc = 0.25;
026 % hc = 3;
027 % lowerTail = 0;
028 % upperTail = ~lowerTail
029 % thscdf(hc,sc,Hm0,Tp)           % Prob(Hd<Hc,Scf<Ec)
030 % thscdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>Hc,Scf>Ec)
031 %
032 %  % Conditional probability of steep and high waves given seastates
033 %  % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz)
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 = thscdf(hc,sc,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, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp);
069 else
070   [errorcode, Scf,Hd] = comnsize(Scf,Hd);
071 end
072 if errorcode > 0
073   error('Requires non-scalar arguments to match in size.');
074 end
075
076 global THSPAR
077 if isempty(THSPAR)
080 end
081
082 Tpp  = THSPAR.Tp;
083 Hm00 = THSPAR.Hm0;
084 Tm020 = THSPAR.Tm02;
085 % Interpolation method
086 method = '*cubic';% Faster interpolation
087
088 [Tp1,Hs1] = meshgrid(Tpp,Hm00);
089 if 1,
090   Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method);
091 else
092   Tm02 = Tp;
093   for ix= 1:100
094     dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method));
095     Tp = Tp+dTp;
096     if all(abs(dTp)<0.01)
097       %dTp
098       ix
099       break
100     end
101   end
102 end
103 %  w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing
104   %w    = linspace(0,10,2*1024+1).';
105 %  St = torsethaugen(w,[Hm0,Tp]);
106 %  ch   = spec2char(St,{'Tm02','eps2'});
107 %  Tm02 = ch(1);
108 %  eps2 = ch(2);
109 Hrms = Hm0/sqrt(2);
110 Erms = 1.25*Hm0./(Tm02.^2); % Erms
111
112 s = Scf./Erms;
113 hMax = 5;
114 h = min(Hd./Hrms,hMax);
115
116 eps2 = 1e-6;
117
118 f = zeros(size(Hd));
119
120 % Only compute within valid range
121 %k0 = find((2<=Tp) & (Tp<=21) & (Hm0<=(Tp-2)*12/11) & (Hm0<=12));
122 upLimit = 6.5;
123 loLimit = 2.5;
124 k0 = find((2<=Tp) & (Tp<=21) & (loLimit*sqrt(Hm0)<Tp) & (Hm0<=12));
125 if any(k0)
126   if multipleSeaStates
127     h = h(k0);
128     s = s(k0);
129     Hm0 = Hm0(k0);
130     Tp = Tp(k0);
131   else
132     k0 = 1:prod(size(Hd));
133   end
134
135   hlim    = h;
136
137   normalizedInput = 1;
138   lowerTail = 0;
139   if 0
140     % This is a trick to get the html documentation correct.
141     k = thspdf(1,1,2,3);
142   end
143
144   if tail==lowerTail
145     k       = find(h>2.5);
146     hlim(k) = 2.5;
147     f(k0) = gaussq('thspdf',0,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,5)...
148     + gaussq('thspdf',hlim,h,eps2/2,[],s,Hm0,Tp,normalizedInput,5);
149   else % upper tail
150     k       = find(h<2.5);
151     hlim(k) = 2.5;
152     f(k0) = gaussq('thspdf',h,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,7)...
153     + gaussq('thspdf',hlim,hMax,eps2/2,[],s,Hm0,Tp,normalizedInput,7);
154   end
155 end
156 return
157
158```

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