pro  rbf_eval,w,x,xx,yy

;+
;
;  Radial basis interpolation in nd dimensions
;
;	IN:   
;          w --  1D float array             Array of weights
;          x --  2D float array 	    Coordinates for n known data points
;					                        (nd columns, n rows)
;          xx --  2D float array	    Coordinates for interpolation
;					                        (nd columns, m rows)
;
;	OUT: yy -- float array		      Interpolated values
;
;
;
;-

;sanity checks
if N_params() lt 3 then begin
  print,'% RBF_eval: use -- rbf_eval,w,xx,yy'
  return
endif

sw=size(w)
if sw[0] gt 1 then begin
   print,'% RBF_eval: error -- input w must be a 1D array'
   return
endif

sx=size(x)
if sx[0] ne 2 then begin
  print,'% RBF_eval: error -- input x must be a 2D array'
  return
endif
n=sx[2]

sxx=size(xx)
if sxx[0] gt 2 then begin
   print,'% RBF_eval: error -- input xx must be a 1D or 2D array'
   return
endif
if sxx[0] eq 2 then m=sxx[2]
if sxx[0] eq 1 then m=1

;evaluate interpolations for xx
yy=fltarr(m)
for i=0l,m-1 do begin
  r=fltarr(n)
  for j=0,n-1 do begin 
    r[j]=sqrt(total((xx[*,i]-x[*,j])^2))
  endfor
  f=fltarr(n)
  wg=where(r gt 0.0)
  if max(wg) gt -1 then f[wg]=r[wg]^2*alog10(r[wg])
  yy[i]=total(w*f)
endfor


end