Home > wafo > wavemodels > ohhsspdf.m

ohhsspdf

PURPOSE

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

SYNOPSIS

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

DESCRIPTION

``` OHHSSPDF Joint (Scf,Hd) PDF linear waves in space with Ochi-Hubble spectra.

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

f   = pdf evaluated at (Scf,Hd).
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

OHHSSPDF approximates the joint distribution of (Scf, Hd), i.e., crest
front steepness (Ac/Lcf) and wave height in space, for a Gaussian
process with a bimodal Ochi-Hubble spectral density. The empirical
parameters of the model is fitted by least squares to simulated
(Scf,Hd) data for 24 classes of Hm0.
Between 50000 and 300000 zero-downcrossing waves were simulated for
each class of Hm0.
OHHSSPDF 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,4*2*Hm0/Tp)';
[V,H] = meshgrid(v,h);
f = ohhsspdf(H,V,Hm0,def);
w = linspace(0,10,2*1024+1).';
S = ohspec2(w,[Hm0 def]);
Sk = spec2spec(specinterp(S,.55),'k1d');
dk = 1;
x = spec2sdat(Sk,80000,dk); rate = 8;
[vi,hi] = dat2steep(x,rate,1);
fk = kdebin([vi,hi],'epan',[],[],.5,128);
plot(vi,hi,'.'), hold on
contour(v,h,f,fk.cl,'g'),
pdfplot(fk,'r'), hold off

CROSS-REFERENCE INFORMATION

This function calls:
 createpdf PDF class constructor ohhgparfun Wave height, Hd, distribution parameters for Ochi-Hubble spectra. 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. 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:
 ohhsscdf Joint (Scf,Hd) CDF for linear waves in space with Ochi-Hubble spectra. ohhsspdf2 Joint (Scf,Hd) PDF linear waves in space with Ochi-Hubble spectra.

SOURCE CODE

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

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