Home > wafo > wavemodels > cav76pdf.m

# cav76pdf

## PURPOSE

Cavanie et al. (1976) approximation of the density (Tc,Ac)

## SYNOPSIS

f = cav76pdf(t,h,mom,g)

## DESCRIPTION

 CAV76PDF Cavanie et al. (1976) approximation of the density  (Tc,Ac)
in a stationary Gaussian transform process X(t) where
Y(t) = g(X(t)) (Y zero-mean Gaussian, X  non-Gaussian).

CALL:  f = cav76pdf(t,h,[m0,m2,m4],g);

f    = density of wave characteristics of half-wavelength
in a stationary Gaussian transformed process X(t),
where Y(t) = g(X(t)) (Y zero-mean Gaussian)
t,h   = vectors of periods and amplitudes, respectively.
default depending on the spectral moments
m0,m2,m4    = the 0'th, 2'nd and 4'th moment of the spectral density
with angular frequency.
g    = space transformation, Y(t)=g(X(t)), default: g is identity
transformation, i.e. X(t) = Y(t)  is Gaussian,
The transformation, g, can be estimated using lc2tr
or dat2tr or given a priori by ochi.
[]    = default values are used.

Example:
mom = spec2mom(jonswap,4);
f   = cav76pdf([],[],mom);

See also  lh83pdf, lc2tr, dat2tr

## CROSS-REFERENCE INFORMATION

This function calls:
 createpdf PDF class constructor gaus2dat Transforms xx using the inverse of transformation g. levels Calculates discrete levels given the parameter matrix. pdfplot Plot contents of pdf structures qlevels Calculates quantile levels which encloses P% of PDF tranproc Transforms process X and up to four derivatives error Display message and abort function. meshgrid X and Y arrays for 3-D plots.
This function is called by:
 Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics

## SOURCE CODE

