Home > wafo > wavemodels > lh83pdf.m

# lh83pdf

## PURPOSE

Longuet-Higgins (1983) approximation of the density (Tc,Ac)

## SYNOPSIS

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

## DESCRIPTION

``` LH83PDF Longuet-Higgins (1983) 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   = lh83pdf(t,h,[m0,m1,m2],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,m1,m2 = the 0'th,1'st and 2'nd 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 apriori by ochi.
[]    = default values are used.

Example:
even = 0; nr =2;
mom  = spec2mom(jonswap,nr,[],even);
f    = lh83pdf([],[],mom);
pdfplot(f)

## 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 = lh83pdf(t,h,mom,g)
002 %LH83PDF Longuet-Higgins (1983) 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   = lh83pdf(t,h,[m0,m1,m2],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,m1,m2 = the 0'th,1'st and 2'nd 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 apriori by ochi.
019 %       []    = default values are used.
020 %
021 % Example:
022 %  even = 0; nr =2;
023 %  mom  = spec2mom(jonswap,nr,[],even);
024 %  f    = lh83pdf([],[],mom);
025 %  pdfplot(f)
026 %
028
029 % References
030 % Longuet-Higgins,  M.S. (1983)
031 %"On the joint distribution wave periods and amplitudes in a
032 % random wave field", Proc. R. Soc. A389, pp 24--258
033 %
034 % Longuet-Higgins,  M.S. (1975)
035 %"On the joint distribution wave periods and amplitudes of sea waves",
036 % J. geophys. Res. 80, pp 2688--2694
037
038 % tested on: matlab 5.3
039 % History:
040 % Revised pab 01.04.2001
042 % - Better automatic scaling for h,t
043 % revised by IR 18.06.2000, fixing transformation and transposing t and h to fit simpson req.
044 % revised by pab 28.09.1999
045 %   made more efficient calculation of f
046 % by Igor Rychlik
047
048
049 if nargin<3|isempty(mom)
050   error('requires the moments')
051 elseif length(mom)~=3
052   error('not enough moments')
053 else
054   m0=mom(1);
055   m1=mom(2);
056   m2=mom(3);
057 end
058
059 if nargin<4|isempty(g)
060    g=[(-5:0.02:5)' (-5:0.02:5)'];
061    g(:,1)=g(:,1)*sqrt(m0);
062 end
063
064 L0=m0;
065 L1=m1/(2*pi);
066 L2=m2/(2*pi)^2;
067
068
069 if nargin<1|isempty(t)
070     tt1=2*pi*sqrt(m0/m2);
071     paramt=[0 1.7*tt1 51];
072     t=levels(paramt);
073 end
074
075
076 if nargin<2|isempty(h)
077     px=gaus2dat([0. 0.;1 4.],g);
078     px=abs(px(2,2)-px(1,2));
079     paramh=[0 1.3*px 41];
080     h=levels(paramh);
081 end
082
083
084 eps2 = sqrt((L2*L0)/(L1^2)-1);
085
086 if ~isreal(eps2)
087   error('input moments are not correct')
088 end
089 const=4/sqrt(pi)/eps2/(1+1/sqrt(1+eps2^2));
090
091
092 a   = length(h);
093 b   = length(t);
094 der = ones(a,1);
095
096 h_lh = tranproc([h' der],g);
097 der  = abs(h_lh(:,2));
098 h_lh = h_lh(:,1);
099
100 % Normalization + transformation of t and h ???????
101 % Without any transformation
102
103 t_lh  = t/(L0/L1);
104 %h_lh = h_lh/sqrt(2*L0);
105 h_lh  = h_lh/sqrt(2);
106
107 t_lh  = 2*t_lh;
108
109
110 % Computation of the distribution
111 if 0, %old call
112   for i=1:a,
113     for j=2:b,
114       f_th(i,j)=const*der(i)*(h_lh(i)/t_lh(j))^2*....
115         exp(-h_lh(i)^2*(1+((1-1/t_lh(j))/eps2)^2))/((L0/L1)*sqrt(2)/2);
116     end
117   end
118 else %new call
119   [T,H] = meshgrid(t_lh(2:b),h_lh);
120   f_th=zeros(a,b);
121   der=der(:);
122   f_th(:,2:b)=const*der(:,ones(1,b-1)).*(H./T).^2.*....
123       exp(-H.^2.*(1+((1-1./T)/eps2).^2))/((L0/L1)*sqrt(2)/2);
124 end
125
126 f      = createpdf(2);
127 f.f    = f_th;
128 f.x{1} = t(:);
129 f.x{2} = h(:);
130 f.labx = {'Tc','Ac'};
131 f.title = 'Joint density of (Tc,Ac) - Longuet-Higgins (1983)';
132 f.pl   = [10:20:90 95 99];
133 f.cl   = qlevels(f.f,f.pl);
134 pdfplot(f)
135```

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