# testing the Shepard interpolation method
# Last edited on 2008-07-21 12:04:15 by stolfi

set xrange [-2.1:+4.1]
set yrange [-0.1:+7.1]

# Test data points:
a = 0.3;  b = 0.6;  c = 0.9;  d = 2.5
fa = 1.0; fb = 3.0; fc = 5.0; fd = 2.0

# Linear interpolant (to visualize the data):
L(x) = (x <= a ? fa : (x <= b ? fa + (x-a)/(b-a)*(fb-fa) : (x <= c ? fb + (x-b)/(c-b)*(fc-fb) : (x <= d ? fc + (x-c)/(d-c)*(fd-fc) : fd))))

# Element support radii
# Ra = 0.7; Rb = 0.7; Rc = 0.7; Rd = 0.7
Ra = 1.5; Rb = 1.5; Rc = 1.5; Rd = 1.5

# The Gauss-weighted Shepard mother function with denominator exponent M:
M = 1;
w(r,R) = exp(-(abs(r)/(R/3))**2/2)/(abs(r)/(R/3) + 0.00000000000001)**M

# The Shepard elements:
wsum(x) = w(x - a, Ra) + w(x - b, Rb) + w(x - c, Rc) + w(x - d, Rd)
wa(x) = w(x - a, Ra)/wsum(x)
wb(x) = w(x - b, Rb)/wsum(x)
wc(x) = w(x - c, Rc)/wsum(x)
wd(x) = w(x - d, Rd)/wsum(x)

# The Shepard interpolant:
f(x) = fa*wa(x) + fb*wb(x) + fc*wc(x) + fd*wd(x)

# ----------------------------------------------------------------------
M = 0; set title "Shepard weight functions (exponent = 0)"
plot w(x-a,Ra) with lines, w(x-b,Rb) with lines, w(x-c,Rc) with lines , w(x-d,Rd) with lines 
pause 10

M = 0; set title "Shepard elements (exponent = 0)"
plot wa(x) with lines, wb(x) with lines, wc(x) with lines , wd(x) with lines 
pause 20

M = 0; set title "Shepard interpolant (exponent = 0)"
plot f(x) with lines, L(x) with lines
pause 20

# ----------------------------------------------------------------------
M = 1; set title "Shepard weight functions (exponent = 1)"
plot w(x-a,Ra) with lines, w(x-b,Rb) with lines, w(x-c,Rc) with lines, w(x-d,Rd) with lines 
pause 10

M = 1; set title "Shepard elements (exponent = 1)"
plot wa(x) with lines, wb(x) with lines, wc(x) with lines , wd(x) with lines 
pause 20

M = 1; set title "Shepard interpolant (exponent = 1)"
plot f(x) with lines, L(x) with lines
pause 20

# ----------------------------------------------------------------------
M = 2; set title "Shepard weight functions (exponent = 2)"
plot w(x-a,Ra) with lines, w(x-b,Rb) with lines, w(x-c,Rc) with lines, w(x-d,Rd) with lines  
pause 10

M = 2; set title "Shepard elements (exponent = 2)"
plot wa(x) with lines, wb(x) with lines, wc(x) with lines , wd(x) with lines 
pause 20

M = 2; set title "Shepard interpolant (exponent = 2)"
plot f(x) with lines, L(x) with lines
pause 20

quit