update docs and add error

Author: mileslucas

README.md

@@ -74,7 +74,13 @@ m[0, 0]      # "index" the model at [y, x]
 m[:, 0]
 m(0.3, 1.0)  # directly query value at (x, y)
 m([1.2, 0.4])
+```
+
+**Important:**:
 
+It is important to recognize the difference in the order of the dimensions between indexing and calling. Indexing is reverse of the cartesian order, which is the natural way of indexing a multi-dimensional array. If you try calling a model as a function with an index, an error will be thrown.
+
+```julia
 # scalar multiplication or division will create a ScaledPSFModel
 20 * m # or `m * 20`
 m / 20

src/PSFModels.jl

@@ -23,8 +23,15 @@ julia> m = PSFModels.Gaussian(5); # fwhm of 5 pixels, centered at (0, 0)
 
 julia> m[0, 0] # [y, x] for indexing
 1.0
+
+julia> m(0, 0) # (x, y) for evaluating
+1.0
 ```
 
+!!! note "axis order"
+
+    It's important to note the difference in the axis ordering between the index-style calls and the function-style calls. The index-style calls are reverse cartesian order (e.g., `(z, y, x)`), while function calls are the typical cartesian order `(x, y, z)`. Regardless, the constructors are always in cartesian order (`(x, y, z)`).
+
 To control the amplitude, the best method is using scalar multiplication or division. These operations create another lazy object ([`ScaledPSFModel`](@ref)) that scales the original model without having to broadcast and potentially allocate.
 
 ```jldoctest model
@@ -55,7 +62,7 @@ if we want to collect the model into a dense matrix, regardless of the indexing
 julia> stamp = collect(m);
 ```
 
-these axes are merely a convenience for bounding the model, since they accept any real number as input. 
+these axes are merely a convenience for bounding the model, since they accept any real number as input.
 
 ```jldoctest model
 julia> m[100, 10000] # index-like inputs [y, x]
@@ -153,11 +160,11 @@ Base.checkbounds(::Type{Bool}, ::PSFModel, idx::CartesianIndex) = true
 # in general, parse to static vector
 (model::PSFModel)(point...) = model(SVector(point))
 (model::PSFModel)(point::Tuple) = model(SVector(point))
-# make sure to reverse indices
-(model::PSFModel)(idx::CartesianIndex) = model(SVector(reverse(idx.I)))
+# disallow index to avoid confusion 
+(model::PSFModel)(::CartesianIndex) = error("PSF models should be indexed using `getindex`, (equivalently `[]`)")
 
 # getindex just calls model with reversed indices
-Base.getindex(model::PSFModel, idx::Vararg{Int,2}) = model(reverse(idx))
+Base.getindex(model::PSFModel, idx::Vararg{<:Number,2}) = model(reverse(idx))
 
 # broadcasting hack to slurp other axes (doesn't work for numbers)
 Broadcast.combine_axes(kern::PSFModel, other) = axes(other)

test/runtests.jl

@@ -11,9 +11,10 @@ function test_model_interface(K)
     @test m.pos ≈ SA[0, 0]
     @test eltype(m) == Float64
 
-    @test m[0, 0] ≈ m(0, 0) ≈ m(SA[0, 0]) ≈ m(CartesianIndex(0, 0)) ≈ 1
+    @test m[0, 0] ≈ m(0, 0) ≈ m(SA[0, 0]) ≈ 1
+    @test_throws ErrorException m(CartesianIndex(0, 0))
     @test m[-100, -10] ≈ m(-10, -100)
-    @test_throws ArgumentError m[1.0, 1.0]
+    @test m(1.0, 1.0) ≈ m[1.0, 1.0]
     @test maximum(m) ≈ 1
     @test 0 ≤ minimum(m) ≤ 1