diff --git a/README.md b/README.md
index acda37d6..045d21af 100644
--- a/README.md
+++ b/README.md
@@ -39,6 +39,8 @@ from the Julia REPL.
 
 ## General usage
 
+Note: the current version of `Interpolations` supports interpolation evaluation using index calls `[]`, but this feature will be deprecated in future. We highly recommend function calls with `()` as follows.
+
 Given an `AbstractArray` `A`, construct an "interpolation object" `itp` as
 ```jl
 itp = interpolate(A, options...)
@@ -49,7 +51,7 @@ samples in `A` are equally-spaced.
 
 To evaluate the interpolation at position `(x, y, ...)`, simply do
 ```jl
-v = itp[x, y, ...]
+v = itp(x, y, ...)
 ```
 
 Some interpolation objects support computation of the gradient, which
@@ -89,14 +91,14 @@ Julia's iterator objects, e.g.,
 ```jl
 function ongrid!(dest, itp)
     for I in CartesianRange(size(itp))
-        dest[I] = itp[I]
+        dest[I] = itp(I)
     end
 end
 ```
 would store the on-grid value at each grid point of `itp` in the output `dest`.
 Finally, courtesy of Julia's indexing rules, you can also use
 ```jl
-fine = itp[linspace(1,10,1001), linspace(1,15,201)]
+fine = itp(linspace(1,10,1001), linspace(1,15,201))
 ```
 
 
@@ -114,11 +116,11 @@ Some examples:
 ```jl
 # Nearest-neighbor interpolation
 itp = interpolate(a, BSpline(Constant()), OnCell())
-v = itp[5.4]   # returns a[5]
+v = itp(5.4)   # returns a[5]
 
 # (Multi)linear interpolation
 itp = interpolate(A, BSpline(Linear()), OnGrid())
-v = itp[3.2, 4.1]  # returns 0.9*(0.8*A[3,4]+0.2*A[4,4]) + 0.1*(0.8*A[3,5]+0.2*A[4,5])
+v = itp(3.2, 4.1)  # returns 0.9*(0.8*A[3,4]+0.2*A[4,4]) + 0.1*(0.8*A[3,5]+0.2*A[4,5])
 
 # Quadratic interpolation with reflecting boundary conditions
 # Quadratic is the lowest order that has continuous gradient
@@ -126,7 +128,7 @@ itp = interpolate(A, BSpline(Quadratic(Reflect())), OnCell())
 
 # Linear interpolation in the first dimension, and no interpolation (just lookup) in the second
 itp = interpolate(A, (BSpline(Linear()), NoInterp()), OnGrid())
-v = itp[3.65, 5]  # returns  0.35*A[3,5] + 0.65*A[4,5]
+v = itp(3.65, 5)  # returns  0.35*A[3,5] + 0.65*A[4,5]
 ```
 There are more options available, for example:
 ```jl
@@ -145,8 +147,8 @@ A_x = 1.:2.:40.
 A = [log(x) for x in A_x]
 itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())
 sitp = scale(itp, A_x)
-sitp[3.] # exactly log(3.)
-sitp[3.5] # approximately log(3.5)
+sitp(3.) # exactly log(3.)
+sitp(3.5) # approximately log(3.5)
 ```
 
 For multidimensional uniformly spaced grids
@@ -157,8 +159,8 @@ f(x1, x2) = log(x1+x2)
 A = [f(x1,x2) for x1 in A_x1, x2 in A_x2]
 itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())
 sitp = scale(itp, A_x1, A_x2)
-sitp[5., 10.] # exactly log(5 + 10)
-sitp[5.6, 7.1] # approximately log(5.6 + 7.1)
+sitp(5., 10.) # exactly log(5 + 10)
+sitp(5.6, 7.1) # approximately log(5.6 + 7.1)
 ```
 ### Gridded interpolation
 
@@ -173,7 +175,7 @@ A = rand(20)
 A_x = collect(1.0:2.0:40.0)
 knots = (A_x,)
 itp = interpolate(knots, A, Gridded(Linear()))
