Home > wafo > wavemodels > thsnlcdf.m

thsnlcdf

PURPOSE ^

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

SYNOPSIS ^

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

DESCRIPTION ^

 THSNLCDF Joint (Scf,Hd) CDF for nonlinear waves with Torsethaugen spectra. 
  
   CALL: f = thsnlcdf(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) 
    
  THSNLCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 
  steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order 
  nonlinear waves 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. 
  THSNLCDF 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   
  thsnlcdf(hc,sc,Hm0,Tp)           % Prob(Hd<hc,Scf<sc) 
  thsnlcdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>hc,Scf>sc)   
  
   % 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 = thsnlcdf(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  thsnlpdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function f = thsnlcdf(Hd,Scf,Hm0,Tp,tail) 
002 %THSNLCDF Joint (Scf,Hd) CDF for nonlinear waves with Torsethaugen spectra. 
003 % 
004 %  CALL: f = thsnlcdf(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 % THSNLCDF approximates the joint CDF of (Scf, Hd), i.e., crest front 
015 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for 2nd order 
016 % nonlinear waves with a Torsethaugen spectral density. 
017 %  The empirical parameters of the 
018 % model is fitted by least squares to simulated (Scf,Hd) data for 600 
019 % classes of Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves 
020 % were simulated for each class of Hm0 and Tp. 
021 % THSNLCDF is restricted to the following range for Hm0 and Tp:  
022 %  0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11. 
023 % 
024 % Example: 
025 % Hm0 = 6;Tp = 8; 
026 % sc = 0.25; 
027 % hc = 3; 
028 % lowerTail = 0; 
029 % upperTail = ~lowerTail   
030 % thsnlcdf(hc,sc,Hm0,Tp)           % Prob(Hd<hc,Scf<sc) 
031 % thsnlcdf(hc,sc,Hm0,Tp,upperTail) % Prob(Hd>hc,Scf>sc)   
032 % 
033 %  % Conditional probability of steep and high waves given seastates 
034 %  % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz)   
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 = thsnlcdf(hc,sc,H,T,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  thsnlpdf 
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 % revised pab 09.08.2003 
056 % changed input   
057 % validated 20.11.2002 
058 % By pab 20.12.2000 
059  
060  
061 error(nargchk(3,5,nargin))   
062  
063 if (nargin < 5|isempty(tail)),tail = 0; end 
064 if (nargin < 4|isempty(Tp)),Tp = 8; end 
065 if (nargin < 3|isempty(Hm0)), Hm0 = 6; end 
066  
067 multipleSeaStates = any(prod(size(Hm0))>1|prod(size(Tp))>1); 
068 if multipleSeaStates 
069   [errorcode, Scf,Hd,Hm0,Tp] = comnsize(Scf,Hd,Hm0,Tp); 
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  
077 global THSNLPAR 
078 if isempty(THSNLPAR) 
079   %THSNLPAR = load('thsnlpar.mat'); 
080   THSNLPAR = load('thsnlpar06-Aug-2004.mat')   
081 end 
082  
083 Tpp  = THSNLPAR.Tp; 
084 Hm00 = THSNLPAR.Hm0; 
085 Tm020 = THSNLPAR.Tm02; 
086 % Interpolation method 
087 method = '*cubic';% Faster interpolation 
088  
089 [Tp1,Hs1] = meshgrid(Tpp,Hm00); 
090 if 1, 
091   Tm02 = interp2(Tp1,Hs1,Tm020,Tp,Hm0,method); 
092 else 
093      
094   Tm02 = Tp; 
095   for ix= 1:100 
096     dTp = (Tm02-interp2(Tp1,Hs1,Tm020,Tp,Hm0,method)); 
097     Tp = Tp+dTp; 
098     if all(abs(dTp)<0.01) 
099       %dTp 
100       ix 
101       break 
102     end 
103   end 
104 end 
105  
106  
107 %  w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
108   %w    = linspace(0,10,2*1024+1).';  
109 %  St = torsethaugen(w,[Hm0,Tp]); 
110 %  ch   = spec2char(St,{'Tm02','eps2'}); 
111 %  Tm02 = ch(1); 
112 %  eps2 = ch(2); 
113 Hrms = Hm0/sqrt(2); 
114 Erms = 1.25*Hm0./(Tm02.^2); % Erms 
115  
116 s = Scf./Erms; 
117 hMax = 5; 
118 h = min(Hd./Hrms,hMax); 
119  
120 eps2 = 1e-6; 
121  
122 f = zeros(size(Hd)); 
123  
124 % Only compute within valid range 
125 %k0 = find((2<=Tp) & (Tp<=21) & (Hm0<=(Tp-2)*12/11) & (Hm0<=12)); 
126 upLimit = 6.5; 
127 loLimit = 2.5; 
128 k0 = find((2<=Tp) & (Tp<=21) & (loLimit*sqrt(Hm0)<Tp) & (Hm0<=12)); 
129 if any(k0) 
130   if multipleSeaStates 
131     h = h(k0); 
132     s = s(k0); 
133     Hm0 = Hm0(k0); 
134     Tp = Tp(k0); 
135   else 
136     k0 = 1:prod(size(Hd)); 
137   end 
138  
139   hlim    = h; 
140  
141   normalizedInput = 1; 
142   lowerTail = 0; 
143    
144   if 0 
145     % This is a trick to get the html documentation correct. 
146     k = thsnlpdf(1,1,2,3); 
147   end 
148    
149   if tail==lowerTail 
150     k       = find(h>2.5); 
151     hlim(k) = 2.5; 
152     f(k0) = gaussq('thsnlpdf',0,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,5)... 
153     + gaussq('thsnlpdf',hlim,h,eps2/2,[],s,Hm0,Tp,normalizedInput,5);  
154   else % upper tail 
155     k       = find(h<2.5); 
156     hlim(k) = 2.5; 
157     f(k0) = gaussq('thsnlpdf',h,hlim,eps2/2,[],s,Hm0,Tp,normalizedInput,7)... 
158     + gaussq('thsnlpdf',hlim,hMax,eps2/2,[],s,Hm0,Tp,normalizedInput,7);  
159   end 
160 end 
161 return 
162  
163

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