Home > wafo > spec > w2k.m

w2k

PURPOSE ^

Translates from frequency to wave number

SYNOPSIS ^

[k,k2,ind]=w2k(w,th,h,g),

DESCRIPTION ^

  W2K Translates from frequency to wave number
      using the dispersion relation
 
  CALL:  [k,k2,ind]=w2k(w,theta,h,g)
 
     k,k2 = wave numbers
      ind = index list of frequencies for which the routine failed (rare)
      w   = angular frequency   
    theta = direction           (default 0)
      h   = water depth         (default Inf)
      g   = constant of gravity (default see gravity)
 
  Uses Newton Raphson method to find the wave number k
  in the dispersion relation w^2= g*k*tanh(k*h).
  The solution k(w) => k = k(w)*cos(theta)
                       k2= k(w)*sin(theta)
  The size of k,k2 is the common size of w and theta if they are matrices
  OR length(theta) x length(w) if they are vectors. If w or theta is scalar
  it functions as a constant matrix of the same size as the other. 
 
  See also  k2w

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [k,k2,ind]=w2k(w,th,h,g),
002 % W2K Translates from frequency to wave number
003 %     using the dispersion relation
004 %
005 % CALL:  [k,k2,ind]=w2k(w,theta,h,g)
006 %
007 %    k,k2 = wave numbers
008 %     ind = index list of frequencies for which the routine failed (rare)
009 %     w   = angular frequency   
010 %   theta = direction           (default 0)
011 %     h   = water depth         (default Inf)
012 %     g   = constant of gravity (default see gravity)
013 %
014 % Uses Newton Raphson method to find the wave number k
015 % in the dispersion relation w^2= g*k*tanh(k*h).
016 % The solution k(w) => k = k(w)*cos(theta)
017 %                      k2= k(w)*sin(theta)
018 % The size of k,k2 is the common size of w and theta if they are matrices
019 % OR length(theta) x length(w) if they are vectors. If w or theta is scalar
020 % it functions as a constant matrix of the same size as the other. 
021 %
022 % See also  k2w
023 
024 % Tested on: Matlab 5.3  
025 % History: 
026 %  revised by es 23.05.00, made 'If w or theta...' in help really true  
027 %  revised by es 21.01.2000  allow length(g)=2 for h=inf (3D normalization)
028 %  by es, pab 13.08.99
029 
030 if nargin<1|isempty(w)
031   k=[];k2=[];
032   return
033 end
034 if nargin<2|isempty(th)
035   th=0;
036 end
037 if nargin<3|isempty(h)
038   h=inf;
039 end
040 if nargin<4|isempty(g)
041   g=gravity;
042 end
043 
044 % Control of dimension of w and theta
045 thtype=0;
046 if prod(size(th))>1 % non-scalar
047   if size(th,1)==1|size(th,2)==1 % theta vector
048     th=th(:);
049     thtype=1;
050   end
051   if size(w,1)==1|size(w,2)==1 % w vector or scalar
052     if thtype==1 % theta also vector
053       th=th*ones(1,length(w));
054 %      w=ones(size(th,1),1)*w(:)';
055     else
056       error('Input dimensions do not match')
057     end
058   else % both matrices
059     if any(size(w)~=size(th))
060       error('Input dimensions do not match')
061     end    
062     thtype=2;
063   end
064 end
065 
066 k=sign(w).*w.^2;  % deep water
067 if h==inf,
068   if thtype==1 % vector
069     k=ones(size(th,1),1)*k(:)';
070   end
071   k2=k.*sin(th)/g(end); %size np x nf
072   k=k.*cos(th)/g(1);
073   return
074 end
075 
076 if length(g)>1
077   error('Finite depth in combination with 3D normalization (=> length(g)=2) is not implemented yet.')
078 end
079 % Newton's Method
080 % Permit no more than count_limit interations.
081 count_limit = 100;
082 count = 0;
083 %disp('Finding the wave numbers')
084 
085 %k = find(w);
086 hn=zeros(size(k));
087 ix=find(w>0 | w<0);
088 
089 % Break out of the iteration loop for three reasons:
090 %  1) the last update is very small (compared to x)
091 %  2) the last update is very small (compared to sqrt(eps))
092 %  3) There are more than 100 iterations. This should NEVER happen. 
093 while(any(ix) & count < count_limit), 
094 
095   count = count + 1;
096   ki=k(ix);
097   hn(ix)  = (ki.*tanh(ki.*h)-w(ix).^2/g)./(tanh(ki.*h)+ki.*h./(cosh(ki.*h).^2));
098   knew = ki - hn(ix);
099   % Make sure that the current guess is not zero.
100   % When Newton's Method suggests steps that lead to zero guesses
101   % take a step 9/10ths of the way to zero:
102   ksmall = find(abs(knew) == 0);
103   if any(ksmall),
104     knew(ksmall) = ki(ksmall) / 10;
105     hn(ix(ksmall)) = ki(ksmall)-knew(ksmall);
106   end
107   
108   k(ix) = knew;
109 %   disp(['Iteration ',num2str(count),'  Number of points left:  ' num2str(length(ix)) ]),
110   
111   ix=find((abs(hn) > sqrt(eps)*abs(k))  &  abs(hn) > sqrt(eps));   
112 end
113 
114 if count == count_limit, 
115   fprintf('\nWarning: w2k did not converge.\n');
116   str = 'The last step was:  ';
117   outstr = sprintf([str,'%13.8f'],hn(ix));
118   fprintf(outstr);
119   ind=ix;
120 else
121   ind=[];
122   %disp('All wave numbers found')
123 end
124 
125 if thtype==1
126   k=ones(size(th,1),1)*k(:)';
127 end
128 k2=k.*sin(th);
129 k =k.*cos(th);
130

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