Home > wafo > wavemodels > jhvpdf.m

jhvpdf

PURPOSE ^

Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum.

SYNOPSIS ^

[f,Hrms,Vrms,fA,fB] = jhvpdf(Hd,Vcf,Hm0,Tp,gam,normalizedInput,condon)

DESCRIPTION ^

 JHVPDF Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. 
  
   CALL: f = jhvpdf(Hd,Vcf,Hm0,Tp,gamma) 
   
   f     = pdf evaluated at (Vcf,Hd) 
   Hd    = zero down crossing wave height [m] 
   Vcf   = crest front velocity    [m/s] 
   Hm0   = significant wave height [m] 
   Tp    = Spectral peak period    [s] 
   Gamma = Peakedness parameter of the JONSWAP spectrum 
  
  JHVPDF approximates the joint distribution of (Vcf, Hd), i.e., crest 
  front velocity (Ac/Tcf) and wave height, for a Gaussian process with a 
  JONSWAP spectral density. The empirical parameters of the model is 
  fitted by least squares to simulated (Vcf,Hd) data for 13 classes of 
  GAMMA between 1 and 7. About 100000 zero-downcrossing waves 
  were simulated for each class. 
  JHVPDF is restricted to the following range for GAMMA:  
   1 <= GAMMA <= 7  
  
  NOTE:  The size of f is the common size of the input arguments.  
  
  Example: 
  Hm0 = 6;Tp = 9; gam=3.5 
  h = linspace(0,4*Hm0/sqrt(2))';  
  v = linspace(0,4*2*Hm0/Tp)'; 
  [V,H] = meshgrid(v,h);   
  f = jhvpdf(H,V,Hm0,Tp,gam); 
  w = linspace(0,40,5*1024+1).'; 
  S = jonswap(w,[Hm0, Tp, gam]); 
  dt = .3; 
  x = spec2sdat(S,80000,dt); rate = 4; 
  [vi,hi] = dat2steep(x,rate,1); 
  fk = kdebin([vi,hi],'epan',[],[],.5,128);  
  plot(vi,hi,'.'), hold on 
  contour(v,h,f,fk.cl), 
  pdfplot(fk,'r'),hold off 
  
  See also  thvpdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [f,Hrms,Vrms,fA,fB] = jhvpdf(Hd,Vcf,Hm0,Tp,gam,normalizedInput,condon) 
