Home > wafo > wavemodels > thsnlpdf.m

thsnlpdf

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 THSNLPDF Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. 
  
   CALL: f = thsnlpdf(Hd,Scf,Hm0,Tp) 
   
     f  = pdf 
    Hd  = zero down crossing wave height 
    Scf = crest front steepness 
    Hm0 = significant wave height 
    Tp  = Spectral peak period  
  
  THSNLPDF approximates the joint PDF of (Scf, Hd), i.e., crest 
  steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for for 2nd order 
  non-linear waves 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. 
  THSNLPDF 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 = thsnlpdf(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] = thsnlpdf(Hd,Scf,Hm0,Tp,normalizedInput,condon) 
002 %THSNLPDF Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. 
003 % 
004 %  CALL: f = thsnlpdf(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 % THSNLPDF approximates the joint PDF of (Scf, Hd), i.e., crest 
013 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for for 2nd order 
014 % non-linear waves with a Torsethaugen spectral density. The empirical 
015 % parameters of the model is 
016 % fitted by least squares to simulated (Scf,Hd) data for 600 classes of 
017 % Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves were 
018 % simulated for each class of Hm0 and Tp. 
019 % THSNLPDF is restricted to the following range for Hm0 and Tp:  
020 %  0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11. 
021 %  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 % s = linspace(0,6*1.25*Hm0/Tp^2,101); 
027 % [S,H] = meshgrid(s,h);  
028 % f = thsnlpdf(H,S,Hm0,Tp); 
029 % contourf(s,h,f)   
030 % 
031 % See also  thspdf2, 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    
040 % History 
041 % revised pab 09.08.2003 
042 %  changed input and help header   
043 % validated 20.11.2002 
044 % By pab 20.12.2000 
045  
046 error(nargchk(3,6,nargin)) 
047  
048  
049 if (nargin < 6|isempty(condon)),  condon  = 0; end 
050 if (nargin < 5|isempty(normalizedInput)),  normalizedInput  = 0;end 
051 if (nargin < 4|isempty(Tp)),  Tp  = 8;end 
052 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end 
053  
054 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 
055 if multipleSeaStates 
056   [errorcode, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp); 
057 else 
058   [errorcode, Scf,Hd] = comnsize(Scf,Hd); 
059 end 
060 if errorcode > 0 
061     error('Requires non-scalar arguments to match in size.'); 
062 end 
063 displayWarning = 0; 
064 if displayWarning, 
065   if any(Hm0>11| Hm0>(Tp-2)*12/11)  
066     disp('Warning: Hm0 is outside the valid range') 
067     disp('The validity of the Joint (Hd,Scf) distribution is questionable') 
068   end 
069   if any(Tp>20|Tp<3) 
070     disp('Warning: Tp is outside the valid range') 
071     disp('The validity of the Joint (Hd,Scf) distribution is questionable') 
072   end 
073 end 
074    
075 global THSNLPAR 
076 if isempty(THSNLPAR) 
077   %THSNLPAR = load('thsnlpar.mat'); 
078   THSNLPAR = load('thsnlpar06-Aug-2004.mat')   
079 end 
080  
081 Tpp  = THSNLPAR.Tp; 
082 Hm00 = THSNLPAR.Hm0; 
083 Tm020 = THSNLPAR.Tm02; 
084 h2   = THSNLPAR.h2; 
085 % Interpolation method 
086 method = '*cubic';% Faster interpolation 
087  
088 if normalizedInput, 
089   Hrms = 1; 
090   Vrms = 1; 
091 else 
092   [Tp1,Hs1] = meshgrid(Tpp,Hm00); 
093   
094   if 1, 
095     Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 
096     %Tp/Tm02 
097   else 
098      
099     Tm02 = Tp; 
100     for ix= 1:100 
101       dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method)); 
102       Tp = Tp+dTp; 
103       if all(abs(dTp)<0.01) 
104     dTp 
105     ix 
106     break 
107       end 
108     end 
109     %Tp,    Tp/Tm02 
110   end 
111    
112 %  w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
113   %w    = linspace(0,10,2*1024+1).';  
114 %  St = torsethaugen(w,[Hm0,Tp]); 
115 %  ch   = spec2char(St,{'Tm02','eps2'}); 
116 %  Tm02 = ch(1); 
117 %  eps2 = ch(2); 
118   Hrms = Hm0/sqrt(2); 
119   Vrms = 1.25*Hm0./(Tm02.^2); % Erms 
120 end 
121  
122 h = Hd./Hrms; 
123 v = Scf./Vrms; 
124 cSize = size(h); % common size 
125  
126 % Gamma distribution parameters as a function of Tp Hm0 and h2 
127 A11 = THSNLPAR.A11s; 
128 B11 = THSNLPAR.B11s; 
129  
130 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2); 
131 Nh2 = length(h2); 
132 if multipleSeaStates 
133   h   = h(:); 
134   v   = v(:); 
135   Tp  = Tp(:); 
136   Hm0 = Hm0(:); 
137   A1 = zeros(length(h),1); 
138   B1 = A1; 
139   [TpHm0,ix,jx] = unique([Tp,Hm0],'rows'); 
140   numSeaStates = length(ix); 
141   Tpi = zeros(Nh2,1); 
142   Hm0i = zeros(Nh2,1); 
143   for iz=1:numSeaStates 
144     k = find(jx==iz); 
145     Tpi(:)  = TpHm0(iz,1); 
146     Hm0i(:) = TpHm0(iz,2); 
147     A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),... 
148                1,h(k),1)); 
149     B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),... 
150                1,h(k),1)); 
151   end 
152 else 
153   Tpi  = repmat(Tp,[Nh2,1]); 
154   Hm0i = repmat(Hm0,[Nh2,1]); 
155   A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1)); 
156   B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1)); 
157 end 
158 % Waveheight distribution in time 
159 % Truncated Weibull  distribution parameters as a function of Tp, Hm0  
160 [A0, B0, C0] = thwparfun(Hm0,Tp,'time'); 
161 switch condon, 
162  case 0, % regular pdf is returned  
163   f = wtweibpdf(h,A0,B0,C0).*wgampdf(v,A1,B1); 
164  case 1, %pdf conditioned on x1 ie. p(x2|x1)  
165   f = wgampdf(v,A1,B1); 
166  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
167   f = v.*wgampdf(v,A1,B1); 
168  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
169   f = v.^2.*wgampdf(v,A1,B1); 
170  case 5, % p(h)*P(V|h) is returned special case used by thscdf 
171   f = wtweibpdf(h,A0,B0,C0).*wgamcdf(v,A1,B1); 
172  case 6, % P(V|h) is returned special case used by thscdf 
173   f = wgamcdf(v,A1,B1); 
174  case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf 
175   f = wtweibpdf(h,A0,B0,C0).*(1-wgamcdf(v,A1,B1)); 
176  otherwise error('unknown option') 
177 end 
178 if multipleSeaStates 
179   f = reshape(f,cSize); 
180 end 
181  
182 if condon~=6 
183   f = f./Hrms./Vrms; 
184 end 
185 f(find(isnan(f)|isinf(f) ))=0; 
186 if any(size(f)~=cSize) 
187   disp('Wrong size') 
188 end 
189  
190 return 
191  
192  
193  
194  
195

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