-itp[2.0]
+itp(2.0)
 ```
 
 The spacing between adjacent samples need not be constant, you can use the syntax
@@ -188,7 +190,7 @@ For example:
 A = rand(8,20)
 knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
 itp = interpolate(knots, A, Gridded(Linear()))
-itp[4,1.2]  # approximately A[2,6]
+itp(4,1.2)  # approximately A[2,6]
 ```
 One may also mix modes, by specifying a mode vector in the form of an explicit tuple:
 ```jl
@@ -225,7 +227,7 @@ A = hcat(x,y)
 itp = scale(interpolate(A, (BSpline(Cubic(Natural())), NoInterp()), OnGrid()), t, 1:2)
 
 tfine = 0:.01:1
-xs, ys = [itp[t,1] for t in tfine], [itp[t,2] for t in tfine]
+xs, ys = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
 ```
 
 We can then plot the spline with:
diff --git a/src/b-splines/indexing.jl b/src/b-splines/indexing.jl
index d34719b6..20eaae04 100644
--- a/src/b-splines/indexing.jl
+++ b/src/b-splines/indexing.jl
@@ -50,6 +50,11 @@ end
     getindex_impl(itp)
 end
 
+function (itp::BSplineInterpolation{T,N,TCoefs,IT,GT,pad})(args...) where {T,N,TCoefs,IT,GT,pad}
+    # support function calls
+    itp[args...]
+end
+
 function gradient_impl(itp::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}}) where {T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}
     meta = Expr(:meta, :inline)
     # For each component of the gradient, alternately calculate
diff --git a/src/extrapolation/extrapolation.jl b/src/extrapolation/extrapolation.jl
index 88164bae..66fb302a 100644
--- a/src/extrapolation/extrapolation.jl
+++ b/src/extrapolation/extrapolation.jl
@@ -63,6 +63,11 @@ end
     getindex_impl(etp, xs...)
 end
 
+function (etp::Extrapolation{T,N,ITPT,IT,GT,ET})(args...) where {T,N,ITPT,IT,GT,ET}
+    # support function calls
+    etp[args...]
+end
+
 checkbounds(::AbstractExtrapolation,I...) = nothing
 
 
diff --git a/src/extrapolation/filled.jl b/src/extrapolation/filled.jl
index dcf3b026..a0bbc277 100644
--- a/src/extrapolation/filled.jl
+++ b/src/extrapolation/filled.jl
@@ -25,6 +25,11 @@ extrapolate(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) where {T,N,IT,GT}
     convert(Tret, fitp.fillvalue)
 end
 
+function (fitp::FilledExtrapolation{T,N,ITP,IT,GT,FT})(args...) where {T,N,ITP,IT,GT,FT}
+    # support function calls
+    fitp[args...]
+end
+
 @inline Base.checkbounds(::Type{Bool}, A::FilledExtrapolation, I...) = _checkbounds(A, 1, indices(A), I)
 @inline _checkbounds(A, d::Int, IA::TT1, I::TT2) where {TT1,TT2} =
     (I[1] >= lbound(A, d, IA[1])) & (I[1] <= ubound(A, d, IA[1])) & _checkbounds(A, d+1, Base.tail(IA), Base.tail(I))
diff --git a/src/gridded/indexing.jl b/src/gridded/indexing.jl
index cef66e39..184d99db 100644
--- a/src/gridded/indexing.jl
+++ b/src/gridded/indexing.jl
@@ -36,6 +36,11 @@ end
     :(getindex(itp, $(args...)))
 end
 
+function (itp::GriddedInterpolation{T,N,TCoefs,IT,K,pad})(args...) where {T,N,TCoefs,IT,K,pad}
+    # support function calls
+    itp[args...]
+end
+
 # Indexing with vector inputs. Here, it pays to pre-process the input indexes,
 # because N*n is much smaller than n^N.
 # TODO: special-case N=1, because there is no reason to separately cache the indexes.
