Home > wafo > wavemodels > thvcdf.m

thvcdf

PURPOSE ^

Joint (Vcf,Hd) CDF for linear waves with Torsethaugen spectra.

SYNOPSIS ^

f = thvcdf(Hd,Vcf,Hm0,Tp,tail)

DESCRIPTION ^

 THVCDF Joint (Vcf,Hd) CDF for linear waves with Torsethaugen spectra. 
  
   CALL: f = thvcdf(Hd,Vcf,Hm0,Tp,tail) 
   
    f   = CDF evaluated at (Vcf,Hd) 
    Hd  = zero down crossing wave height  
    Vcf = crest front velocity 
    Hm0 = significant wave height [m] 
    Tp  = Spectral peak period    [s] 
   tail = 1 if upper tail is calculated    
          0 if lower tail is calulated (default) 
    
  THVCDF approximates the joint CDF of (Vcf, Hd), i.e., crest front 
  velocity (Ac/Tcf) and wave height, for a Gaussian process with a 
  Torsethaugen spectral density. The empirical parameters of the model is 
  fitted by least squares to simulated (Vcf,Hd) data for 600 classes of 
  Hm0 and Tp. Between 50000 and 150000 zero-downcrossing waves were 
  simulated for each class of Hm0 and Tp. 
  THVCDF 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. 
  
  Example: 
  Hm0 = 6;Tp = 8; 
  vc = 3; 
  hc = 3; 
  lowerTail = 0; 
  upperTail = ~lowerTail   
  thvcdf(hc,vc,Hm0,Tp)           % Prob(Hd<Hc,Vcf<Vc) 
  thvcdf(hc,vc,Hm0,Tp,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 = thvcdf(hc,vc,H,T,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 = thvcdf(Hd,Vcf,Hm0,Tp,tail) 
002 %THVCDF Joint (Vcf,Hd) CDF for linear waves with Torsethaugen spectra. 
003 % 
004 %  CALL: f = thvcdf(Hd,Vcf,Hm0,Tp,tail) 
005 %  
006 %   f   = CDF evaluated at (Vcf,Hd) 
007 %   Hd  = zero down crossing wave height  
008 %   Vcf = crest front velocity 
009 %   Hm0 = significant wave height [m] 
010 %   Tp  = Spectral peak period    [s] 
011 %  tail = 1 if upper tail is calculated    
012 %         0 if lower tail is calulated (default) 
013 %   
014 % THVCDF approximates the joint CDF of (Vcf, Hd), i.e., crest front 
015 % velocity (Ac/Tcf) and wave height, for a Gaussian process with a 
016 % Torsethaugen spectral density. The empirical parameters of the model is 
017 % fitted by least squares to simulated (Vcf,Hd) data for 600 classes of 
018 % Hm0 and Tp. Between 50000 and 150000 zero-downcrossing waves were 
019 % simulated for each class of Hm0 and Tp. 
020 % THVCDF is restricted to the following range for Hm0 and Tp:  
021 %  0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11. 
022 % 
023 % Example: 
024 % Hm0 = 6;Tp = 8; 
025 % vc = 3; 
026 % hc = 3; 
027 % lowerTail = 0; 
028 % upperTail = ~lowerTail   
029 % thvcdf(hc,vc,Hm0,Tp)           % Prob(Hd<Hc,Vcf<Vc) 
030 % thvcdf(hc,vc,Hm0,Tp,upperTail) % Prob(Hd>Hc,Vcf>Vc)   
031 %   
032 %  % Conditional probability of steep and high waves given seastates 
033 %  % i.e., Prob(Hd>hc,Vcf>vc|Hs,Tp)   
034 %  upperTail = 1; 
035 %  Hs = linspace(2.5,11.5,10); 
036 %  Tp = linspace(4.5,19.5,16); 
037 %  [T,H] = meshgrid(Tp,Hs);  
038 %  p = thvcdf(hc,vc,H,T,upperTail); 
039 %  v = 10.^(-6:-1);   
040 %  contourf(Tp,Hs,log10(p),log10(v)) 
041 %  xlabel('Tp') 
042 %  ylabel('Hs')   
043 %  fcolorbar(log10(v))   
044 %   
045 % See also  thvpdf 
046  
047 % Reference  
048 % P. A. Brodtkorb (2004),   
049 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
050 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
051 % Trondheim, Norway.     
052    
053 % History 
054 % revised pab 09.08.2003 
055 % changed input   
056 % validated 20.11.2002 
057 % By pab 20.12.2000 
058  
059  
060 error(nargchk(3,5,nargin))   
061  
062 if (nargin < 5|isempty(tail)),tail = 0; end 
063 if (nargin < 4|isempty(Tp)),Tp = 8; end 
064 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end 
065  
066 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 
067 if multipleSeaStates 
068   [errorcode, Vcf,Hd,Hm0,Tp] = comnsize(Vcf,Hd,Hm0,Tp); 
069 else 
070   [errorcode, Vcf,Hd] = comnsize(Vcf,Hd); 
071 end 
072 if errorcode > 0 
073   error('Requires non-scalar arguments to match in size.'); 
074 end 
075  
076 global THVPAR 
077 if isempty(THVPAR) 
078   THVPAR = load('thvpar.mat'); 
079 end 
080  
081 Tpp  = THVPAR.Tp; 
082 Hm00 = THVPAR.Hm0; 
083 Tm020 = THVPAR.Tm02; 
084 % Interpolation method 
085 method = '*cubic';% Faster interpolation 
086  
087 [Tp1,Hs1] = meshgrid(Tpp,Hm00); 
088 Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 
089 %  w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
090   %w    = linspace(0,10,2*1024+1).';  
091 %  St = torsethaugen(w,[Hm0,Tp]); 
092 %  ch   = spec2char(St,{'Tm02','eps2'}); 
093 %  Tm02 = ch(1); 
094 %  eps2 = ch(2); 
095 Hrms = Hm0/sqrt(2); 
096 Vrms = 2*Hm0./Tm02; % Erms 
097  
098 v = Vcf./Vrms; 
099 f = zeros(size(Hd)); 
100  
101 % Only compute within valid range 
102 k0 = find((2<=Tp) & (Tp<=21) & (Hm0<=(Tp-2)*12/11) & (Hm0<=12)); 
103 if any(k0) 
104   hMax = 5; 
105   eps2 = 1e-6; 
106   h = min(Hd./Hrms,hMax); 
107    
108   if multipleSeaStates 
109     h = h(k0); 
110     v = v(k0); 
111     Hm0 = Hm0(k0); 
112     Tp = Tp(k0); 
113   else 
114     k0 = 1:prod(size(Hd)); 
115   end 
116   if 0 
117     % This is a trick to get the html documentation correct. 
118     k = thvpdf(1,1,2,3); 
119   end 
120   normalizedInput = 1; 
121   utprb = gaussq('thvpdf',hMax,2*hMax,eps2/2,[],mean(v(:)),mean(Hm0(:)),mean(Tp(:)),normalizedInput,7); 
122   if eps2<utprb 
123     warning('Check the accuracy of integration!') 
124   end 
125  
126   
127  
128   hlim    = h; 
129  
130  
131   lowerTail = 0; 
132   if tail==lowerTail, 
133     k       = find(h>2.5); 
134     hlim(k) = 2.5; 
135     k       = find(h>1.3*v); 
136     hlim(k) = 1.3*v(k); 
137     f(k0) = gaussq('thvpdf',0,hlim,eps2/2,[],v,Hm0,Tp,normalizedInput,5)... 
138     + gaussq('thvpdf',hlim,h,eps2/2,[],v,Hm0,Tp,normalizedInput,5);  
139   else % upper tail 
140     k       = find(h<1.3*v); 
141     hlim(k) = 1.3*v(k); 
142     f(k0) = gaussq('thvpdf',h,hlim,eps2/2,[],v,Hm0,Tp,normalizedInput,7)... 
143     + gaussq('thvpdf',hlim,hMax,eps2/2,[],v,Hm0,Tp,normalizedInput,7);  
144   end 
145 end 
146 return 
147  
148

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