Home > wafo > wavemodels > thsspdf.m

thsspdf

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 THSSPDF Joint (Scf,Hd) PDF for linear waves in space with Torsethaugen spectra. 
  
   CALL: f = thsspdf(Hd,Scf,Hm0,Tp) 
   
   f   = pdf  
   Hd  = zero down crossing wave height 
   Scf = crest front steepness 
   Hm0 = significant wave height [m] 
   Tp  = Spectral peak period    [s] 
  
  THSSPDF 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 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 100000 and 1000000 zero-downcrossing 
  waves were simulated for each class of Hm0 and Tp. 
  THSSPDF 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. 
  
  Example: 
  Hm0 = 6;Tp = 8; 
  h = linspace(0,4*Hm0/sqrt(2));  
  v = linspace(0,6*1.25*Hm0/Tp^2); 
  [V,H] = meshgrid(v,h);   
  f = thsspdf(H,V,Hm0,Tp); 
  contourf(v,h,f)   
    
  w = linspace(0,10,2*1024+1).';  
  S = torsethaugen(w,[Hm0 Tp]); 
  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); 
  fk.title = f.title; fk.labx = f.labx;  
  plot(vi,hi,'.'), hold on 
  pdfplot(f) 
  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] = thsspdf(Hd,Scf,Hm0,Tp,normalizedInput,condon) 
