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
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 = 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
Comments or corrections to the WAFO group