--- a/src/Grids/equidistant_grid.jl	Fri Feb 10 08:36:56 2023 +0100
+++ b/src/Grids/equidistant_grid.jl	Sat Feb 18 11:37:35 2023 +0100
@@ -1,26 +1,77 @@
-"""
-    EquidistantGrid{Dim,T<:Real} <: Grid
+struct EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1,1}
+    points::R
+end
 
-`Dim`-dimensional equidistant grid with coordinates of type `T`.
+Base.eltype(g::EquidistantGrid{T}) where T = T
+Base.getindex(g::EquidistantGrid, i) = g.points[i]
+Base.size(g::EquidistantGrid) = size(g.points)
+Base.length(g::EquidistantGrid) = length(g.points)
+Base.eachindex(g::EquidistantGrid) = eachindex(g.points)
+
+# TODO: Make sure collect works!
+
+
 """
-struct EquidistantGrid{Dim,T<:Real} <: Grid
-    size::NTuple{Dim, Int}
-    limit_lower::NTuple{Dim, T}
-    limit_upper::NTuple{Dim, T}
+    spacing(grid::EquidistantGrid)
+
+The spacing between grid points.
+"""
+spacing(g::EquidistantGrid) = step(g.points)
+
 
-    function EquidistantGrid{Dim,T}(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where {Dim,T}
-        if any(size .<= 0)
-            throw(DomainError("all components of size must be postive"))
-        end
-        if any(limit_upper.-limit_lower .<= 0)
-            throw(DomainError("all side lengths must be postive"))
-        end
-        return new{Dim,T}(size, limit_lower, limit_upper)
-    end
+"""
+    inverse_spacing(grid::EquidistantGrid)
+
+The reciprocal of the spacing between grid points.
+"""
+inverse_spacing(g::EquidistantGrid) = 1/step(g.points)
+
+
+boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
+boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
+boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
+
+
+"""
+    refine(g::EquidistantGrid, r::Int)
+
+Refines `grid` by a factor `r`. The factor is applied to the number of
+intervals which is 1 less than the size of the grid.
+
+See also: [`coarsen`](@ref)
+"""
+function refine(g::EquidistantGrid, r::Int)
+    new_sz = (length(g) - 1)*r + 1
+    return EquidistantGrid(change_length(g.points, new_sz))
 end
 
 """
-    EquidistantGrid(size, limit_lower, limit_upper)
+    coarsen(grid::EquidistantGrid, r::Int)
+
+Coarsens `grid` by a factor `r`. The factor is applied to the number of
+intervals which is 1 less than the size of the grid. If the number of
+intervals are not divisible by `r` an error is raised.
+
+See also: [`refine`](@ref)
+"""
+function coarsen(g::EquidistantGrid, r::Int)
+    if (length(g)-1)%r != 0
+        throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
+    end
+
+    new_sz = (length(g) - 1)÷r + 1
+
+    return EquidistantGrid(change_length(g.points), new_sz)
+end
+
+
+
+
+
+
+
+"""
+    equidistant_grid(size::Dims, limit_lower, limit_upper)
 
 Construct an equidistant grid with corners at the coordinates `limit_lower` and
 `limit_upper`.
