ohhspdf

 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) 
  
  See also  ohhpdf, thspdf

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 % 
035 % See also  ohhpdf, thspdf 
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) 
074   OHHSPAR = load('ohhspar.mat'); 
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