diff --git a/src/scaling/scaling.jl b/src/scaling/scaling.jl
index d2c30d62..a95e421f 100644
--- a/src/scaling/scaling.jl
+++ b/src/scaling/scaling.jl
@@ -45,6 +45,10 @@ end
 
 getindex(sitp::ScaledInterpolation{T,1}, x::Number, y::Int) where {T} = y == 1 ? sitp[x] : throw(BoundsError())
 
+function (sitp::ScaledInterpolation{T,N,ITPT,IT})(args...) where {T,N,ITPT,IT<:DimSpec}
+    sitp[args...]
+end
+
 size(sitp::ScaledInterpolation, d) = size(sitp.itp, d)
 lbound(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) where {T,N,ITPT,IT} = 1 <= d <= N ? sitp.ranges[d][1] : throw(BoundsError())
 lbound(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) where {T,N,ITPT,IT} = 1 <= d <= N ? sitp.ranges[d][1] - boundstep(sitp.ranges[d]) : throw(BoundsError())
diff --git a/test/b-splines/function-call-syntax.jl b/test/b-splines/function-call-syntax.jl
new file mode 100644
index 00000000..adba1bf2
--- /dev/null
+++ b/test/b-splines/function-call-syntax.jl
@@ -0,0 +1,24 @@
+module ExtrapFunctionCallSyntax
+
+using Base.Test, Interpolations, DualNumbers
+
+# Test if b-spline interpolation by function syntax yields identical results
+f(x) = sin((x-3)*2pi/9 - 1)
+xmax = 10
+A = Float64[f(x) for x in 1:xmax]
+itpg = interpolate(A, BSpline(Linear()), OnGrid())
+schemes = (Flat,Line,Free)
+
+for T in (Cubic, Quadratic), GC in (OnGrid, OnCell)
+    for etp in map(S -> @inferred(interpolate(A, BSpline(T(S())), GC())), schemes),
+        x in linspace(1,xmax,100)
+        @test (getindex(etp, x)) == etp(x)
+    end
+end
+
+for T in (Constant, Linear), GC in (OnGrid, OnCell), x in linspace(1,xmax,100)
+    etp = interpolate(A, BSpline(T()), GC())
+    @test (getindex(etp, x)) == etp(x)
+end
+
+end
diff --git a/test/b-splines/runtests.jl b/test/b-splines/runtests.jl
index 9ea4c231..981a6c9e 100644
--- a/test/b-splines/runtests.jl
+++ b/test/b-splines/runtests.jl
@@ -7,5 +7,6 @@ include("cubic.jl")
 include("mixed.jl")
 include("multivalued.jl")
 include("non1.jl")
+include("function-call-syntax.jl")
 
 end
