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
PDF class constructor | |
Transforms xx using the inverse of transformation g. | |
Calculates discrete levels given the parameter matrix. | |
Plot contents of pdf structures | |
Calculates quantile levels which encloses P% of PDF | |
Transforms process X and up to four derivatives | |
Display message and abort function. | |
X and Y arrays for 3-D plots. |
% CHAPTER3 Demonstrates distributions of wave characteristics |
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
Comments or corrections to the WAFO group