Home > wafo > wavemodels > jhspdf.m

jhspdf

PURPOSE ^

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

SYNOPSIS ^

[f,Hrms,Vrms] = jhspdf(Hd,Scf,Hm0,Tp,gam,normalizedInput,condon)

DESCRIPTION ^

 JHSPDF Joint (Scf,Hd) PDF for linear waves with JONSWAP spectra. 
  
   CALL: f = jhspdf(Hd,Scf,Hm0,Tp,gam) 
   
     f  = pdf 
    Hd  = zero down crossing wave height 
    Scf = crest front steepness 
    Hm0 = significant wave height 
    Tp  = Spectral peak period  
  
  JHSPDF approximates the joint PDF of (Scf, Hd), i.e., crest front 
  steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process with a 
  JONSWAP spectral density. The empirical parameters of the model is 
  fitted by least squares to simulated (Scf,Hd) data for 13 classes of 
  GAMMA. Between 47000 and 55000 zero-downcrossing waves were 
  simulated for each class of GAMMA. 
  JHSPDF is restricted to the following range for GAMMA:  
   1 <= GAMMA <= 7  
    
  NOTE: The size of f is the common size of the input arguments. 
    
  Example: 
  Hm0 = 6;Tp = 8; 
  h = linspace(0,4*Hm0/sqrt(2));  
  s = linspace(0,6*1.25*Hm0/Tp^2,101); 
  [S,H] = meshgrid(s,h);  
  f = jhspdf(H,S,Hm0,Tp); 
  contourf(s,h,f)   
  
  See also  jhspdf2

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [f,Hrms,Vrms] = jhspdf(Hd,Scf,Hm0,Tp,gam,normalizedInput,condon) 
002 %JHSPDF Joint (Scf,Hd) PDF for linear waves with JONSWAP spectra. 
003 % 
004 %  CALL: f = jhspdf(Hd,Scf,Hm0,Tp,gam) 
005 %  
006 %    f  = pdf 
007 %   Hd  = zero down crossing wave height 
008 %   Scf = crest front steepness 
009 %   Hm0 = significant wave height 
010 %   Tp  = Spectral peak period  
011 % 
012 % JHSPDF approximates the joint PDF of (Scf, Hd), i.e., crest front 
013 % steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, for a Gaussian process with a 
014 % JONSWAP spectral density. The empirical parameters of the model is 
015 % fitted by least squares to simulated (Scf,Hd) data for 13 classes of 
016 % GAMMA. Between 47000 and 55000 zero-downcrossing waves were 
017 % simulated for each class of GAMMA. 
018 % JHSPDF is restricted to the following range for GAMMA:  
019 %  1 <= GAMMA <= 7  
020 %   
021 % NOTE: The size of f is the common size of the input arguments. 
022 %   
023 % Example: 
024 % Hm0 = 6;Tp = 8; 
025 % h = linspace(0,4*Hm0/sqrt(2));  
026 % s = linspace(0,6*1.25*Hm0/Tp^2,101); 
027 % [S,H] = meshgrid(s,h);  
028 % f = jhspdf(H,S,Hm0,Tp); 
029 % contourf(s,h,f)   
030 % 
031 % See also  jhspdf2 
032  
033    
034 % Reference   
035 % P. A. Brodtkorb (2004),   
036 % The Probability of Occurrence of Dangerous Wave Situations at Sea. 
037 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 
038 % Trondheim, Norway.   
039    
040 % History 
041 % By pab 17.01.2003 
042  
043 error(nargchk(3,7,nargin)) 
044  
045  
046 if (nargin < 7|isempty(condon)),  condon  = 0; end 
047 if (nargin < 6|isempty(normalizedInput)),  normalizedInput  = 0;end 
048  
049 if (nargin < 4|isempty(Tp)),  Tp  = 8;end 
050 if (nargin < 3|isempty(Hm0)), Hm0 = 6;end 
051 if (nargin < 5|isempty(gam)) 
052    gam = getjonswappeakedness(Hm0,Tp); 
053 end 
054  
055 multipleSeaStates = any(prod(size(Hm0))>1|... 
056             prod(size(Tp)) >1|... 
057             prod(size(gam))>1); 
058 if multipleSeaStates 
059   [errorcode, Scf,Hd,Hm0,Tp,gam] = comnsize(Scf,Hd,Hm0,Tp,gam); 
060 else 
061   [errorcode, Scf,Hd] = comnsize(Scf,Hd); 
062 end 
063 if errorcode > 0 
064     error('Requires non-scalar arguments to match in size.'); 
065 end 
066 displayWarning = 0; 
067 if displayWarning 
068   if any(any(Tp>5*sqrt(Hm0) | Tp<3.6*sqrt(Hm0))) 
069     disp('Warning: Hm0,Tp is outside the JONSWAP range') 
070     disp('The validity of the parameters returned are questionable') 
071   end 
072 end 
073 k = find(gam<1); 
074 if any(k)  
075   if displayWarning, 
076     warning('Peakedness parameter less than 1. Must be larger than 1.') 
077   end 
078   gam(k)=1; 
079 end 
080 k1 = find(7<gam); 
081 if any(k1) 
082   if displayWarning 
083   warning('Peakedness parameter larger than 7. The pdf returned is questionable') 
084   end 
085   gam(k1) = 7; 
086 end 
087  
088 global JHSPAR 
089 if isempty(JHSPAR) 
090    JHSPAR = load('jhspar13-Jul-2004.mat'); 
091 %  JHSPAR = load('jhspar.mat'); 
092 end 
093 %Generalized Gamma distribution parameters as a function of e2 and h2 
094 A11 = JHSPAR.A11s; 
095 B11 = JHSPAR.B11s; 
096 C11 = 1.0; 
097 %C11 = 1.5; 
098 e2  = JHSPAR.gam(:); % gamma 
099 h2  = JHSPAR.h2(:);  % Hd/Hrms 
100 [E2 H2] = meshgrid(e2,h2); 
101  
102  
103 if normalizedInput, 
104   Hrms = 1; 
105   Vrms = 1; 
106 else 
107   %Tm02 = Tp./(1.30301-0.01698*gam+0.12102./gam); 
108   %dev = 2e-5; 
109   c1 =[ 0.16183666835624   1.53691936441548   1.55852759524555]; 
110   c2 =[ 0.15659478203944   1.15736959972513   1]; 
111   Tm02 = Tp.*(polyval(c2,gam)./polyval(c1,gam)); 
112   Hrms = Hm0/sqrt(2); 
113   Vrms = 1.25*Hm0./(Tm02.^2); % Erms 
114 end 
115  
116 h = Hd./Hrms; 
117 v = Scf./Vrms; 
118 cSize = size(h); % common size 
119  
120 Nh2 = length(h2); 
121 method ='*cubic'; 
122 if multipleSeaStates 
123    h   = h(:); 
124   v   = v(:); 
125   Tp  = Tp(:); 
126   Hm0 = Hm0(:); 
127   gam = gam(:); 
128 %  eps2 = eps2(:); 
129   A1 = zeros(length(h),1); 
130   B1 = A1; 
131   [gamu,ix,jx] = unique(gam); 
132   numSeaStates = length(gamu); 
133   gami = zeros(Nh2,1); 
134    
135   for iz=1:numSeaStates 
136     k = find(jx==iz); 
137     gami(:) = gamu(iz); 
138     A1(k) = exp(smooth(h2,interp2(E2,H2,log(A11),gami,h2,method),1,h(k),1)); 
139     B1(k) = exp(smooth(h2,interp2(E2,H2,log(B11),gami,h2,method),1,h(k),1)); 
140   end 
141 else 
142   A1 = exp(smooth(h2,interp2(E2,H2,log(A11),... 
143                  gam(ones(size(h2))),h2,method),1,h,1)); 
144   B1 = exp(smooth(h2,interp2(E2,H2,log(B11),... 
145                  gam(ones(size(h2))),h2,method),1,h,1)); 
146 end 
147 % Waveheight distribution in time 
148 % Truncated Weibull  distribution parameters as a function of Tp, Hm0, gam  
149 [A0, B0, C0] = jhwparfun(Hm0,Tp,gam,'time'); 
150 switch condon, 
151  case 0, % regular pdf is returned  
152   f = wtweibpdf(h,A0,B0,C0).*wggampdf(v,A1,B1,C11); 
153  case 1, %pdf conditioned on x1 ie. p(x2|x1)  
154   f = wggampdf(v,A1,B1,C11); 
155  case 3, % secret option  used by XXstat: returns x2*p(x2|x1)  
156   f = v.*wggampdf(v,A1,B1,C11); 
157  case 4, % secret option  used by XXstat: returns x2.^2*p(x2|x1)  
158   f = v.^2.*wggampdf(v,A1,B1,C11); 
159  case 5, % p(h)*P(V|h) is returned special case used by thscdf 
160   f = wtweibpdf(h,A0,B0,C0).*wggamcdf(v,A1,B1,C11); 
161  case 6, % P(V|h) is returned special case used by thscdf 
162   f = wggamcdf(v,A1,B1,C11); 
163  case 7,% p(h)*(1-P(V|h)) is returned special case used by thscdf 
164   f = wtweibpdf(h,A0,B0,C0).*(1-wggamcdf(v,A1,B1,C11)); 
165  otherwise error('unknown option') 
166 end 
167 if multipleSeaStates 
168   f = reshape(f,cSize); 
169 end 
170  
171 if condon~=6 
172   f = f./Hrms./Vrms; 
173 end 
174 f(find(isnan(f)|isinf(f) ))=0; 
175 if any(size(f)~=cSize) 
176   disp('Wrong size') 
177 end 
178  
179 return 
180  
181  
182  
183

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