diff --git a/docs/lit/examples/02-ellipse.jl b/docs/lit/examples/02-ellipse.jl
index b0abe64..15e93af 100644
--- a/docs/lit/examples/02-ellipse.jl
+++ b/docs/lit/examples/02-ellipse.jl
@@ -114,7 +114,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
 err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
 @assert err < 2e-2
 p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact),
-    "Difference × 10³"; xlabel, ylabel)
+    "|Difference| × 10³"; xlabel, ylabel)
 jim(p1, p4, p2, p3)
 
 
diff --git a/docs/lit/examples/03-rect.jl b/docs/lit/examples/03-rect.jl
index 7ee3d1a..63a495e 100644
--- a/docs/lit/examples/03-rect.jl
+++ b/docs/lit/examples/03-rect.jl
@@ -112,7 +112,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
 err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
 @assert err < 2e-2
 p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact),
-    "Difference × 10³"; xlabel, ylabel)
+    "|Difference| × 10³"; xlabel, ylabel)
 jim(p1, p4, p2, p3)
 
 
@@ -125,7 +125,7 @@ dr = 0.2mm # radial sample spacing
 nr = 2^10 # radial sinogram bins
 r = (-nr÷2:nr÷2-1) * dr # radial samples
 fr = (-nr÷2:nr÷2-1) / nr / dr # corresponding spectral axis
-ϕ = deg2rad.(0:180) # * Unitful.rad # todo round unitful Unitful.°
+ϕ = deg2rad.(0:180)
 sino = radon(ob).(r, ϕ') # sample Radon transform of a single shape object
 smax = ob.value * sqrt(sum(abs2, maximum(width)))
 p5 = jim(r, rad2deg.(ϕ), sino; title="sinogram",
diff --git a/docs/lit/examples/04-gauss.jl b/docs/lit/examples/04-gauss.jl
index 0f70162..6b00975 100644
--- a/docs/lit/examples/04-gauss.jl
+++ b/docs/lit/examples/04-gauss.jl
@@ -114,7 +114,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
 err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
 @assert err < 4e-4
 p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact),
-    "Difference × 10³"; xlabel, ylabel)
+    "|Difference| × 10³"; xlabel, ylabel)
 jim(p1, p4, p2, p3)
 
 
@@ -127,7 +127,7 @@ dr = 0.2mm # radial sample spacing
 nr = 2^10 # radial sinogram bins
 r = (-nr÷2:nr÷2-1) * dr # radial samples
 fr = (-nr÷2:nr÷2-1) / nr / dr # corresponding spectral axis
-ϕ = deg2rad.(0:180) # * Unitful.rad # todo round unitful Unitful.°
+ϕ = deg2rad.(0:180)
 sino = radon(ob).(r, ϕ') # sample Radon transform of a single shape object
 smax = ob.value * IP.fwhm2spread(50mm)
 p5 = jim(r, rad2deg.(ϕ), sino; title="sinogram",
diff --git a/docs/lit/examples/05-triangle.jl b/docs/lit/examples/05-triangle.jl
index 042637f..8832e05 100644
--- a/docs/lit/examples/05-triangle.jl
+++ b/docs/lit/examples/05-triangle.jl
@@ -117,7 +117,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
 err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
 @assert err < 2e-2
 p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact),
-    "Difference × 10³"; xlabel, ylabel)
+    "|Difference| × 10³"; xlabel, ylabel)
 jim(p1, p4, p2, p3)
 
 
@@ -130,7 +130,7 @@ dr = 0.2mm # radial sample spacing
 nr = 2^10 # radial sinogram bins
 r = (-nr÷2:nr÷2-1) * dr # radial samples
 fr = (-nr÷2:nr÷2-1) / nr / dr # corresponding spectral axis
