Home > wafo > wavemodels > thpdf.m

thpdf

PURPOSE ^

Marginal wave height, Hd, PDF for Torsethaugen spectra.

SYNOPSIS ^

f = thpdf(h,Hm0,Tp,dim)

DESCRIPTION ^

 THPDF Marginal wave height, Hd, PDF for Torsethaugen spectra. 
  
   CALL: f = thpdf(h,Hm0,Tp,dim) 
   
   f   = pdf evaluated at h. 
   h   = vectors of evaluation points. 
   Hm0 = significant wave height [m]. 
   Tp  = Spectral peak period    [s]. 
   dim = 'time'  : Hd distribution in time (default) 
         'space' : Hd distribution in space 
  
  THPDF approximates the marginal distribution of Hd, i.e., 
  zero-downcrossing wave height, for a Gaussian process with a Torsethaugen 
  spectral density. The empirical parameters of the model is fitted by 
  least squares to simulated 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. 
  THPDF 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; 
  h = linspace(0,4*Hm0/sqrt(2))';  
  f = thpdf(h,Hm0,Tp); 
  dt = 0.4; w = linspace(0,2*pi/dt,256)'; 
  S = torsethaugen(w,[Hm0 Tp]); 
  xs = spec2sdat(S,20000,dt); rate=8; method=1; 
  [S,H] = dat2steep(xs,rate,method); 
  fk = kdebin(H,'epan',[],[],.5,128); 
  plot(h,f), hold on, pdfplot(fk,'r'), hold off 
  
  See also  thvpdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function f = thpdf(h,Hm0,Tp,dim) 