002 %THSSPDF Joint (Scf,Hd) PDF for linear waves in space with Torsethaugen spectra. 
003 % 
004 %  CALL: f = thsspdf(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 [m] 
010 %  Tp  = Spectral peak period    [s] 
011 % 
012 % THSSPDF approximates the joint distribution of (Scf, Hd), i.e., crest 
013 % front steepness (Ac/Lcf) and wave height in space, for a Gaussian 
014 % process with a Torsethaugen spectral density. The empirical parameters 
015 % of the model is fitted by least squares to simulated (Scf,Hd) data for 
016 % 600 classes of Hm0 and Tp. Between 100000 and 1000000 zero-downcrossing 
017 % waves were simulated for each class of Hm0 and Tp. 
018 % THSSPDF 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 % Example: 
022 % Hm0 = 6;Tp = 8; 
023 % h = linspace(0,4*Hm0/sqrt(2));  
024 % v = linspace(0,6*1.25*Hm0/Tp^2); 
025 % [V,H] = meshgrid(v,h);   
026 % f = thsspdf(H,V,Hm0,Tp); 
027 % contourf(v,h,f)   
028 %   
029 % w = linspace(0,10,2*1024+1).';  
030 % S = torsethaugen(w,[Hm0 Tp]); 
031 % Sk = spec2spec(specinterp(S,.55),'k1d'); 
032 % dk = 1; 
033 % x = spec2sdat(Sk,80000,dk); rate = 8; 
034 % [vi,hi] = dat2steep(x,rate,1); 
035 % fk = kdebin([vi,hi],'epan',[],[],.5,128); 
036 % fk.title = f.title; fk.labx = f.labx;  
037 % plot(vi,hi,'.'), hold on 
038 % pdfplot(f) 
039 % pdfplot(fk,'r'), hold off 
040 % 
041 % See also  thvpdf 
042  
043    
044 % Reference   
045 % P. A. Brodtkorb (2004),   
046 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
047 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
048 % Trondheim, Norway.   
049    
050 % History 
051 % revised pab 09.08.2003   
052 % By pab 20.12.2000 
053  
054 error(nargchk(3,6,nargin)) 
055 if (nargin < 6|isempty(condon)),  condon  = 0;end 
056 if (nargin < 5|isempty(normalizedInput)),  normalizedInput  = 0;end 
057 if (nargin < 4|isempty(Tp)),  Tp  = 8;end 
058 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end 
059  
060  
061 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 
062 if multipleSeaStates 
063   [errorcode, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp); 
064 else 
065   [errorcode, Scf,Hd] = comnsize(Scf,Hd); 
066 end 
067 if errorcode > 0, 
068   error('Requires non-scalar arguments to match in size.'); 
069 end 
070 displayWarning = 0; 
071 if displayWarning, 
072   if any(Hm0>11| Hm0>(Tp-2)*12/11)  
073     disp('Warning: Hm0 is outside the valid range') 
074     disp('The validity of the Joint (Hd,Scf) distribution in space is questionable') 
075   end 
076   if any(Tp>20|Tp<3 ) 
077     disp('Warning: Tp is outside the valid range') 
078     disp('The validity of the Joint (Hd,Scf) distribution in space is questionable') 
079   end 
080 end 
081 useWeibull = 1; 
082 if useWeibull 
083   global THSSPARW 
084   if isempty(THSSPARW) 
085     THSSPARW = load('thsspar27-Jul-2004.mat'); 
086   end 
087  
088   % Gamma distribution parameters as a function of Tp Hm0 and h2 
089   A11 = THSSPARW.A11s; 
090   B11 = THSSPARW.B11s; 
091    
092   % Truncated weibull distribution parameters as a function of Tp and Hm0  
093   A00 = THSSPARW.A00s; 
094   B00 = THSSPARW.B00s; 
095   C00 = THSSPARW.C00s; 
096  
097   Tpp  = THSSPARW.Tp; 
098   Hm00 = THSSPARW.Hm0; 
099   h2   = THSSPARW.h2(:); 
100   Tm020 = THSSPARW.Tm02; 
101 else 
102   global THSSPARG 
103   if isempty(THSSPARG) 
104     THSSPARG = load('thsspar.mat'); 
105   end 
106  
107   % Gamma distribution parameters as a function of Tp Hm0 and h2 
108   A11 = THSSPARG.A11s; 
109   B11 = THSSPARG.B11s; 
110    
111   % Generalized gamma  distribution parameters as a function of Tp and Hm0  
112   A00 = THSSPARG.A00s; 
113   B00 = THSSPARG.B00s; 
114   C00 = THSSPARG.C00s; 
115  
116   Tpp  = THSSPARG.Tp; 
117   Hm00 = THSSPARG.Hm0; 
118   h2   = THSSPARG.h2; 
119   Tm020 = THSSPARG.Tm02; 
120 end 
121 if normalizedInput, 
122   Hrms = 1; 
123   Vrms = 1; 
124 else 
125   method = 'cubic'; 
126   [Tp1,Hs1] = meshgrid(Tpp,Hm00); 
127   Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 
128   %w    = linspace(0,10,2*1024+1).';  
129   %S = spec2spec(specinterp(torsethaugen(w,[Hm0,Tp]),.55),'k1d'); 
130   %ch   = spec2char(S,{'Tm02'}) 
131   %Tm02 = ch(1); 
132   Hrms = Hm0/sqrt(2); 
133   Vrms = 2*Hm0./Tm02; % Srms 
134 end 
135  
136 %Hrms = Hm0/sqrt(2); 
137 %%w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
138 %w    = linspace(0,10,2*1024+1).';  
139 %S = spec2spec(specinterp(torsethaugen(w,[Hm0,Tp]),.55),'k1d'); 
140  
141 %ch   = spec2char(S,{'Tm02'}); 
142 %Tm02 = ch(1); 
143 %Vrms = 2*Hm0/Tm02; % Actually  Srms 
144  
145  
146 h = Hd./Hrms; 
147 v = Scf./Vrms; 
148 cSize = size(h); % common size 
149  
150 method = '*cubic';% Faster interpolation 
151 %method ='spline'; 
152  
153 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2); 
154 Nh2 = length(h2); 
155  
156 if multipleSeaStates 
157   h   = h(:); 
158   v   = v(:); 
159   Tp  = Tp(:); 
160   Hm0 = Hm0(:); 
161   A1 = zeros(length(h),1); 
162   B1 = A1; 
163   [TpHm0,ix,jx] = unique([Tp,Hm0],'rows'); 
164   numSeaStates = length(ix); 
165   Tpi = zeros(Nh2,1); 
166   Hm0i = zeros(Nh2,1); 
167   for iz=1:numSeaStates 
168     k = find(jx==iz); 
169     Tpi(:)  = TpHm0(iz,1); 
170     Hm0i(:) = TpHm0(iz,2); 
171     A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),... 
172                1,h(k),1)); 
173     B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),... 
174                1,h(k),1)); 
175   end 
176 else 
177   Tpi  = repmat(Tp,[Nh2,1]); 
178   Hm0i = repmat(Hm0,[Nh2,1]); 
179   A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1)); 
180   B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1)); 
181 end 
182  
183 [E1, H1] = meshgrid(Tpp,Hm00); 
184 A0 = interp2(E1,H1,A00,Tp,Hm0,method); 
185 B0 = interp2(E1,H1,B00,Tp,Hm0,method); 
186 C0 = interp2(E1,H1,C00,Tp,Hm0,method); 
187  
188 if useWeibull 
189   switch condon, 
190    case 0, % regular pdf is returned  
191     f = wtweibpdf(h,A0,B0,C0).*wgampdf(v,A1,B1); 
192    case 1, %pdf conditioned on x1 ie. p(x2|x1)  
193     f = wgampdf(v,A1,B1); 
194    case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
195     f = v.*wgampdf(v,A1,B1); 
196    case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
197     f = v.^2.*wgampdf(v,A1,B1); 
198    case 5, % p(h)*P(V|h) is returned special case used by thscdf 
199     f = wtweibpdf(h,A0,B0,C0).*wgamcdf(v,A1,B1); 
200    case 6, % P(V|h) is returned special case used by thscdf 
201     f = wgamcdf(v,A1,B1); 
202    case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf 
203     f = wtweibpdf(h,A0,B0,C0).*(1-wgamcdf(v,A1,B1)); 
204    otherwise error('unknown option') 
205   end 
206 else 
207    
208   switch condon, 
209    case 0, % regular pdf is returned  
210     f = wggampdf(h,A0,B0,C0).*wgampdf(v,A1,B1); 
211    case 1, %pdf conditioned on x1 ie. p(x2|x1)  
212     f = wgampdf(v,A1,B1); 
213    case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
214     f = v.*wgampdf(v,A1,B1); 
215    case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
216     f = v.^2.*wgampdf(v,A1,B1); 
217    case 5, % p(h)*P(V|h) is returned special case used by thscdf 
218     f = wggampdf(h,A0,B0,C0).*wgamcdf(v,A1,B1); 
219    case 6, % P(V|h) is returned special case used by thscdf 
220     f = wgamcdf(v,A1,B1); 
221    case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf 
222     f = wggampdf(h,A0,B0,C0).*(1-wgamcdf(v,A1,B1)); 
223    otherwise error('unknown option') 
224   end 
225 end 
226 if multipleSeaStates 
227   f = reshape(f,cSize); 
228 end 
229  
230 if condon~=6 
231   f = f./Hrms./Vrms; 
232 end 
233 f(find(isnan(f)|isinf(f) ))=0; 
234 if any(size(f)~=cSize) 
235   disp('Wrong size') 
236 end 
237  
238 if nargout>3, 
239   fA      = createpdf(2); 
240   fA.x    = {Tpp,Hm00}; 
241   fA.labx = {'Tp', 'Hm0'}; 
242   fA(3)   = fA(1); 
243   fA(2)   = fA(1); 
244    
245   fA(1).f    = A00; 
246   fA(2).f    = B00; 
247   fA(3).f    = C00; 
248    
249   fA(1).title = 'wggampdf parameter A'; 
250   fA(2).title = 'wggampdf parameter B'; 
251   fA(3).title = 'wggampdf parameter C'; 
252    
253   txt1 = ['The Wggampdf distribution Parameter ']; 
254   txt2=[' for wave heigth in space as a function of Tp and Hm0 for' ... 
255     ' the Torsethaugen spectrum']; 
256   fA(1).note =[txt1 'A' txt2]; 
257   fA(2).note =[txt1 'B' txt2]; 
258   fA(3).note =[txt1 'C' txt2]; 
259    
260   tmp= [A00(:) B00(:) C00(:)]; 
261   ra = range(tmp); 
262   st = round(min(tmp)*100)/100; 
263   en = max(tmp); 
264   for ix=1:3 
265     fA(ix).cl   = st(ix):ra(ix)/20:en(ix); 
266   end 
267 end 
268 if nargout>4, 
269   fB      = createpdf(3); 
270   fB.x    = {Tpp,Hm00,h2}; 
271   fB.labx = {'Tp','Hm0', 'h'}; 
272   fB(2) = fB(1); 
273    
274   fB(1).f = A11; 
275   fB(2).f = B11; 
276    
277   txt11 = 'The conditonal Wgampdf distribution Parameter '; 
278   txt22 = [' for Scf given h=Hd/Hrms in space as function of Tp' ... 
279     ' and Hm0 for the Torsethaugen spectrum']; 
280   fB(1).title = 'wgampdf parameter A'; 
281   fB(2).title = 'wgampdf parameter B'; 
282    
283   fB(1).note = [txt11,'A',txt22]; 
284   fB(2).note = [txt11,'B',txt22]; 
285   %fB(2).note = ['The conditonal Wggampdf distribution Parameter B(h)/Hrms' ... 
286 %    'for wave heigth in space as a function of h=Hd/Hrms and eps2 for' ... 
287 %    'the Torsethaugen spectrum']; 
288  tmp= [A11(:) B11(:)]; 
289   ra = range(tmp); 
290   st = round(min(tmp)*100)/100; 
291   en = max(tmp); 
292   for ix=1:2 
293     fB(ix).cl   = st(ix):ra(ix)/20:en(ix); 
294   end 
295 end 
296 return 
297  
298  
299  
300  
301  
302  
303

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