Home > wafo > wavemodels > thspdf.m

thspdf

PURPOSE ^

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

SYNOPSIS ^

[f,Hrms,Vrms] = thspdf(Hd,Scf,Hm0,Tp,normalizedInput,condon)

DESCRIPTION ^

 THSPDF Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. 
  
   CALL: f = thspdf(Hd,Scf,Hm0,Tp) 
   
     f  = pdf 
    Hd  = zero down crossing wave height 
    Scf = crest front steepness 
    Hm0 = significant wave height 
    Tp  = Spectral peak period  
  
  THSPDF approximates the joint PDF of (Scf, Hd), i.e., crest 
  steepness (2*pi*Ac/(g*Td*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 (Scf,Hd) data for 600 classes of 
  Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves were 
  simulated for each class of Hm0 and Tp. 
  THSPDF 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. 
   The size of f is the common size of the input arguments. 
    
  Example: 
  Hm0 = 6;Tp = 8; 
  h = linspace(0,4*Hm0/sqrt(2));  
  s = linspace(0,6*1.25*Hm0/Tp^2,101); 
  [S,H] = meshgrid(s,h);  
  f = thspdf(H,S,Hm0,Tp); 
  contourf(s,h,f)   
  
  See also  thspdf2, thsspdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [f,Hrms,Vrms] = thspdf(Hd,Scf,Hm0,Tp,normalizedInput,condon) 
002 %THSPDF Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. 
003 % 
004 %  CALL: f = thspdf(Hd,Scf,Hm0,Tp) 
005 %  
006 %    f  = pdf 
007 %   Hd  = zero down crossing wave height 
008 %   Scf = crest front steepness 
009 %   Hm0 = significant wave height 
010 %   Tp  = Spectral peak period  
011 % 
012 % THSPDF approximates the joint PDF of (Scf, Hd), i.e., crest 
013 % steepness (2*pi*Ac/(g*Td*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 (Scf,Hd) data for 600 classes of 
016 % Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves were 
017 % simulated for each class of Hm0 and Tp. 
018 % THSPDF 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 %  The size of f is the common size of the input arguments. 
021 %   
022 % Example: 
023 % Hm0 = 6;Tp = 8; 
024 % h = linspace(0,4*Hm0/sqrt(2));  
025 % s = linspace(0,6*1.25*Hm0/Tp^2,101); 
026 % [S,H] = meshgrid(s,h);  
027 % f = thspdf(H,S,Hm0,Tp); 
028 % contourf(s,h,f)   
029 % 
030 % See also  thspdf2, thsspdf 
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 % revised pab 09.08.2003 
040 %  changed input and help header   
041 % validated 20.11.2002 
042 % By pab 20.12.2000 
043  
044 error(nargchk(3,6,nargin)) 
045  
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, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp); 
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 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,Scf) 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,Scf) distribution is questionable') 
070   end 
071 end 
072    
073 global THSPAR 
074 if isempty(THSPAR) 
075   %THSPAR = load('thspar.mat'); 
076   THSPAR = load('thspar19-Jul-2004.mat'); 
077 end 
078  
079 Tpp  = THSPAR.Tp; 
080 Hm00 = THSPAR.Hm0; 
081 Tm020 = THSPAR.Tm02; 
082 h2   = THSPAR.h2(:); 
083  
084 % Interpolation method 
085 method = '*cubic';% Faster interpolation 
086  
087 if normalizedInput, 
088   Hrms = 1; 
089   Vrms = 1; 
090 else 
091   [Tp1,Hs1] = meshgrid(Tpp,Hm00); 
092   if 1, 
093     Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 
094     %Tp/Tm02 
095   else 
096      
097     Tm02 = Tp; 
098     for ix= 1:100 
099       dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method)); 
100       Tp = Tp+dTp; 
101       if all(abs(dTp)<0.01) 
102     dTp 
103     ix 
104     break 
105       end 
106     end 
107     %Tp,    Tp/Tm02 
108   end 
109    
110 %  w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
111   %w    = linspace(0,10,2*1024+1).';  
112 %  St = torsethaugen(w,[Hm0,Tp]); 
113 %  ch   = spec2char(St,{'Tm02','eps2'}); 
114 %  Tm02 = ch(1); 
115 %  eps2 = ch(2); 
116   Hrms = Hm0/sqrt(2); 
117   Vrms = 1.25*Hm0./(Tm02.^2); % Erms 
118 end 
119  
120 h = Hd./Hrms; 
121 v = Scf./Vrms; 
122 cSize = size(h); % common size 
123  
124 % Gamma distribution parameters as a function of Tp Hm0 and h2 
125 A11 = THSPAR.A11s; 
126 B11 = THSPAR.B11s; 
127  
128 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2); 
129 Nh2 = length(h2); 
130  
131 if multipleSeaStates 
132   h   = h(:); 
133   v   = v(:); 
134   Tp  = Tp(:); 
135   Hm0 = Hm0(:); 
136   A1 = zeros(length(h),1); 
137   B1 = A1; 
138   [TpHm0,ix,jx] = unique([Tp,Hm0],'rows'); 
139   numSeaStates = length(ix); 
140   Tpi = zeros(Nh2,1); 
141   Hm0i = zeros(Nh2,1); 
142   for iz=1:numSeaStates 
143     k = find(jx==iz); 
144     Tpi(:)  = TpHm0(iz,1); 
145     Hm0i(:) = TpHm0(iz,2); 
146     A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),... 
147                1,h(k),1)); 
148     B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),... 
149                1,h(k),1)); 
150   end 
151 else 
152   Tpi  = repmat(Tp,[Nh2,1]); 
153   Hm0i = repmat(Hm0,[Nh2,1]); 
154   A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1)); 
155   B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1)); 
156 end 
157 % Waveheight distribution in time 
158 % Truncated Weibull  distribution parameters as a function of Tp, Hm0  
159 [A0, B0, C0] = thwparfun(Hm0,Tp,'time'); 
160 switch condon, 
161  case 0, % regular pdf is returned  
162   f = wtweibpdf(h,A0,B0,C0).*wgampdf(v,A1,B1); 
163  case 1, %pdf conditioned on x1 ie. p(x2|x1)  
164   f = wgampdf(v,A1,B1); 
165  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
166   f = v.*wgampdf(v,A1,B1); 
167  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
168   f = v.^2.*wgampdf(v,A1,B1); 
169  case 5, % p(h)*P(V|h) is returned special case used by thscdf 
170   f = wtweibpdf(h,A0,B0,C0).*wgamcdf(v,A1,B1); 
171  case 6, % P(V|h) is returned special case used by thscdf 
172   f = wgamcdf(v,A1,B1); 
173  case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf 
174   f = wtweibpdf(h,A0,B0,C0).*(1-wgamcdf(v,A1,B1)); 
175  otherwise error('unknown option') 
176 end 
177 if multipleSeaStates 
178   f = reshape(f,cSize); 
179 end 
180  
181 if condon~=6 
182   f = f./Hrms./Vrms; 
183 end 
184 f(find(isnan(f)|isinf(f) ))=0; 
185 if any(size(f)~=cSize) 
186   disp('Wrong size') 
187 end 
188  
189 if nargout>3, 
190   fA      = createpdf(2); 
191   fA.x    = {Tpp,Hm00}; 
192   fA.labx = {'Tp', 'Hm0'}; 
193   fA(3)   = fA(1); 
194   fA(2)   = fA(1); 
195    
196   % Truncated Weibull  distribution parameters as a function of Tp, Hm0  
197   A00 = THSPAR.A00s; 
198   B00 = THSPAR.B00s; 
199   C00 = THSPAR.C00s; 
200    
201   fA(1).f    = A00; 
202   fA(2).f    = B00; 
203   fA(3).f    = C00; 
204    
205   fA(1).title = 'wtweibpdf parameter A'; 
206   fA(2).title = 'wtweibpdf parameter B'; 
207   fA(3).title = 'wtweibpdf parameter C'; 
208    
209   txt1 = ['The Wtweibpdf  distribution Parameter ']; 
210   txt2=[' for wave heigth in time as a function of Tp and Hm0 for' ... 
211     'the Torsethaugen spectrum']; 
212   fA(1).note =[txt1 'A' txt2]; 
213   fA(2).note =[txt1 'B' txt2]; 
214   fA(3).note =[txt1 'C' txt2]; 
215    
216   tmp = [A00(:) B00(:) C00(:)]; 
217   ra  = range(tmp); 
218   st  = round(min(tmp)*100)/100; 
219   en  = max(tmp); 
220   for ix = 1:3, 
221     fA(ix).cl   = st(ix):ra(ix)/20:en(ix); 
222   end 
223 end 
224 if nargout>4, 
225   fB      = createpdf(3); 
226   fB.x    = {Tpp,Hm00,h2}; 
227   fB.labx = {'Tp','Hm0', 'h'}; 
228   fB(2)   = fB(1); 
229    
230   fB(1).f = A11; 
231   fB(2).f = B11; 
232    
233   txt11 = 'The conditonal Wgampdf distribution Parameter '; 
234   txt22 = [' for Scf given h=Hd/Hrms in time as function of Tp' ... 
235     ' and Hm0 for the Torsethaugen spectrum']; 
236   fB(1).title = 'wgampdf parameter A'; 
237   fB(2).title = 'wgampdf parameter B'; 
238   fB(1).note = [txt11,'A',txt22]; 
239   fB(2).note = [txt11,'B',txt22]; 
240    
241   %fB(2).note = ['The conditonal Wggampdf distribution Parameter B(h)/Hrms', ...%    ' for crest front steepness as a function of Tp,Hm0 and',... 
242   %    ' h=Hd/Hrms for the Torsethaugen spectrum']; 
243   tmp= [A11(:) B11(:)]; 
244   ra = range(tmp); 
245   st = round(min(tmp)*100)/100; 
246   en = max(tmp); 
247   for ix=1:2 
248     fB(ix).cl   = st(ix):ra(ix)/20:en(ix); 
249   end 
250 end 
251 return 
252  
253  
254  
255  
256

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