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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 % 
043 % See also  ohhpdf, thvpdf 
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) 
085   OHHSSPAR = load('ohhsspar.mat'); 
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