Home > wafo > wavemodels > thvpdf.m

thvpdf

PURPOSE ^

Joint (Vcf,Hd) PDF for linear waves with Torsethaugen spectra.

SYNOPSIS ^

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

DESCRIPTION ^

 THVPDF Joint (Vcf,Hd) PDF for linear waves with Torsethaugen spectra. 
  
   CALL: f = thvpdf(Hd,Vcf,Hm0,Tp) 
   
   f   = pdf evaluated at (Hd,Vcf). 
   Hd  = zero down crossing wave height 
   Vcf = crest front velocity 
   Hm0 = significant wave height [m]. (default Hm0 = 6) 
   Tp  = Spectral peak period    [s]. (default  Tp = 8) 
  
  THVPDF approximates the joint PDF of (Vcf, Hd), i.e., crest 
  front velocity (Ac/Tcf) and wave height, for a Gaussian process with a 
  Torsethaugen spectral density. The empirical parameters of the model is 
  fitted by least squares to simulated (Vcf,Hd) data for 600 classes of 
  Hm0 and Tp. Between 50000 and 150000 zero-downcrossing waves were 
  simulated for each class of Hm0 and Tp. 
  THVPDF is restricted to the following range for Hm0 and Tp:  
  0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11. 
    
  NOTE: The size of f is the common size of the input arguments. 
  
  Example: 
  Hm0 = 6;Tp = 8; 
  h = linspace(0,4*Hm0/sqrt(2))';  
  v = linspace(0,4*2*Hm0/Tp)'; 
  [V,H] = meshgrid(v,h);   
  f = thvpdf(H,V,Hm0,Tp); 
  contourf(v,h,f)   
  
  See also  thspdf, thsspdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [f,Hrms,Vrms,fA,fB] = thvpdf(Hd,Vcf,Hm0,Tp,normalizedInput,condon) 