diff --git a/test/extrapolation/function-call-syntax.jl b/test/extrapolation/function-call-syntax.jl
new file mode 100644
index 00000000..466f2191
--- /dev/null
+++ b/test/extrapolation/function-call-syntax.jl
@@ -0,0 +1,104 @@
+module ExtrapFunctionCallSyntax
+
+using Base.Test, Interpolations, DualNumbers
+
+# Test if extrapolation by function syntax yields identical results
+f(x) = sin((x-3)*2pi/9 - 1)
+xmax = 10
+A = Float64[f(x) for x in 1:xmax]
+itpg = interpolate(A, BSpline(Linear()), OnGrid())
+
+schemes = (
+    Flat,
+    Linear,
+    Reflect,
+    Periodic
+)
+
+for etp in map(E -> @inferred(extrapolate(itpg, E())), schemes),
+    x in [
+        # In-bounds evaluation
+        3.4, 3, dual(3.1),
+        # Out-of-bounds evaluation
+        -3.4, -3, dual(-3,1),
+        13.4, 13, dual(13,1)
+    ]
+    @test (getindex(etp, x)) == etp(x)
+end
+
+etpg = extrapolate(itpg, Flat())
+@test typeof(etpg) <: AbstractExtrapolation
+
+@test etpg(-3) == etpg(-4.5) == etpg(0.9) == etpg(1.0) == A[1]
+@test etpg(10.1) == etpg(11) == etpg(148.298452) == A[end]
+
+etpf = @inferred(extrapolate(itpg, NaN))
+@test typeof(etpf) <: Interpolations.FilledExtrapolation
+@test parent(etpf) === itpg
+
+@test @inferred(size(etpf)) == (xmax,)
+@test isnan(@inferred(etpf(-2.5)))
+@test isnan(etpf(0.999))
+@test @inferred(etpf(1)) ≈ f(1)
+@test etpf(10) ≈ f(10)
+@test isnan(@inferred(etpf(10.001)))
+
+@test etpf(2.5,1) == etpf(2.5)   # for show method
+@test_throws BoundsError etpf(2.5,2)
+@test_throws BoundsError etpf(2.5,2,1)
+
+x =  @inferred(etpf(dual(-2.5,1)))
+@test isa(x, Dual)
+
+etpl = extrapolate(itpg, Linear())
+k_lo = A[2] - A[1]
+x_lo = -3.2
+@test etpl(x_lo) ≈ A[1] + k_lo * (x_lo - 1)
+k_hi = A[end] - A[end-1]
+x_hi = xmax + 5.7
+@test etpl(x_hi) ≈ A[end] + k_hi * (x_hi - xmax)
+
+xmax, ymax = 8,8
+g(x, y) = (x^2 + 3x - 8) * (-2y^2 + y + 1)
+
+itp2g = interpolate(Float64[g(x,y) for x in 1:xmax, y in 1:ymax], (BSpline(Quadratic(Free())), BSpline(Linear())), OnGrid())
+etp2g = extrapolate(itp2g, (Linear(), Flat()))
+
+@test @inferred(etp2g(-0.5,4)) ≈ itp2g(1,4) - 1.5 * epsilon(etp2g(dual(1,1),4))
+@test @inferred(etp2g(5,100)) ≈ itp2g(5,ymax)
+
+etp2ud = extrapolate(itp2g, ((Linear(), Flat()), Flat()))
+@test @inferred(etp2ud(-0.5,4)) ≈ itp2g(1,4) - 1.5 * epsilon(etp2g(dual(1,1),4))
+@test @inferred(etp2ud(5, -4)) == etp2ud(5,1)
+@test @inferred(etp2ud(100, 4)) == etp2ud(8,4)
+@test @inferred(etp2ud(-.5, 100)) == itp2g(1,8) - 1.5 * epsilon(etp2g(dual(1,1),8))
+
+etp2ll = extrapolate(itp2g, Linear())
+@test @inferred(etp2ll(-0.5,100)) ≈ (itp2g(1,8) - 1.5 * epsilon(etp2ll(dual(1,1),8))) + (100 - 8) * epsilon(etp2ll(1,dual(8,1)))
+
+# Allow element types that don't support conversion to Int (#87):
+etp87g = extrapolate(interpolate([1.0im, 2.0im, 3.0im], BSpline(Linear()), OnGrid()), 0.0im)
+@test @inferred(etp87g(1)) == 1.0im
+@test @inferred(etp87g(1.5)) == 1.5im
+@test @inferred(etp87g(0.75)) == 0.0im
+@test @inferred(etp87g(3.25)) == 0.0im
+
+etp87c = extrapolate(interpolate([1.0im, 2.0im, 3.0im], BSpline(Linear()), OnCell()), 0.0im)
+@test @inferred(etp87c(1)) == 1.0im
+@test @inferred(etp87c(1.5)) == 1.5im
+@test @inferred(etp87c(0.75)) == 0.75im
+@test @inferred(etp87c(3.25)) == 3.25im
+@test @inferred(etp87g(0)) == 0.0im
+@test @inferred(etp87g(3.7)) == 0.0im
+
+# Make sure it works with Gridded too
+etp100g = extrapolate(interpolate(([10;20],),[100;110], Gridded(Linear())), Flat())
+@test @inferred(etp100g(5)) == 100
+@test @inferred(etp100g(15)) == 105
+@test @inferred(etp100g(25)) == 110
+# issue #178
+a = randn(10,10) + im*rand(10,10)
+etp = @inferred(extrapolate(interpolate((1:10, 1:10), a, Gridded(Linear())), 0.0))
+@test @inferred(etp(-1,0)) === 0.0+0.0im
+
+end
diff --git a/test/extrapolation/runtests.jl b/test/extrapolation/runtests.jl
index 888a6a50..95608195 100644
--- a/test/extrapolation/runtests.jl
+++ b/test/extrapolation/runtests.jl
@@ -103,3 +103,4 @@ end
 
 include("type-stability.jl")
 include("non1.jl")
