home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / POLYFUN.DI$ / INTERP3.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.1 KB  |  96 lines
  1. function [zi,weights] = interp3(x,y,z,xi,yi,weights)
  2. %INTERP3 2-D biharmonic data interpolation and gridding.
  3. %    ZI = INTERP3(X,Y,Z,XI,YI) uses biharmonic interpolation to
  4. %    find ZI, the values of the underlying 2-D function in Z at the points
  5. %    in matrices XI and YI.  Matrices X and Y specify the points at which 
  6. %    the data Z is given.  X and Y can also be vectors specifying the 
  7. %    abscissae for the matrix Z as for MESHGRID. In both cases, X
  8. %    and Y must be equally spaced and monotonic.
  9. %
  10. %    If X, Y, and Z are all vectors, INTERP3 interpolates within the
  11. %    data given by the (possibly irregularly spaced) ordered triplets
  12. %    (X,Y,Z) to return matrix ZI.
  13. %
  14. %    ZI = INTERP3(Z,M,N) returns an M-by-N matrix ZI whose values
  15. %    are interpolated within matrix Z.
  16. %
  17. %    [ZI,WEIGHTS]=INTERP3(X,Y,Z,XI,YI) also returns the computed
  18. %    interpolation WEIGHTS, which can be reused as follows:
  19. %    INTERP3(X,Y,Z,XI,YI,WEIGHTS) uses the given WEIGHTS vector,
  20. %    instead of recalculating them from X, Y, and Z.
  21. %
  22. %    INTERP4 and INTERP5 are much faster than INTERP3 for
  23. %    data interpolation and should be used unless irregularly
  24. %    spaced vector triplet gridding is required.
  25. %
  26. %    See also INTERP4, INTERP5, INTERP2, SPLINE, INTERP1.
  27.  
  28. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  29. %    Charles R. Denham, Revised by Clay M. Thompson 7-29-91.
  30.  
  31. %     Reference:  David T. Sandwell, Biharmonic spline
  32. %    interpolation of GEOS-3 and SEASAT altimeter
  33. %    data, Geophysical Research Letters, 2, 139-142,
  34. %    1987.  Describes interpolation using value or
  35. %    gradient of value in any dimension.
  36.  
  37. if nargin < 3
  38.   error('Not enough input arguments.');
  39.  
  40. elseif nargin == 3
  41.    rows = y; cols = z;
  42.    z = x;
  43.    [xi,yi] = meshgrid([0:rows-1]./(rows-1),[0:cols-1]./(cols-1));
  44.    [m,n] = size(z);
  45.    [x, y] = meshgrid([0:m-1]./(m-1),[0:n-1]./(n-1));
  46.  
  47. elseif nargin == 4
  48.   error('Not enough or too many input arguments.')
  49.  
  50. end
  51. disp('This usage of interp3 is obsolete and will be eliminated')
  52. disp('in future versions. Please use griddata(x,y,z,xi,yi) instead.')
  53.  
  54. if prod(size(x)) ~= prod(size(z))
  55.    if length(x(:)) .* length(y(:)) == prod(size(z))
  56.       [x, y] = meshgrid(x, y);
  57.      else
  58.       error('Sizes of x, y, and z must match.')
  59.    end
  60. end
  61.  
  62. if size(xi) ~= size(yi)
  63.   error('The size of XI and YI must be the same.');
  64. end
  65.  
  66. xy = x(:) + y(:)*sqrt(-1);
  67.  
  68. % Determine weights for interpolation
  69. if nargin < 6
  70.    d = xy * ones(1,length(xy));
  71.    d = abs(d - d.');
  72.    mask = find(d == 0);
  73.    d(mask) = ones(length(mask),1);
  74.    g = (d.^2) .* (log(d)-1);   % Green's function.
  75.    g(mask) = zeros(length(mask),1); % Value of Green's function at zero 
  76.    weights = g \ z(:);
  77. end
  78.  
  79. [m,n] = size(xi);
  80. zi = zeros(m,n);
  81. jay = sqrt(-1);
  82. xy = xy.';
  83.  
  84. % Evaluate at requested points (xi,yi).  Loop to save memory.
  85. for i=1:m
  86.   for j=1:n
  87.     d = abs(xi(i,j)+jay*yi(i,j) - xy);
  88.     mask = find(d == 0);
  89.     if length(mask)>0, d(mask) = ones(length(mask),1); end
  90.     g = (d.^2) .* (log(d)-1);   % Green's function.
  91.     % Value of Green's function at zero 
  92.     if length(mask)>0, g(mask) = zeros(length(mask),1); end 
  93.     zi(i,j) = g * weights;
  94.   end
  95. end
  96.