002 %THVPDF Joint (Vcf,Hd) PDF for linear waves with Torsethaugen spectra. 
003 % 
004 %  CALL: f = thvpdf(Hd,Vcf,Hm0,Tp) 
005 %  
006 %  f   = pdf evaluated at (Hd,Vcf). 
007 %  Hd  = zero down crossing wave height 
008 %  Vcf = crest front velocity 
009 %  Hm0 = significant wave height [m]. (default Hm0 = 6) 
010 %  Tp  = Spectral peak period    [s]. (default  Tp = 8) 
011 % 
012 % THVPDF approximates the joint PDF of (Vcf, Hd), i.e., crest 
013 % front velocity (Ac/Tcf) and wave height, for a Gaussian process with a 
014 % Torsethaugen spectral density. The empirical parameters of the model is 
015 % fitted by least squares to simulated (Vcf,Hd) data for 600 classes of 
016 % Hm0 and Tp. Between 50000 and 150000 zero-downcrossing waves were 
017 % simulated for each class of Hm0 and Tp. 
018 % THVPDF is restricted to the following range for Hm0 and Tp:  
019 % 0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11. 
020 %   
021 % NOTE: The size of f is the common size of the input arguments. 
022 % 
023 % Example: 
024 % Hm0 = 6;Tp = 8; 
025 % h = linspace(0,4*Hm0/sqrt(2))';  
026 % v = linspace(0,4*2*Hm0/Tp)'; 
027 % [V,H] = meshgrid(v,h);   
028 % f = thvpdf(H,V,Hm0,Tp); 
029 % contourf(v,h,f)   
030 % 
031 % See also  thspdf, thsspdf 
032  
033 % Reference   
034 % P. A. Brodtkorb (2004),   
035 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
036 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
037 % Trondheim, Norway.    
038    
039 % History 
040 % revised pab 10.08.2003   
041 % -revised pab 28.11.2002 
042 %   extended the example 
043 % By pab 20.12.2000 
044  
045 error(nargchk(3,6,nargin)) 
046  
047 if (nargin < 6|isempty(condon)),  condon  = 0; end 
048 if (nargin < 5|isempty(normalizedInput)),  normalizedInput  = 0;end 
049 if (nargin < 4|isempty(Tp)),  Tp  = 8;end 
050 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end 
051  
052 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 
053 if multipleSeaStates 
054   [errorcode, Vcf,Hd,Hm0,Tp] = comnsize(Vcf,Hd,Hm0,Tp); 
055 else 
056   [errorcode, Vcf,Hd] = comnsize(Vcf,Hd); 
057 end 
058 if errorcode > 0 
059   error('Requires non-scalar arguments to match in size.'); 
060 end 
061 displayWarning = 0; 
062 if displayWarning, 
063   if any(Hm0>11| Hm0>(Tp-2)*12/11)  
064     disp('Warning: Hm0 is outside the valid range') 
065     disp('The validity of the Joint (Hd,Vcf) distribution is questionable') 
066   end 
067   if any(Tp>20|Tp<3) 
068     disp('Warning: Tp is outside the valid range') 
069     disp('The validity of the Joint (Hd,Vcf) distribution is questionable') 
070   end 
071 end 
072  
073 global THVPAR 
074 if isempty(THVPAR) 
075   THVPAR = load('thvpar.mat'); 
076 end 
077 % Weibull distribution parameters as a function of e2 and h2 
078 A11 = THVPAR.A11s; 
079 B11 = THVPAR.B11s; 
080 e2  = THVPAR.e2; 
081 h2  = THVPAR.h2; 
082 Tpp   = THVPAR.Tp; 
083 Hm00  = THVPAR.Hm0; 
084 Tm020 = THVPAR.Tm02; 
085 eps20 = THVPAR.eps2; 
086  
087 [Tp1,Hs1] = meshgrid(Tpp,Hm00); 
088 method = '*cubic'; %'spline' 
089  
090 eps2 = interp2(Tp1,Hs1,eps20,Tp,Hm0,method); 
091  
092 if normalizedInput, 
093   Hrms = 1; 
094   Vrms = 1; 
095 else   
096   Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 
097 %  w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
098 %  St = torsethaugen(w,[Hm0,Tp]); 
099 %  ch   = spec2char(St,{'Tm02','eps2'}); 
100 %  Tm02 = ch(1); 
101 %  eps2 = ch(2); 
102   Hrms = Hm0/sqrt(2); 
103   Vrms = 2*Hm0./Tm02; % Erms 
104 end 
105    
106 %Fh = thpdf(h(:)/Hrms,Hm0,Tp,eps2,1); 
107 if displayWarning 
108   if eps2<0.4  
109     disp('Warning: eps2 is less than 0.4. The pdf returned is questionable') 
110   elseif eps2>1.3 
111     disp('Warning: eps2 is larger than 1.3. The pdf returned is questionable') 
112   end 
113 end 
114 h = Hd/Hrms; 
115 v = Vcf/Vrms; 
116 cSize = size(h); % common size of input 
117  
118 [E2 H2] = meshgrid(e2,h2); 
119 Nh2 = length(h2); 
120 if multipleSeaStates 
121   h   = h(:); 
122   v   = v(:); 
123   Tp  = Tp(:); 
124   Hm0 = Hm0(:); 
125   eps2 = eps2(:); 
126   A1 = zeros(length(h),1); 
127   B1 = A1; 
128   [eps2u,ix,jx] = unique(eps2); 
129   numSeaStates = length(ix); 
130   eps2i = zeros(Nh2,1); 
131   for iz=1:numSeaStates 
132     k = find(jx==iz); 
133     eps2i(:)  = eps2u(iz); 
134     A1(k) = smooth(h2,interp2(E2,H2,A11,eps2i,h2,method),1,h(k),1); 
135     B1(k) = smooth(h2,interp2(E2,H2,B11,eps2i,h2,method),1,h(k),1); 
136   end 
137 else 
138   eps2i = repmat(eps2,[Nh2,1]); 
139   A1 = smooth(h2,interp2(E2,H2,A11,eps2i,h2,method),1,h,1); 
140   B1 = smooth(h2,interp2(E2,H2,B11,eps2i,h2,method),1,h,1); 
141 end 
142 % Note if eps2<0.4 then B1 is questionable 
143  
144 % Waveheight distribution in time 
145 % Truncated Weibull  distribution parameters as a function of Tp, Hm0  
146 [A0, B0, C0] = thwparfun(Hm0,Tp,'time'); 
147  
148 switch condon, 
149  case 0, % regular pdf is returned  
150   f = wtweibpdf(h,A0,B0,C0).*wweibpdf(v,A1,B1); 
151  case 1, %pdf conditioned on x1 ie. p(x2|x1)  
152   f = wweibpdf(v,A1,B1); 
153  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
154   f = v.*wweibpdf(v,A1,B1); 
155  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
156   f = v.^2.*wweibpdf(v,A1,B1); 
157  case 5, % p(h)*P(V|h) is returned special case used by thvcdf 
158   f = wtweibpdf(h,A0,B0,C0).*wweibcdf(v,A1,B1); 
159  case 6, % P(V|h) is returned special case used by thvcdf 
160   f = wweibcdf(v,A1,B1); 
161  case 7,% p(h)*(1-P(V|h)) is returned special case used by thvcdf 
162   f = wtweibpdf(h,A0,B0,C0).*(1-wweibcdf(v,A1,B1)); 
163  otherwise error('unknown option') 
164 end 
165 if multipleSeaStates 
166   f = reshape(f,cSize); 
167 end 
168  
169 if condon~=6 
170   f = f./Hrms./Vrms; 
171 end 
172 f(find(isnan(f)|isinf(f) ))=0; 
173 if any(size(f)~=cSize) 
174   disp('Wrong size') 
175 end 
176  
177 if nargout>3, 
178   fA      = createpdf(2); 
179   fA.f    = A11; 
180   fA.x    = {e2,h2}; 
181   fA.labx = {'eps2', 'h'}; 
182   fA.title = 'wweibpdf parameter A'; 
183   fA.note = ['The conditonal Weibull distribution Parameter A(h,eps2)/Hrms' ... 
184     'for Vcf given h=Hd/Hrms and eps2 for' ... 
185     'the Torsethaugen spectrum']; 
186      
187   ra = range(A11(:)); 
188   st = round(min(A11(:))*100)/100; 
189   en = max(A11(:)); 
190   fA.cl   = st:ra/20:en; 
191 end 
192 if nargout>4, 
193   fB      = createpdf(2); 
194   fB.f    = B11; 
195   fB.x    = {e2,h2}; 
196   fB.labx = {'eps2', 'h'}; 
197   fB.title = 'wweibpdf parameter B'; 
198   fB.note = ['The conditonal Weibull distribution Parameter B(h,eps2)/Hrms' ... 
199     'for Vcf given h=Hd/Hrms and eps2 for' ... 
200     'the Torsethaugen spectrum']; 
201   ra = range(B11(:)); 
202   st = round(min(B11(:))*100)/100; 
203   en = max(B11(:)); 
204   fB.cl   = st:ra/20:en; 
205 end 
206 return 
207  
208  
209  
210  
211  
212  
213

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