Home > wafo > wavemodels > mk87pdf.m

# mk87pdf

## PURPOSE

Myrhaug and Kjeldsen (1987) joint (Scf,Hd) PDF.

## SYNOPSIS

p = mk87pdf(Hd,Scf,Hs,Tz, condon,normalizedInput)

## DESCRIPTION

``` MK87PDF Myrhaug and Kjeldsen (1987) joint (Scf,Hd) PDF.

CALL:  f = mk87pdf(Hd,Scf,Hs,Tz)

f  = density
Hd  = zero down crossing wave height
Scf = crest front steepness
Hs  = significant wave height.
Tz  = average zero down crossing period.

MK87PDF returns the joint PDF of (Scf, Hd) given Hs and Tz,
i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, given
the seastate. The root mean square values of Hd and Scf (Hrms,Erms) are
related to the significant waveheight and the average zero down
crossing period by:
Hrms = 0.715*Hs;
Erms = 0.0202+0.826*Hs/(Tz^2);

The size of f is the common size of  the input arguments

Example:
Hs = 7;Tz=10;
h  = linspace(0,3*Hs)';
s  = linspace(0,4*Hs/Tz^2)';
[S,H] = meshgrid(s,h);
contour(s,h,mk87pdf(H,S,Hs,Tz))

## CROSS-REFERENCE INFORMATION

This function calls:
 wlogncdf Lognormal cumulative distribution function wlognpdf Lognormal probability density function wweibpdf Weibull probability density function comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function.
This function is called by:
 mk87cdf Myrhaug and Kjeldsen (1987) joint (Scf,Hd) CDF. mk87pdf2 Myrhaug and Kjeldsen (1987) joint (Scf,Hd) PDF.

## SOURCE CODE

```001 function p = mk87pdf(Hd,Scf,Hs,Tz, condon,normalizedInput)
002 %MK87PDF Myrhaug and Kjeldsen (1987) joint (Scf,Hd) PDF.
003 %
004 % CALL:  f = mk87pdf(Hd,Scf,Hs,Tz)
005 %
006 %    f  = density
007 %   Hd  = zero down crossing wave height
008 %   Scf = crest front steepness
009 %   Hs  = significant wave height.
010 %   Tz  = average zero down crossing period.
011 %
012 % MK87PDF returns the joint PDF of (Scf, Hd) given Hs and Tz,
013 % i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, given
014 % the seastate. The root mean square values of Hd and Scf (Hrms,Erms) are
015 % related to the significant waveheight and the average zero down
016 % crossing period by:
017 %             Hrms = 0.715*Hs;
018 %             Erms = 0.0202+0.826*Hs/(Tz^2);
019 %
020 %   The size of f is the common size of  the input arguments
021 %
022 % Example:
023 %  Hs = 7;Tz=10;
024 %  h  = linspace(0,3*Hs)';
025 %  s  = linspace(0,4*Hs/Tz^2)';
026 %  [S,H] = meshgrid(s,h);
027 %  contour(s,h,mk87pdf(H,S,Hs,Tz))
028 %
030
031 %   References:
032 %   Myrhaug, D. and Kjelsen S.P. (1987)
033 %  'Prediction of occurences of steep and high waves in deep water'.
034 %   Journal of waterway, Port, Coastal and Ocean Engineers,
035 %   Vol. 113, pp 122--138
036
037 %
038 %   Myrhaug & Dahle (1984)
039 %   'Parametric modelling of joint probability density
040 %   distributions for steepness and asymmetry in deep water waves'
041 %
042
043 % tested on: matlab 5.1
044 % history:
045 % revised pab 09.08.2003
046 % Changed input + updated help header
047 % revised pab 23.01.2001
048 % - no longer dependent on stats toolbox only wstats
049 % revised pab 10.02.2000
050 % - fixed normalization
051 % by  Per A. Brodtkorb 1998
052
053
054 error(nargchk(3,6,nargin))
055
056 if nargin < 6|isempty(normalizedInput),  normalizedInput=0; end
057 if nargin < 5|isempty(condon),  condon=0; end % regular pdf is returned
058 if nargin < 4|isempty(Tz),  Tz=8;end
059 if nargin < 3|isempty(Hs),  Hs=6;end
060
061 [errorcode, Scf,Hd,Tz,Hs] = comnsize(Scf,Hd,Tz,Hs);
062 if errorcode > 0
063     error('Requires non-scalar arguments to match in size.');
064 end
065 if (normalizedInput>0)
066   Hrms = 1;
067   Erms = 1;
068 else
069   Hrms = 0.715*Hs;
070   Erms = 0.0202 + 0.826*Hs./(Tz.^2);
071 end
072
073 s = Scf./Erms;
074 h = Hd./Hrms;
075
076 sig = (-0.21*atan(2*(h-1.4))+0.325);
077 sig = max(sig,eps); % pab fix: make sure it is positive
078 p   = zeros(size(h));
079
080 % NB! weibpdf must be modified to correspond to
081 % pdf=x^(b-1)/a^b*exp(-(x/a)^b) or else insert
082 % weibpdf=2.39.*h.^1.39/(1.05^2.39).*exp(-(h./1.05).^2.39);
083 switch condon,
084  case 0, % regular pdf is returned
085   p = wweibpdf(h,1.05,2.39).*wlognpdf(s,my(h),sig);
086  case 1, %pdf conditioned on x1 ie. p(x2|x1)
087   p = wlognpdf(s,my(h),sig);
088  case 3, % secret option  used by mk87stat: returns x2*p(x2|x1)
089   p = s.*wlognpdf(s,my(h),sig);
090  case 4, % secret option  used by mk87stat: returns x2.^2*p(x2|x1)
091   p = s.^2.*wlognpdf(s,my(h),sig);
092  case 5, % p(h)*P(V|h) is returned special case used by mk87cdf
093   p = wweibpdf(h,1.05,2.39).*wlogncdf(s,my(h),sig);
094  case 6, % P(V|h) is returned special case used by mk87cdf
095   p = wlogncdf(s,my(h),sig);
096  case 7, % p(h)*(1-P(V|h)) is returned special case used by mk87cdf
097   p = wweibpdf(h,1.05,2.39).*(1-wlogncdf(s,my(h),sig));
098  otherwise error('unknown option')
099 end
100 if condon~=6
101   p = p./Hrms./Erms;
102 end
103 return
104
105 function y = my(h)
106 y   = zeros(size(h));
107 ind = (h <= 1.7);
108 h1  = h(ind);
109 y(ind) = 0.024-1.065.*h1+0.585.*h1.^2;
110
111 %ind=find(h > 1.7);
112 %h2=h(~ind);
113 y(~ind) = 0.32*atan(3.14*(h(~ind)-1.7)) - 0.096;
114 return
115```

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