Home > wafo > wavemodels > jhsnlcdf.m

jhsnlcdf

PURPOSE ^

Joint (Scf,Hd) CDF for non-linear waves with JONSWAP spectrum.

SYNOPSIS ^

f = jhsnlcdf(Hd,Scf,Hm0,Tp,gam,tail)

DESCRIPTION ^

 JHSNLCDF Joint (Scf,Hd) CDF for non-linear waves with JONSWAP spectrum. 
  
   CALL: f = jhsnlcdf(Hd,Scf,Hm0,Tp,Gamma,tail) 
   
    f   = CDF evaluated at (Scf,Hd) 
    Hd  = zero down crossing wave height [m]  
    Scf = crest front steepness 
    Hm0 = significant wave height [m] 
    Tp  = Spectral peak period    [s] 
  Gamma = Peakedness parameter of the JONSWAP spectrum   
   tail = 1 if upper tail is calculated    
          0 if lower tail is calulated (default) 
    
  JHSNLCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 
  steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order Stokes 
  waves with a JONSWAP spectral density. The empirical parameters of the 
  model is fitted by least squares to simulated (Scf,Hd) data for 13 
  classes of GAMMA between 1 and 7. Between 47000 and 55000 
  zero-downcrossing waves were simulated for each class of GAMMA. 
  JHSNLCDF is restricted to the following range for GAMMA:  
   1 <= GAMMA <= 7  
  
  Example: 
  Hm0 = 6;Tp = 9; gam = 3.5; 
  ec = 0.25; 
  hc = 3; 
  lowerTail = 0; 
  upperTail = ~lowerTail   
  jhsnlcdf(hc,ec,Hm0,Tp,gam)           % Prob(Hd<Hc,Scf<ec) 
  jhsnlcdf(hc,ec,Hm0,Tp,gam,upperTail) % Prob(Hd>Hc,Scf>ec)   
    
   % Conditional probability of steep and high waves given seastates 
   % i.e., Prob(Hd>hc,Scf>ec|Hs,Tp)   
   upperTail = 1; 
   Hs = linspace(2.5,11.5,10); 
   Tp = linspace(4.5,19.5,16); 
   [T,H] = meshgrid(Tp,Hs);  
   p = jhsnlcdf(hc,ec,H,T,gam,upperTail); 
   v = 10.^(-6:-1);   
   contourf(Tp,Hs,log10(p),log10(v)) 
   xlabel('Tp') 
   ylabel('Hs')   
   fcolorbar(log10(v))   
    
  See also  jhscdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function f = jhsnlcdf(Hd,Scf,Hm0,Tp,gam,tail) 