001 function f = cav76pdf(t,h,mom,g)
002 %CAV76PDF Cavanie et al. (1976) approximation of the density  (Tc,Ac)
003 %         in a stationary Gaussian transform process X(t) where
004 %         Y(t) = g(X(t)) (Y zero-mean Gaussian, X  non-Gaussian).
005 %
006 % CALL:  f = cav76pdf(t,h,[m0,m2,m4],g);
007 %
008 %        f    = density of wave characteristics of half-wavelength
009 %               in a stationary Gaussian transformed process X(t),
010 %               where Y(t) = g(X(t)) (Y zero-mean Gaussian)
011 %       t,h   = vectors of periods and amplitudes, respectively.
012 %               default depending on the spectral moments
013 % m0,m2,m4    = the 0'th, 2'nd and 4'th moment of the spectral density
014 %               with angular frequency.
015 %        g    = space transformation, Y(t)=g(X(t)), default: g is identity
016 %               transformation, i.e. X(t) = Y(t)  is Gaussian,
017 %               The transformation, g, can be estimated using lc2tr
018 %               or dat2tr or given a priori by ochi.
019 %       []    = default values are used.
020 %
021 % Example:
022 %  mom = spec2mom(jonswap,4);
023 %  f   = cav76pdf([],[],mom);
024 %
026
027 % References:
028 % Cavanie, A., Arhan, M. and Ezraty, R. (1976)
029 % "A statistical relationship between individual heights and periods of
030 %  storm waves".
031 % In Proceedings Conference on Behaviour of Offshore Structures,
032 % Trondheim, pp. 354--360
033 % Norwegian Institute of Technology, Trondheim, Norway
034 %
035 % Lindgren, G. and Rychlik, I. (1982)
036 % Wave Characteristics Distributions for Gaussian Waves --
037 % Wave-lenght, Amplitude and Steepness, Ocean Engng vol 9, pp. 411-432.
038
039 % tested on: matlab 5.3 NB! note
040 % History:
041 % revised pab 04.11.2000
042 % - fixed xlabels i.e. f.labx={'Tc','Ac'}
043 % revised by IR 4 X 2000. fixed transform and normalisation
044 % using Lindgren & Rychlik (1982) paper.
045 % At the end of the function there is a text with derivation of the density.
046 %
047 % revised by jr 21.02.2000
048 % - Introduced cell array for f.x for use with pdfplot
049 % by pab 28.09.1999
050
051 if nargin<3|isempty(mom)
052   error('requires the moments')
053 elseif length(mom)<3
054   error('not enough moments')
055 else
056   m0=mom(1);
057   m2=mom(2);
058   m4=mom(3);
059 end
060
061 if nargin<4|isempty(g)
062   g=[(-5:0.02:5)' (-5:0.02:5)'];
063   g(:,1)=sqrt(m0)*g(:,1);
064 end
065
066 if nargin<1|isempty(t)
067     tt1=2*pi*sqrt(m0/m2);
068     paramt=[0 1.7*tt1 51];
069     t=levels(paramt);
070 end
071
072
073 if nargin<2|isempty(h)
074     px=gaus2dat([0. 0.;1 4.],g);
075     px=abs(px(2,2)-px(1,2));
076     paramh=[0 1.3*px 41];
077     h=levels(paramh);
078 end
079
080 eps4 = 1-m2^2/(m0*m4);
081 alfa = m2/sqrt(m0*m4);
082 if ~isreal(sqrt(eps4))
083   error('input moments are not correct')
084 end
085
086
087
088 a   = length(h);
089 b   = length(t);
090 der = ones(a,1);
091
092 h_lh = tranproc([h' der],g);
093 der  = abs(h_lh(:,2));
094 h_lh = h_lh(:,1);
095
096 % Normalization + transformation of t and h
097
098 pos  = 2/(1+alfa);    % inverse of a fraction of positive maxima
099 cons = 2*pi^4*pos/sqrt(2*pi)/m4/sqrt((1-alfa^2));
100                                         %Tm=2*pi*sqrt(m0/m2)/alpha; %mean period between positive maxima
101
102 t_lh = t;
103 h_lh = sqrt(m0)*h_lh;
104
105
106 % Computation of the distribution
107 [T,H] = meshgrid(t_lh(2:b),h_lh);
108 f_th  = zeros(a,b);
109 der   = der(:);
110 %f_th(:,2:b)=const*der(:,ones(1,b-1)).*(H.^2./(T.^5)).*....
111 %    exp(-(H./(eps4*T.^2)).^2/8.*((T.^2-alpha^2).^2+....
112 %    alpha^4*beta^2 ))/(Tm/4*sqrt(m0));
113 f_th(:,2:b)=cons*der(:,ones(1,b-1)).*(H.^2./(T.^5)).*....
114     exp(-0.5*(H./T.^2).^2.*((T.^2-pi^2*m2/m4).^2/(m0*(1-alfa^2))+....
115     pi^4/m4));
116
117 %pdfplot
118 f   = createpdf(2);
119 f.f = f_th;
120
121 f.x{1} = t;
122 f.x{2} = h;
123 f.labx = {'Tc','Ac'};
124 f.title='Joint density of (Tc,Ac) - Cavanie et al. (1976)';
125 %f.eps2=eps2;
126 f.eps4=eps4;
127 f.pl = [10:20:90 95 99 99.9];
128 f.cl = qlevels(f.f,f.pl);
129 pdfplot(f)
130
131 % Let U,Z be the hight and second deivative (curvature) at a local maximum in a Gaussian proces
132 % with spectral moments m0,m2,m4. The conditional density ($U>0$) has the following form
133 %$$134 % f(z,u)=c \frac{1}{\sqrt{2\pi}}\frac{1}{\sqrt{m0(1-\alpha^2)}}\exp(-0.5\left(\frac{u-z(m2/m4)} 135 % {\sqrt{m0(1-\alpha^2)}}\right)^2)\frac{|z|}{m4}\exp(-0.5z^2/m4), \quad z<0, 136 %$$
137 % where $c=2/(1+\alpha)$, $\alpha=m2/\sqrt{m0\cdot m4}$.
138 %
139 % The cavanie approximation is based on the model $X(t)=U \cos(\pi t/T)$, consequently
140 % we have $U=H$ and by twice differentiation $Z=-U(\pi^2/T)^2\cos(0)$. The variable change has Jacobian
141 % $2\pi^2 H/T^3$ giving the final formula for the density of $T,H$
142 %$$143 % f(t,h)=c \frac{2\pi^4}{\sqrt{2\pi}}\frac{1}{m4\sqrt{m0(1-\alpha^2)}}\frac{h^2}{t^5} 144 % \exp(-0.5\frac{h^2}{t^4}\left(\left(\frac{t^2-\pi^2(m2/m4)} 145 % {\sqrt{m0(1-\alpha^2)}}\right)^2+\frac{\pi^4}{m4}\right)). 146 %$$
147 %
148 %
149

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