+include("function-call-syntax.jl")
diff --git a/test/gridded/function-call-syntax.jl b/test/gridded/function-call-syntax.jl
new file mode 100644
index 00000000..a8b79f7e
--- /dev/null
+++ b/test/gridded/function-call-syntax.jl
@@ -0,0 +1,74 @@
+module GriddedFunctionCallSyntax
+
+using Interpolations, Base.Test
+
+for D in (Constant, Linear), G in (OnCell, OnGrid)
+    ## 1D
+    a = rand(5)
+    knots = (collect(linspace(1,length(a),length(a))),)
+    itp = @inferred(interpolate(knots, a, Gridded(D())))
+    @inferred(getindex(itp, 2))
+    @inferred(getindex(itp, CartesianIndex((2,))))
+    for i = 2:length(a)-1
+        @test itp(i) ≈ a[i]
+        @test itp(CartesianIndex((i,))) ≈ a[i]
+    end
+    @inferred(getindex(itp, knots...))
+    @test itp[knots...] ≈ a
+    # compare scalar indexing and vector indexing
+    x = knots[1]+0.1
+    v = itp(x)
+    for i = 1:length(x)
+        @test v[i] ≈ itp(x[i])
+    end
+    # check the fallback vector indexing
+    x = [2.3,2.2]   # non-increasing order
+    v = itp[x]
+    for i = 1:length(x)
+        @test v[i] ≈ itp(x[i])
+    end
+    # compare against BSpline
+    itpb = @inferred(interpolate(a, BSpline(D()), G()))
+    for x in linspace(1.1,4.9,101)
+        @test itp(x) ≈ itpb(x)
+    end
+
+    ## 2D
+    A = rand(6,5)
+    knots = (collect(linspace(1,size(A,1),size(A,1))),collect(linspace(1,size(A,2),size(A,2))))
+    itp = @inferred(interpolate(knots, A, Gridded(D())))
+    @test parent(itp) === A
+    @inferred(getindex(itp, 2, 2))
+    @inferred(getindex(itp, CartesianIndex((2,2))))
+    for j = 2:size(A,2)-1, i = 2:size(A,1)-1
+        @test itp(i,j) ≈ A[i,j]
+        @test itp(CartesianIndex((i,j))) ≈ A[i,j]
+    end
+    @test itp[knots...] ≈ A
+    @inferred(getindex(itp, knots...))
+    # compare scalar indexing and vector indexing
+    x, y = knots[1]+0.1, knots[2]+0.6
+    v = itp(x,y)
+    for j = 1:length(y), i = 1:length(x)
+        @test v[i,j] ≈ itp(x[i],y[j])
+    end
+    # check the fallback vector indexing
+    x = [2.3,2.2]   # non-increasing order
+    y = [3.5,2.8]
+    v = itp[x,y]
+    for j = 1:length(y), i = 1:length(x)
+        @test v[i,j] ≈ itp(x[i],y[j])
+    end
+    # compare against BSpline
+    itpb = @inferred(interpolate(A, BSpline(D()), G()))
+    for y in linspace(1.1,5.9,101), x in linspace(1.1,4.9,101)
+        @test itp(x,y) ≈ itpb(x,y)
+    end
+
+    A = rand(8,20)
+    knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
+    itp = interpolate(knots, A, Gridded(D()))
+    @test itp(4,1.2) ≈ A[2,6]
+end
+
+end
diff --git a/test/gridded/runtests.jl b/test/gridded/runtests.jl
index e42bbb32..4491d388 100644
--- a/test/gridded/runtests.jl
+++ b/test/gridded/runtests.jl
@@ -2,5 +2,6 @@ module GriddedTests
 
 include("gridded.jl")
 include("mixed.jl")