-ϕ = deg2rad.(0:180) # * Unitful.rad # todo round unitful Unitful.°
+ϕ = deg2rad.(0:180)
 sino = radon(ob).(r, ϕ') # sample Radon transform of a single shape object
 smax = ob.value * sqrt((width[1]/2)^2 + (width[2] * sqrt(3) / 2)^2)
 p5 = jim(r, rad2deg.(ϕ), sino; title="sinogram",
diff --git a/docs/lit/examples/32-ellipsoid.jl b/docs/lit/examples/32-ellipsoid.jl
index 89f8360..48795db 100644
--- a/docs/lit/examples/32-ellipsoid.jl
+++ b/docs/lit/examples/32-ellipsoid.jl
@@ -167,7 +167,7 @@ vols  = round.(((p -> sum(p)*prod(pg.deltas)).(proj2)..., volume) ./ 1mm^3; digi
 
 # Look at a set of projections as the views orbit around the object.
 ϕd = 0:6:360
-ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
+ϕs = deg2rad.(ϕd)
 θs = :(π/7)
 θ = eval(θs)
 projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
@@ -221,7 +221,7 @@ p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
 err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
 @assert err < 1e-3
 p9 = jim(axesf(pg), 1e3*abs.(proj_fft - spectrum_slice);
-    title="Difference × 10³", xlabel, ylabel, prompt=false)
+    title="|Difference| × 10³", xlabel, ylabel, prompt=false)
 jim(p6, p7, p8, p9)
 
 
diff --git a/docs/lit/examples/33-cuboid.jl b/docs/lit/examples/33-cuboid.jl
index 881bcb9..95a013e 100644
--- a/docs/lit/examples/33-cuboid.jl
+++ b/docs/lit/examples/33-cuboid.jl
@@ -170,7 +170,7 @@ vols = round.(((p -> sum(p)*prod(pg.deltas)).(proj3)..., volume) ./ 1mm^3; digit
 
 # Look at a set of projections as the views orbit around the object.
 ϕd = 0:6:360
-ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
+ϕs = deg2rad.(ϕd)
 θs = :(π/7)
 θ = eval(θs)
 projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
@@ -224,7 +224,7 @@ p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
 err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
 @assert err < 2e-3
 p9 = jim(axesf(pg), 1e6*abs.(proj_fft - spectrum_slice);
-    title="Difference × 10⁶", xlabel, ylabel, prompt=false)
+    title="|Difference| × 10⁶", xlabel, ylabel, prompt=false)
 jim(p6, p7, p8, p9)
 
 
diff --git a/docs/lit/examples/34-gauss3.jl b/docs/lit/examples/34-gauss3.jl
index 2c75d6d..07dd8fb 100644
--- a/docs/lit/examples/34-gauss3.jl
+++ b/docs/lit/examples/34-gauss3.jl
@@ -128,7 +128,7 @@ p3 = jim(axesf(ig), sp.(spectrum_fft), "log10|DFT|"; clim, xlabel, ylabel)
 err = maximum(abs, spectrum_exact - spectrum_fft) / maximum(abs, spectrum_exact)
 @assert err < 1e-3
 p4 = jim(axesf(ig), 1e3*abs.(spectrum_fft - spectrum_exact);
-   title="Difference × 10³", xlabel, ylabel)
+   title="|Difference| × 10³", xlabel, ylabel)
 jim(p1, p4, p2, p3)
 
 
@@ -167,7 +167,7 @@ vols = round.(((p -> sum(p)*prod(pg.deltas)).(proj2)..., volume) ./ 1mm^3; digit
 
 # Look at a set of projections as the views orbit around the object.
 ϕd = 0:6:360
-ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
+ϕs = deg2rad.(ϕd)
 θs = :(π/7)
 θ = eval(θs)
 projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
@@ -221,7 +221,7 @@ p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
 err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
 @assert err < 1e-5
 p9 = jim(axesf(pg), 1e6*abs.(proj_fft - spectrum_slice);
-    title="Difference × 10⁶", xlabel, ylabel, prompt=false)
+    title="|Difference| × 10⁶", xlabel, ylabel, prompt=false)
 jim(p6, p7, p8, p9)
 
 
diff --git a/docs/lit/examples/35-cylinder.jl b/docs/lit/examples/35-cylinder.jl
index b449125..718941f 100644
--- a/docs/lit/examples/35-cylinder.jl
+++ b/docs/lit/examples/35-cylinder.jl
@@ -170,7 +170,7 @@ vols = round.(((p -> sum(p)*prod(pg.deltas)).(proj3)..., volume) ./ 1mm^3; digit
 
 # Look at a set of projections as the views orbit around the object.
 ϕd = 0:6:360
-ϕs = deg2rad.(ϕd) # * Unitful.rad # todo round unitful Unitful.°
+ϕs = deg2rad.(ϕd)
 θs = :(π/7)
 θ = eval(θs)
 projs = radon(axes(pg)..., ϕs, [θ], [ob]) # many projection views
@@ -224,7 +224,7 @@ p8 = jim(axesf(pg), sp.(proj_fft); prompt=false,
 err = maximum(abs, spectrum_slice - proj_fft) / maximum(abs, spectrum_slice)
 @assert err < 1e-3
 p9 = jim(axesf(pg), 1e3*abs.(proj_fft - spectrum_slice);
-    title="Difference × 10³", xlabel, ylabel, prompt=false)
+    title="|Difference| × 10³", xlabel, ylabel, prompt=false)
 jim(p6, p7, p8, p9)
 
 
diff --git a/src/cuboid.jl b/src/cuboid.jl
index f7f82c0..a3b6c2a 100644
--- a/src/cuboid.jl
+++ b/src/cuboid.jl
@@ -53,7 +53,7 @@ end
 # methods
 
 
-volume(ob::Object3d{Cuboid}) = prod(ob.width)
+volume1(::Cuboid) = 1 # volume of unit cube
 
 
 """
diff --git a/src/cylinder.jl b/src/cylinder.jl
index 6ea1a39..2aaed1f 100644
--- a/src/cylinder.jl
+++ b/src/cylinder.jl
@@ -29,7 +29,7 @@ cylinder(args... ; kwargs...) = Object(Cylinder(), args...; kwargs...)
 # methods
 
 
-volume(ob::Object3d{Cylinder}) = π * prod(ob.width[1:2]) * ob.width[3]
+volume1(::Cylinder) = π
 
 
 """
diff --git a/src/ellipse.jl b/src/ellipse.jl
index f653170..cba1650 100644
--- a/src/ellipse.jl
+++ b/src/ellipse.jl
@@ -51,6 +51,9 @@ end
 # methods
 
 
+area1(::Ellipse) = π # area of unit circle
+
+
 """
     phantom1(ob::Object2d{Ellipse}, (x,y))
 Evaluate unit circle at `(x,y)`,
diff --git a/src/ellipsoid.jl b/src/ellipsoid.jl
index 43e5e85..1b02fc6 100644
--- a/src/ellipsoid.jl
+++ b/src/ellipsoid.jl
@@ -51,7 +51,7 @@ end
 # methods
 
 
-volume(ob::Object3d{Ellipsoid}) = 4/3 * π * prod(ob.width)
+volume1(::Ellipsoid) = 4/3 * π # volume of unit sphere
 
 
 """
diff --git a/src/gauss2.jl b/src/gauss2.jl
index 2a56e04..6233730 100644
--- a/src/gauss2.jl
+++ b/src/gauss2.jl
@@ -53,6 +53,9 @@ Convert FWHM `w` to equivalent Gaussian spread `s` for ``\\exp(-π (x/s)^2)``.
 # methods
 
 
+area1(::Gauss2) = fwhm2spread(1)^2
+
+
 """
     phantom1(ob::Object2d{Gauss2}, (x,y))
 Evaluate unit gauss2 at `(x,y)`,
diff --git a/src/gauss3.jl b/src/gauss3.jl
index f8219b4..47c28b2 100644
--- a/src/gauss3.jl
+++ b/src/gauss3.jl
@@ -32,7 +32,7 @@ gauss3(args... ; kwargs...) = Object(Gauss3(), args...; kwargs...)
 # methods
 
 
-volume(ob::Object3d{Gauss3}) = fwhm2spread(1)^3 * prod(ob.width)
+volume1(::Gauss3) = fwhm2spread(1)^3
 
 
 """
diff --git a/src/object2.jl b/src/object2.jl
index ea19cf3..b59250c 100644
--- a/src/object2.jl
+++ b/src/object2.jl
@@ -7,6 +7,9 @@ Utilities for 2D objects
 export phantom, radon, spectrum
 
 
+area(ob::Object2d{S}) where S = area1(S()) * prod(ob.width)
+
+
 # rotate
 
 """
diff --git a/src/object3.jl b/src/object3.jl
index 61f8e12..00f467e 100644
--- a/src/object3.jl
+++ b/src/object3.jl
@@ -7,6 +7,9 @@ using LazyGrids: ndgrid
 export phantom, radon, spectrum
 
 
+volume(ob::Object3d{S}) where S = volume1(S()) * prod(ob.width)
+
+
 # rotate
 
 """
diff --git a/src/rect.jl b/src/rect.jl
index 930afbb..2535795 100644
--- a/src/rect.jl
+++ b/src/rect.jl
@@ -73,6 +73,9 @@ end
 # methods
 
 
+area1(::Rect) = 1 # area of unit square
+
+
 """
     phantom1(ob::Object2d{Rect}, (x,y))
 Evaluate unit square at `(x,y)`,
diff --git a/src/triangle.jl b/src/triangle.jl
index 5ee4904..ce09384 100644
--- a/src/triangle.jl
+++ b/src/triangle.jl
@@ -125,6 +125,9 @@ end
 # methods
 
 
+area1(::Triangle) = sqrt(3)/4 # area of unit base equilateral
+
+
 """
     phantom1(ob::Object2d{Triangle}, (x,y))
 Evaluate unit triangle at `(x,y)`,
diff --git a/test/cuboid.jl b/test/cuboid.jl
index d48f6d3..3d51100 100644
--- a/test/cuboid.jl
+++ b/test/cuboid.jl
@@ -158,7 +158,8 @@ end
     v = (-nv÷2:nv÷2-1) * dv
     fu = (-nu÷2:nu÷2-1) / nu / du
     fv = (-nv÷2:nv÷2-1) / nv / dv
-    ϕ = deg2rad.(0:6:180) # * Unitful.rad # todo round unitful Unitful.°
+#   ϕ = (0:6:180) * deg2rad(1)
+    ϕ = deg2rad.(0:6:180) # helps coverage of "abs(e2) < eps(T)"
     θ = [π/7]
     sino = @inferred radon(u, v, ϕ, θ, [ob])
 
diff --git a/test/cylinder.jl b/test/cylinder.jl
index b3b0425..e23bbba 100644
--- a/test/cylinder.jl
+++ b/test/cylinder.jl
@@ -141,12 +141,9 @@ end
     v = (-nv÷2:nv÷2-1) * dv
     fu = (-nu÷2:nu÷2-1) / nu / du
     fv = (-nv÷2:nv÷2-1) / nv / dv
-    ϕ = deg2rad.(0:6:180) # * Unitful.rad # todo round unitful Unitful.°
+    ϕ = (0:30:180) * deg2rad(1)
     θ = [π/7]
     sino = @inferred radon(u, v, ϕ, θ, [ob])
 
-#=
-todo projection slice
-=#
-
+# todo projection slice
 end
diff --git a/test/ellipse.jl b/test/ellipse.jl
index 898e832..48ee0f6 100644
--- a/test/ellipse.jl
+++ b/test/ellipse.jl
@@ -99,7 +99,7 @@ end
     for ob in obs
         nr, dr = 2^4, 0.02m
         r = (-nr÷2:nr÷2-1) * dr .+ ob.center[1]
-        ϕ = deg2rad.(0:20:360)
+        ϕ = (0:30:360) * deg2rad(1)
         @inferred IP._radon(ob, r[1], ϕ[1])
         sino1 = @inferred radon([r[1]], [ϕ[1]], [ob])
         sino = @inferred radon(r, ϕ, [ob])
@@ -118,7 +118,7 @@ end
     ob = shape((4m, 3m), radii, π/6, 1.0f0)
     img = @inferred phantom(x, y, [ob])
 
-    zscale = 1 / π / prod(radii) # normalize spectra by area
+    zscale = 1 / (ob.value * IP.area(ob)) # normalize spectra by area
     fx = (-M÷2:M÷2-1) / M / dx
     fy = (-N÷2:N÷2-1) / N / dy
     X = myfft(img) * dx * dy * zscale
@@ -145,8 +145,7 @@ end
     nr = 2^10
     r = (-nr÷2:nr÷2-1) * dr
     fr = (-nr÷2:nr÷2-1) / nr / dr
-    ϕ = deg2rad.(0:360) # * Unitful.rad # todo round unitful Unitful.°
-#   ϕ = deg2rad.((0:360)°) # not yet due to Unitful issue
+    ϕ = (0:30:360) * deg2rad(1)
     sino = @inferred radon(r, ϕ, [ob])
 
     ia = argmin(abs.(ϕ .- deg2rad(55)))
diff --git a/test/ellipsoid.jl b/test/ellipsoid.jl
index c78f326..31ee528 100644
--- a/test/ellipsoid.jl
+++ b/test/ellipsoid.jl
@@ -144,7 +144,7 @@ end
     v = (-nv÷2:nv÷2-1) * dv
     fu = (-nu÷2:nu÷2-1) / nu / du
     fv = (-nv÷2:nv÷2-1) / nv / dv
-    ϕ = deg2rad.(0:6:180) # * Unitful.rad # todo round unitful Unitful.°
+    ϕ = (0:30:180) * deg2rad(1)
     θ = [π/7]
     sino = @inferred radon(u, v, ϕ, θ, [ob])
 
diff --git a/test/gauss2.jl b/test/gauss2.jl
index 6aeb570..4fafa2f 100644
--- a/test/gauss2.jl
+++ b/test/gauss2.jl
@@ -108,7 +108,7 @@ end
     ob = shape((2m, 3m), width, π/6, 1.0f0)
     img = @inferred phantom(x, y, [ob])
 
-    zscale = 1 / IP.fwhm2spread(1)^2 / prod(width) # normalize spectra by area
+    zscale = 1 / (ob.value * IP.area(ob)) # normalize spectra by area
     fx = (-M÷2:M÷2-1) / M / dx
     fy = (-N÷2:N÷2-1) / N / dy
     X = myfft(img) * dx * dy * zscale
@@ -135,8 +135,7 @@ end
     nr = 2^10
     r = (-nr÷2:nr÷2-1) * dr
     fr = (-nr÷2:nr÷2-1) / nr / dr
-    ϕ = deg2rad.(0:360) # * Unitful.rad # todo round unitful Unitful.°
-#   ϕ = deg2rad.((0:180)°) # not yet due to Unitful issue
+    ϕ = (0:30:360) * deg2rad(1)
     sino = @NOTinferred radon(r, ϕ, [ob])
 
     ia = argmin(abs.(ϕ .- deg2rad(55)))
diff --git a/test/gauss3.jl b/test/gauss3.jl
index f88b46f..4266486 100644
--- a/test/gauss3.jl
+++ b/test/gauss3.jl
@@ -138,11 +138,11 @@ end
     v = (-nv÷2:nv÷2-1) * dv
     fu = (-nu÷2:nu÷2-1) / nu / du
     fv = (-nv÷2:nv÷2-1) / nv / dv
-    ϕ = deg2rad.(0:6:180) # * Unitful.rad # todo round unitful Unitful.°
+    ϕ = (0:30:180) * deg2rad(1)
     θ = [π/7]
     sino = @inferred radon(u, v, ϕ, θ, [ob])
-end
 
 #=
 todo projection slice
 =#
+end
diff --git a/test/rect.jl b/test/rect.jl
index fdb0696..64b80d0 100644
--- a/test/rect.jl
+++ b/test/rect.jl
@@ -98,7 +98,7 @@ end
     ob = shape((4m, 3m), width, π/6, 1.0f0)
     img = @inferred phantom(x, y, [ob])
 
-    zscale = 1 / prod(width) # normalize spectra by area
+    zscale = 1 / (ob.value * IP.area(ob)) # normalize spectra by area
     fx = (-M÷2:M÷2-1) / M / dx
     fy = (-N÷2:N÷2-1) / N / dy
     X = myfft(img) * dx * dy * zscale
@@ -125,8 +125,7 @@ end
     nr = 2^10
     r = (-nr÷2:nr÷2-1) * dr
     fr = (-nr÷2:nr÷2-1) / nr / dr
-    ϕ = deg2rad.(0:360) # * Unitful.rad # todo round unitful Unitful.°
-#   ϕ = deg2rad.((0:360)°) # not yet due to Unitful issue
+    ϕ = (0:30:360) * deg2rad(1)
     sino = @inferred radon(r, ϕ, [ob])
 
     ia = argmin(abs.(ϕ .- deg2rad(55)))
diff --git a/test/shepplogan.jl b/test/shepplogan.jl
index ce87811..b2a825c 100644
--- a/test/shepplogan.jl
+++ b/test/shepplogan.jl
@@ -2,7 +2,6 @@
 shepplogan.jl
 =#
 
-#using MIRTjim: jim, prompt
 using ImagePhantoms # many
 using Test: @test, @testset, @test_throws, @inferred
 
@@ -14,38 +13,33 @@ using Test: @test, @testset, @test_throws, @inferred
     end
     @test ellipse_parameters(SouthPark()) isa Matrix
 
-    image0 = @NOTinferred shepp_logan(256, SheppLoganEmis())
+    image0 = @NOTinferred shepp_logan(2^6, SheppLoganEmis())
     @test image0 isa Matrix
-    image1 = @NOTinferred shepp_logan(256, SheppLoganEmis(); oversample=1)
-    image2 = @NOTinferred shepp_logan(256, SheppLoganEmis(); oversample=3)
+    image1 = @NOTinferred shepp_logan(2^6, SheppLoganEmis(); oversample=1)
+    image2 = @NOTinferred shepp_logan(2^6, SheppLoganEmis(); oversample=3)
     @test image0 == image2
-#   jim(jim(image1), jim(image2))
 
     ob = shepp_logan(SheppLoganEmis())
-    x = LinRange(-1,1,201) * 0.5
-    y = LinRange(-1,1,200) * 0.5
+    x = LinRange(-1,1,2^6+1) * 0.5
+    y = LinRange(-1,1,2^6) * 0.5
     image = phantom(x, y, ob)
     @test image isa Matrix
     image = phantom(ob).(x,y')
     @test image isa Matrix
-#   jim(x, y, image)
 
-    r = LinRange(-0.5,0.5,101)
-    ϕ = deg2rad.(0:180)
+    r = LinRange(-0.5,0.5,2^5)
+    ϕ = (0:30:180) * deg2rad(1)
     sino = radon(ob).(r,ϕ')
     @test sino isa Matrix
     sino = radon(r, ϕ, ob)
     @test sino isa Matrix
-#   jim(r, ϕ, sino; aspect_ratio=:none, yflip=false)
 
-    kx = LinRange(-1,1,100) * 9
-    ky = LinRange(-1,1,101) * 9
+    kx = LinRange(-1,1,2^6) * 9
+    ky = LinRange(-1,1,2^6+1) * 9
     kspace = spectrum(ob).(kx, ky')
     @test kspace isa Matrix
     kspace = spectrum(kx, ky, ob)
     @test kspace isa Matrix
-#   jim(kx, ky, kspace)
 
-    image = shepp_logan(80, 100, SouthPark(), fovs=(1,1))
-#   jim(image)
+    image = shepp_logan(40, 50, SouthPark(), fovs=(1,1))
 end
diff --git a/test/triangle.jl b/test/triangle.jl
index f02f691..490b55c 100644
--- a/test/triangle.jl
+++ b/test/triangle.jl
@@ -100,7 +100,7 @@ end
     ob = shape((-3m, -7m), width, π/6, 1.0f0)
     img = @inferred phantom(x, y, [ob])
 
-    zscale = 1 / (sqrt(3)/4 * prod(width)) # normalize spectra by area
+    zscale = 1 / (ob.value * IP.area(ob)) # normalize spectra by area
     fx = (-M÷2:M÷2-1) / M / dx
     fy = (-N÷2:N÷2-1) / N / dy
     X = myfft(img) * dx * dy * zscale
@@ -127,8 +127,7 @@ end
     nr = 2^12
     r = (-nr÷2:nr÷2-1) * dr
     fr = (-nr÷2:nr÷2-1) / nr / dr
-    ϕ = deg2rad.(0:360) # * Unitful.rad # todo round unitful Unitful.°
-#   ϕ = deg2rad.((0:360)°) # not yet due to Unitful issue
+    ϕ = (0:30:360) * deg2rad(1)
     sino = @NOTinferred radon(r, ϕ, [ob])
 
     ia = argmin(abs.(ϕ .- deg2rad(35)))
@@ -155,5 +154,5 @@ if DEBUG
     plot(p1, p2, p3, p4); gui()
 end
 
-    @test maximum(abs, ideal - Slice) / maximum(abs, ideal) < 4e-6
+    @test maximum(abs, ideal - Slice) / maximum(abs, ideal) < 5e-6
 end