002 %THPDF Marginal wave height, Hd, PDF for Torsethaugen spectra. 
003 % 
004 %  CALL: f = thpdf(h,Hm0,Tp,dim) 
005 %  
006 %  f   = pdf evaluated at h. 
007 %  h   = vectors of evaluation points. 
008 %  Hm0 = significant wave height [m]. 
009 %  Tp  = Spectral peak period    [s]. 
010 %  dim = 'time'  : Hd distribution in time (default) 
011 %        'space' : Hd distribution in space 
012 % 
013 % THPDF approximates the marginal distribution of Hd, i.e., 
014 % zero-downcrossing wave height, for a Gaussian process with a Torsethaugen 
015 % spectral density. The empirical parameters of the model is fitted by 
016 % least squares to simulated Hd data for 600 classes of Hm0 and 
017 % Tp. Between 50000 and 150000 zero-downcrossing waves were simulated for 
018 % each class of Hm0 and Tp. 
019 % THPDF is restricted to the following range for Hm0 and Tp:  
020 %  0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11. 
021 % 
022 % Example: 
023 % Hm0 = 6;Tp = 8; 
024 % h = linspace(0,4*Hm0/sqrt(2))';  
025 % f = thpdf(h,Hm0,Tp); 
026 % dt = 0.4; w = linspace(0,2*pi/dt,256)'; 
027 % S = torsethaugen(w,[Hm0 Tp]); 
028 % xs = spec2sdat(S,20000,dt); rate=8; method=1; 
029 % [S,H] = dat2steep(xs,rate,method); 
030 % fk = kdebin(H,'epan',[],[],.5,128); 
031 % plot(h,f), hold on, pdfplot(fk,'r'), hold off 
032 % 
033 % See also  thvpdf 
034  
035  
036 % Reference   
037 % P. A. Brodtkorb (2004),   
038 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
039 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
040 % Trondheim, Norway.   
041    
042 % History 
043 % revised pab jan2004   
044 % By pab 20.12.2000 
045  
046  
047  
048 error(nargchk(3,4,nargin)) 
049  
050 if nargin<4|isempty(dim), 
051   dim = 'time';% dim='time'->wtweibpdf, dim='space'->wggampdf 
052 end 
053  
054 if Hm0>12| Hm0>(Tp-2)*12/11  
055   disp('Warning: Hm0 is outside the valid range') 
056   disp('The validity of the Hd distribution is questionable') 
057 end 
058 if Tp>20|Tp<3  
059   disp('Warning: Tp is outside the valid range') 
060   disp('The validity of the Hd distribution is questionable') 
061 end 
062 Hrms = Hm0/sqrt(2); 
063  
064 [a b c] = thwparfun(Hm0,Tp,dim); 
065   f = wtweibpdf(h/Hrms,a,b,c)/Hrms; 
066    
067   return 
068   % old code kept just in case 
069 if strncmpi(dim,'t',1)     
070   [a b c] = thwparfun(Hm0,Tp,dim); 
071   f = wtweibpdf(h/Hrms,a,b,c)/Hrms; 
072 else  
073   [a b c] = thgparfun(Hm0,Tp,dim); 
074   f = wggampdf(h/Hrms,a,b,c)/Hrms; 
075 end 
076  
077 return 
078  
079 % old code: saved just in case 
080 % Weibull distribution parameters as a function of e2 and h2 
081 % Hrms = Hm0/sqrt(2); 
082 % if nargin<1 |isempty(h), h=linspace(0,4*Hrms).'; end 
083 % if nargin>4, 
084 %   h = h*Hrms; 
085 %   if pdef==1 
086 %     f = hggam(h,Hm0,Tp,eps2,1); 
087 %   else 
088 %     f = htweib(h,Hm0,Tp,eps2,1); 
089 %   end 
090 % elseif pdef==1 
091 %   f = hggam(h,Hm0,Tp,eps2); 
092 % else 
093 %   f = htweib(h,Hm0,Tp,eps2);  
094 % end 
095  
096 % f=f.f; 
097 % return 
098  
099 % function f = hggam(h,Hm0,Tp,eps2,norm) 
100 %   pardef =2; 
101 %   if pardef == 1, 
102 %     if nargin<4|isempty(eps2), 
103 %       w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
104 %       %w    = linspace(0,1000/Tp,8*1024+1).'; % adaptiv spacing 
105 %       eps2   = spec2char(torsethaugen(w,[Hm0,Tp]),{'eps2'}); 
106 %     end 
107 %     if eps2<0.4  
108 %       disp('Warning: eps2 is less than 0.4. The pdf returned is questionable') 
109 %     elseif eps2>1.3 
110 %       disp('Warning: eps2 is larger than 1.3. The pdf returned is questionable') 
111 %     end  
112 %   end 
113  
114 % switch pardef 
115 %   case 1,% Wggampdf distribution parameters as a function of eps2 
116 %     % Best fit by smoothing spline forcing a=1, b=1 and c=2 for eps2=0.  
117 %     % Then approximate the spline with a rational polynomial 
118      
119 %     da1 = [0.3579    0.9601   -1.8152    1.0009]; 
120 %     da2 = [1.5625   -0.1537   -1.1390    1.0000]; 
121 %     db1 = [-1.8555   3.5375   -3.4056    1.0026];  
122 %     db2 = [-2.0855   4.1624   -3.6399    1.0000]; 
123 %     dc1 = [6.9572   -6.7932    2.1760    1.9974]; 
124 %     dc2 = [3.0979   -3.0781    0.5369    1.0000]; 
125 %     A0 = polyval(da1,eps2)./polyval(da2,eps2); 
126 %     B0 = polyval(db1,eps2)./polyval(db2,eps2); 
127 %     C0 = polyval(dc1,eps2)./polyval(dc2,eps2); 
128 %   case 2, % wave height distribution in time 
129 %     global THGPAR 
130 %     if isempty(THGPAR) 
131 %       THGPAR = load('thgpar.mat'); 
132 %     end 
133 %     % Generalized Gamma  distribution parameters as a function of Tp, Hm0  
134 %     A00 = THGPAR.A00s; 
135 %     B00 = THGPAR.B00s; 
136 %     C00 = THGPAR.C00s; 
137  
138 %     Tpp  = THGPAR.Tp; 
139 %     Hm00 = THGPAR.Hm0; 
140 %     [E1, H1] = meshgrid(Tpp,Hm00); 
141 %     method = '*cubic'; 
142 %     A0 = interp2(E1,H1,A00,Tp,Hm0,method); 
143 %     B0 = interp2(E1,H1,B00,Tp,Hm0,method); 
144 %     C0 = interp2(E1,H1,C00,Tp,Hm0,method); 
145 %   case 3,% Waveheight distribution in space 
146 %      global THSSPAR 
147 %     if isempty(THSSPAR) 
148 %       THSSPAR = load('thsspar.mat'); 
149 %     end 
150 %     % Generalized Gamma  distribution parameters as a function of Tp, Hm0  
151 %     A00 = THSSPAR.A00s; 
152 %     B00 = THSSPAR.B00s; 
153 %     C00 = THSSPAR.C00s; 
154  
155 %     Tpp  = THSSPAR.Tp; 
156 %     Hm00 = THSSPAR.Hm0; 
157 %     [E1, H1] = meshgrid(Tpp,Hm00); 
158 %     method = '*cubic'; 
159 %     A0 = interp2(E1,H1,A00,Tp,Hm0,method); 
160 %     B0 = interp2(E1,H1,B00,Tp,Hm0,method); 
161 %     C0 = interp2(E1,H1,C00,Tp,Hm0,method); 
162 % end 
163 % Hrms = Hm0/sqrt(2); 
164 % f.f    = wggampdf(h/Hrms,A0,B0,C0)/Hrms; 
165 % return 
166  
167 % function f = htweib(h,Hm0,Tp,eps2,norm) 
168 %  pardef =7; 
169 %   if pardef < 7, 
170 %     if nargin<4|isempty(eps2), 
171 %       w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
172 %       %w    = linspace(0,1000/Tp,8*1024+1).'; % adaptiv spacing 
173 %       S = torsethaugen(w,[Hm0,Tp]); 
174 %       eps2   = spec2char(S,{'eps2'}); 
175 %     end 
176 %     if eps2<0.4  
177 %       disp('Warning: eps2 is less than 0.4. The pdf returned is questionable') 
178 %     elseif eps2>1.3 
179 %       disp('Warning: eps2 is larger than 1.3. The pdf returned is questionable') 
180 %     end  
181 %   end 
182 % switch pardef, 
183 %   case 1,% Wtweibpdf distribution parameters as a function of eps2 
184 %     % Best fit by smoothing spline forcing a=1, b=2 and c=0 for eps2=0.  
185 %     brks = [0 .1 .2 .4 ,.6, .8, 1, 1.1 1.2]'; 
186 %     coefa = [0                  0   0.02260415153596   0.99807186986167; ... 
187 %    2.19065400617385                  0   0.02260415153596  1.00033228501527; ... 
188 %    4.34015195156053   0.65719620185215   0.03709199393185    1.00478335417504; ... 
189 %   -1.59533089716870   3.26128737278847   0.80983543882910   1.07321081664798; ... 
190 %   -6.81273221810880   2.30408883448726   1.92291068028425   1.35286675214799; ... 
191 %   -3.69498826658975  -1.78355049637802   2.09369407217829   1.77511058383946; ... 
192 %   13.33514485443956  -4.00054345633187   0.94471547491809   2.09294747228728; ... 
193 %                   0                  0   0.54466112928490   2.16074873007021]; 
194          
195 %   coefb = [ 0                  0   0.32503235228616   1.99054481866418; ... 
196 %    3.28321899128157                  0   0.32503235228616    2.02304805389280; ... 
197 %    5.67672309005450   0.98496569738447   0.37964649056830    2.05883450811270; ... 
198 %   -5.29907238080822   4.39099955141717   1.43842344537222   2.21957621884217; .... 
199 %   -5.89663569823287   1.21155612293224   2.55893458024211   2.64050831092684; ... 
200 %   -6.21824739906323  -2.32642529600749   2.43691697455115   3.15358438630669; ... 
201 %   20.19124578481806  -6.05737373544542   0.77134599291113   3.49816479018411; ... 
202 %                   0                  0   0.16560861936659   3.53491689790559]; 
203 % coefc =[                0                  0   0.04818579357214       -0.00817761487085; ... 
204 %    2.94432030165157                  0   0.04818579357214    -0.00335903551363; ... 
205 %    4.77660844045250   0.88329609049547   0.09917317190900    0.00440386414523; ... 
206 %   -1.24578770271258   3.74926115476697   1.01096301945323   0.09778320967047; ... 
207 %   -7.70868155645400   3.00178853313943   2.36117295703451   0.43997995813009; ... 
208 %   -3.98346578867600  -1.62342040073298   2.70373824808144   0.97061663841094; ... 
209 %   13.37833291312857  -4.01349987393858   1.57933014446810   1.41455974568850; ... 
210 %                   0                  0   1.17798015707424   1.54573609430905]; 
211 %   pa = mkpp(brks,coefa); 
212 %   pb = mkpp(brks,coefb); 
213 %   pc = mkpp(brks,coefc); 
214 %   A0 = ppval(pa,eps2);         
215 %   B0 = ppval(pb,eps2);        
216 %   C0 = ppval(pc,eps2);        
217 %   case 2,% Wtweibpdf distribution parameters as a function of eps2 
218 %     % Best fit by smoothing spline forcing a=1, b=2 and c=0 for eps2=0.  
219 %     % Then approximate the spline with a rational polynomial. 
220      
221 %     da1 = [1.8262  -2.1141   1.0034]; 
222 %     da2 = [-0.1371  0.1626   1.3189  -2.0016   1.0000]; 
223 %     db1 = [2.4360  -1.6331   1.9932]; 
224 %     db2 = [-1.3266  4.3755  -4.6795   2.4763  -1.0450  1.0000]; 
225 %     dc1 = [0.3162  -0.1647   0.0803  -0.0104]; 
226 %     dc2 = [-1.4325  6.3773 -11.1967  10.1490  -4.7401  1.0000]; 
227 %     A0 = polyval(da1,eps2)./polyval(da2,eps2); 
228 %     B0 = polyval(db1,eps2)./polyval(db2,eps2); 
229 %     C0 = polyval(dc1,eps2)./polyval(dc2,eps2); 
230 %   case 3, % Wtweibpdf distribution parameters as a function of eps2  
231 %     % LS fit to data 
232 %     % best fit to torsethaugen for eps2 = 0.4 to 1.2; 
233 %     % and forcing through A0=1 B0=2 C0=0 for eps2==0)  
234 %     A0 = sqrt(1+2.0094*eps2.^1.8492); 
235 %     B0 = 1+sqrt(1+3.4117*eps2.^1.4972); 
236 %     C0 = 0.9994*eps2.^1.6728;  
237 %   case 4,% Wtweibpdf distribution parameters as a function of eps2 
238 %     % best fit to torsethaugen eps2 = 0.4 to 1.2; 
239 %     A0 = 0.7315+eps2.^0.9916; 
240 %     B0 = 2.0628+eps2.^1.1927; 
241 %     C0 = 0.9994*eps2.^1.6728; 
242 %   case 5,% Naess (1985) 
243 %     R  = spec2cov(S); 
244 % %    A0 = sqrt((1-min(R.R)/R.R(1))/2);% Naess (1985) 
245 %     A0 = sqrt((1-min(R.R)/R.R(1))/2)+0.03;% Modified approach broadbanded time 
246 % %    A0 = sqrt((1-min(R.R)/R.R(1))/2)+0.1;% Modified approach broadbanded space                                  
247                        
248 %     B0 = 2; 
249 %     C0 = 0; 
250 %   case 6,% Naess (1985) wave height in space 
251 %     R  = spec2cov(spec2spec(specinterp(S,0.55),'k1d')); 
252 %     %A0 = sqrt((1-min(R.R)/R.R(1))/2); % Naess (1985) 
253 %     A0 = sqrt((1-min(R.R)/R.R(1))/2)+0.03; % modified approach     
254 %     B0 = 2; 
255 %     C0 = 0;  
256 %   case 7,% wave height distribution in time 
257 %     global THWPAR 
258 %     if isempty(THWPAR) 
259 %       THWPAR = load('thwpar.mat'); 
260 %     end 
261 %     % Truncated Weibull  distribution parameters as a function of Tp, Hm0  
262 %     A00 = THWPAR.A00s; 
263 %     B00 = THWPAR.B00s; 
264 %     C00 = THWPAR.C00s; 
265  
266 %     Tpp  = THWPAR.Tp; 
267 %     Hm00 = THWPAR.Hm0; 
268 %     [E1, H1] = meshgrid(Tpp,Hm00); 
269 %     method = '*cubic'; 
270 %     A0 = interp2(E1,H1,A00,Tp,Hm0,method); 
271 %     B0 = interp2(E1,H1,B00,Tp,Hm0,method); 
272 %     C0 = interp2(E1,H1,C00,Tp,Hm0,method); 
273   
274 % end 
275  
276 % Hrms = Hm0/sqrt(2); 
277 % f = wtweibpdf(h/Hrms,A0,B0,C0)/Hrms; 
278  
279 % return 
280

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