Home > wafo > spec > private > time2spa.m

time2spa

PURPOSE ^

Transform of spectrum from frequency to wave no. (used in spec2spec)

SYNOPSIS ^

Sk=time2spa(S,k,k2,g,rate)

DESCRIPTION ^

 TIME2SPA Transform of spectrum from frequency to wave no. (used in spec2spec)
 
  CALL:  Sk = time2spa(S,k,k2,g)
 
    Sk = wave number spectrum, 1- or 2D depending on input
    S  = spectrum struct
   k,k2= wave number lag of output (default depending on input spectrum)
    g  = constant of gravity (default see gravity)
  rate = defining the ratio between the wave numbers and angular
         frequencies. (default 2)
 
  Transforms a frequency or directional spectrum to 
  a wave number spectrum.
  Input 1D gives output 1D, input 2D gives output 2D
  Rotation, transformation etc is preserved
 
  See also definitions, spec2spec

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function Sk=time2spa(S,k,k2,g,rate)
002 %TIME2SPA Transform of spectrum from frequency to wave no. (used in spec2spec)
003 %
004 % CALL:  Sk = time2spa(S,k,k2,g)
005 %
006 %   Sk = wave number spectrum, 1- or 2D depending on input
007 %   S  = spectrum struct
008 %  k,k2= wave number lag of output (default depending on input spectrum)
009 %   g  = constant of gravity (default see gravity)
010 % rate = defining the ratio between the wave numbers and angular
011 %        frequencies. (default 2)
012 %
013 % Transforms a frequency or directional spectrum to 
014 % a wave number spectrum.
015 % Input 1D gives output 1D, input 2D gives output 2D
016 % Rotation, transformation etc is preserved
017 %
018 % See also definitions, spec2spec
019 
020 % Tested on Matlab 5.3
021 % History:  revised by så 03.10.2005: corrected bug when input k=dk
022 %           revised by es 29.11.1999: norm. in both x and y (length(g) 1 OR 2)
023 %           by es 99.08.17
024 
025 error(nargchk(1,5,nargin))
026 if nargin<5|isempty(rate), rate = 2; end
027 rate = max(round(abs(rate)),1);
028 
029 if nargin<3
030   if isfield(S,'g')
031     g=S.g;
032   else
033     g=gravity;
034   end
035 end
036 if strcmpi(S.type(end-2:end),'k1d')|strcmpi(S.type(end-2:end),'k2d')
037   error('Spectrum already in space domain')
038 end
039 Sk=S;
040 if isfield(S,'f')
041   w=2*pi*S.f(:);
042   Sf=S.S/2/pi;
043   Sk=rmfield(Sk,'f');
044 else
045   w=S.w(:);
046   Sf=S.S;
047   Sk=rmfield(Sk,'w');
048 end
049 if isfield(Sk,'v')
050   Sk=rmfield(Sk,{'v','phi'});
051   % enc.velocity  and direction make no sense for wave no spectrum
052 end
053 Sk.g=g;
054 
055 if strcmpi(S.type,'freq')|strcmpi(S.type,'enc')
056   if nargin<2|isempty(k)
057     k=linspace(0,w2k(w(end),0,S.h,g(1)),rate*length(w))';
058   else
059     if length(k)==1 % then interpret k as dk (step length)
060       k=((0:length(w)-1)*k)';
061     end
062   end
063   length(k);
064   wk = k2w(k,0,Sk.h,g(1));
065   in = (wk>min(w))&(wk<=w(end));
066   Gk = zeros(size(w));
067   Gk(in)=interp1(w,Sf,wk(in));
068   Sk.S=zeros(size(k'));
069   Sk.S=(Gk.*dwdk(wk,k,Sk.h,g(1)))';
070   Sk.type='k1D';
071   Sk.k=k';
072 else % directional spectrum
073   if S.h<inf
074     disp('This transformation for finite depth is not available yet')
075     return
076   end
077   if nargin<2|(isempty(k)&isempty(k2)) % no arg-in for wave-numbers
078     nw=max(rate-1,1)*length(w);
079     k=linspace(0,w(end)^2/g(1),nw);
080     if g(end)>0
081       k2=linspace(-w(end)^2/g(end), w(end)^2/g(end), 2*nw-1)';
082     else
083        k2=linspace(-w(end)^2/g(1), w(end)^2/g(1), 2*nw-1)';
084     end     
085   else % Some arg-in for k
086     if nargin<3|isempty(k2) % no arg-in for k2
087       k2=k;
088     end
089     if length(k)==1 % k scalar, ie dk
090       k=(0:(nw-1))*k;
091     end
092     if length(k2)==1 % ditto for k2
093       k2=(-(nw-1):(nw-1))'*k2;
094     end
095   end
096   % k,k2 vectors. Make matrices
097   k=k(:)'; % make k row vector
098   if (S.theta(1)<-pi/2)|(S.theta(end)>pi/2);
099     %thetas on the left half plane => k's<0
100     k=[-fliplr(k(:,2:end)) k];
101   end
102   k2=k2(:); % make k2 column vector
103   K1=ones(size(k2))*k;
104   K2=k2*ones(size(k));
105   Sk.k=k;
106   Sk.k2=k2;
107 
108   % Matrices K1 and K2 are equidistant
109   W=sqrt(sqrt(g(1)^2*K1.^2+g(end)^2*K2.^2)); % Non-equidistant W derived from (K1,K2)
110   TH=atan2(g(end)*K2,g(1)*K1); % => -pi < TH <= pi, Non-eq.dist. TH ----"----
111   Gk=zeros(size(W));
112   in=(W>w(1))&(W<=w(end))&(TH>S.theta(1))&(TH<=S.theta(end));
113   Gk(in)=interp2(w,S.theta,Sf,W(in),TH(in),'*linear');
114   Sk.S=zeros(size(Gk));
115   Sk.S(W>0)=Gk(W>0)./abs(W(W>0)).^3/2*g(1)*g(end);  % Wave-number spectrum
116   Sk=rmfield(Sk,'theta');
117   Sk.type='k2D';
118 end
119 return
120 
121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
122 
123 function dw=dwdk(w,k,h,g)
124 
125 dw=zeros(size(w));
126 if h==inf
127   dw(w>0)=g/2./w(w>0);
128 else
129   dw(w>0)=(g*tanh(k(w>0)*h)+g*h*k(w>0).*(1-tanh(k(w>0)*h).^2))./w(w>0)/2;
130 end
131 return
132 
133 
134 
135 
136 
137 
138 
139 
140 
141

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