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)

## CROSS-REFERENCE INFORMATION

This function calls:
 getjonswappeakedness Peakedness factor Gamma given Hm0 and Tp for JONSWAP jhwparfun Wave height, Hd, distribution parameters for Jonswap spectrum. smooth Calculates a smoothing spline. wggamcdf Generalized Gamma cumulative distribution function wggampdf Generalized Gamma probability density function wtweibpdf Truncated Weibull probability density function comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function. interp2 2-D interpolation (table lookup). meshgrid X and Y arrays for 3-D plots. polyval Evaluate polynomial. unique Set unique. warning Display warning message; disable or enable warning messages.
This function is called by:
 jhscdf Joint (Scf,Hd) CDF for linear waves with a JONSWAP spectrum. jhspdf2 Joint (Scf,Hd) PDF for linear waves with a JONSWAP spectrum.

## 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 %
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)
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