Home > wafo > wavemodels > mk87cdf.m

# mk87cdf

## PURPOSE

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

## SYNOPSIS

[p, eps2, a, b] = mk87cdf(Hd,Scf,Hs,Tz,tail)

## DESCRIPTION

``` MK87CDF  Myrhaug and Kjeldsen (1987) joint (Scf,Hd) CDF.

CALL: F = mk87cdf(Hd,Scf,Hs,Tz,tail)

F  = CDF
Hd  = zero down crossing wave height.
Scf = crest front steepness.
Hs  = significant wave height.
Tz  = average zero down crossing period.
tail = 1 if upper tail is calculated
0 if lower tail is calulated (default)

MK87CDF returns the joint CDF 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

Examples:
Hs = 5.5;
Tz = 8.5;
sc = 0.25;
hc = 4;
p = mk87cdf(hc,sc,Hs,Tz) % Prob(Hd<=hc,Scf<=sc|Hs,Tz) = 0.59

% Conditional probability of steep and high waves given seastates
% i.e., Prob(Hd>hc,Scf>sc|Hs,Tz)
upperTail = 1;
Hs = linspace(2.5,18.5,17);
Tz = linspace(4.5,19.5,16);
[T,H] = meshgrid(Tz,Hs);
p = mk87cdf(hc,sc,H,T,upperTail);
v = 10.^(-6:-1);
contourf(Tz,Hs,log10(p),log10(v))
xlabel('Tz')
ylabel('Hs')
fcolorbar(log10(v))

## CROSS-REFERENCE INFORMATION

This function calls:
 gaussq Numerically evaluates a integral using a Gauss quadrature. mk87pdf Myrhaug and Kjeldsen (1987) joint (Scf,Hd) PDF. comnsize Check if all input arguments are either scalar or of common size. error Display message and abort function.
This function is called by:

## SOURCE CODE

```001 function [p, eps2, a, b] = mk87cdf(Hd,Scf,Hs,Tz,tail)
002 %MK87CDF  Myrhaug and Kjeldsen (1987) joint (Scf,Hd) CDF.
003 %
004 % CALL: F = mk87cdf(Hd,Scf,Hs,Tz,tail)
005 %
006 %    F  = CDF
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 %  tail = 1 if upper tail is calculated
012 %         0 if lower tail is calulated (default)
013 %
014 % MK87CDF returns the joint CDF of (Scf, Hd) given Hs and Tz,
015 % i.e., crest front steepness (2*pi*Ac/(g*Td*Tcf)) and wave height, given
016 % the seastate. The root mean square values of Hd and Scf (Hrms,Erms) are
017 % related to the significant waveheight and the average zero down
018 % crossing period by:
019 %             Hrms = 0.715*Hs;
020 %             Erms = 0.0202+0.826*Hs/(Tz^2);
021 %
022 %   The size of F is the common size of  the input arguments
023 %
024 % Examples:
025 %  Hs = 5.5;
026 %  Tz = 8.5;
027 %  sc = 0.25;
028 %  hc = 4;
029 %  p = mk87cdf(hc,sc,Hs,Tz) % Prob(Hd<=hc,Scf<=sc|Hs,Tz) = 0.59
030 %
031 %  % Conditional probability of steep and high waves given seastates
032 %  % i.e., Prob(Hd>hc,Scf>sc|Hs,Tz)
033 %  upperTail = 1;
034 %  Hs = linspace(2.5,18.5,17);
035 %  Tz = linspace(4.5,19.5,16);
036 %  [T,H] = meshgrid(Tz,Hs);
037 %  p = mk87cdf(hc,sc,H,T,upperTail);
038 %  v = 10.^(-6:-1);
039 %  contourf(Tz,Hs,log10(p),log10(v))
040 %  xlabel('Tz')
041 %  ylabel('Hs')
042 %  fcolorbar(log10(v))
043 %
045
046 %   References:
047 %   Myrhaug, D. and Kjelsen S.P. (1987)
048 %  'Prediction of occurences of steep and high waves in deep water'.
049 %   Journal of waterway, Port, Coastal and Ocean Engineers,
050 %   Vol. 113, pp 122--138
051 %
052 %   Myrhaug & Dahle (1984) Parametric modelling of joint probability
053 %   density distributions for steepness and asymmetry in deep water
054 %   waves
055 %
056
057 %tested on matlab 5.2
058 %history:
059 % revised pab 09.08.2003
060 % Changed input + updated help header
061 % revised pab 04.11.2000
062 % - no dependency on stats toolbox anymore.
063 % by  Per A. brodtkorb 1998
064
065
066 error(nargchk(3,5,nargin))
067
068 if (nargin < 5|isempty(tail)),tail = 0; end
069 if (nargin < 4|isempty(Tz)),Tz = 8; end
070 if (nargin < 3|isempty(Hs)), Hs = 6; end
071
072 multipleSeaStates = any(prod(size(Hs))>1|prod(size(Tz))>1);
073 if multipleSeaStates
074 [errorcode, Scf,Hd,Tz,Hs] = comnsize(Scf,Hd,Tz,Hs);
075 else
076   [errorcode, Scf,Hd] = comnsize(Scf,Hd);
077 end
078 if errorcode > 0
079   error('Requires non-scalar arguments to match in size.');
080 end
081
082 Hrms = 0.715*Hs;
083 Erms = 0.0202 + 0.826*Hs./(Tz.^2);
084
085 s = Scf./Erms;
086 hMax = 20;
087 h = min(Hd./Hrms,hMax);
088
089 eps2 = 1e-6;
090
091
092 p = zeros(size(Hd));
093 %k0 = find((Hs<=(Tz-4)*13/6+4));
094 upLimit = 6.5/1.4;
095 loLimit = 2.5/1.26;;
096 k0 = find((loLimit*sqrt(Hs)<Tz));
097 if any(k0)
098   if multipleSeaStates
099     h = h(k0);
100     s = s(k0);
101     Hs = Hs(k0);
102     Tz = Tz(k0);
103   else
104     k0 = 1:prod(size(Hd));
105   end
106   hlim    = h;
107
108   normalizedInput = 1;
109   lowerTail = 0;
110
111
112
113   if 0
114     % This is a trick to get the html documentation correct.
115     k = mk87pdf(1,1,2,3);
116   end
117
118   if (tail == lowerTail)
119     k       = find(h>2.5);
120     hlim(k) = 2.5;
121     p(k0) = gaussq('mk87pdf',0,hlim,eps2/2,[],s,Hs,Tz,5,normalizedInput)...
122     + gaussq('mk87pdf',hlim,h,eps2/2,[],s,Hs,Tz,5,normalizedInput);
123   else
124     k       = find(h<2.5);
125     hlim(k) = 2.5;
126     p(k0) = gaussq('mk87pdf',h,hlim,eps2/2,[],s,Hs,Tz,7,normalizedInput)...
127     + gaussq('mk87pdf',hlim,hMax,eps2/2,[],s,Hs,Tz,7,normalizedInput);
128   end
129 end
130 return
131
132```

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