Home > wafo > wavemodels > ohhspdf.m

# ohhspdf

## PURPOSE

Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra.

## SYNOPSIS

[f,Hrms,Vrms,fA,fB] = ohhspdf(Hd,Scf,Hm0,def,normalizedInput, condon)

## DESCRIPTION

``` OHHSPDF Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra.

CALL: f = ohhspdf(Hd,Scf,Hm0,def)

f   = pdf
Hd  = zero down crossing wave height
Scf = crest front steepness
Hm0 = significant wave height [m].
def = defines the parametrization of the spectral density (default 1)
1 : The most probable spectrum  (default)
2,3,...11 : gives 95% Confidence spectra

OHHSPDF approximates the joint distribution of (Scf, Hd) in time,
i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height,
for a Gaussian process with a bimodal Ochi-Hubble spectral density
(ohspec2). The empirical
parameters of the model is fitted by least squares to simulated
(Scf,Hd) data for 24 classes of Hm0. Between 50000 and 150000
zero-downcrossing waves were simulated for each class of Hm0.
OHHSPDF is restricted to the following range for Hm0:
0.5 < Hm0 [m] < 12
The size of f is the common size of the input arguments, Hd, Scf and
Hm0.

Example:
Hm0 = 6;Tp = 8;def= 2;
h = linspace(0,4*Hm0/sqrt(2))';
v = linspace(0,6*1.25*Hm0/Tp^2)';
[V,H]  = meshgrid(v,h);
f = ohhspdf(H,V,Hm0,def);
contourf(v,h,f)

## CROSS-REFERENCE INFORMATION

This function calls:
 createpdf PDF class constructor ohhgparfun Wave height, Hd, distribution parameters for Ochi-Hubble spectra. qlevels Calculates quantile levels which encloses P% of PDF range Calculates the difference between the maximum and minimum values. smooth Calculates a smoothing spline. wggampdf Generalized Gamma probability density function wweibcdf Weibull cumulative distribution function wweibpdf Weibull probability density function comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. interp1 1-D interpolation (table lookup) interp2 2-D interpolation (table lookup). meshgrid X and Y arrays for 3-D plots. num2str Convert number to string. (Fast version) sprintf Write formatted data to string. squeeze Remove singleton dimensions. unique Set unique. warning Display warning message; disable or enable warning messages.
This function is called by:
 ohhscdf Joint (Scf,Hd) CDF for linear waves with Ochi-Hubble spectra. ohhspdf2 Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra.

## SOURCE CODE

```001 function [f,Hrms,Vrms,fA,fB] = ohhspdf(Hd,Scf,Hm0,def,normalizedInput, ...
002                        condon)
003 %OHHSPDF Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra.
004 %
005 %  CALL: f = ohhspdf(Hd,Scf,Hm0,def)
006 %
007 %  f   = pdf
008 %  Hd  = zero down crossing wave height
009 %  Scf = crest front steepness
010 %  Hm0 = significant wave height [m].
011 %  def = defines the parametrization of the spectral density (default 1)
012 %        1 : The most probable spectrum  (default)
013 %        2,3,...11 : gives 95% Confidence spectra
014 %
015 % OHHSPDF approximates the joint distribution of (Scf, Hd) in time,
016 % i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height,
017 %  for a Gaussian process with a bimodal Ochi-Hubble spectral density
018 % (ohspec2). The empirical
019 % parameters of the model is fitted by least squares to simulated
020 % (Scf,Hd) data for 24 classes of Hm0. Between 50000 and 150000
021 % zero-downcrossing waves were simulated for each class of Hm0.
022 % OHHSPDF is restricted to the following range for Hm0:
023 % 0.5 < Hm0 [m] < 12
024 % The size of f is the common size of the input arguments, Hd, Scf and
025 % Hm0.
026 %
027 % Example:
028 % Hm0 = 6;Tp = 8;def= 2;
029 % h = linspace(0,4*Hm0/sqrt(2))';
030 % v = linspace(0,6*1.25*Hm0/Tp^2)';
031 % [V,H]  = meshgrid(v,h);
032 % f = ohhspdf(H,V,Hm0,def);
033 % contourf(v,h,f)
034 %
036
037 % Reference
038 % P. A. Brodtkorb (2004),
039 % The Probability of Occurrence of Dangerous Wave Situations at Sea.
040 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU,
041 % Trondheim, Norway.
042
043 % History
044 % revised pab 09.09.2003
045 % By pab 06.02.2001
046
047 error(nargchk(3,6,nargin))
048 if (nargin < 6|isempty(condon)),  condon  = 0;end
049 if (nargin < 5|isempty(normalizedInput)),  normalizedInput  = 0;end
050 if (nargin < 4|isempty(def)), def=1;end
051
052 multipleSeaStates = any(prod(size(Hm0))>1);
053 if multipleSeaStates
054   [errorcode, Scf,Hd,Hm0] = comnsize(Scf,Hd,Hm0);
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
062 if any(Hm0>12| Hm0<=0.5)
063   disp('Warning: Hm0 is outside the valid range')
064   disp('The validity of the Hd distribution is questionable')
065 end
066
067 if def>11|def<1
068   Warning('DEF is outside the valid range')
069   def = mod(def-1,11)+1;
070 end
071
072 global OHHSPAR
073 if isempty(OHHSPAR)
075 end
076
077
078 %w    = linspace(0,100,16*1024+1).'; %original spacing
079 %ch   = spec2char(ohspec2(w,[Hm0,def]),{'Tm02','eps2'});
080 %Tm02 = ch(1);
081 %Vrms = 2*Hm0/Tm02;
082 Hm00  = OHHSPAR.Hm0;
083
084 if normalizedInput,
085   Hrms = 1;
086   Vrms = 1;
087 else
088   method = 'cubic'; %'spline'
089   Tm020 = OHHSPAR.Tm02;
090   Hrms = Hm0/sqrt(2);
091   Tm02 = interp1(Hm00,Tm020(:,def),Hm0,method);
092   Vrms = 1.25*Hm0./(Tm02.^2); % Erms
093 end
094
095
096 h = Hd./Hrms;
097 v = Scf./Vrms;
098 cSize = size(h); % common size
099
100 % Logarithm of Weibull distribution parameters as a function of Hm0 and h2
101 A11  = squeeze(OHHSPAR.A11s(:,def,:));
102 B11  = squeeze(OHHSPAR.B11s(:,def,:));
103
104 h2    = OHHSPAR.h2; %Hd/Hrms
105 [E2 H2] = meshgrid(Hm00,h2);
106 method = '*cubic'; %'spline'
107 if multipleSeaStates
108   h   = h(:);
109   v   = v(:);
110   Hm0 = Hm0(:);
111   A1 = zeros(length(h),1);
112   B1 = A1;
113   [Hm0u,ix,jx] = unique(Hm0);
114   numSeaStates = length(ix);
115   Nh2 = length(h2);
116   Hm0i = zeros(Nh2,1);
117   for iz=1:numSeaStates
118     k = find(jx==iz);
119     Hm0i(:) = Hm0u(iz);
120     A1(k) = exp(smooth(h2,interp2(E2,H2,A11.',Hm0i,h2,method),...
121                1,h(k),1));
122     B1(k) = (smooth(h2,interp2(E2,H2,exp(B11).',Hm0i,h2,method),...
123                1,h(k),1));
124   end
125 else
126   Hm0i = Hm0(ones(size(h2)));
127   A1 = exp(smooth(h2,interp2(E2,H2,A11.', ...
128                  Hm0i,h2,method),1,h,1));
129   B1 = (smooth(h2,interp2(E2,H2,exp(B11).', ...
130               Hm0i,h2,method),1,h,1));
131 end
132 %fh = ohhpdf(h(:)/Hrms,Hm0,def,'time',1);
133 % Fh = fh.f;
134 % Waveheight distribution in time
135 % Generalized gamma distribution parameters as a function of Hm0
136 [A0 B0 C0] = ohhgparfun(Hm0,def,'time');
137
138
139 switch condon,
140  case 0, % regular pdf is returned
141   f = wggampdf(h,A0,B0,C0).*wweibpdf(v,A1,B1);
142  case 1, %pdf conditioned on x1 ie. p(x2|x1)
143   f = wweibpdf(v,A1,B1);
144  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)
145   f = v.*wweibpdf(v,A1,B1);
146  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)
147   f = v.^2.*wweibpdf(v,A1,B1);
148  case 5, % p(h)*P(V|h) is returned special case used by thscdf
149   f = wggampdf(h,A0,B0,C0).*wweibcdf(v,A1,B1);
150  case 6, % P(V|h) is returned special case used by thscdf
151   f = wweibcdf(v,A1,B1);
152  case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf
153   f = wggampdf(h,A0,B0,C0).*(1-wweibcdf(v,A1,B1));
154  otherwise error('unknown option')
155 end
156 if multipleSeaStates
157   f = reshape(f,cSize);
158 end
159
160 if condon~=6
161   f = f./Hrms./Vrms;
162 end
163 f(find(isnan(f)|isinf(f) ))=0;
164 if any(size(f)~=cSize)
165   disp('Wrong size')
166 end
167 if nargout>3,
168   fA      = createpdf(2);
169   fA.f    = A11;
170   fA.x    = {Hm00, h2};
171   fA.labx = {'Hm0', 'h'};
172   fA.note = sprintf('%s \n',...
173             'The conditonal Weibull distribution Parameter',...
174             'A(h)/Hrms for wave heigth as a function of',...
175             'h=Hd/Hrms and Hm0 for the ohspec2 spectrum, ',...
176             sprintf('given def = %d.',def));
177
178   ra = range(A11(:));
179   st = round(min(A11(:))*100)/100;
180   en = max(A11(:));
181   fA.cl   = st:ra/20:en;
182 end
183 if nargout>4,
184   fB      = createpdf(2);
185   fB.f    = B11;
186   fB.x    = {Hm00,h2};
187   fB.labx = {'Hm0','h'};
188   fB.note =  sprintf('%s \n',...
189             'The conditonal Weibull distribution Parameter',...
190             'B(h)/Hrms for wave heigth as a function of',...
191             'h=Hd/Hrms and Hm0 for the ohspec2 spectrum, ',...
192             sprintf('given def = %d.',def));
193   ra = range(B11(:));
194   st = round(min(B11(:))*100)/100;
195   en = max(B11(:));
196   fB.cl   = st:ra/20:en;
197 end
198 return
199
200
201
202
203
204
205
206 return % old call
207 Fh = wggampdf(h(:)/Hrms,A0,B0,C0);
208
209
210
211 %C1 = exp(smooth(h2,interp2(E2,H2,log(squeeze(C11(:,def,:))).',Hm0(ones(size(h2))),h2,method),1,h/Hrms,1));
212 %figure(1), plot(h/Hrms, A1,h2,squeeze(exp(A11(10,def,:))),'r')
213 %figure(2), plot(h/Hrms,B1,h2,squeeze(exp(B11(10,def,:))),'r')
214 [V1 A1] = meshgrid(v/Vrms,A1);
215 [V1 B1] = meshgrid(v/Vrms,B1);
216 %[V1 C1] = meshgrid(v/Vrms,C1);
217
218 f = createpdf(2);
219 f.title = 'Joint distribution of (Hd,Scf) in time';
220 if norm==0
221   f.f = repmat(Fh/Hrms,[1 length(v)]).*wweibpdf(V1,A1,B1)/Vrms;
222   f.x = {v,h};
223   f.norm=0;
224   f.labx={'Scf', 'Hd [m]'};
225 else
226   f.f = repmat(Fh,[1 length(v)]).*wweibpdf(V1,A1,B1);
227   f.x = {v/Vrms,h/Hrms};
228   f.labx={'Scf', 'Hd'};
229   f.norm = 1;
230 end
231 f.note = ['Ochi-Hubble Hm0=' num2str(Hm0) ' def = ' num2str(def)];
232 if 1,
233   f.f(isinf(f.f)|isnan(f.f))=0;
234   % Some places numerical problems occur
235   %f.f(f.f>30)=0; % Needed to calculate the contour levels
236   [f.cl,f.pl] = qlevels(f.f);%,[10:20:90 95 99 99.9 99.99 99.999 99.9999]);
237 end
238```

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