+include("function-call-syntax.jl")
 
-end
\ No newline at end of file
+end
diff --git a/test/readme-examples.jl b/test/readme-examples.jl
new file mode 100644
index 00000000..336e9f36
--- /dev/null
+++ b/test/readme-examples.jl
@@ -0,0 +1,53 @@
+module ReadmeExampleTests
+# verify examples from README.md run 
+using Interpolations, Base.Test
+
+## Bsplines
+a = randn(5)
+A = randn(5, 5)
+
+# Nearest-neighbor interpolation
+itp = interpolate(a, BSpline(Constant()), OnCell())
+v = itp(5.4)   # returns a[5]
+@test v ≈ a[5]
+
+# (Multi)linear interpolation
+itp = interpolate(A, BSpline(Linear()), OnGrid())
+v = itp(3.2, 4.1)  # returns 0.9*(0.8*A[3,4]+0.2*A[4,4]) + 0.1*(0.8*A[3,5]+0.2*A[4,5])
+@test v ≈ (0.9*(0.8*A[3,4]+0.2*A[4,4]) + 0.1*(0.8*A[3,5]+0.2*A[4,5]))
+
+# Quadratic interpolation with reflecting boundary conditions
+# Quadratic is the lowest order that has continuous gradien
+itp = interpolate(A, BSpline(Quadratic(Reflect())), OnCell())
+
+# Linear interpolation in the first dimension, and no interpolation (just lookup) in the second
+itp = interpolate(A, (BSpline(Linear()), NoInterp()), OnGrid())
+v = itp(3.65, 5)  # returns  0.35*A[3,5] + 0.65*A[4,5]
+@test v ≈ (0.35*A[3,5] + 0.65*A[4,5])
+
+
+## Scaled Bsplines
+A_x = 1.:2.:40.
+A = [log(x) for x in A_x]
+itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())
+sitp = scale(itp, A_x)
+@test sitp(3.) ≈ log(3.) # exactly log(3.)
+@test sitp(3.5) ≈ log(3.5) atol=.1 # approximately log(3.5)
+
+# For multidimensional uniformly spaced grids
+A_x1 = 1:.1:10
+A_x2 = 1:.5:20
+f(x1, x2) = log(x1+x2)
+A = [f(x1,x2) for x1 in A_x1, x2 in A_x2]
+itp = interpolate(A, BSpline(Cubic(Line())), OnGrid())
+sitp = scale(itp, A_x1, A_x2)
+@test sitp(5., 10.) ≈ log(5 + 10) # exactly log(5 + 10)
+@test sitp(5.6, 7.1) ≈ log(5.6 + 7.1) atol=.1 # approximately log(5.6 + 7.1)
+
+## Gridded interpolation
+A = rand(8,20)
+knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
+itp = interpolate(knots, A, Gridded(Linear()))
+@test itp[4,1.2] ≈ A[2,6] atol=.1  # approximately A[2,6]
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 87fc5fb9..0f52b74d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -28,6 +28,7 @@ include("typing.jl")
 include("issues/runtests.jl")
 
 include("io.jl")
+include("readme-examples.jl")
 
 end
 
