Home > wafo > wavemodels > thspdf.m

# thspdf

## PURPOSE

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

## SYNOPSIS

[f,Hrms,Vrms] = thspdf(Hd,Scf,Hm0,Tp,normalizedInput,condon)

## DESCRIPTION

``` THSPDF Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra.

CALL: f = thspdf(Hd,Scf,Hm0,Tp)

f  = pdf
Hd  = zero down crossing wave height
Scf = crest front steepness
Hm0 = significant wave height
Tp  = Spectral peak period

THSPDF approximates the joint PDF of (Scf, Hd), i.e., crest
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.
THSPDF 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.
The size of f is the common size of the input arguments.

Example:
Hm0 = 6;Tp = 8;
h = linspace(0,4*Hm0/sqrt(2));
s = linspace(0,6*1.25*Hm0/Tp^2,101);
[S,H] = meshgrid(s,h);
f = thspdf(H,S,Hm0,Tp);
contourf(s,h,f)

## CROSS-REFERENCE INFORMATION

This function calls:
 createpdf PDF class constructor range Calculates the difference between the maximum and minimum values. smooth Calculates a smoothing spline. thwparfun Wave height, Hd, distribution parameters for Torsethaugen spectra. wgamcdf Gamma cumulative distribution function wgampdf Gamma probability density function wtweibpdf Truncated 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). interp3 3-D interpolation (table lookup). meshgrid X and Y arrays for 3-D plots. unique Set unique.
This function is called by:
 thscdf Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra. thspdf2 Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra.

## SOURCE CODE

```001 function [f,Hrms,Vrms] = thspdf(Hd,Scf,Hm0,Tp,normalizedInput,condon)
002 %THSPDF Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra.
003 %
004 %  CALL: f = thspdf(Hd,Scf,Hm0,Tp)
005 %
006 %    f  = pdf
007 %   Hd  = zero down crossing wave height
008 %   Scf = crest front steepness
009 %   Hm0 = significant wave height
010 %   Tp  = Spectral peak period
011 %
012 % THSPDF approximates the joint PDF of (Scf, Hd), i.e., crest
013 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process with a
014 % Torsethaugen spectral density. The empirical parameters of the model is
015 % fitted by least squares to simulated (Scf,Hd) data for 600 classes of
016 % Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves were
017 % simulated for each class of Hm0 and Tp.
018 % THSPDF is restricted to the following range for Hm0 and Tp:
019 %  0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11.
020 %  The size of f is the common size of the input arguments.
021 %
022 % Example:
023 % Hm0 = 6;Tp = 8;
024 % h = linspace(0,4*Hm0/sqrt(2));
025 % s = linspace(0,6*1.25*Hm0/Tp^2,101);
026 % [S,H] = meshgrid(s,h);
027 % f = thspdf(H,S,Hm0,Tp);
028 % contourf(s,h,f)
029 %
031
032 % Reference
033 % P. A. Brodtkorb (2004),
034 % The Probability of Occurrence of Dangerous Wave Situations at Sea.
035 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
036 % Trondheim, Norway.
037
038 % History
039 % revised pab 09.08.2003
040 %  changed input and help header
041 % validated 20.11.2002
042 % By pab 20.12.2000
043
044 error(nargchk(3,6,nargin))
045
046
047 if (nargin < 6|isempty(condon)),  condon  = 0; end
048 if (nargin < 5|isempty(normalizedInput)),  normalizedInput  = 0;end
049 if (nargin < 4|isempty(Tp)),  Tp  = 8;end
050 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end
051
052 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1);
053 if multipleSeaStates
054   [errorcode, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp);
055 else
056   [errorcode, Scf,Hd] = comnsize(Scf,Hd);
057 end
058 if errorcode > 0
059     error('Requires non-scalar arguments to match in size.');
060 end
061 displayWarning = 0;
062 if displayWarning,
063   if any(Hm0>11| Hm0>(Tp-2)*12/11)
064     disp('Warning: Hm0 is outside the valid range')
065     disp('The validity of the Joint (Hd,Scf) distribution is questionable')
066   end
067   if any(Tp>20|Tp<3)
068     disp('Warning: Tp is outside the valid range')
069     disp('The validity of the Joint (Hd,Scf) distribution is questionable')
070   end
071 end
072
073 global THSPAR
074 if isempty(THSPAR)
077 end
078
079 Tpp  = THSPAR.Tp;
080 Hm00 = THSPAR.Hm0;
081 Tm020 = THSPAR.Tm02;
082 h2   = THSPAR.h2(:);
083
084 % Interpolation method
085 method = '*cubic';% Faster interpolation
086
087 if normalizedInput,
088   Hrms = 1;
089   Vrms = 1;
090 else
091   [Tp1,Hs1] = meshgrid(Tpp,Hm00);
092   if 1,
093     Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method);
094     %Tp/Tm02
095   else
096
097     Tm02 = Tp;
098     for ix= 1:100
099       dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method));
100       Tp = Tp+dTp;
101       if all(abs(dTp)<0.01)
102     dTp
103     ix
104     break
105       end
106     end
107     %Tp,    Tp/Tm02
108   end
109
110 %  w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing
111   %w    = linspace(0,10,2*1024+1).';
112 %  St = torsethaugen(w,[Hm0,Tp]);
113 %  ch   = spec2char(St,{'Tm02','eps2'});
114 %  Tm02 = ch(1);
115 %  eps2 = ch(2);
116   Hrms = Hm0/sqrt(2);
117   Vrms = 1.25*Hm0./(Tm02.^2); % Erms
118 end
119
120 h = Hd./Hrms;
121 v = Scf./Vrms;
122 cSize = size(h); % common size
123
124 % Gamma distribution parameters as a function of Tp Hm0 and h2
125 A11 = THSPAR.A11s;
126 B11 = THSPAR.B11s;
127
128 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2);
129 Nh2 = length(h2);
130
131 if multipleSeaStates
132   h   = h(:);
133   v   = v(:);
134   Tp  = Tp(:);
135   Hm0 = Hm0(:);
136   A1 = zeros(length(h),1);
137   B1 = A1;
138   [TpHm0,ix,jx] = unique([Tp,Hm0],'rows');
139   numSeaStates = length(ix);
140   Tpi = zeros(Nh2,1);
141   Hm0i = zeros(Nh2,1);
142   for iz=1:numSeaStates
143     k = find(jx==iz);
144     Tpi(:)  = TpHm0(iz,1);
145     Hm0i(:) = TpHm0(iz,2);
146     A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),...
147                1,h(k),1));
148     B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),...
149                1,h(k),1));
150   end
151 else
152   Tpi  = repmat(Tp,[Nh2,1]);
153   Hm0i = repmat(Hm0,[Nh2,1]);
154   A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1));
155   B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1));
156 end
157 % Waveheight distribution in time
158 % Truncated Weibull  distribution parameters as a function of Tp, Hm0
159 [A0, B0, C0] = thwparfun(Hm0,Tp,'time');
160 switch condon,
161  case 0, % regular pdf is returned
162   f = wtweibpdf(h,A0,B0,C0).*wgampdf(v,A1,B1);
163  case 1, %pdf conditioned on x1 ie. p(x2|x1)
164   f = wgampdf(v,A1,B1);
165  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)
166   f = v.*wgampdf(v,A1,B1);
167  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)
168   f = v.^2.*wgampdf(v,A1,B1);
169  case 5, % p(h)*P(V|h) is returned special case used by thscdf
170   f = wtweibpdf(h,A0,B0,C0).*wgamcdf(v,A1,B1);
171  case 6, % P(V|h) is returned special case used by thscdf
172   f = wgamcdf(v,A1,B1);
173  case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf
174   f = wtweibpdf(h,A0,B0,C0).*(1-wgamcdf(v,A1,B1));
175  otherwise error('unknown option')
176 end
177 if multipleSeaStates
178   f = reshape(f,cSize);
179 end
180
181 if condon~=6
182   f = f./Hrms./Vrms;
183 end
184 f(find(isnan(f)|isinf(f) ))=0;
185 if any(size(f)~=cSize)
186   disp('Wrong size')
187 end
188
189 if nargout>3,
190   fA      = createpdf(2);
191   fA.x    = {Tpp,Hm00};
192   fA.labx = {'Tp', 'Hm0'};
193   fA(3)   = fA(1);
194   fA(2)   = fA(1);
195
196   % Truncated Weibull  distribution parameters as a function of Tp, Hm0
197   A00 = THSPAR.A00s;
198   B00 = THSPAR.B00s;
199   C00 = THSPAR.C00s;
200
201   fA(1).f    = A00;
202   fA(2).f    = B00;
203   fA(3).f    = C00;
204
205   fA(1).title = 'wtweibpdf parameter A';
206   fA(2).title = 'wtweibpdf parameter B';
207   fA(3).title = 'wtweibpdf parameter C';
208
209   txt1 = ['The Wtweibpdf  distribution Parameter '];
210   txt2=[' for wave heigth in time as a function of Tp and Hm0 for' ...
211     'the Torsethaugen spectrum'];
212   fA(1).note =[txt1 'A' txt2];
213   fA(2).note =[txt1 'B' txt2];
214   fA(3).note =[txt1 'C' txt2];
215
216   tmp = [A00(:) B00(:) C00(:)];
217   ra  = range(tmp);
218   st  = round(min(tmp)*100)/100;
219   en  = max(tmp);
220   for ix = 1:3,
221     fA(ix).cl   = st(ix):ra(ix)/20:en(ix);
222   end
223 end
224 if nargout>4,
225   fB      = createpdf(3);
226   fB.x    = {Tpp,Hm00,h2};
227   fB.labx = {'Tp','Hm0', 'h'};
228   fB(2)   = fB(1);
229
230   fB(1).f = A11;
231   fB(2).f = B11;
232
233   txt11 = 'The conditonal Wgampdf distribution Parameter ';
234   txt22 = [' for Scf given h=Hd/Hrms in time as function of Tp' ...
235     ' and Hm0 for the Torsethaugen spectrum'];
236   fB(1).title = 'wgampdf parameter A';
237   fB(2).title = 'wgampdf parameter B';
238   fB(1).note = [txt11,'A',txt22];
239   fB(2).note = [txt11,'B',txt22];
240
241   %fB(2).note = ['The conditonal Wggampdf distribution Parameter B(h)/Hrms', ...%    ' for crest front steepness as a function of Tp,Hm0 and',...
242   %    ' h=Hd/Hrms for the Torsethaugen spectrum'];
243   tmp= [A11(:) B11(:)];
244   ra = range(tmp);
245   st = round(min(tmp)*100)/100;
246   en = max(tmp);
247   for ix=1:2
248     fB(ix).cl   = st(ix):ra(ix)/20:en(ix);
249   end
250 end
251 return
252
253
254
255
256```

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