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: This function is called by:

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 %
025 % See also  lh83pdf, lc2tr, dat2tr
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