Home > wafo > wavemodels > jhvcdf.m

jhvcdf

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 JHVCDF Joint (Vcf,Hd) CDF for linear waves with JONSWAP spectrum. 
  
   CALL: f = jhvcdf(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) 
    
  JHVCDF 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. About 100000 zero-downcrossing waves were 
  simulated for each class of GAMMA. 
  JHVCDF 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   
  jhvcdf(hc,vc,Hm0,Tp,gam)           % Prob(Hd<Hc,Vcf<Vc) 
  jhvcdf(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);  
   p = jhvcdf(hc,vc,H,T,gam,upperTail); 
   v = 10.^(-6:-1);   
   contourf(Tp,Hs,log10(p),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 = jhvcdf(Hd,Vcf,Hm0,Tp,gam,tail) 
002 %JHVCDF Joint (Vcf,Hd) CDF for linear waves with JONSWAP spectrum. 
003 % 
004 %  CALL: f = jhvcdf(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 % JHVCDF 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. About 100000 zero-downcrossing waves were 
020 % simulated for each class of GAMMA. 
021 % JHVCDF 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 % jhvcdf(hc,vc,Hm0,Tp,gam)           % Prob(Hd<Hc,Vcf<Vc) 
031 % jhvcdf(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 %  p = jhvcdf(hc,vc,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  thvpdf 
047  
048    
049 % Reference  
050 % P. A. Brodtkorb (2004),   
051 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
052 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
053 % Trondheim, Norway.   
054    
055 % History 
056 % revised pab 09.08.2003 
057 % changed input   
058 % validated 20.11.2002 
059 % By pab 20.12.2000 
060  
061  
062 error(nargchk(2,6,nargin))   
063  
064 if (nargin < 6|isempty(tail)),tail = 0; end 
065 if (nargin < 4|isempty(Tp)),Tp = 8; end 
066 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end 
067 if (nargin < 5|isempty(gam)), 
068    gam = getjonswappeakedness(Hm0,Tp); 
069 end 
070  
071 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 
072 if multipleSeaStates 
073   [errorcode, Vcf,Hd,Hm0,Tp,gam] = comnsize(Vcf,Hd,Hm0,Tp,gam); 
074 else 
075   [errorcode, Vcf,Hd] = comnsize(Vcf,Hd); 
076 end 
077 if errorcode > 0 
078   error('Requires non-scalar arguments to match in size.'); 
079 end 
080 displayWarning = 0; 
081 if displayWarning 
082   if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0))) 
083     disp('Warning: Hm0,Tp is outside the JONSWAP range') 
084     disp('The validity of the parameters returned are questionable') 
085   end 
086 end 
087 %dev = 2e-5; 
088 c1 =[ 0.16183666835624   1.53691936441548   1.55852759524555]; 
089 c2 =[ 0.15659478203944   1.15736959972513   1]; 
090 Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam)); 
091  
092 Hrms = Hm0/sqrt(2); 
093 Vrms = 2*Hm0./Tm02; % Erms 
094  
095 v = Vcf./Vrms; 
096  
097 hMax = 6; 
098 eps2 = 1e-6; 
099 normalizedInput = 1; 
100  
101  
102 h = min(Hd./Hrms,hMax); 
103 f = zeros(size(Hd)); 
104 % Only compute 
105 k0 = find((Tp<5*sqrt(Hm0)) & (3.6*sqrt(Hm0)<Tp)); 
106 if any(k0) 
107   if multipleSeaStates 
108     h = h(k0); 
109     v = v(k0); 
110     Hm0 = Hm0(k0); 
111     Tp = Tp(k0); 
112     gam = gam(k0);     
113   else 
114     k0 = 1:prod(size(Hd)); 
115   end 
116   utprb = gaussq('jhvpdf',hMax,2*hMax,eps2/2,[],mean(v(:)),mean(Hm0(:)),mean(Tp(:)),mean(gam(:)),normalizedInput,7); 
117   if eps2<utprb 
118     warning('Check the accuracy of integration!') 
119   end 
120    
121   if 0 
122     % This is a trick to get the html documentation correct. 
123     k = jhvpdf(1,1,2,3); 
124   end 
125    
126   hlim  = h; 
127  
128  
129   lowerTail = 0; 
130   if tail==lowerTail, 
131     k       = find(h>2*v); 
132     hlim(k) = 2*v(k); 
133     f(k0) = gaussq('jhvpdf',0,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5)... 
134         + gaussq('jhvpdf',hlim,h,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,5);  
135   else % upper tail 
136     k       = find(h<2*v); 
137     hlim(k) = 2*v(k); 
138     f(k0) = gaussq('jhvpdf',h,hlim,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7)... 
139         + gaussq('jhvpdf',hlim,hMax,eps2/2,[],v,Hm0,Tp,gam,normalizedInput,7);  
140   end 
141 end 
142 return 
143  
144

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