@@ -32,168 +83,33 @@
 The number of equidistantly spaced points in each coordinate direction are given
 by the tuple `size`.
 """
-function EquidistantGrid(size, limit_lower, limit_upper)
-    return EquidistantGrid{length(size), eltype(limit_lower)}(size, limit_lower, limit_upper)
-end
-
-"""
-    EquidistantGrid{T}()
-
-Constructs a 0-dimensional grid.
-"""
-EquidistantGrid{T}() where T = EquidistantGrid{0,T}((),(),()) # Convenience constructor for 0-dim grid
-
-
-"""
-    EquidistantGrid(size::Int, limit_lower::T, limit_upper::T)
-
-Convenience constructor for 1D grids.
-"""
-function EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) where T
-	return EquidistantGrid((size,),(limit_lower,),(limit_upper,))
-end
-
-Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T
-
-Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size)
-
-Base.size(g::EquidistantGrid) = g.size
-
-Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim
-
-function Base.getindex(g::EquidistantGrid, I::Vararg{Int})
-    h = spacing(g)
-    return g.limit_lower .+ (I.-1).*h
-end
-
-Base.getindex(g::EquidistantGrid, I::CartesianIndex) = g[Tuple(I)...]
-
-# Review: 
-# Is it not strange that evalOn(::Grid) is non-lazy while evalOn(::EquidistantGrid) is?
-# Also: Change name to evalon or eval_on!!!!!!
-function evalOn(grid::EquidistantGrid, f::Function)
-    F(I...) = f(grid[I...]...)
+function equidistant_grid(size::Dims, limit_lower, limit_upper)
+    gs = map(size, limit_lower, limit_upper) do s,l,u
+        EquidistantGrid(range(l, u, length=s)) # TBD: Should it use LinRange instead?
+    end
 
-    return LazyFunctionArray(F, size(grid))
-end
-
-"""
-    spacing(grid::EquidistantGrid)
-
-The spacing between grid points.
-"""
-spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1)
-
-
-"""
-    inverse_spacing(grid::EquidistantGrid)
-
-The reciprocal of the spacing between grid points.
-"""
-inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid)
-
-
-"""
-    points(grid::EquidistantGrid)
-
-The point of the grid as an array of tuples with the same dimension as the grid.
-The points are stored as [(x1,y1), (x1,y2), … (x1,yn);
-						  (x2,y1), (x2,y2), … (x2,yn);
-						  	⋮		 ⋮            ⋮
-						  (xm,y1), (xm,y2), … (xm,yn)]
-"""
-function points(grid::EquidistantGrid)
-    indices = Tuple.(CartesianIndices(grid.size))
-    h = spacing(grid)
-    return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
-end
-
-"""
-    restrict(::EquidistantGrid, dim)
-
-Pick out given dimensions from the grid and return a grid for them.
-"""
-function restrict(grid::EquidistantGrid, dim)
-    size = grid.size[dim]
-    limit_lower = grid.limit_lower[dim]
-    limit_upper = grid.limit_upper[dim]
-
-    return EquidistantGrid(size, limit_lower, limit_upper)
+    return TensorGrid(gs...)
 end
 
 
 """
-    orthogonal_dims(grid::EquidistantGrid,dim)
+    equidistant_grid(size::Int, limit_lower::T, limit_upper::T)
 
-Returns the dimensions of grid orthogonal to that of dim.
+Constructs a 1D equidistant grid.
 """
-function orthogonal_dims(grid::EquidistantGrid, dim)
-    orth_dims = filter(i -> i != dim, dims(grid))
-	if orth_dims == dims(grid)
-		throw(DomainError(string("dimension ",string(dim)," not matching grid")))
-	end
-    return orth_dims
+function equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T
+	return equidistant_grid((size,),(limit_lower,),(limit_upper,))
 end
 
 
-"""
-    boundary_identifiers(::EquidistantGrid)
-
-Returns a tuple containing the boundary identifiers for the grid, stored as
-	(CartesianBoundary(1,Lower),
-	 CartesianBoundary(1,Upper),
-	 CartesianBoundary(2,Lower),
-	 ...)
-"""
-boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,)
-
 
 """
-    boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
+    change_length(::AbstractRange, n)
 
-Creates the lower-dimensional restriciton of `grid` spanned by the dimensions
-orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional
-grid is a zero-dimensional grid.
-"""
-function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
-    orth_dims = orthogonal_dims(grid, dim(id))
-    return restrict(grid, orth_dims)
-end
-boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
-
-
-"""
-    refine(grid::EquidistantGrid, r::Int)
-
-Refines `grid` by a factor `r`. The factor is applied to the number of
-intervals which is 1 less than the size of the grid.
-
-See also: [`coarsen`](@ref)
+Change the length of a range to `n`, keeping the same start and stop.
 """
-function refine(grid::EquidistantGrid, r::Int)
-    sz = size(grid)
-    new_sz = (sz .- 1).*r .+ 1
-    return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
-end
-
-
-"""
-    coarsen(grid::EquidistantGrid, r::Int)
-
-Coarsens `grid` by a factor `r`. The factor is applied to the number of
-intervals which is 1 less than the size of the grid. If the number of
-intervals are not divisible by `r` an error is raised.
+function change_length(::AbstractRange, n) end
 
-See also: [`refine`](@ref)
-"""
-function coarsen(grid::EquidistantGrid, r::Int)
-    sz = size(grid)
-
-    if !all(n -> (n % r == 0), sz.-1)
-        throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
-    end
-
-    new_sz = (sz .- 1).÷r .+ 1
-
-    return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
-end
+change_length(r::LinRange, n) = LinRange(r[begin], r[end], n)
+change_length(r::StepRangeLen, n) = range(r[begin], r[end], n)
+# TODO: Test the above