Home > wafo > wavemodels > jhvnlcdf.m

jhvnlcdf

PURPOSE ^

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

SYNOPSIS ^

f = jhvnlcdf(Hd,Vcf,Hm0,Tp,gam,tail)

DESCRIPTION ^

 JHVNLCDF Joint (Vcf,Hd) CDF for non-linear waves with JONSWAP spectrum. 
  
   CALL: f = jhvnlcdf(Hd,Vcf,Hm0,Tp,Gamma,tail) 
   
    f   = CDF 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   
   tail = 1 if upper tail is calculated    
          0 if lower tail is calulated (default) 
    
  JHVNLCDF approximates the joint CDF of (Vcf, Hd), i.e., crest front 
  velocity (Ac/Tcf) and wave height, for a Gaussian process with a 
  JONSWAP spectral density. The empirical parameters of the model is 
  fitted by least squares to simulated (Vcf,Hd) data for 13 classes of 
  GAMMA between 1 and 7. Between 50000 and 150000 zero-downcrossing waves were 
  simulated for each class of GAMMA. 
  JHVNLCDF is restricted to the following range for GAMMA:  
   1 <= GAMMA <= 7  
  
  Example: 
  Hm0 = 6;Tp = 8; gam = 3.5; 
  vc = 3; 
  hc = 3; 
  lowerTail = 0; 
  upperTail = ~lowerTail   
  jhvnlcdf(hc,vc,Hm0,Tp,gam)           % Prob(Hd<Hc,Vcf<Vc) 
  jhvnlcdf(hc,vc,Hm0,Tp,gam,upperTail) % Prob(Hd>Hc,Vcf>Vc)   
    
   % Conditional probability of steep and high waves given seastates 
   % i.e., Prob(Hd>hc,Vcf>vc|Hs,Tp)   
   upperTail = 1; 
   Hs = linspace(2.5,11.5,10); 
   Tp = linspace(4.5,19.5,16); 
   [T,H] = meshgrid(Tp,Hs);  
   pl = jhvcdf(hc,vc,H,T,gam,upperTail);  
   p = jhvnlcdf(hc,vc,H,T,gam,upperTail); 
   v = 10.^(-6:-1);   
   contourf(Tp,Hs,log10(p),log10(v)) 
   xlabel('Tp'), ylabel('Hs'),  fcolorbar(log10(v))   
   figure(2)   
   contourf(Tp,Hs,log10(pl),log10(v)) 
   xlabel('Tp'), ylabel('Hs'),  fcolorbar(log10(v))   
      
    
  See also  thvpdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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