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)
 
  See also  cav76pdf,  lc2tr, dat2tr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

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 %
027 % See also  cav76pdf,  lc2tr, dat2tr
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
041 % - Added example
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