Home > wafo > wavemodels > jhvnlpdf.m

jhvnlpdf

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 JHVNLPDF Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. 
  
   CALL: f = jhvnlpdf(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 
  
  JHVNLPDF approximates the joint distribution of (Vcf, Hd), i.e., crest 
  front velocity (Ac/Tcf) and wave height, for 2nd order Stokes waves 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. Between 140000 and 150000 zero-downcrossing waves 
  were simulated for each class. 
  JHVNLPDF is restricted to the following range for GAMMA:  
   1 <= GAMMA <= 7  
  
  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 = jhvnlpdf(H,V,Hm0,Tp,gam); 
  contourf(v,h,f)     
  
  See also  thvpdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [f,Hrms,Vrms,fA,fB] = jhvnlpdf(Hd,Vcf,Hm0,Tp,gam,normalizedInput,condon) 
002 %JHVNLPDF Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. 
003 % 
004 %  CALL: f = jhvnlpdf(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 % JHVNLPDF approximates the joint distribution of (Vcf, Hd), i.e., crest 
014 % front velocity (Ac/Tcf) and wave height, for 2nd order Stokes waves 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. Between 140000 and 150000 zero-downcrossing waves 
018 % were simulated for each class. 
019 % JHVNLPDF is restricted to the following range for GAMMA:  
020 %  1 <= GAMMA <= 7  
021 % 
022 % Example: 
023 % Hm0 = 6;Tp = 9; gam=3.5 
024 % h = linspace(0,4*Hm0/sqrt(2))';  
025 % v = linspace(0,4*2*Hm0/Tp)'; 
026 % [V,H] = meshgrid(v,h);  
027 % f = jhvnlpdf(H,V,Hm0,Tp,gam); 
028 % contourf(v,h,f)     
029 % 
030 % See also  thvpdf 
031  
032 % Reference   
033 % P. A. Brodtkorb (2004),   
034 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
035 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
036 % Trondheim, Norway.  
037    
038 % History 
039 % By pab 20.12.2000 
040  
041 error(nargchk(4,7,nargin)) 
042 if (nargin < 7|isempty(condon)),  condon  = 0; end 
043 if (nargin < 6|isempty(normalizedInput)),  normalizedInput  = 0;end 
044 if (nargin < 5|isempty(gam)) 
045    gam = getjonswappeakedness(Hm0,Tp); 
046 end 
047  
048 multipleSeaStates = any(prod(size(Hm0))>1|... 
049             prod(size(Tp)) >1|... 
050             prod(size(gam))>1); 
051 if multipleSeaStates 
052   [errorcode, Vcf,Hd,Hm0,Tp,gam] = comnsize(Vcf,Hd,Hm0,Tp,gam); 
053 else 
054   [errorcode, Vcf,Hd] = comnsize(Vcf,Hd); 
055 end 
056 if errorcode > 0 
057   error('Requires non-scalar arguments to match in size.'); 
058 end 
059  
060 displayWarning = 0; 
061 if displayWarning 
062   if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0))) 
063     disp('Warning: Hm0,Tp is outside the JONSWAP range') 
064     disp('The validity of the parameters returned are questionable') 
065   end 
066 end 
067 k = find(gam<1); 
068 if any(k)  
069   if displayWarning, 
070     warning('Peakedness parameter less than 1. Must be larger than 1.') 
071   end 
072   gam(k)=1; 
073 end 
074 k1 = find(7<gam); 
075 if any(k1) 
076   if displayWarning 
077   warning('Peakedness parameter larger than 7. The pdf returned is questionable') 
078   end 
079   gam(k1) = 7; 
080 end 
081  
082 global JHVNLPAR 
083 if isempty(JHVNLPAR) 
084   JHVNLPAR = load('jhvnlpar.mat'); 
085 end 
086 % Weibull distribution parameters as a function of e2 and h2 
087 A11 = JHVNLPAR.A11s; 
088 B11 = JHVNLPAR.B11s; 
089  
090  
091  
092  
093 %dev = 2e-5; 
094 c1 =[ 0.16183666835624   1.53691936441548   1.55852759524555]; 
095 c2 =[ 0.15659478203944   1.15736959972513   1]; 
096 Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam)); 
097  
098 %w    = linspace(0,10.5,16*1024+1).'; % jonswap original spacing 
099 %ch   = spec2char(jonswap(w,[Hm0,Tp gam]),{'Tm02','eps2'}); 
100 %Tm02 = ch(1); 
101 %eps2 = ch(2); 
102  
103  
104 if normalizedInput 
105   Hrms = 1; 
106   Vrms = 1; 
107 else 
108   Hrms = Hm0/sqrt(2); 
109   Vrms = 2*Hm0./Tm02; % Erms 
110 end 
111  
112 v = Vcf./Vrms; 
113 h = Hd./Hrms; 
114 cSize = size(h);  
115   
116 if gam<1  
117   disp('Warning: Peakedness parameter less than 1. Must be larger than 1.') 
118     gam=1; 
119 elseif gam>7 
120   disp('Warning: Peakedness parameter larger than 7. The pdf returned is questionable') 
121 end  
122  
123  
124 h2  = JHVNLPAR.h2(:);  % Hd/Hrms 
125 method ='*cubic'; 
126 Nh2 = length(h2); 
127 if 1, 
128    
129   Tpp  = JHVNLPAR.Tp; % gamma 
130   Hm00 = JHVNLPAR.Hm0; 
131   [E1, H1, H2] = meshgrid(Tpp,Hm00,h2); 
132   Nh2 = length(h2); 
133   if multipleSeaStates 
134     h   = h(:); 
135     v   = v(:); 
136     Tp  = Tp(:); 
137     Hm0 = Hm0(:); 
138     gam = gam(:); 
139     A1 = zeros(length(h),1); 
140     B1 = A1; 
141     [TpHm0,ix,jx] = unique([Tp,Hm0],'rows'); 
142     numSeaStates = length(ix); 
143     Tpi = zeros(Nh2,1); 
144     Hm0i = zeros(Nh2,1); 
145     for iz=1:numSeaStates 
146       k = find(jx==iz); 
147       Tpi(:)  = TpHm0(iz,1); 
148       Hm0i(:) = TpHm0(iz,2); 
149       A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),... 
150              1,h(k),1)); 
151       B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),... 
152              1,h(k),1)); 
153     end 
154   else 
155     Tpi  = repmat(Tp,[Nh2,1]); 
156     Hm0i = repmat(Hm0,[Nh2,1]); 
157     A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1)); 
158     B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1)); 
159   end 
160 else 
161   e2  = JHVNLPAR.gam; % gamma 
162    
163   [E2 H2] = meshgrid(e2,h2); 
164 if multipleSeaStates 
165   h   = h(:); 
166   v   = v(:); 
167   Tp  = Tp(:); 
168   Hm0 = Hm0(:); 
169   gam = gam(:); 
170 %  eps2 = eps2(:); 
171   A1 = zeros(length(h),1); 
172   B1 = A1; 
173   [gamu,ix,jx] = unique(gam); 
174   numSeaStates = length(gamu); 
175   gami = zeros(Nh2,1); 
176   for iz=1:numSeaStates 
177     k = find(jx==iz); 
178     gami(:) = gamu(iz); 
179     A1(k) = smooth(h2,interp2(E2,H2,A11,gami,h2,method),1,h(k),1); 
180     B1(k) = smooth(h2,interp2(E2,H2,B11,gami,h2,method),1,h(k),1); 
181   end 
182 else 
183   A1 = smooth(h2,interp2(E2,H2,A11,gam(ones(size(h2))),h2,method),1,h,1); 
184   B1 = smooth(h2,interp2(E2,H2,B11,gam(ones(size(h2))),h2,method),1,h,1); 
185 end 
186 end 
187 [A0,B0,C0] = jhnlwparfun(Hm0,Tp,gam); 
188  
189 switch condon, 
190  case 0, % regular pdf is returned  
191   f = wtweibpdf(h,A0,B0,C0).*wweibpdf(v,A1,B1); 
192  case 1, %pdf conditioned on x1 ie. p(x2|x1)  
193   f = wweibpdf(v,A1,B1); 
194  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
195   f = v.*wweibpdf(v,A1,B1); 
196  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
197   f = v.^2.*wweibpdf(v,A1,B1); 
198  case 5, % p(h)*P(V|h) is returned special case used by jhvnlcdf2 
199   f = wtweibpdf(h,A0,B0,C0).*wweibcdf(v,A1,B1); 
200  case 6, % P(V|h) is returned special case used by jhvnlcdf 
201   f = wweibcdf(v,A1,B1); 
202  case 7,% p(h)*(1-P(V|h)) is returned special case used by jhvnlcdf 
203   f = wtweibpdf(h,A0,B0,C0).*(1-wweibcdf(v,A1,B1)); 
204  otherwise error('unknown option') 
205 end 
206 if multipleSeaStates 
207   f = reshape(f,cSize); 
208 end 
209  
210 if condon~=6 
211   f = f./Hrms./Vrms; 
212 end 
213 f(find(isnan(f)|isinf(f) ))=0; 
214 if any(size(f)~=cSize) 
215   disp('Wrong size') 
216 end 
217  
218 if nargout>3, 
219   fA      = createpdf(2); 
220   fA.f    = A11; 
221   fA.x    = {e2,h2}; 
222   fA.labx = {'Gamma', 'h'}; 
223   fA.note = ['The conditonal Weibull distribution Parameter A(h,gamma)/Hrms' ... 
224     'for wave heigth as a function of h=Hd/Hrms and gamma for' ... 
225     'the Jonswap spectrum']; 
226      
227   ra = range(A11(:)); 
228   st = round(min(A11(:))*100)/100; 
229   en = max(A11(:)); 
230   fA.cl   = st:ra/20:en; 
231 end 
232 if nargout>4, 
233   fB      = createpdf(2); 
234   fB.f    = B11; 
235   fB.x    = {e2,h2}; 
236   fB.labx = {'Gamma', 'h'}; 
237   fB.note = ['The conditonal Weibull distribution Parameter B(h,gamma)/Hrms' ... 
238     'for wave heigth as a function of h=Hd/Hrms and gamma for' ... 
239     'the Jonswap spectrum']; 
240   ra = range(B11(:)); 
241   st = round(min(B11(:))*100)/100; 
242   en = max(B11(:)); 
243   fB.cl   = st:ra/20:en; 
244 end 
245 return 
246  
247  
248  
249  
250  
251  
252

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