Home > wafo > wavemodels > jhsnlpdf.m

jhsnlpdf

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 JHSNLPDF Joint (Scf,Hd) PDF for nonlinear waves with a JONSWAP spectra. 
  
   CALL: f = jhsnlpdf(Hd,Scf,Hm0,Tp,gam) 
   
     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 
  
  JHSNLPDF approximates the joint PDF of (Scf, Hd), i.e., crest front 
  steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order nonlinear 
  waves with a JONSWAP spectral density. The empirical parameters of the 
  model is fitted by least squares to simulated (Scf,Hd) data for 70 
  classes of Hm0 and Tp. Between 40000 and 220000 zero-downcrossing 
  waves were simulated for each class. 
  JHSNLPDF is restricted to the following range:  
  0.5 < Hm0 < 12 and 3 < Tp < 20 
    
  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));  
  s = linspace(0,6*1.25*Hm0/Tp^2,101); 
  [S,H] = meshgrid(s,h);  
  f = jhsnlpdf(H,S,Hm0,Tp); 
  contourf(s,h,f)   
  
  See also  jhsnlpdf2

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [f,Hrms,Vrms] = jhsnlpdf(Hd,Scf,Hm0,Tp,gam,normalizedInput,condon) 
002 %JHSNLPDF Joint (Scf,Hd) PDF for nonlinear waves with a JONSWAP spectra. 
003 % 
004 %  CALL: f = jhsnlpdf(Hd,Scf,Hm0,Tp,gam) 
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 % JHSNLPDF approximates the joint PDF of (Scf, Hd), i.e., crest front 
014 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order nonlinear 
015 % waves with a JONSWAP spectral density. The empirical parameters of the 
016 % model is fitted by least squares to simulated (Scf,Hd) data for 70 
017 % classes of Hm0 and Tp. Between 40000 and 220000 zero-downcrossing 
018 % waves were simulated for each class. 
019 % JHSNLPDF is restricted to the following range:  
020 % 0.5 < Hm0 < 12 and 3 < Tp < 20 
021 %   
022 % NOTE: The size of f is the common size of the input arguments. 
023 %   
024 % Example: 
025 % Hm0 = 6;Tp = 8; 
026 % h = linspace(0,4*Hm0/sqrt(2));  
027 % s = linspace(0,6*1.25*Hm0/Tp^2,101); 
028 % [S,H] = meshgrid(s,h);  
029 % f = jhsnlpdf(H,S,Hm0,Tp); 
030 % contourf(s,h,f)   
031 % 
032 % See also  jhsnlpdf2 
033  
034 % Reference   
035 % P. A. Brodtkorb (2004),   
036 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
037 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
038 % Trondheim, Norway.  
039    
040 % History 
041 % Validated pab 7 april 2004 
042 % By pab 17.01.2003 
043  
044 error(nargchk(3,7,nargin)) 
045  
046  
047 if (nargin < 7|isempty(condon)),  condon  = 0; end 
048 if (nargin < 6|isempty(normalizedInput)),  normalizedInput  = 0;end 
049  
050 if (nargin < 4|isempty(Tp)),  Tp  = 8;end 
051 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end 
052 if (nargin < 5|isempty(gam)) 
053 %   gam = getjonswappeakedness(Hm0,Tp); 
054 end 
055  
056 gam = getjonswappeakedness(Hm0,Tp); 
057 multipleSeaStates = any(prod(size(Hm0))>1|... 
058             prod(size(Tp)) >1|... 
059             prod(size(gam))>1); 
060 if multipleSeaStates 
061   [errorcode, Scf,Hd,Hm0,Tp,gam] = comnsize(Scf,Hd,Hm0,Tp,gam); 
062 else 
063   [errorcode, Scf,Hd] = comnsize(Scf,Hd); 
064 end 
065 if errorcode > 0 
066     error('Requires non-scalar arguments to match in size.'); 
067 end 
068 displayWarning = 0; 
069 if displayWarning 
070   if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0))) 
071     disp('Warning: Hm0,Tp is outside the JONSWAP range') 
072     disp('The validity of the parameters returned are questionable') 
073   end 
074 end 
075 k = find(gam<1); 
076 if any(k)  
077   if displayWarning, 
078     warning('Peakedness parameter less than 1. Must be larger than 1.') 
079   end 
080   gam(k)=1; 
081 end 
082 k1 = find(7<gam); 
083 if any(k1) 
084   if displayWarning 
085   warning('Peakedness parameter larger than 7. The pdf returned is questionable') 
086   end 
087   gam(k1) = 7; 
088 end 
089  
090 global JHSNLPAR 
091 if isempty(JHSNLPAR) 
092   %JHSNLPAR = load('jhsnlpar.mat'); 
093   JHSNLPAR = load('jhsnlpar21-Jul-2004.mat'); 
094 end 
095 % Gamma distribution parameters as a function of e2 and h2 
096 A11 = JHSNLPAR.A11s; 
097 B11 = JHSNLPAR.B11s; 
098 C1 = 1; 
099  
100  
101  
102 if normalizedInput, 
103   Hrms = 1; 
104   Vrms = 1; 
105 else 
106   %Tm02 = Tp./(1.30301-0.01698*gam+0.12102./gam); 
107   %dev = 2e-5; 
108   c1 =[ 0.16183666835624   1.53691936441548   1.55852759524555]; 
109   c2 =[ 0.15659478203944   1.15736959972513   1]; 
110   Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam)); 
111   Hrms = Hm0/sqrt(2); 
112   Vrms = 1.25*Hm0./(Tm02.^2); % Erms 
113 end 
114  
115 h = Hd./Hrms; 
116 v = Scf./Vrms; 
117 cSize = size(h); % common size 
118  
119 h2  = JHSNLPAR.h2(:);  % Hd/Hrms 
120 Nh2 = length(h2); 
121 method ='*cubic'; 
122 if 1, 
123   Tpp  = JHSNLPAR.Tp; % gamma 
124   Hm00 = JHSNLPAR.Hm0; 
125   [E1, H1, H2] = meshgrid(Tpp,Hm00,h2); 
126   Nh2 = length(h2); 
127   if multipleSeaStates 
128     h   = h(:); 
129     v   = v(:); 
130     Tp  = Tp(:); 
131     Hm0 = Hm0(:); 
132     gam = gam(:); 
133     A1 = zeros(length(h),1); 
134     B1 = A1; 
135     [TpHm0,ix,jx] = unique([Tp,Hm0],'rows'); 
136     numSeaStates = length(ix); 
137     Tpi = zeros(Nh2,1); 
138     Hm0i = zeros(Nh2,1); 
139     for iz=1:numSeaStates 
140       k = find(jx==iz); 
141       Tpi(:)  = TpHm0(iz,1); 
142       Hm0i(:) = TpHm0(iz,2); 
143       A1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),... 
144              1,h(k),1)); 
145       B1(k) = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),... 
146              1,h(k),1)); 
147     end 
148   else 
149     Tpi  = repmat(Tp,[Nh2,1]); 
150     Hm0i = repmat(Hm0,[Nh2,1]); 
151     A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tpi,Hm0i,h2,method),1,h,1)); 
152     B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tpi,Hm0i,h2,method),1,h,1)); 
153   end 
154 else % old call 
155   e2  = JHSNLPAR.gam(:); % gamma 
156   h2  = JHSNLPAR.h2(:);  % Hd/Hrms 
157   [E2 H2] = meshgrid(e2,h2); 
158   if multipleSeaStates 
159     h   = h(:); 
160     v   = v(:); 
161     Tp  = Tp(:); 
162     Hm0 = Hm0(:); 
163     gam = gam(:); 
164     %  eps2 = eps2(:); 
165     A1 = zeros(length(h),1); 
166     B1 = A1; 
167     [gamu,ix,jx] = unique(gam); 
168     numSeaStates = length(gamu); 
169     gami = zeros(Nh2,1); 
170     for iz=1:numSeaStates 
171       k = find(jx==iz); 
172       gami(:) = gamu(iz); 
173       A1(k) = exp(smooth(h2,interp2(E2,H2,log(A11),gami,h2,method),1,h(k),1)); 
174       B1(k) = exp(smooth(h2,interp2(E2,H2,log(B11),gami,h2,method),1,h(k),1)); 
175     end 
176   else 
177     A1 = exp(smooth(h2,interp2(E2,H2,log(A11),... 
178                    gam(ones(size(h2))),h2,method),1,h,1)); 
179     B1 = exp(smooth(h2,interp2(E2,H2,log(B11),... 
180                    gam(ones(size(h2))),h2,method),1,h,1)); 
181   end 
182 end 
183 % Waveheight distribution in time 
184 % Truncated Weibull  distribution parameters as a function of Tp, Hm0, gam  
185 [A0, B0, C0] = jhnlwparfun(Hm0,Tp,gam,'time'); 
186 switch condon, 
187  case 0, % regular pdf is returned  
188   f = wtweibpdf(h,A0,B0,C0).*wggampdf(v,A1,B1,C1); 
189  case 1, %pdf conditioned on x1 ie. p(x2|x1)  
190   f = wggampdf(v,A1,B1,C1); 
191  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
192   f = v.*wggampdf(v,A1,B1,C1); 
193  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
194   f = v.^2.*wggampdf(v,A1,B1,C1); 
195  case 5, % p(h)*P(V|h) is returned special case used by thscdf 
196   f = wtweibpdf(h,A0,B0,C0).*wggamcdf(v,A1,B1,C1); 
197  case 6, % P(V|h) is returned special case used by thscdf 
198   f = wggamcdf(v,A1,B1,C1); 
199  case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf 
200   f = wtweibpdf(h,A0,B0,C0).*(1-wggamcdf(v,A1,B1,C1)); 
201  otherwise error('unknown option') 
202 end 
203 if multipleSeaStates 
204   f = reshape(f,cSize); 
205 end 
206  
207 if condon~=6 
208   f = f./Hrms./Vrms; 
209 end 
210 f(find(isnan(f)|isinf(f) ))=0; 
211 if any(size(f)~=cSize) 
212   disp('Wrong size') 
213 end 
214  
215 return 
216  
217  
218  
219  
220

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