;**************************
;interpolation matrices: [R+(var/q)I]A=S
;where R is the correlation matrix between data points.
;R=e^-(r^2/c^2), where r is distance between points,
;and c is correlation length (both in km).

;S is the correlation matrix between any one data point and particular grid
;point. 
;var is the variance in the error field of the measured values (vector).
;q is the variance in the (known) measured value field (scalar).
;I is identity matrix
;A is the array of weighting coefficients (alpha) =(R^-1)S

;output oi array of values
p=fltarr(3,nx*ny)

;for each grid point
b=0
xg=fltarr(nx)
yg=fltarr(ny)
for i=0,nx-1 do begin
for ii=0,ny-1 do begin

;convert each lat lon point to x-y one at a time so that c,md can remain in 
;kilometers and xg,yg fall on the correct 1degx1deg (varying size) grid.


LLtoXY,la(i),lo(ii),h,v,55.
xg(i)=h
yg(ii)=v

;distance between all data points and each grid point.
d=sqrt((data(0,*)-xg(i))^2 + (data(1,*)-yg(ii))^2)

;consider only nuse points closest to grid point.
;l returns array indices of nuse nearest points only if there are that many
;within c of grid point, otherwise it returns only the number of points
;within c. Nuse is variable.

;h is indices of all points within md of grid point.
h=where(d le md)

if h(0) ge 0 then begin
;u is ordered index of h indices
u=sort(d(h))
;l is h indices in order of increasing d
l=intarr(n_elements(u))
l=h(u)

if n_elements(l) gt n then begin
l=l(0:n-1)
nuse=n
endif else nuse=n_elements(l)
endif else nuse=0

if nuse ne 0 then begin

r=fltarr(nuse,nuse)
s=fltarr(nuse)
alpha=fltarr(nuse)


	for j=0,nuse-1 do begin
	for k=0,nuse-1 do begin

		r(k,j)=exp(-e*( (data(0,l(k))-data(0,l(j)))^2 + $
			        (data(1,l(k))-data(1,l(j)))^2 )/(c^2))
	endfor
	s(j)=exp(-d(l(j))^2/(c^2))
	r(j,j)=r(j,j)+error(2,l(j))/var
	endfor
	
	rinv=invert(r,status)
;just make sure the matrices are not ill-conditioned
	if status ne 0 then begin
	b=b+1
	endif

; for every grid point (i,ii), calculate the weighting coefficients (alpha), 
; and the expectation value (p_hat).
	p_hat=0
	for m=0,nuse-1 do begin
		alpha(m)=total(rinv(*,m)*s)
;weighted sum of expectation deviation from background bg which resides in data(3,*) in 
;order to be at same location as input data.  bg contains the background values at
;output grid points.
		p_hat=p_hat + alpha(m)*(data(2,l(m))-data(3,l(m)))
	endfor
endif else p_hat=0
	
	p(0,i*360.+ii)=la(i)
	p(1,i*360.+ii)=lo(ii)
	p(2,i*360.+ii)=p_hat+bg(2,i*360.+ii)

endfor
endfor

Return to Table of Contents