diff --git a/test/scaling/function-call-syntax.jl b/test/scaling/function-call-syntax.jl
new file mode 100644
index 00000000..cb26658b
--- /dev/null
+++ b/test/scaling/function-call-syntax.jl
@@ -0,0 +1,100 @@
+module ScalingFunctionCallTests
+
+using Interpolations
+using Base.Test
+
+# Model linear interpolation of y = -3 + .5x by interpolating y=x
+# and then scaling to the new x range
+
+itp = interpolate(1:1.0:10, BSpline(Linear()), OnGrid())
+
+sitp = @inferred(scale(itp, -3:.5:1.5))
+@test typeof(sitp) <: Interpolations.ScaledInterpolation
+@test parent(sitp) === itp
+
+for (x,y) in zip(-3:.05:1.5, 1:.1:10)
+    @test sitp(x) ≈ y
+end
+
+# Verify that it works in >1D, with different types of ranges
+
+gauss(phi, mu, sigma) = exp(-(phi-mu)^2 / (2sigma)^2)
+testfunction(x,y) = gauss(x, 0.5, 4) * gauss(y, -.5, 2)
+
+xs = -5:.5:5
+ys = -4:.2:4
+zs = Float64[testfunction(x,y) for x in xs, y in ys]
+
+itp2 = interpolate(zs, BSpline(Quadratic(Flat())), OnGrid())
+sitp2 = @inferred scale(itp2, xs, ys)
+
+for x in xs, y in ys
+    @test testfunction(x,y) ≈ sitp2(x,y)
+end
+
+# Iteration
+itp = interpolate(rand(3,3,3), BSpline(Quadratic(Flat())), OnCell())
+knots = map(d->1:10:21, 1:3)
+sitp = @inferred scale(itp, knots...)
+
+iter = @inferred(eachvalue(sitp))
+state = @inferred(start(iter))
+@test !(@inferred(done(iter, state)))
+val, state = @inferred(next(iter, state))
+
+function foo!(dest, sitp)
+    i = 0
+    for s in eachvalue(sitp)
+        dest[i+=1] = s
+    end
+    dest
+end
+function bar!(dest, sitp)
+    for I in CartesianRange(size(dest))
+        dest[I] = sitp(I)
+    end
+    dest
+end
+rfoo = Array{Float64}( Interpolations.ssize(sitp))
+rbar = similar(rfoo)
+foo!(rfoo, sitp)
+bar!(rbar, sitp)
+@test rfoo ≈ rbar
+
+# with extrapolation
+END = 10
+xs = linspace(-5, 5, END)
+ys = map(sin, xs)
+
+function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT,OnGrid}, itp) where {T,N,IT}
+    for x in xs
+        @test ≈(sut[x],sin(x),atol=sqrt(eps(sin(x))))
+    end
+    @test sut(-5) == sut(-5.1) == sut(-15.8) == sut(-Inf) == itp(1)
+    @test sut(5) == sut(5.1) == sut(15.8) == sut(Inf) == itp(END)
+end
+
+function run_tests(sut::Interpolations.AbstractInterpolation{T,N,IT,OnCell}, itp) where {T,N,IT}
+    halfcell = (xs[2] - xs[1]) / 2
+
+    for x in (5 + halfcell, 5 + 1.1halfcell, 15.8, Inf)
+        @test sut(-x) == itp(.5)
+        @test sut(x) == itp(END+.5)
+    end
+end
+
+for GT in (OnGrid, OnCell)
+    itp = interpolate(ys, BSpline(Quadratic(Flat())), GT())
+
+    # Test extrapolating, then scaling
+    eitp = extrapolate(itp, Flat())
+    seitp = scale(eitp, xs)
+    run_tests(seitp, itp)
+
+    # Test scaling, then extrapolating
+    sitp = scale(itp, xs)
+    esitp = extrapolate(sitp, Flat())
+    run_tests(esitp, itp)
+end
+
+end
diff --git a/test/scaling/runtests.jl b/test/scaling/runtests.jl
index f1176d71..512c294d 100644
--- a/test/scaling/runtests.jl
+++ b/test/scaling/runtests.jl
@@ -1,4 +1,5 @@
 include("scaling.jl")
 include("dimspecs.jl")
 include("nointerp.jl")
-include("withextrap.jl")
\ No newline at end of file
+include("withextrap.jl")
+include("function-call-syntax.jl")