Home > wafo > wavemodels > thscdf.m

thscdf

PURPOSE ^

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

SYNOPSIS ^

f = thscdf(Hd,Scf,Hm0,Tp,tail)

DESCRIPTION ^

 THSCDF Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra. 
  
   CALL: f = thscdf(Hd,Scf,Hm0,Tp,tail) 
   
    f   = CDF evaluated at (Scf,Hd) 
    Hd  = zero down crossing wave height (vector) 
    Scf = crest front steepness (vector)  
    Hm0 = significant wave height [m] 
    Tp  = Spectral peak period    [s] 
   tail = 1 if upper tail is calculated    
          0 if lower tail is calulated (default) 
    
  THSCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 
  steepness (2*pi*Ac/(g*Td*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 (Scf,Hd) data for 600 
  classes of Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves 
  were simulated for each class of Hm0 and Tp. 
  THSCDF 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; 
  sc = 0.25; 
  hc = 3; 
  lowerTail = 0; 
  upperTail = ~lowerTail   
  thscdf(hc,sc,Hm0,Tp)           % Prob(Hd<Hc,Scf<Ec) 
  thscdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>Hc,Scf>Ec)   
  
   % Conditional probability of steep and high waves given seastates 
   % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz)   
   upperTail = 1; 
   Hs = linspace(2.5,11.5,10); 
   Tp = linspace(4.5,19.5,16); 
   [T,H] = meshgrid(Tp,Hs);  
   p = thscdf(hc,sc,H,T,upperTail); 
   v = 10.^(-6:-1);   
   contourf(Tp,Hs,log10(p),log10(v)) 
   xlabel('Tp') 
   ylabel('Hs')   
   fcolorbar(log10(v))   
    
  See also  thspdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function f = thscdf(Hd,Scf,Hm0,Tp,tail) 
002 %THSCDF Joint (Scf,Hd) CDF for linear waves with Torsethaugen spectra. 
003 % 
004 %  CALL: f = thscdf(Hd,Scf,Hm0,Tp,tail) 
005 %  
006 %   f   = CDF evaluated at (Scf,Hd) 
007 %   Hd  = zero down crossing wave height (vector) 
008 %   Scf = crest front steepness (vector)  
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 % THSCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 
015 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process 
016 % with a Torsethaugen spectral density. The empirical parameters of the 
017 % model is fitted by least squares to simulated (Scf,Hd) data for 600 
018 % classes of Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves 
019 % were simulated for each class of Hm0 and Tp. 
020 % THSCDF 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 % sc = 0.25; 
026 % hc = 3; 
027 % lowerTail = 0; 
028 % upperTail = ~lowerTail   
029 % thscdf(hc,sc,Hm0,Tp)           % Prob(Hd<Hc,Scf<Ec) 
030 % thscdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>Hc,Scf>Ec)   
031 % 
032 %  % Conditional probability of steep and high waves given seastates 
033 %  % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz)   
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 = thscdf(hc,sc,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  thspdf 
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, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp); 
069 else 
070   [errorcode, Scf,Hd] = comnsize(Scf,Hd); 
071 end 
072 if errorcode > 0 
073   error('Requires non-scalar arguments to match in size.'); 
074 end 
075  
076 global THSPAR 
077 if isempty(THSPAR) 
078   %THSPAR = load('thspar.mat'); 
079   THSPAR = load('thspar19-Jul-2004.mat'); 
080 end 
081  
082 Tpp  = THSPAR.Tp; 
083 Hm00 = THSPAR.Hm0; 
084 Tm020 = THSPAR.Tm02; 
085 % Interpolation method 
086 method = '*cubic';% Faster interpolation 
087  
088 [Tp1,Hs1] = meshgrid(Tpp,Hm00); 
089 if 1, 
090   Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 
091 else 
092   Tm02 = Tp; 
093   for ix= 1:100 
094     dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method)); 
095     Tp = Tp+dTp; 
096     if all(abs(dTp)<0.01) 
097       %dTp 
098       ix 
099       break 
100     end 
101   end 
102 end 
103 %  w = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
104   %w    = linspace(0,10,2*1024+1).';  
105 %  St = torsethaugen(w,[Hm0,Tp]); 
106 %  ch   = spec2char(St,{'Tm02','eps2'}); 
107 %  Tm02 = ch(1); 
108 %  eps2 = ch(2); 
109 Hrms = Hm0/sqrt(2); 
110 Erms = 1.25*Hm0./(Tm02.^2); % Erms 
111  
112 s = Scf./Erms; 
113 hMax = 5; 
114 h = min(Hd./Hrms,hMax); 
115  
116 eps2 = 1e-6; 
117  
118 f = zeros(size(Hd)); 
119  
120 % Only compute within valid range 
121 %k0 = find((2<=Tp) & (Tp<=21) & (Hm0<=(Tp-2)*12/11) & (Hm0<=12)); 
122 upLimit = 6.5; 
123 loLimit = 2.5; 
124 k0 = find((2<=Tp) & (Tp<=21) & (loLimit*sqrt(Hm0)<Tp) & (Hm0<=12)); 
125 if any(k0) 
126   if multipleSeaStates 
127     h = h(k0); 
128     s = s(k0); 
129     Hm0 = Hm0(k0); 
130     Tp = Tp(k0); 
131   else 
132     k0 = 1:prod(size(Hd)); 
133   end 
134  
135   hlim    = h; 
136  
137   normalizedInput = 1; 
138   lowerTail = 0; 
139   if 0 
140     % This is a trick to get the html documentation correct. 
141     k = thspdf(1,1,2,3); 
142   end 
143    
144   if tail==lowerTail 
145     k       = find(h>2.5); 
146     hlim(k) = 2.5; 
147     f(k0) = gaussq('thspdf',0,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,5)... 
148     + gaussq('thspdf',hlim,h,eps2/2,[],s,Hm0,Tp,normalizedInput,5);  
149   else % upper tail 
150     k       = find(h<2.5); 
151     hlim(k) = 2.5; 
152     f(k0) = gaussq('thspdf',h,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,7)... 
153     + gaussq('thspdf',hlim,hMax,eps2/2,[],s,Hm0,Tp,normalizedInput,7);  
154   end 
155 end 
156 return 
157  
158

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