002 %JHSNLCDF Joint (Scf,Hd) CDF for non-linear waves with JONSWAP spectrum. 
003 % 
004 %  CALL: f = jhsnlcdf(Hd,Scf,Hm0,Tp,Gamma,tail) 
005 %  
006 %   f   = CDF evaluated at (Scf,Hd) 
007 %   Hd  = zero down crossing wave height [m]  
008 %   Scf = crest front steepness 
009 %   Hm0 = significant wave height [m] 
010 %   Tp  = Spectral peak period    [s] 
011 % Gamma = Peakedness parameter of the JONSWAP spectrum   
012 %  tail = 1 if upper tail is calculated    
013 %         0 if lower tail is calulated (default) 
014 %   
015 % JHSNLCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 
016 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order Stokes 
017 % waves with a JONSWAP spectral density. The empirical parameters of the 
018 % model is fitted by least squares to simulated (Scf,Hd) data for 13 
019 % classes of GAMMA between 1 and 7. Between 47000 and 55000 
020 % zero-downcrossing waves were simulated for each class of GAMMA. 
021 % JHSNLCDF is restricted to the following range for GAMMA:  
022 %  1 <= GAMMA <= 7  
023 % 
024 % Example: 
025 % Hm0 = 6;Tp = 9; gam = 3.5; 
026 % ec = 0.25; 
027 % hc = 3; 
028 % lowerTail = 0; 
029 % upperTail = ~lowerTail   
030 % jhsnlcdf(hc,ec,Hm0,Tp,gam)           % Prob(Hd<Hc,Scf<ec) 
031 % jhsnlcdf(hc,ec,Hm0,Tp,gam,upperTail) % Prob(Hd>Hc,Scf>ec)   
032 %   
033 %  % Conditional probability of steep and high waves given seastates 
034 %  % i.e., Prob(Hd>hc,Scf>ec|Hs,Tp)   
035 %  upperTail = 1; 
036 %  Hs = linspace(2.5,11.5,10); 
037 %  Tp = linspace(4.5,19.5,16); 
038 %  [T,H] = meshgrid(Tp,Hs);  
039 %  p = jhsnlcdf(hc,ec,H,T,gam,upperTail); 
040 %  v = 10.^(-6:-1);   
041 %  contourf(Tp,Hs,log10(p),log10(v)) 
042 %  xlabel('Tp') 
043 %  ylabel('Hs')   
044 %  fcolorbar(log10(v))   
045 %   
046 % See also  jhscdf 
047  
048 % Reference  
049 % P. A. Brodtkorb (2004),   
050 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
051 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
052 % Trondheim, Norway.   
053    
054 % History 
055 % By pab 20.01.2003 
056  
057  
058 error(nargchk(2,6,nargin))   
059  
060 if (nargin < 6|isempty(tail)),tail = 0; end 
061 if (nargin < 4|isempty(Tp)),Tp = 8; end 
062 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end 
063 if (nargin < 5|isempty(gam)), 
064    gam = getjonswappeakedness(Hm0,Tp); 
065 end 
066  
067 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 
068 if multipleSeaStates 
069   [errorcode, Scf,Hd,Hm0,Tp,gam] = comnsize(Scf,Hd,Hm0,Tp,gam); 
070 else 
071   [errorcode, Scf,Hd] = comnsize(Scf,Hd); 
072 end 
073 if errorcode > 0 
074   error('Requires non-scalar arguments to match in size.'); 
075 end 
076 displayWarning = 0; 
077 if displayWarning 
078   if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0))) 
079     disp('Warning: Hm0,Tp is outside the JONSWAP range') 
080     disp('The validity of the parameters returned are questionable') 
081   end 
082 end 
083 %dev = 2e-5; 
084 c1 =[ 0.16183666835624   1.53691936441548   1.55852759524555]; 
085 c2 =[ 0.15659478203944   1.15736959972513   1]; 
086 if 0, 
087   Tm02 = Tp 
088   Tp = Tm02.*(polyval(c1,gam)./polyval(c2,gam)); 
089 else 
090   Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam)); 
091 end 
092 Hrms = Hm0/sqrt(2); 
093 Vrms = 1.25*Hm0./(Tm02.^2); % Erms 
094  
095 v = Scf./Vrms; 
096  
097 hMax = 6; 
098 eps2 = 1e-6; 
099 normalizedInput = 1; 
100 utprb = gaussq('jhsnlpdf',hMax,2*hMax,eps2/2,[],mean(v(:)),mean(Hm0(:)),mean(Tp(:)),mean(gam(:)),normalizedInput,7); 
101 if eps2<utprb 
102   warning('Check the accuracy of integration!') 
103 end 
104  
105 h = min(Hd./Hrms,hMax); 
106  
107 f = zeros(size(Hd)); 
108 % Only compute 
109 if 0, % haver parametrization 
110   loLimit = 3.6; 
111   upLimit = 5; 
112 else 
113   loLimit = 2.5; 
114   upLimit = 6.5; 
115 end 
116  
117 k0 = find( (loLimit*sqrt(Hm0)<Tp) & (Tp<upLimit*sqrt(Hm0)) ); 
118 if any(k0) 
119   if multipleSeaStates 
120     h = h(k0); 
121     v = v(k0); 
122     Hm0 = Hm0(k0); 
123     Tp = Tp(k0); 
124     gam = gam(k0);     
125   else 
126     k0 = 1:prod(size(Hd)); 
127   end 
128 else 
129   return 
130 end 
131  
132 if 0 
133   % This is a trick to get the html documentation correct. 
134   k = jhsnlpdf(1,1,2,3); 
135 end 
136  
137 hlim  = h; 
138  
139 lowerTail = 0; 
140 if tail==lowerTail, 
141   k       = find(h>2.5); 
142   hlim(k) = 2.5; 
143   f(k0) = gaussq('jhsnlpdf',0,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5)... 
144       + gaussq('jhsnlpdf',hlim,h,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5);  
145 else % upper tail 
146   k       = find(h<2.5); 
147   hlim(k) = 2.5; 
148   f(k0) = gaussq('jhsnlpdf',h,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7)... 
149       + gaussq('jhsnlpdf',hlim,hMax,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7);  
150 end 
151 return 
152  
153

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