002 %JHVPDF Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. 
003 % 
004 %  CALL: f = jhvpdf(Hd,Vcf,Hm0,Tp,gamma) 
005 %  
006 %  f     = pdf evaluated at (Vcf,Hd) 
007 %  Hd    = zero down crossing wave height [m] 
008 %  Vcf   = crest front velocity    [m/s] 
009 %  Hm0   = significant wave height [m] 
010 %  Tp    = Spectral peak period    [s] 
011 %  Gamma = Peakedness parameter of the JONSWAP spectrum 
012 % 
013 % JHVPDF approximates the joint distribution of (Vcf, Hd), i.e., crest 
014 % front velocity (Ac/Tcf) and wave height, for a Gaussian process with a 
015 % JONSWAP spectral density. The empirical parameters of the model is 
016 % fitted by least squares to simulated (Vcf,Hd) data for 13 classes of 
017 % GAMMA between 1 and 7. About 100000 zero-downcrossing waves 
018 % were simulated for each class. 
019 % JHVPDF is restricted to the following range for GAMMA:  
020 %  1 <= GAMMA <= 7  
021 % 
022 % NOTE:  The size of f is the common size of the input arguments.  
023 % 
024 % Example: 
025 % Hm0 = 6;Tp = 9; gam=3.5 
026 % h = linspace(0,4*Hm0/sqrt(2))';  
027 % v = linspace(0,4*2*Hm0/Tp)'; 
028 % [V,H] = meshgrid(v,h);   
029 % f = jhvpdf(H,V,Hm0,Tp,gam); 
030 % w = linspace(0,40,5*1024+1).'; 
031 % S = jonswap(w,[Hm0, Tp, gam]); 
032 % dt = .3; 
033 % x = spec2sdat(S,80000,dt); rate = 4; 
034 % [vi,hi] = dat2steep(x,rate,1); 
035 % fk = kdebin([vi,hi],'epan',[],[],.5,128);  
036 % plot(vi,hi,'.'), hold on 
037 % contour(v,h,f,fk.cl), 
038 % pdfplot(fk,'r'),hold off 
039 % 
040 % See also  thvpdf 
041  
042 % Reference   
043 % P. A. Brodtkorb (2004),   
044 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
045 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
046 % Trondheim, Norway.  
047    
048 % History 
049 % revised pab 10 Jan 2004   
050 % By pab 20.12.2000 
051  
052 error(nargchk(4,7,nargin)) 
053 if (nargin < 7|isempty(condon)),  condon  = 0; end 
054 if (nargin < 6|isempty(normalizedInput)),  normalizedInput  = 0;end 
055 if (nargin < 5|isempty(gam)) 
056    gam = getjonswappeakedness(Hm0,Tp); 
057 end 
058  
059 multipleSeaStates = any(prod(size(Hm0))>1|... 
060             prod(size(Tp)) >1|... 
061             prod(size(gam))>1); 
062 if multipleSeaStates 
063   [errorcode, Vcf,Hd,Hm0,Tp,gam] = comnsize(Vcf,Hd,Hm0,Tp,gam); 
064 else 
065   [errorcode, Vcf,Hd] = comnsize(Vcf,Hd); 
066 end 
067 if errorcode > 0 
068   error('Requires non-scalar arguments to match in size.'); 
069 end 
070  
071 displayWarning = 0; 
072 if displayWarning 
073   if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0))) 
074     disp('Warning: Hm0,Tp is outside the JONSWAP range') 
075     disp('The validity of the parameters returned are questionable') 
076   end 
077 end 
078 k = find(gam<1); 
079 if any(k)  
080   if displayWarning, 
081     warning('Peakedness parameter less than 1. Must be larger than 1.') 
082   end 
083   gam(k)=1; 
084 end 
085 k1 = find(7<gam); 
086 if any(k1) 
087   if displayWarning 
088   warning('Peakedness parameter larger than 7. The pdf returned is questionable') 
089   end 
090   gam(k1) = 7; 
091 end 
092  
093 global JHVPAR 
094 if isempty(JHVPAR) 
095   JHVPAR = load('jhvpar.mat'); 
096 end 
097 % Weibull distribution parameters as a function of e2 and h2 
098 A11 = JHVPAR.A11s; 
099 B11 = JHVPAR.B11s; 
100 e2  = JHVPAR.gam; % gamma 
101 h2  = JHVPAR.h2;  % Hd/Hrms 
102 [E2 H2] = meshgrid(e2,h2); 
103  
104 if 1, 
105    
106   %Tm02 = Tp./(1.30301-0.01698*gam+0.12102./gam); 
107   %dev = 2e-5; 
108   c1 =[ 0.16183666835624   1.53691936441548   1.55852759524555]; 
109   c2 =[ 0.15659478203944   1.15736959972513   1]; 
110   Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam)); 
111    
112    %dev = 2e-4; 
113   %EPS2cof = [0.00068263671017  -0.01802256231624   0.44176198490431]; 
114   %eps2 = polyval(EPS2cof,gam); 
115 else 
116   w    = linspace(0,100,16*1024+1).'; % jonswap original spacing 
117   %Hm0 = 6; 
118   %gam = linspace(1,7,32); 
119   Tm02 = zeros(size(gam)); 
120   eps2 = Tm02; 
121   for ix=1:length(gam(:)) 
122     ch   = spec2char(jonswap(w,[Hm0(ix),Tp(ix) gam(ix)]),{'Tm02','eps2'}); 
123     Tm02(ix) = ch(1); 
124     eps2(ix) = ch(2); 
125   end 
126 end 
127  
128  
129 if normalizedInput 
130   Hrms = 1; 
131   Vrms = 1; 
132 else 
133   Hrms = Hm0/sqrt(2); 
134   Vrms = 2*Hm0./Tm02; % Erms 
135 end 
136  
137 v = Vcf./Vrms; 
138 h = Hd./Hrms; 
139 cSize = size(h); % common size of input 
140  
141  
142  
143 method ='*cubic'; 
144 Nh2 = length(h2); 
145 if multipleSeaStates 
146   h   = h(:); 
147   v   = v(:); 
148   Tp  = Tp(:); 
149   Hm0 = Hm0(:); 
150   gam = gam(:); 
151 %  eps2 = eps2(:); 
152   A1 = zeros(length(h),1); 
153   B1 = A1; 
154   [gamu,ix,jx] = unique(gam); 
155   numSeaStates = length(gamu); 
156   gami = zeros(Nh2,1); 
157   for iz=1:numSeaStates 
158     k = find(jx==iz); 
159     gami(:) = gamu(iz); 
160     A1(k) = smooth(h2,interp2(E2,H2,A11,gami,h2,method),1,h(k),1); 
161     B1(k) = smooth(h2,interp2(E2,H2,B11,gami,h2,method),1,h(k),1); 
162   end 
163 else 
164   A1 = smooth(h2,interp2(E2,H2,A11,gam(ones(size(h2))),h2,method),1,h,1); 
165   B1 = smooth(h2,interp2(E2,H2,B11,gam(ones(size(h2))),h2,method),1,h,1); 
166 end 
167  
168 % Waveheight distribution in time 
169 % Truncated Weibull  distribution parameters as a function of Tp, Hm0, gam  
170 [A0, B0, C0] = jhwparfun(Hm0,Tp,gam,'time'); 
171  
172 switch condon, 
173  case 0, % regular pdf is returned  
174   f = wtweibpdf(h,A0,B0,C0).*wweibpdf(v,A1,B1); 
175  case 1, %pdf conditioned on x1 ie. p(x2|x1)  
176   f = wweibpdf(v,A1,B1); 
177  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
178   f = v.*wweibpdf(v,A1,B1); 
179  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
180   f = v.^2.*wweibpdf(v,A1,B1); 
181  case 5, % p(h)*P(V|h) is returned special case used by jhvcdf2 
182   f = wtweibpdf(h,A0,B0,C0).*wweibcdf(v,A1,B1); 
183  case 6, % P(V|h) is returned special case used by jhvcdf2 
184   f = wweibcdf(v,A1,B1); 
185  case 7,% p(h)*(1-P(V|h)) is returned special case used by jhvcdf2 
186   f = wtweibpdf(h,A0,B0,C0).*(1-wweibcdf(v,A1,B1)); 
187  otherwise error('unknown option') 
188 end 
189  
190 if multipleSeaStates 
191   f = reshape(f,cSize); 
192 end 
193  
194 if condon~=6 
195   f = f./Hrms./Vrms; 
196 end 
197 f(find(isnan(f)|isinf(f) ))=0; 
198 if any(size(f)~=cSize) 
199   disp('Wrong size') 
200 end 
201  
202  
203 if nargout>3, 
204   fA      = createpdf(2); 
205   fA.f    = A11; 
206   fA.x    = {e2,h2}; 
207   fA.labx = {'Gamma', 'h'}; 
208   fA.note = ['The conditonal Weibull distribution Parameter A(h,gamma)/Hrms' ... 
209     'for Vcf as a function of h=Hd/Hrms and gamma for' ... 
210     'the Jonswap spectrum']; 
211      
212   ra = range(A11(:)); 
213   st = round(min(A11(:))*100)/100; 
214   en = max(A11(:)); 
215   fA.cl   = st:ra/20:en; 
216 end 
217 if nargout>4, 
218   fB      = createpdf(2); 
219   fB.f    = B11; 
220   fB.x    = {e2,h2}; 
221   fB.labx = {'Gamma', 'h'}; 
222   fB.note = ['The conditonal Weibull distribution Parameter B(h,gamma)/Hrms' ... 
223     'for Vcf as a function of h=Hd/Hrms and gamma for' ... 
224     'the Jonswap spectrum']; 
225   ra = range(B11(:)); 
226   st = round(min(B11(:))*100)/100; 
227   en = max(B11(:)); 
228   fB.cl   = st:ra/20:en; 
229 end 
230 return 
231  
232  
233  
234  
235  
236  
237

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