; 
; Copyright (c) 1993, Research Systems, Inc. All rights reserved.  
;	Unauthorized reproduction prohibited.  
;
FUNCTION Krig_expon, d, t       ;Return Exponential Covariance Fcn
r = t(2) * exp((-3./t(0)) * d)
z = where(d eq 0.0, count)
if count gt 0 then r(z) = t(1) + t(2)
return, r
end

FUNCTION Krig_sphere, d, t      ;Return Spherical Covariance Fcn
r = d/t(0) < 1.0                ;Normalized distance
r = t(2) * (1. - 1.5 * r + 0.5 * r^3)
z = where(d eq 0.0, count)
if count gt 0 then r(z) = t(1) + t(2)
return, r
end


FUNCTION krig2d, z, x, y, REGULAR = regular, XGRID=xgrid, $
	XVALUES = xvalues, YGRID = ygrid, YVALUES = yvalues, $
	GS = gs, BOUNDS = bounds, NX = nx0, NY = ny0, EXPONENTIAL = ex, $
	SPHERICAL = sp, NESTED = nest, C0 = C0
;+ 
; NAME: 
;	KRIG_2D 
; 
; PURPOSE: 
;	Interpolate a regular or irregularly gridded set of points 
;	z = F(x,y) using kriging.
; CATEGORY: 
;	Interpolation, Surface Fitting 
; CALLING SEQUENCE: 
;	Result = KRIG2D(z [, x, y]) 
; INPUTS: 
;	X, Y, Z = arrays containing the x, y, and z locations of the 
;		data points on the surface.  Need not necessarily be 
;		regularly gridded.  For regularly gridded input data, 
;		X and Y are not used: the grid spacing is specified 
;		via the XGRID or XVALUES, and YGRID or YVALUES, 
;		keywords, and Z must be a two dimensional array.  
;		For irregular grids, all three parameters must be present 
;		and have the same number of elements.  
; KEYWORD PARAMETERS: 
;  Model Parameters:
;	EXPONENTIAL = if set, use an exponential semivariogram model, with 
;		parameters of [ A, C0, C1], which are explained below.
;		If only a two element vector is supplied,
;		C1 is set to the variance of Z.
;	SPHERICAL = if set, use a spherical semivariogram model, with 
;		parameters of [ A, C0, C1], with meanings as follows: 
;		C0 = the "nugget effect", which provides a discontinuity 
;		 at the origin.  
;		(C0 + C1) = the "sill", which is the variogram value for 
;		 very large distances. Also, C1 is the covariance value 
;		 for a zero distance, and the variance of the random 
;		 sample  Z variable.  If only a two element vector is supplied,
;		 C1 is set to the sample variance.
;		A = the range. At distances beyond A, the semivariogram 
;		 or covariance remains essentialy constant.  
;		See the definition of the functions below.  
;  Input grid description:
;	REGULAR = if set, the Z parameter is a two dimensional array, 
;		of dimensions (N,M), containing measurements over a 
;		regular grid.  If any of XGRID, YGRID, XVALUES, YVALUES 
;		are specified, REGULAR is implied. REGULAR is also 
;		implied if there is only one parameter, Z.  If REGULAR is 
;		set, and no grid (_VALUE or _GRID) specifications are present, 
;		the respective grid is set to (0, 1, 2, ...).  
;	XGRID = contains a two element array, [xstart, xspacing], 
;		defining the input grid in the X direction.  Do not specify 
;		both XGRID and XVALUES. 
;	XVALUES = if present, Xvalues(i) contains the X location 
;		of Z(i,j).  Xvalues must be dimensioned with N elements.  
;	YGRID, YVALUES = same meaning as XGRID, XVALUES except for the 
;		Y locations of the input points.  
;  Output grid description:
;	GS = spacing of output grid.  
;		If present, GS a two-element vector 
;		[XS, YS], where XS is the horizontal spacing between 
;	 	grid points and YS is the vertical spacing. The 
;		default is based on the extents of X and Y. If the 
;	 	grid starts at X value Xmin and ends at Xmax, then the 
; 		default horizontal spacing is (Xmax - Xmin)/(NX-1).  
;		YS is computed in the same way. The default grid 
;		size, if neither NX or NY are specified, is 26 by 26.  
;	BOUNDS = a four element array containing the grid limits in X and 
;		Y of the output grid: [ Xmin, Ymin, Xmax, Ymax].  
;		If not specified, the grid limits are set to the extend 
;		of X and Y.  
;	NX = Output grid size in the X direction.  Default = 26, need 
;		not be specified if the size can be inferred by 
;		GS and BOUNDS.  
;	NY = Output grid size in the Y direction.  See NX.  
; 
; OUTPUTS: 
;	A two dimensional floating point array containing the interpolated 
;	surface, sampled at the grid points.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;	None.
; RESTRICTIONS:
;	Limited by the single precision floating point accuracy of the
;	machine.
;		SAMPLE EXECUTION TIMES (measured on a Sun IPX)
;	# of input points	# of output points	Seconds
;	10			676			1.1
;	20			676			1.5
;	40			676			2.6
;	80			676			7.8
;	10			1024			1.6
;	10			4096			5.9
;	10			16384			23
;
; PROCEDURE: 
;	Ordinary kriging is used to fit the surface described by the
;	data points X,Y, and Z.  See: Isaaks and Srivastava,
;	"An Introduction to Applied Geostatistics", Oxford Univ
;	Press, 1989, Chapter 12.
;
;	The parameters of the data model, the range, nugget, and
;	sill, are highly dependent upon the degree and type of spatial
;	variation of your data, and should be determined statistically.
;	Experimentation, or preferrably rigorus analysis, is required.
;
;	For N data points, a system of N+1 simultaneous
;	equations are solved for the coefficients of the 
;	surface.  For any interpolation point, the interpolated value 
;	is: 
;          F(x,y) = Sum( w(i) * C(x(i),y(i), x, y) 
;	Formulas used to model the variogram functions:
;		d(i,j) = distance from point i to point j.
;		V = variance of samples.
;		C(i,j) = Covariance of sample i with sample j.
;               C(x0,y0,x1,y1) = Covariance of point (x0,y0) with (x1,y1).
;       Exponential covar: C(d) = C1 * EXP(-3*d/A)   if d ne 0.
;                               = C1 + C0          if d eq 0.
;       Spherical covar:   C(d) = (1.0 - 1.5 * d/a + 0.5 * (d/a)^3)
;                               = C1 + C0           if d eq 0.
;                               = 0                 if d > a.
; EXAMPLE: IRREGULARLY GRIDDED CASES 
;	Make a random set of points that lie on a gaussian: 
;	n = 15		;# random points
;	x = RANDOMU(seed, n) 
;	y = RANDOMU(seed, n) 
;	z = exp(-2 * ((x-.5)^2 + (y-.5)^2))  ;The gaussian 
; get a 26 by 26 grid over the rectangle bounding x and y: 
;	e = [ 0.25, 0.0]	;Range and nugget are 0.25, and 0.
;			;(These numbers are dependent upon your data model.)
;	r = krig2d(z, x, y, EXPON = e)	;Get the surface.  
; Or: get a surface over the unit square, with spacing of 0.05: 
;	r = krig2d(z, x, y, EXPON=e, GS=[0.05, 0.05], BOUNDS=[0,0,1,1])
; Or: get a 10 by 10 surface over the rectangle bounding x and y: 
;	r = krig2d(z, x, y, EXPON=e, NX=10, NY=10) 
; 
;		REGULARLY GRIDDED CASES 
;	s = [ 10., 0.2]			;Range and sill, data dependent.
;	z = randomu(seed, 5, 6)		;Make some random data
;		interpolate to a 26 x 26 grid: 
;	CONTOUR, krig2d(z, /REGULAR, SPH = s)
; 
; MODIFICATION HISTORY: 
;	DMS, RSI, March, 1993.  Written.  
;-

on_error, 2

s = size(z)		;Assume 2D
nx = s(1)
ny = s(2)

reg = keyword_set(regular) or (n_params() eq 1)

if n_elements(xgrid) eq 2 then begin
	x = findgen(nx) * xgrid(1) + xgrid(0)
	reg = 1
endif else if n_elements(xvalues) gt 0 then begin
	if n_elements(xvalues) ne nx then $
		message,'Xvalues must have '+string(nx)+' elements.'
	x = xvalues
	reg = 1
endif

if n_elements(ygrid) eq 2 then begin
	y = findgen(ny) * ygrid(1) + ygrid(0)
	reg = 1
endif else if n_elements(yvalues) gt 0 then begin
	if n_elements(yvalues) ne ny then $
		message,'Yvalues must have '+string(ny)+' elements.'
	y = yvalues
	reg = 1
endif

if reg then begin
	if s(0) ne 2 then message,'Z array must be 2D for regular grids'
	if n_elements(x) ne nx then x = findgen(nx)
	if n_elements(y) ne ny then y = findgen(ny)
	x = x # replicate(1., ny)	;Expand to full arrays.
	y = replicate(1.,nx) # y
	endif

n = n_elements(x)
if n ne n_elements(y) or n ne n_elements(z) then $
	message,'x, y, and z must have same number of elements.'

if keyword_set(ex) then begin       ;Get model params
    t = ex
    fname = 'KRIG_EXPON'
endif else if keyword_set(sp) then begin
    t = sp
    fname = 'KRIG_SPHERE'
endif else MESSAGE,'Either EXPONENTIAL or SPHERICAL model must be selected.'

if n_elements(t) eq 2 then begin    ;Default value for variance?
    mz = total(z) / n		;Mean of z
    var = total((z - mz)^2)/n	;Variance of Z
    t = [t, var-t(1)]	;Default value for C1
    endif

m = n + 1			;# of eqns to solve
a = fltarr(m, m)

for i=0, n-2 do for j=i,n-1 do begin  ;Only upper diagonal elements
    d = (x(i)-x(j))^2 + (y(i)-y(j))^2  ;Distance squared
    a(i,j) = d & a(j,i) = d             ;symmetric
    endfor

a = call_function(fname, sqrt(a), t)        ;Get coefficient matrix
a(n,*) = 1.0            ;Fill edges
a(*,n) = 1.0
a(n,n) = 0.0

; c = invert(a)           ;Solution using inverse
ludcmp, a, indx, even_odd       ;Solution using LU decomposition

if n_elements(nx0) le 0 then nx0 = 26	;Defaults for nx and ny
if n_elements(ny0) le 0 then ny0 = 26

xmin = min(x, max = xmax)		;Make the grid...
ymin = min(y, max = ymax)

if n_elements(gs) lt 2 then $
	gs = [(xmax-xmin)/(nx0-1.), (ymax-ymin)/(ny0-1.)]
if n_elements(bounds) lt 4 then bounds = [xmin, ymin, xmax, ymax]

nx = ceil((bounds(2)-bounds(0))/gs(0))+1	;# of elements
ny = ceil((bounds(3)-bounds(1))/gs(1))+1

d = fltarr(m)                       ;One extra for lagrange constranint
r = fltarr(nx,ny,/nozero)           ;Result

for j=0,ny-1 do begin               ;Each output point
  y0 = bounds(1) + gs(1) * j
  for i=0,nx-1 do begin
    x0 = bounds(0) + gs(0) * i
    d(0) = sqrt((x-x0)^2 + (y-y0)^2) ;distance
    d = call_function(fname, d, t)          ;Get rhs
    d(n) = 1.0                              ;lagrange constr
    lubksb, a, indx, d
    r(i,j) = total(d * z)
  endfor
endfor
return, r
end