Home > wafo > wavemodels > thspdf2.m

thspdf2

PURPOSE ^

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

SYNOPSIS ^

[f,varargout] = thspdf2(Hd,Scf,Hm0,Tp,normalizedInput)

DESCRIPTION ^

 THSPDF2 Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. 
  
   CALL: f = thspdf2(Hd,Scf,Hm0,Tp) 
   
    f   = pdf struct evaluated at meshgrid(Scf,Hd) 
    Hd  = zero down crossing wave height (vector) 
    Scf = crest front steepness (vector)  
    Hm0 = significant wave height [m] 
    Tp  = Spectral peak period    [s] 
  
  THSPDF2 approximates the joint distribution of (Scf, Hd), i.e., crest 
  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. 
  THSPDF 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));  
  s = linspace(0,6*1.25*Hm0/Tp^2); 
  f = thspdf2(h,s,Hm0,Tp); 
  w = linspace(0,40,5*1024+1).'; 
  S = torsethaugen(w,[Hm0 Tp]); 
  dt = 0.3; 
  x = spec2sdat(S,80000,.2); rate = 8; 
  [si,hi] = dat2steep(x,rate,2); 
  fk = kdebin([si,hi],'epan',[],[],.5,128); 
   fk.title = f.title; fk.labx = f.labx;  
  plot(si,hi,'.'), hold on 
  pdfplot(f),pdfplot(fk,'r'),hold off 
  
  See also  thsspdf

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [f,varargout] = thspdf2(Hd,Scf,Hm0,Tp,normalizedInput) 
002 %THSPDF2 Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. 
003 % 
004 %  CALL: f = thspdf2(Hd,Scf,Hm0,Tp) 
005 %  
006 %   f   = pdf struct evaluated at meshgrid(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 % 
012 % THSPDF2 approximates the joint distribution of (Scf, Hd), i.e., crest 
013 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process with a 
014 % Torsethaugen spectral density. The empirical parameters of the model is 
015 % fitted by least squares to simulated (Scf,Hd) data for 600 classes of 
016 % Hm0 and Tp. Between 40000 and 200000 zero-downcrossing waves were 
017 % simulated for each class of Hm0 and Tp. 
018 % THSPDF is restricted to the following range for Hm0 and Tp:  
019 %  0.5 < Hm0 [m] < 12,  3.5 < Tp [s] < 20,  and  Hm0 < (Tp-2)*12/11. 
020 % 
021 % Example: 
022 % Hm0 = 6;Tp = 8; 
023 % h = linspace(0,4*Hm0/sqrt(2));  
024 % s = linspace(0,6*1.25*Hm0/Tp^2); 
025 % f = thspdf2(h,s,Hm0,Tp); 
026 % w = linspace(0,40,5*1024+1).'; 
027 % S = torsethaugen(w,[Hm0 Tp]); 
028 % dt = 0.3; 
029 % x = spec2sdat(S,80000,.2); rate = 8; 
030 % [si,hi] = dat2steep(x,rate,2); 
031 % fk = kdebin([si,hi],'epan',[],[],.5,128); 
032 %  fk.title = f.title; fk.labx = f.labx;  
033 % plot(si,hi,'.'), hold on 
034 % pdfplot(f),pdfplot(fk,'r'),hold off 
035 % 
036 % See also  thsspdf 
037  
038    
039 % Reference   
040 % P. A. Brodtkorb (2004),   
041 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
042 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
043 % Trondheim, Norway.   
044    
045 % History 
046 % revised pab 09.08.2003 
047 % changed input   
048 % validated 20.11.2002 
049 % By pab 20.12.2000 
050  
051 error(nargchk(4,5,nargin)) 
052  
053 if (nargin < 5|isempty(normalizedInput)),  normalizedInput  = 0;end 
054 if (nargin < 4|isempty(Tp)),  Tp  = 8;end 
055 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end 
056  
057  
058 [V,H] = meshgrid(Scf,Hd); 
059  
060 f = createpdf(2); 
061 [f.f,Hrms,Vrms,varargout{1:nargout-1}]  = thspdf(H,V,Hm0,Tp,normalizedInput); 
062  
063  f.x = {Scf(:),Hd(:)}; 
064   
065 if (normalizedInput) 
066   f.labx={'Scf', 'Hd'}; 
067   f.norm = 1; 
068 else 
069   f.norm=0; 
070   f.labx={'Scf', 'Hd [m]'}; 
071 end 
072 f.title = 'Joint distribution of (Hd,Scf) in time'; 
073 f.note = ['Torsethaugen Hm0=' num2str(Hm0) ' Tp = ' num2str(Tp)]; 
074 [f.cl,f.pl] = qlevels(f.f); 
075  
076 return  
077 % old call 
078 if Hm0>11| Hm0>(Tp-2)*12/11  
079   disp('Warning: Hm0 is outside the valid range') 
080   disp('The validity of the Joint (Hd,Scf) distribution is questionable') 
081 end 
082 if Tp>20|Tp<3  
083   disp('Warning: Tp is outside the valid range') 
084   disp('The validity of the Joint (Hd,Scf) distribution is questionable') 
085 end 
086  
087 global THSPAR 
088 if isempty(THSPAR) 
089   THSPAR = load('thspar.mat'); 
090 end 
091 % Gamma distribution parameters as a function of Tp Hm0 and h2 
092 A11 = THSPAR.A11s; 
093 B11 = THSPAR.B11s; 
094  
095 % Waveheight distribution in time 
096  
097 if 0, 
098   % Truncated Weibull  distribution parameters as a function of Tp, Hm0  
099   global THWPAR 
100   if isempty(THWPAR) 
101     THWPAR = load('thwpar.mat'); 
102   end 
103   A00 = THWPAR.A00s; 
104   B00 = THWPAR.B00s; 
105   C00 = THWPAR.C00s; 
106 else 
107  
108   % Truncated Weibull  distribution parameters as a function of Tp, Hm0  
109   A00 = THSPAR.A00s; 
110   B00 = THSPAR.B00s; 
111   C00 = THSPAR.C00s; 
112 end 
113  
114 Tpp  = THSPAR.Tp; 
115 Hm00 = THSPAR.Hm0; 
116 h2   = THSPAR.h2; 
117  
118  
119  
120  
121 w    = linspace(0,100,16*1024+1).'; % torsethaugen original spacing 
122 %w    = linspace(0,10,2*1024+1).';  
123 S = torsethaugen(w,[Hm0,Tp]); 
124 ch   = spec2char(S,{'Tm02','eps2'}); 
125 Tm02 = ch(1); 
126 eps2 = ch(2); 
127  
128 Hrms = Hm0/sqrt(2); 
129 Vrms = 1.25*Hm0/(Tm02^2); % Erms 
130  
131 if nargin<1 |isempty(v), v=linspace(0,4*Vrms); end 
132 if nargin<2 |isempty(h), h=linspace(0,4*Hrms); end 
133  
134 if nargin>4, 
135   h = h*Hrms; 
136   v = v*Vrms; 
137 end 
138  
139 %Fh = thpdf(h(:)/Hrms,Hm0,Tp,eps2,1); 
140  
141 [A0, B0, C0] = thwparfun(Hm0,Tp,'time'); 
142  
143 method = '*cubic';% Faster interpolation 
144 [E1, H1, H2] = meshgrid(Tpp,Hm00,h2); 
145 A1 = exp(smooth(h2,interp3(E1,H1,H2,log(A11),Tp(ones(size(h2))),... 
146     Hm0(ones(size(h2))) ,h2,method),1,h/Hrms,1)); 
147 B1 = exp(smooth(h2,interp3(E1,H1,H2,log(B11),Tp(ones(size(h2))),... 
148      Hm0(ones(size(h2))),h2,method),1,h/Hrms,1)); 
149  
150 [V1 H1] = meshgrid(v/Vrms,h/Hrms); 
151 [V1 A1] = meshgrid(v/Vrms,A1); 
152 [V1 B1] = meshgrid(v/Vrms,B1); 
153  
154 f = createpdf(2); 
155 f.title = 'Joint distribution of (Hd,Scf) in time'; 
156  
157 if nargin<5  
158    
159    f.f = wtweibpdf(H1,A0,B0,C0).*wgampdf(V1,A1,B1)/Vrms/Hrms; 
160  % f.f = repmat(Fh/Hrms,[1 length(v)]).*wgampdf(V1,A1,B1)/Vrms; 
161   f.x = {v,h}; 
162   f.norm=0; 
163   f.labx={'Scf', 'Hd [m]'}; 
164 else 
165   f.f = wtweibpdf(H1,A0,B0,C0).*wgampdf(V1,A1,B1); 
166   %f.f = repmat(Fh,[1 length(v)]).*wgampdf(V1,A1,B1); 
167   f.x = {v/Vrms,h/Hrms}; 
168   f.labx={'Scf', 'Hd'}; 
169   f.norm = 1; 
170 end 
171 f.f(find(isnan(f.f)|isinf(f.f) ))=0; 
172  
173 f.note = ['Torsethaugen Hm0=' num2str(Hm0) ' Tp = ' num2str(Tp)]; 
174 [f.cl,f.pl] = qlevels(f.f); 
175 if nargout>1, 
176   fA      = createpdf(2); 
177   fA.x    = {Tpp,Hm00}; 
178   fA.labx = {'Tp', 'Hm0'}; 
179   fA(3)   = fA(1); 
180   fA(2)   = fA(1); 
181    
182   fA(1).f    = A00; 
183   fA(2).f    = B00; 
184   fA(3).f    = C00; 
185    
186   fA(1).title = 'wtweibpdf parameter A'; 
187   fA(2).title = 'wtweibpdf parameter B'; 
188   fA(3).title = 'wtweibpdf parameter C'; 
189    
190   txt1 = ['The Wtweibpdf  distribution Parameter ']; 
191   txt2=[' for wave heigth in time as a function of Tp and Hm0 for' ... 
192     'the Torsethaugen spectrum']; 
193   fA(1).note =[txt1 'A' txt2]; 
194   fA(2).note =[txt1 'B' txt2]; 
195   fA(3).note =[txt1 'C' txt2]; 
196    
197   tmp = [A00(:) B00(:) C00(:)]; 
198   ra  = range(tmp); 
199   st  = round(min(tmp)*100)/100; 
200   en  = max(tmp); 
201   for ix = 1:3, 
202     fA(ix).cl   = st(ix):ra(ix)/20:en(ix); 
203   end 
204 end 
205 if nargout>2, 
206   fB      = createpdf(3); 
207   fB.x    = {Tpp,Hm00,h2}; 
208   fB.labx = {'Tp','Hm0', 'h'}; 
209   fB(2)   = fB(1); 
210    
211   fB(1).f = A11; 
212   fB(2).f = B11; 
213    
214   txt11 = 'The conditonal Wgampdf distribution Parameter '; 
215   txt22 = [' for Scf given h=Hd/Hrms in time as function of Tp' ... 
216     ' and Hm0 for the Torsethaugen spectrum']; 
217   fB(1).title = 'wgampdf parameter A'; 
218   fB(2).title = 'wgampdf parameter B'; 
219   fB(1).note = [txt11,'A',txt22]; 
220   fB(2).note = [txt11,'B',txt22]; 
221    
222   %fB(2).note = ['The conditonal Wggampdf distribution Parameter B(h)/Hrms', ...%    ' for crest front steepness as a function of Tp,Hm0 and',... 
223   %    ' h=Hd/Hrms for the Torsethaugen spectrum']; 
224   tmp= [A11(:) B11(:)]; 
225   ra = range(tmp); 
226   st = round(min(tmp)*100)/100; 
227   en = max(tmp); 
228   for ix=1:2 
229     fB(ix).cl   = st(ix):ra(ix)/20:en(ix); 
230   end 
231 end 
232 return 
233  
234  
235  
236  
237  
238  
239 %Old calls 
240 %method ='spline'; 
241 switch 1 
242  case 3,% Best fit by smoothing spline 
243     brks = [0 .1 .2 .4 ,.6, .8, 1, 1.1 1.2]'; 
244     coefa = [0                  0   0.02260415153596   0.99807186986167; ... 
245    2.19065400617385                  0   0.02260415153596  1.00033228501527; ... 
246    4.34015195156053   0.65719620185215   0.03709199393185    1.00478335417504; ... 
247   -1.59533089716870   3.26128737278847   0.80983543882910   1.07321081664798; ... 
248   -6.81273221810880   2.30408883448726   1.92291068028425   1.35286675214799; ... 
249   -3.69498826658975  -1.78355049637802   2.09369407217829   1.77511058383946; ... 
250   13.33514485443956  -4.00054345633187   0.94471547491809   2.09294747228728; ... 
251                   0                  0   0.54466112928490   2.16074873007021]; 
252          
253   coefb = [ 0                  0   0.32503235228616   1.99054481866418; ... 
254    3.28321899128157                  0   0.32503235228616    2.02304805389280; ... 
255    5.67672309005450   0.98496569738447   0.37964649056830    2.05883450811270; ... 
256   -5.29907238080822   4.39099955141717   1.43842344537222   2.21957621884217; .... 
257   -5.89663569823287   1.21155612293224   2.55893458024211   2.64050831092684; ... 
258   -6.21824739906323  -2.32642529600749   2.43691697455115   3.15358438630669; ... 
259   20.19124578481806  -6.05737373544542   0.77134599291113   3.49816479018411; ... 
260                   0                  0   0.16560861936659   3.53491689790559]; 
261 coefc =[                0                  0   0.04818579357214       -0.00817761487085; ... 
262    2.94432030165157                  0   0.04818579357214    -0.00335903551363; ... 
263    4.77660844045250   0.88329609049547   0.09917317190900    0.00440386414523; ... 
264   -1.24578770271258   3.74926115476697   1.01096301945323   0.09778320967047; ... 
265   -7.70868155645400   3.00178853313943   2.36117295703451   0.43997995813009; ... 
266   -3.98346578867600  -1.62342040073298   2.70373824808144   0.97061663841094; ... 
267   13.37833291312857  -4.01349987393858   1.57933014446810   1.41455974568850; ... 
268                   0                  0   1.17798015707424   1.54573609430905]; 
269   pa = mkpp(brks,coefa); 
270   pb = mkpp(brks,coefb); 
271   pc = mkpp(brks,coefc); 
272   A0 = ppval(pa,eps2);         
273   B0 = ppval(pb,eps2);        
274   C0 = ppval(pc,eps2);        
275 case 1, 
276   
277   [E1, H1] = meshgrid(Tpp,Hm00); 
278   A0 = interp2(E1,H1,A00,Tp,Hm0,method); 
279   B0 = interp2(E1,H1,B00,Tp,Hm0,method); 
280   C0 = interp2(E1,H1,C00,Tp,Hm0,method); 
281 end

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