diff --git a/Manifest.toml b/Manifest.toml
index 7961b15..69b999f 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -95,9 +95,9 @@ version = "3.10.2"
 
 [[ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
-git-tree-sha1 = "4bffea7ed1a9f0f3d1a131bbcd4b925548d75288"
+git-tree-sha1 = "5e9769a17f17b587c951d57ba4319782b40c3513"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.10.9"
+version = "0.10.10"
 
 [[ColorVectorSpace]]
 deps = ["ColorTypes", "Colors", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "StatsBase"]
@@ -203,10 +203,10 @@ uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
 version = "1.1.0"
 
 [[FileIO]]
-deps = ["Pkg"]
-git-tree-sha1 = "fee8955b9dfa7bec67117ef48085fb2b559b9c22"
+deps = ["Pkg", "Requires", "UUIDs"]
+git-tree-sha1 = "8800ec70aee7292931a3d3c10a3be3445b9c6141"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.4.5"
+version = "1.6.2"
 
 [[FixedPointNumbers]]
 deps = ["Statistics"]
@@ -300,9 +300,9 @@ version = "0.8.20"
 
 [[ImageMagick]]
 deps = ["FileIO", "ImageCore", "ImageMagick_jll", "InteractiveUtils", "Libdl", "Pkg", "Random"]
-git-tree-sha1 = "02558f83932fde6ebd3ab007dbff6bd8740a8247"
+git-tree-sha1 = "f35fb0080591e6f271b7e1cf4c8323a4a78d5150"
 uuid = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
-version = "1.1.6"
+version = "1.1.7"
 
 [[ImageMagick_jll]]
 deps = ["JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg", "Zlib_jll", "libpng_jll"]
@@ -380,15 +380,15 @@ uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
 version = "2.10.0+3"
 
 [[LaTeXStrings]]
-git-tree-sha1 = "c7aebfecb1a60d59c0fe023a68ec947a208b1e6b"
+git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-version = "1.2.0"
+version = "1.2.1"
 
 [[Latexify]]
 deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
-git-tree-sha1 = "3a0084cec7bf157edcb45a67fac0647f88fe5eaf"
+git-tree-sha1 = "537df0b32bc0a99620872cc6d06278e454da4533"
 uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
-version = "0.14.7"
+version = "0.14.10"
 
 [[LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
@@ -472,9 +472,9 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[Luxor]]
 deps = ["Base64", "Cairo", "Colors", "Dates", "FileIO", "ImageMagick", "Juno", "QuartzImageIO", "Random", "Rsvg"]
-git-tree-sha1 = "c018cde3ab7b0867d9f9bd49d290cf15dc9fb926"
+git-tree-sha1 = "8796dfbdde81ec4a811a2d972a6d5f8a7c4113d0"
 uuid = "ae8d54c2-7ccd-5906-9d76-62fc9837b5bc"
-version = "2.9.0"
+version = "2.10.0"
 
 [[METIS_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -578,9 +578,9 @@ version = "1.1.1"
 
 [[OffsetArrays]]
 deps = ["Adapt"]
-git-tree-sha1 = "f64fbf703e7f5af5d82ea7b7385b443bb7765ce9"
+git-tree-sha1 = "b3dfef5f2be7d7eb0e782ba9146a5271ee426e90"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.6.1"
+version = "1.6.2"
 
 [[OpenBLAS32_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
@@ -619,9 +619,9 @@ version = "1.42.4+10"
 
 [[Parsers]]
 deps = ["Dates"]
-git-tree-sha1 = "50c9a9ed8c714945e01cd53a21007ed3865ed714"
+git-tree-sha1 = "223a825cccef2228f3fdbf2ecc7ca93363059073"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "1.0.15"
+version = "1.0.16"
 
 [[Pixman_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -667,9 +667,9 @@ version = "1.0.0"
 
 [[Requires]]
 deps = ["UUIDs"]
-git-tree-sha1 = "cfbac6c1ed70c002ec6361e7fd334f02820d6419"
+git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.1.2"
+version = "1.1.3"
 
 [[Roots]]
 deps = ["Printf"]
diff --git a/Project.toml b/Project.toml
index f2b9a1b..d951e3b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MechanicalSketch"
 uuid = "183f2c95-4b91-4df3-9bd2-197e8f9edb13"
 authors = ["hustf <hustf@users.noreply.github.com> and contributors"]
-version = "0.4"
+version = "0.4.1"
 
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
diff --git a/src/MechanicalSketch.jl b/src/MechanicalSketch.jl
index ab1cac5..18d810c 100644
--- a/src/MechanicalSketch.jl
+++ b/src/MechanicalSketch.jl
@@ -1,6 +1,6 @@
 module MechanicalSketch
 import Luxor
-import Luxor: Drawing,Turtle, Pencolor, Penwidth, Forward, Turn, HueShift, Point, SVGimage,
+import Luxor: Drawing,Turtle, Pencolor, Penwidth, Forward, Turn, HueShift, SVGimage,
 paper_sizes,
 Tiler, Partition,
 rescale,
@@ -283,6 +283,7 @@ include("rope.jl")
 include("power.jl")
 include("table.jl")
 include("curves.jl")
+include("labels.jl")
 include("flow.jl")
 include("runge_kutta.jl")
 include("autodiff_unitfu.jl")
@@ -292,5 +293,6 @@ include("matrix_drawing.jl")
 include("colorlegends.jl")
 include("colorlegends_vector.jl")
 include("place_image.jl")
+include("latex.jl")
 include("chart.jl")
 end # module
diff --git a/src/arrow.jl b/src/arrow.jl
index 7f9aa46..36ecb44 100644
--- a/src/arrow.jl
+++ b/src/arrow.jl
@@ -50,8 +50,6 @@ function arrow(p::Point, Fx::F, Fy::F;
                 linewidth = PT,
                 backgroundcolor = colorant"white",
                 labellength::Bool = false) where {F <: Force}
-    # TODO linewidth - drop or implement
-    # TODO backgroundcolor - drop or implement
     ΔFx = Point(Fx * cos(α), Fx *sin(α))
     ΔFy = Point(-Fy * sin(α), Fy *cos(α))
     ΔF = ΔFx + ΔFy
@@ -74,7 +72,7 @@ function arrow(p::Point, Fx::F, Fy::F;
         arrow_nofill(p, p + ΔFx, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT / 2)
         shaftlength = hypot(ΔFy)
         arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3
-        arrow_nofill(p, p + ΔFy, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT / 2)
+        arrow_nofill(p, p + ΔFy, arrowheadlength = arrowheadlength, arrowheadangle = π / 12, linewidth = PT / 2)
         grestore()
     end
 
@@ -177,36 +175,7 @@ function arrow(p::Point, vx::V, vy::V;
     return nothing
 end
 
-"Add a label at the end point, rounded absolute length with unit."
-function setlabel(endpoint, vx, vy, α)
-    Δvx = (vx * cos(α), vx *sin(α))
-    Δvy = (-vy * sin(α), vy *cos(α))
-    Δv = Δvx + Δvy
-    l = hypot(Δv)
-    rounded = round(unit(l), l, digits = 1)
-    direction = atan(Δv[2], Δv[1]) |> °
-    labdir = direction + 90°
-    labpos = if -90° <= labdir < -70°
-        :S
-    elseif labdir < -20°
-        :SE
-    elseif labdir < 20°
-        :E
-    elseif labdir < 70°
-        :NE
-    elseif labdir < 110°
-        :N
-    elseif labdir < 160°
-        :NW
-    elseif labdir < 200°
-        :W
-    elseif labdir < 250°
-        :SW
-    else
-        :S
-    end
-    label(string(rounded) , labpos, endpoint, offset = div(EM, 4))
-end
+
 function curved_point_arrow_with_vanes(p, Δv, endpoint, linewidth)
     arrowheadangle = pi/24
 
@@ -286,8 +255,8 @@ Same as Luxor.arrow, except for leaving the arrow head unfilled
 """
 function arrow_nofill(startpoint::Point, endpoint::Point;
     linewidth=1.0,
-    arrowheadlength=10,
-    arrowheadangle=pi/8,
+    arrowheadlength= EM / 4.4,
+    arrowheadangle = π / 8,
     decoration = 0.5,
     decorate = () -> ())
     gsave()
@@ -297,7 +266,9 @@ function arrow_nofill(startpoint::Point, endpoint::Point;
     isapprox(startpoint, endpoint) && throw(error("can't draw arrow between two identical points"))
     shaftlength = distance(startpoint, endpoint)
     shaftangle = atan(startpoint.y - endpoint.y, startpoint.x - endpoint.x)
+
     arrowheadtopsideangle = shaftangle + arrowheadangle
+
     # shorten the length so that lines
     # stop before we get to the arrow
     # thus wide shafts won't stick out through the head of the arrow.
@@ -312,23 +283,139 @@ function arrow_nofill(startpoint::Point, endpoint::Point;
     newpath()
     line(Point(fromx, fromy), Point(tox, toy), :stroke)
 
+    # prepare to add decorations at point along shaft
+    for decpos in decoration
+        decpoint = between(startpoint, endpoint, decpos)
+        slp = slope(startpoint, endpoint)
+        @layer begin
+            translate(decpoint)
+            rotate(slp)
+            decorate()
+        end
+    end
+    # draw the arrowhead
+    draw_arrowhead(endpoint, shaftangle; 
+        arrowheadangle, 
+        arrowheadlength, filled = false)
+    grestore()
+end
+"""
+    draw_arrowhead(endpoint, shaftangle; 
+        arrowheadangle = π / 8, 
+        arrowheadlength  = EM / 4.4; filled = false)
+"""
+function draw_arrowhead(endpoint, shaftangle; 
+    arrowheadangle = π / 8, 
+    arrowheadlength  = EM / 4.4, filled = false)
+    arrowheadtopsideangle = shaftangle + arrowheadangle
     # draw the arrowhead
     topx = endpoint.x + cos(arrowheadtopsideangle) * arrowheadlength
     topy = endpoint.y + sin(arrowheadtopsideangle) * arrowheadlength
     arrowheadbottomsideangle = shaftangle - arrowheadangle
     botx = endpoint.x + cos(arrowheadbottomsideangle) * arrowheadlength
     boty = endpoint.y + sin(arrowheadbottomsideangle) * arrowheadlength
-    poly([Point(topx, topy), endpoint, Point(botx, boty)], :stroke)
+    poly([Point(topx, topy), endpoint, Point(botx, boty)], filled ? :fill : :stroke)
+end
+
+
+"""
+    arrow(centerpos::Point, radius, α_sta::Angle, α_end::Angle;
+            linewidth=1.0,
+            arrowheadlength = EM / 4.4,
+            arrowheadangle = π / 8,
+            decoration = 0.5,
+            decorate = () -> (),
+            clockwise = false)
+
+Draw a curved arrow, an arc centered at `centerpos` starting at `α_sta` and
+ending at `α_end`.
+
+The arrowhead at the end is optional. Angles are measured counter-clockwise
+from the positive x-axis.
+
+Arrows don't use the current linewidth setting (`setline()`); you can specify
+the linewidth.
+
+The `decorate` keyword argument accepts a function that can execute code at
+locations on the arrow's shaft. The inherited graphic environment is centered at
+points on the curve between 0 and 1 given by scalar or vector `decoration`, and
+the x-axis is aligned with the direction of the curve at that point.
+"""
+function arrow(centerpos::Point, radius, α_sta::Angle, α_end::Angle;
+        linewidth=1.0,
+        arrowheadlength = EM / 4.4,
+        arrowheadangle = π / 8,
+        decoration = 0.5,
+        decorate = () -> (),
+        clockwise = false,
+        filled = false)
+    @assert radius >= zero(radius) " radius < 0 not supported. Flip keyword argument clockwise?"
+    gsave()
+    setlinejoin("butt")
+    setline(linewidth)
+    # Radius in points, unitless
+    r = scale_to_pt(oneunit(radius)) * radius / unit(radius)
+    while α_sta > α_end
+        α_end += unit(α_end)(360°)
+    end
+    # Full arc ccw angle span. Negative values for clockwise pointing arrows
+    Δα = clockwise ? α_sta - α_end : α_end - α_sta
+
+    # End arrow arc ccw angle span. Negative values for clockwise pointing arrows.
+    Δα_head = unit(α_end)(arrowheadlength / r) * (clockwise ? -1 : 1)
+    α_head = α_end - Δα_head
+
+    # Is there room for one arrow head?
+    if sign(Δα_head) != sign(Δα) && abs(Δα_head) < abs(Δα) 
+        @warn "Arrow head too large, try negative keyword argument arrowheadlength"
+        return
+    end
+
+    # In case we are drawing a sharp arrow head point and a wide shaft, we want to end the arc somewhere inside of the head.
+    # Sign convention: Positive values for cutting off ccw arcs.
+    Δα_cutoff = unit(α_end)((linewidth / 2) / tan(arrowheadangle) / r) * sign(Δα)
+    # Draw the arc to
+    α_cutoff = unit(α_end)(α_end - Δα_cutoff)
+
+    # Local origo at centre
+    translate(centerpos)
+    # Start, end, head crossing points and end of drawn arc in local pt coordinates, where +y is down
+    p0 = Point(r∙cos(α_sta), -r∙sin(α_sta))
+    p1 = Point(r∙cos(α_end), -r∙sin(α_end))
+    p1h = Point(r∙cos(α_end - Δα_head), -r∙sin(α_end - Δα_head))
+    p1c = Point(r∙cos(α_end - Δα_cutoff), -r∙sin(α_end - Δα_cutoff))
+
+    # Unitless, clockwise radian angles from origo 
+    β_sta = rad(-α_sta) / 1rad
+    β_head = rad(-α_head) / 1rad
+    β_cutoff = rad(-α_cutoff) / 1rad
+    β_end = rad(-α_end) / 1rad
+
+    # Draw the cut-off shaft arc
+    newpath()
+    move(p0.x, p0.y)
+    if clockwise
+        arc(0, 0, r, β_sta, β_cutoff, :stroke)
+    else
+        carc(0, 0, r,  β_sta, β_cutoff, :stroke)
+    end
+    closepath()
+    # prepare to add decorations at points along shaft
 
-    # prepare to add decorations at point along shaft
     for decpos in decoration
-        decpoint = between(startpoint, endpoint, decpos)
-        slp = slope(startpoint, endpoint)
+        decorationangle = rescale(decpos, 0, 1, β_sta, β_end)
+        decorationpoint = Point(r * cos(decorationangle), r * sin(decorationangle))
+        perp = perpendicular(decorationpoint)
         @layer begin
-            translate(decpoint)
-            rotate(slp)
+            translate(decorationpoint)
+            rotate(slope(decorationpoint, perp))
             decorate()
         end
     end
+
+    # Head
+    draw_arrowhead(p1, β_cutoff + (clockwise ? -π / 2 : π / 2 ); 
+        arrowheadangle, arrowheadlength, filled)
     grestore()
-end
+    return
+end
\ No newline at end of file
diff --git a/src/autodiff_unitfu.jl b/src/autodiff_unitfu.jl
index 5573f6a..68307d8 100644
--- a/src/autodiff_unitfu.jl
+++ b/src/autodiff_unitfu.jl
@@ -69,7 +69,7 @@ function ∇_R²_in_CQ_out(R²_to_R, ulxy, one_q_out::Quantity{T}) where {T}
     # Note could be sped up using in-place evaluation, but would need to drop StaticArrays.
     map(ulxy) do ul
         r² = ForwardDiff.gradient(R²_to_R, ul, cfg, Val{false}())::SVector{2, T}
-        complex(r²[1]*one_q_out, r²[2]∙one_q_out)
+        complex(r²[1] ∙ one_q_out, r²[2]∙one_q_out)
     end
 end
 
@@ -124,15 +124,7 @@ function ∇_rectangle(CQ_to_Q;
     physwidth = 10.0m,
     physheight = 4.0m,
     cutoff = NaN)
-#=
-    # Discretize per pixel
-    nx = round(Int64, scale_to_pt(physwidth))
-    ny = round(Int64, scale_to_pt(physheight))
-
-    # Iterators for each pixel relative to the center, O
-    pixiterx = (1 - div(nx + 1, 2):(nx - div(nx, 2)))
-    pixitery = (1 - div(ny + 1, 2):(ny - div(ny, 2)))
-=#  
+
     xs, ys = x_y_iterators_at_pixels(;physwidth, physheight)
 
     # Matrix of 2-element static vectors, one per pixel
@@ -167,10 +159,19 @@ Differentiates the complex potential function CQ_to_Q. Evaluates at pixels in a
 physical dimensions and current sketch scale. Limits values to cutoff.
 """
 function clamped_velocity_matrix(CQ_to_Q; physwidth = 1.0m, physheight = 1.0m, cutoff = 0.5m/s)
-    unclamped = ∇_rectangle(CQ_to_Q,
-        physwidth = physwidth,
-        physheight = physheight);
+    unclamped = ∇_rectangle(CQ_to_Q; physwidth, physheight)
     map(unclamped) do u
         hypot(u) > cutoff ? cutoff∙u / hypot(u) : u
     end
 end
+
+"""
+    ∇(f; physwidth, physheight)
+    → derivative function based on linear interpolation between pixels.
+
+The values of the output function varies linearly between pixels. Take
+care if differentiating twice!
+"""
+function ∇(f; physwidth = 10.0m, physheight = 4.0m)
+    matrix_to_function(∇_rectangle(f; physwidth, physheight))
+end
\ No newline at end of file
diff --git a/src/chart.jl b/src/chart.jl
index 9e1c5f6..9658fa2 100644
--- a/src/chart.jl
+++ b/src/chart.jl
@@ -41,7 +41,12 @@ end
                 rotatehue_degrees_total = 270°, firstcolor = PALETTE[3],
                 height = 2EM, top_circle_radius = 0.04EM)
 
+    - origo is a Point.
+    - samples is the vector to plot as bars
+    - width is the horizontal distance over which to distribute bar
+
 Draw a bar plot without decorations; arrows and the first sample in the same colour.
+
 """
 function draw_barplot(origo, samples, width; 
                       first_sample_no = 1, 
@@ -55,4 +60,74 @@ function draw_barplot(origo, samples, width;
         sethue(color_from_palette("red"))
         _draw_bars(origo, samples, width; rotatehue_degrees_total, first_sample_no, top_circle_radius)
     end
-end
\ No newline at end of file
+end
+
+"""
+    color_matrix(f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true,
+                 legend = missing, kwargs...) 
+    -> (color_matrix, legend)
+
+Return a matrix of colors suitable for display with 'place_image'. Colors 
+are found by calling 'legend' with single values out of function f_xy.
+
+'legend' can be input as a keyword argument. If not, a legend is generated
+internally, as modified by keyword arguments relevant to the proposed legend type. 
+See: BinLegend or BinLegendVector.
+
+The legend is also returned. Internally generated legends are of type 'AbstractColorLegend',
+and can be displayed with 'draw_legend'.
+"""
+function color_matrix(f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true,
+    legend = missing, kwargs...) 
+
+    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight, centered)
+    valmat = [f_xy(x, y) for y in ys, x in xs];
+    _color_matrix(valmat, f_xy;  legend, kwargs...)
+end
+function _color_matrix(valmat::Matrix{T}, f_xy;  legend = missing, centered = true, kwargs...) where T<:Tuple
+    magnitude_matrix = hypot.(valmat)
+    if ismissing(legend)
+        legend = propose_legend(magnitude_matrix; kwargs...)
+    end
+    streamlinepixels = falses(size(valmat))
+    streamlines_matrix!(streamlinepixels, f_xy)
+    complex_convolution_matrix = convolution_matrix(f_xy, streamlinepixels)
+    curmat = lic_matrix_current(complex_convolution_matrix, 0, zeroval = -0.0)
+    colmat = map(legend.(magnitude_matrix), curmat) do col, lu
+        color_with_lumin(col, 50 + 50 * lu)
+    end
+    colmat, legend
+end
+function _color_matrix(valmat::Matrix, f_xy;  legend = missing, centered = true, kwargs...)
+    @error "Not implemented"
+end
+
+"""
+    propose_legend(scalmat::Matrix{T}, kwargs...) where T <: Number
+    -> AbstractColorLegend
+
+Return a proposed legend based on the unit and values of magmat.
+Keyword arguments set here overrule proposed arguments. See:
+BinLegend or BinLegendVector. Possible values include 'missing',
+which will also overrule proposals from this function.
+"""
+function propose_legend(mat::Matrix{T}; kwargs...) where T <: Number
+    lum = 50
+    mi, ma = lenient_min_max(mat)
+    maxlegend = ma  # kwargs overrule if present
+    noofbins = if !haskey(kwargs, :binwidth) && !haskey(kwargs, :noofbins) && !haskey(kwargs, :binbounds)
+        # None of three present
+        7
+    else
+        # Sufficient arguments are present.
+        # IF noofbins is actually contained in kwargs, it will be placed after in the call, and thus overrule:
+        missing
+    end
+    if T <: Velocity
+        name = :Speed # kwargs overrule if present
+        colorscheme = ColorSchemes.isoluminant_cgo_80_c38_n256.colors .|> co -> color_with_lumin(co, lum)
+        BinLegend(;maxlegend, noofbins, colorscheme, name, kwargs...)
+    else
+        @error "$T scalar quantity not implemented."
+    end
+end
diff --git a/src/curves.jl b/src/curves.jl
index ff08bc8..f1f1ccc 100644
--- a/src/curves.jl
+++ b/src/curves.jl
@@ -97,24 +97,5 @@ function trace_rotate_hue(center, vx, vy; rotatehue_degrees_total = 270°)
     end
     grestore()
 end
-"""
-    box_fill_outline(pt_topleft, colfill; 
-                     height = row_height(), width = EM, luminfac = 2)
 
-Draw a filled box.
-Default height adapts to the current 'Top API' font, which is used for text tables.
-Outline color luminosity is 200% the luminosity of fill color, if possible.
-For dark backgrounds, consider luminfac = 0.5
-"""
-function box_fill_outline(pt_topleft, colfill; 
-                          height = row_height(), width = EM, luminfac = 2)
-    gsave()
-    sethue(colfill)
-    box(pt_topleft, pt_topleft + (width, height), :fill)
-    # outline
-    lu = get_current_lumin()
-    sethue(color_with_lumin(colfill, min(100, luminfac * lu)))
-    box(pt_topleft, pt_topleft + (width, height), :stroke)
-    grestore()
-end
 setline(x::Length) = setline(scale_to_pt(x))
\ No newline at end of file
diff --git a/src/flow.jl b/src/flow.jl
index 294be73..8fa6ee6 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -110,66 +110,6 @@ function normalize_datarange(A)
 end
 
 
-"""
-    hue_from_complex_argument!(color_values, A)
-    Changes the hue in collection color_values by roting based on complex arguments in A
-
-Take a matrix of colors, modify elements of it by changing hue to be according to
-the argument (polar angle) of a complex number
-
-Ref. domain colouring, complex plane https://en.wikipedia.org/wiki/Domain_coloring
-"""
-function hue_from_complex_argument!(color_values, A)
-    @error "Replaced by Legend"
-    @assert length(A) == length(color_values)
-    for i in 1:length(A) # Works for matrices also
-        cqua = A[i]
-        if !isnan(cqua) && !isinf(cqua)
-            ang = angle(complex(cqua))rad
-            col = color_values[i]
-            if col != RGB(0.0, 0.0, 0.0)
-                color_values[i] = rotate_hue(col, ang)
-            end
-        end
-    end
-end
-
-
-"""
-    color_matrix(qua::AbstractArray; normalize_data_range = true) -> RGBA
-
-Map a collection of quantities to colors, transparent
-pixels for NaN and Inf values.
-
-# TODO adapt ColorLegend for this functionality?
-"""
-function color_matrix(qua::AbstractArray; normalize_data_range = true)
-    @error "color_matrix is replaced"
-    # Boolean collection, valid elements which should be opaque
-    valid_element = map(x-> isnan(x) || isinf(x) ? false : true, qua)
-
-     # Map value or magnitude to [0.0, 1.0], invalid elements are 1.0
-    norm_real_values = hypot.(normalize_data_range ? normalize_datarange(qua) : qua)
-
-    # Map [0.0, 1.0] to perceived luminosity color map
-    color_values = if eltype(qua) <: RealQuantity || eltype(qua) <: Real
-        get(absolute_scale(), norm_real_values)
-    else
-        # Complex quanitity
-        get(complex_arg0_scale(), norm_real_values)
-    end
-    # If input is complex, adjust the hue
-    if !(eltype(qua) <: RealQuantity || eltype(qua) <: Real)
-        # Set hue from complex argument (angle)
-        # while not changing perceived luminance
-        hue_from_complex_argument!(color_values, qua)
-    end
-
-    # Add transparency for any invalid elements
-    map(color_values, valid_element) do co, valid
-        RGBA(ColorSchemes.red(co), ColorSchemes.green(co), ColorSchemes.blue(co), Float64(valid))
-    end
-end
 
 
 
diff --git a/src/labels.jl b/src/labels.jl
new file mode 100644
index 0000000..81445b4
--- /dev/null
+++ b/src/labels.jl
@@ -0,0 +1,73 @@
+"""
+    box_fill_outline(pt_topleft, colfill; 
+                     height = row_height(), width = EM, luminfac = 2,
+                     opacity = missing)
+
+Draw a filled box.
+Default height adapts to the current 'Top API' font, which is used for text tables.
+Outline color luminosity is 200% the luminosity of fill color, if possible.
+For dark backgrounds, consider luminfac = 0.5.
+By default, opacity is not changed from the current state. Range 0-1, where 1 is opaque.
+"""
+function box_fill_outline(pt_topleft, colfill; 
+                          height = row_height(), width = EM, luminfac = 2, boxopacity = missing)
+    gsave()
+    !ismissing(boxopacity) && setopacity(boxopacity)
+    sethue(colfill)
+    box(pt_topleft, pt_topleft + (width, height), :fill)
+    # outline
+    lu = get_current_lumin()
+    sethue(color_with_lumin(colfill, min(100, luminfac * lu)))
+    !ismissing(boxopacity) && setopacity(boxopacity)
+    box(pt_topleft, pt_topleft + (width, height), :stroke)
+    grestore()
+end
+"""
+label_boxed(pt, str; colfill = color_with_lumin(get_current_RGBA(), 0.5 * get_current_lumin()),
+               luminfac = 1.0, boxopacity = 0.75, valign = :bottom)
+
+Using the 'Toy API. For keyword arguments, see 'box_fill_outline'. Could be extended with 'text' arguments.
+"""
+function label_boxed(pt, str; colfill = color_with_lumin(get_current_RGBA(), 0.5 * get_current_lumin()),
+               luminfac = 1.0, boxopacity = 0.75, valign = :bottom)
+    @layer begin
+        box_fill_outline(pt, colfill; luminfac, height = -FS / 2, width = pixelwidth(str), boxopacity)
+        sethue(PALETTE[10])
+        text(str, pt; valign)
+    end
+end
+
+"""
+    setlabel(endpoint, vx, vy, α)
+
+Add a label at the end point, rounded absolute length with unit.
+"""
+function setlabel(endpoint, vx, vy, α)
+    Δvx = (vx * cos(α), vx *sin(α))
+    Δvy = (-vy * sin(α), vy *cos(α))
+    Δv = Δvx + Δvy
+    l = hypot(Δv)
+    rounded = round(unit(l), l, digits = 1)
+    direction = atan(Δv[2], Δv[1]) |> °
+    labdir = direction + 90°
+    labpos = if -90° <= labdir < -70°
+        :S
+    elseif labdir < -20°
+        :SE
+    elseif labdir < 20°
+        :E
+    elseif labdir < 70°
+        :NE
+    elseif labdir < 110°
+        :N
+    elseif labdir < 160°
+        :NW
+    elseif labdir < 200°
+        :W
+    elseif labdir < 250°
+        :SW
+    else
+        :S
+    end
+    label(string(rounded) , labpos, endpoint, offset = div(EM, 4))
+end
\ No newline at end of file
diff --git a/src/latex.jl b/src/latex.jl
new file mode 100644
index 0000000..b84d621
--- /dev/null
+++ b/src/latex.jl
@@ -0,0 +1,215 @@
+"""
+    tex2svg_string(formula::String)
+    tex2svg_string(formula::LaTeXString; modifylatex = true, modifysvgcolor = true)
+    --> string
+
+Relies on NodeJS to call the package matjax, tex2svg. Based off DocumenterEpub 'prerender_mathjax'.
+You may need elevated privileges to spawn NodeJS.
+
+modifylatex = true: Calls 'modify_latex', which subscripts bracketed indexes and other things.
+modifysvgcolor = true: Calls Replace currentColor with the last color and opacity from 'setcolor', 'setopacity' or 'sethue'.
+"""
+function tex2svg_string(formula::String)
+    mathjaxfile = escape_string(realpath(joinpath(@__DIR__, "..", "mathjax", "node-main.js")))
+    @assert isfile(mathjaxfile)
+    adapted_formula = escape_string(replace(formula, '\'' => "\\prime"))
+    svgs = cd(dirname(mathjaxfile)) do
+        js_command = """
+        require('$(mathjaxfile)').init({
+            loader: {load: ['input/tex', 'output/svg', '[tex]/boldsymbol']}, svg: {
+                fontCache: 'none',
+                'localID':$(rand(1:100000))}
+        }).then((mathjaxfile) => {
+            const svg = MathJax.tex2svg('$adapted_formula');
+            console.log(MathJax.startup.adaptor.innerHTML(svg));
+        }).catch((err) => console.log("NodeJS was ran and output: " + err.message));
+        """
+        fullcmd = `$(NodeJS.nodejs_cmd()) -e $js_command`
+        try String(read(fullcmd))
+        catch err
+            @warn "If on windows, you may need to run your editor or julia with some administrator privileges."
+            @warn "Program: $(NodeJS.nodejs_cmd())"
+            @warn "Option: -e"
+            @warn "js_command: $js_command"
+            @error err
+        end
+    end
+    @assert !startswith(svgs, "NodeJS was ran and output") svgs
+    svgs
+end
+function tex2svg_string(formula::LaTeXString; modifylatex = true, modifysvgcolor = true)
+    bw = modifylatex ? modify_latex(formula) : String(formula)
+    sv = tex2svg_string(bw)
+    modifysvgcolor ? modify_svg_color(sv) : sv
+end
+
+
+
+
+
+"""
+    modify_latex(formula)
+    --> string
+Iteratively 
+    - get rid of latex outer enclosing dollar and paranthesises.
+    - get rid of latex outer enclosing paranthesises.
+    - drop ∙ used within units: N∙m => Nm
+    - drop {\\vysmblkcircle}, which is used within units
+    - fix display of unit superscripts
+    - remove underscore escaping
+    - reorder √ function to square root
+    - reorder dot multiplication function
+    - replace bold font unicode with
+    - replace square bracket with contents in subscript
+"""
+function modify_latex(formula::String)::String
+    st = String(formula)
+    if startswith(st, "\$") && endswith(st, "\$")
+        # get rid of latex outer enclosing dollar and paranthesises.
+        modify_latex(st[nextind(st, begin, 1):prevind(st, end, 1)])
+    elseif startswith(st, "\\left(") && endswith(st, "\\right)")
+        # get rid of latex outer enclosing paranthesises.
+        output = st[nextind(st, begin, 6):prevind(st, end, 7)]
+        #println("1")
+        #@show st
+        #@show output
+        modify_latex(output)
+    #elseif occursin('∙', st)
+    #    # drop ∙ used within units: N∙m => Nm
+    #    modify_latex(replace(st, '∙' => ""))
+    elseif occursin("{\\vysmblkcircle}", st)
+        # drop {\\vysmblkcircle} used within units (no space)
+        output = replace(st, "{\\vysmblkcircle}" => "")
+        #println("2")
+        #@show st
+        #@show output
+        modify_latex(output)
+    elseif occursin("\\^-{^", st)
+        # fix display of unit superscripts
+        output = replace(st, "\\^-{^" => "{^-}{^")
+        #println("3")
+        #@show st
+        #@show output
+        modify_latex(output)
+    elseif occursin("\\_", st)
+        # remove underscore escaping
+        regex = r"\\_(\w+)"
+        ma = match(regex, st)
+        if !isnothing(ma)
+            captured = ma.captures[1]
+            matched = ma.match
+            replacement = "_{" * captured * "}"
+            output = replace(st, matched => replacement; count = 1)
+            #println("4")
+            #@show st
+            #@show output
+            @assert st != output
+            modify_latex(output)
+        else
+            @error "Bad regex"
+        end
+    elseif occursin("\\sqrt\\left( ", st)
+        # reorder √ function to square root
+        #"\\sqrt\\left( capture \\right)"
+        # => \\sqrt{capture}
+        regex = r"\\sqrt\\left\( (.*?) \\right\)"
+        ma = match(regex, st)
+        if !isnothing(ma)
+            captured = ma.captures[1]
+            matched = ma.match
+            replacement = "\\sqrt{" * captured * "} "
+            output = replace(st, matched => replacement; count = 1)
+            #println("5")
+            #@show st
+            #@show output
+            modify_latex(output)
+        else
+            @error "Bad regex"
+        end
+    elseif occursin("\\vysmblkcircle\\left( ", st)
+        # reorder dot multiplication function:
+        # \\vysmblkcircle\\left(factor1, factor2 \\right)
+        # =>  factor1 \\cdot factor2
+        # A more solid approach is probably defining: #@latexrecipe ∙
+        regex = r"\\vysmblkcircle\\left\( (.*?), (.*?) \\right\)"
+        ma = match(regex, st)
+        if !isnothing(ma)
+            captured1 = ma.captures[1]
+            captured2 = ma.captures[2]
+            matched = ma.match
+            replacement = captured1 * " \\cdot " * captured2
+            output = replace(st, matched => replacement; count = 1)
+            #println("6")
+            #@show st
+            #@show output
+            modify_latex(output)
+        else
+            @error "Bad regex"
+        end
+    elseif occursin(r"\\bf([a-zA-Z])", st)
+        # Replace bold font unicode
+        # \\bfX
+        # =>  \\\boldsymbol{X}
+        # A more solid approach is probably defining: #@latexrecipe ∙
+        regex = r"\\bf([a-zA-Z])"
+        ma = match(regex, st)
+        if !isnothing(ma)
+            captured = ma.captures[1]
+            matched = ma.match
+            replacement = "\\boldsymbol{" * captured * "}"
+            output = replace(st, matched => replacement; count = 1)
+            #println("7")
+            #@show st
+            #@show output
+            modify_latex(output)
+        else
+            @error "Bad regex"
+        end
+    else
+        # replace square bracket with contents in subscript
+        # Also: The only exit point from this iterative function.
+        regex = r"\\left\[(.*?)\\right\]"
+        ma = match(regex, st)
+        if !isnothing(ma)
+            captured = ma.captures[1]
+            matched = ma.match
+            replacement = "_{" * captured * "}"
+            output = replace(st, matched => replacement; count = 1)
+            #println("8")
+            #@show st
+            #@show output
+            @assert st != output
+            modify_latex(output)
+        else
+            st
+        end
+    end
+end
+modify_latex(formula::LaTeXString) = modify_latex(String(formula))
+
+
+"""
+    modify_svg_color(svgstring::String)::String
+
+Remove local color definitions (normally defined by style sheet) with
+an extended style definition.
+Both 'stroke' and 'fill' styles equal the last color and opacity from 'setcolor', 'setopacity' or 'sethue'
+Primarily for svgs output by MathJax via NodeJS.
+"""
+function modify_svg_color(svgstring::String)::String
+    c = get_current_RGBA()
+    r = Int(round(c.r * 255))
+    g = Int(round(c.g * 255))
+    b = Int(round(c.b * 255))
+    a = c.alpha
+    removes = "stroke=\"currentColor\" fill=\"currentColor\" "
+    step1 = replace(svgstring, removes => "")
+    olds = "style=\""
+    cols = "rgb($r, $g, $b)"
+    opacs = "-opacity: $a"
+    sfills = "fill: " * cols * "; " * "fill" * opacs * "; "
+    sstrokes = "stroke: " * cols * "; " * "stroke" * opacs * "; "
+    styles = sfills * sstrokes
+    news = olds * styles
+    replace(step1, olds  => news)
+end
\ No newline at end of file
diff --git a/src/place_image.jl b/src/place_image.jl
index 359003b..6b16930 100644
--- a/src/place_image.jl
+++ b/src/place_image.jl
@@ -1,22 +1,3 @@
-"""
-    _CairoSurfaceBase(colmat::VecOrMat{T}) where T
-
-Convert a matrix or vector to a cairo surface suitable for ´placeimage´.
-One pixel per element.
-
-"""
-function _CairoSurfaceBase(colmat::VecOrMat{T}) where T
-    msg = "Can't display matrices of eltype($(eltype(colmat))). \n\tConsider defining a subtypes(AbstractColorLegend) to convert."
-    @assert T <: Union{RGB, RGBA, Luxor.Cairo.ARGB32, HSV, HSVA,  HSL, HSLA,LCHab, LCHabA, LCHuv, LCHuvA}  msg
-    tempfilename = joinpath(tempname(), "temp_sketch.png")
-    # write to a file using ImageIO. Results for transparent colours seem to differ.
-    save(File(format"PNG", tempfilename), colmat)
-    # Read this image using Cairo's interface to get the right format for placing into other cairo surfaces.
-    png = readpng(tempfilename);
-    rm(tempfilename)
-    png
-end
-
 """
     place_image(p::Point, colmat::VecOrMat; centered = true, alpha = missing)
     -> (upper left point, lower right point)
@@ -139,174 +120,23 @@ Note that the width can be inaccurate, as this is ambiguous for svg input.
 get_width_height(image) =  scale_pt_to_unit(m) .* (image.width, image.height)
 get_width_height(image::Matrix) =  scale_pt_to_unit(m) .* reverse(size(matrix))
 
-"""
-    tex2svg_string(formula::String)
-    tex2svg_string(formula::LaTeXString; modifylatex = true, modifysvgcolor = true)
-    --> string
 
-Relies on NodeJS to call the package matjax, tex2svg. Based off DocumenterEpub 'prerender_mathjax'.
-You may need elevated privileges to spawn NodeJS.
 
-modifylatex = true: Calls 'modify_latex', which subscripts bracketed indexes and other things.
-modifysvgcolor = true: Calls Replace currentColor with the last color and opacity from 'setcolor', 'setopacity' or 'sethue'.
 """
-function tex2svg_string(formula::String)
-    mathjaxfile = escape_string(realpath(joinpath(@__DIR__, "..", "mathjax", "node-main.js")))
-    @assert isfile(mathjaxfile)
-    adapted_formula = escape_string(replace(formula, '\'' => "\\prime"))
-    cd(dirname(mathjaxfile)) do
-        js_command = """
-        require('$(mathjaxfile)').init({
-            loader: {load: ['input/tex', 'output/svg']}, svg: {
-                fontCache: 'none',
-                'localID':$(rand(1:100000))}
-        }).then((mathjaxfile) => {
-            const svg = MathJax.tex2svg('$adapted_formula');
-            console.log(MathJax.startup.adaptor.innerHTML(svg));
-        }).catch((err) => console.log("NodeJS was ran and output: " + err.message));
-        """
-        fullcmd = `$(NodeJS.nodejs_cmd()) -e $js_command`
-        try String(read(fullcmd))
-        catch err
-            @warn "If on windows, you may need to run your editor or julia with some administrator privileges."
-            @warn "Program: $(NodeJS.nodejs_cmd())"
-            @warn "Option: -e"
-            @warn "js_command: $js_command"
-            @error err
-        end
-    end # -> resulting svg string
-end
-function tex2svg_string(formula::LaTeXString; modifylatex = true, modifysvgcolor = true)
-    bw = modifylatex ? modify_latex(formula) : String(formula)
-    sv = tex2svg_string(bw)
-    modifysvgcolor ? modify_svg_color(sv) : sv
-end
+    _CairoSurfaceBase(colmat::VecOrMat{T}) where T
 
+Convert a matrix or vector to a cairo surface suitable for ´placeimage´.
+One pixel per element. Relies on ImageIO / ImageMagick because it seems 
+to work more easily.
 """
-    modify_latex(formula)
-    --> string
-Iteratively 
-    - get rid of latex outer enclosing dollar and paranthesises.
-    - get rid of latex outer enclosing paranthesises.
-    - drop ∙ used within units: N∙m => Nm
-    - drop {\\vysmblkcircle}, which is used within units
-    - fix display of unit superscripts
-    - remove underscore escaping
-    - reorder √ function to square root
-    - reorder dot multiplication function
-    - replace square bracket with contents in subscript
-"""
-function modify_latex(formula::String)::String
-    st = String(formula)
-    if startswith(st, "\$") && endswith(st, "\$")
-        # get rid of latex outer enclosing dollar and paranthesises.
-        modify_latex(st[nextind(st, begin, 1):prevind(st, end, 1)])
-    elseif startswith(formula, "\\left(") && endswith(formula, "\\right)")
-        # get rid of latex outer enclosing paranthesises.
-        modify_latex(st[nextind(st, begin, 6):prevind(st, end, 7)])
-    #elseif occursin('∙', st)
-    #    # drop ∙ used within units: N∙m => Nm
-    #    modify_latex(replace(st, '∙' => ""))
-    elseif occursin("{\\vysmblkcircle}", st)
-        # drop {\\vysmblkcircle} used within units (no space)
-        modify_latex(replace(st, "{\\vysmblkcircle}" => ""))
-    elseif occursin("\\^-{^", st)
-        # fix display of unit superscripts
-        modify_latex(replace(st, "\\^-{^" => "{^-}{^"))
-    elseif occursin("\\_", st)
-        # remove underscore escaping
-        regex = r"\\_(\w+)"
-        ma = match(regex, st)
-        if !isnothing(ma)
-            captured = ma.captures[1]
-            matched = ma.match
-            replacement = "_{" * captured * "}"
-            output = replace(st, matched => replacement; count = 1)
-            @assert st != output
-            modify_latex(output)
-        else
-            @error "Bad regex"
-        end
-    elseif occursin("\\sqrt\\left( ", st)
-        # reorder √ function to square root
-        #"\\sqrt\\left( capture \\right)"
-        # => \\sqrt{capture}
-        regex = r"\\sqrt\\left\( (.*?) \\right\)"
-        ma = match(regex, st)
-        if !isnothing(ma)
-            captured = ma.captures[1]
-            matched = ma.match
-            replacement = "\\sqrt{" * captured * "}"
-            output = replace(st, matched => replacement; count = 1)
-            modify_latex(output)
-        else
-            @error "Bad regex"
-        end
-    elseif occursin("\\vysmblkcircle\\left( ", st)
-        # reorder dot multiplication function:
-        # \\vysmblkcircle\\left(factor1, factor2 \\right)
-        # =>  factor1 \\cdot factor2
-        # A more solid approach is probably defining: #@latexrecipe ∙
-        regex = r"\\vysmblkcircle\\left\( (.*?), (.*?) \\right\)"
-        ma = match(regex, st)
-        if !isnothing(ma)
-            captured1 = ma.captures[1]
-            captured2 = ma.captures[2]
-            matched = ma.match
-            replacement = captured1 * " \\cdot " * captured2
-            output = replace(st, matched => replacement; count = 1)
-#=
-            println("\n# reorder dot multiplication function:")
-            @show st
-            @show captured1
-            @show captured2
-            @show matched
-            @show output
-            @assert st != output
-=#
-            modify_latex(output)
-        else
-            @error "Bad regex"
-        end
-    else
-        # replace square bracket with contents in subscript
-        # Also: The only exit point from this iterative function.
-        regex = r"\\left\[(.*?)\\right\]"
-        ma = match(regex, st)
-        if !isnothing(ma)
-            captured = ma.captures[1]
-            matched = ma.match
-            replacement = "_{" * captured * "}"
-            output = replace(st, matched => replacement; count = 1)
-            @assert st != output
-            modify_latex(output)
-        else
-            st
-        end
-    end
+function _CairoSurfaceBase(colmat::VecOrMat{T}) where T
+    msg = "Can't display matrices of eltype($(eltype(colmat))). \n\tConsider defining a subtypes(AbstractColorLegend) to convert."
+    @assert T <: Union{RGB, RGBA, Luxor.Cairo.ARGB32, HSV, HSVA,  HSL, HSLA,LCHab, LCHabA, LCHuv, LCHuvA}  msg
+    tempfilename = joinpath(tempname(), "temp_sketch.png")
+    # write to a file using ImageIO. Results for transparent colours seem to differ.
+    save(File(format"PNG", tempfilename), colmat)
+    # Read this image using Cairo's interface to get the right format for placing into other cairo surfaces.
+    png = readpng(tempfilename);
+    rm(tempfilename)
+    png
 end
-modify_latex(formula::LaTeXString) = modify_latex(String(formula))
-"""
-    modify_svg_color(svgstring::String)::String
-
-Remove local color definitions (normally defined by style sheet) with
-an extended style definition.
-Both 'stroke' and 'fill' styles equal the last color and opacity from 'setcolor', 'setopacity' or 'sethue'
-"""
-function modify_svg_color(svgstring::String)::String
-    c = get_current_RGBA()
-    r = Int(round(c.r * 255))
-    g = Int(round(c.g * 255))
-    b = Int(round(c.b * 255))
-    a = c.alpha
-    removes = "stroke=\"currentColor\" fill=\"currentColor\" "
-    step1 = replace(svgstring, removes => "")
-    olds = "style=\""
-    cols = "rgb($r, $g, $b)"
-    opacs = "-opacity: $a"
-    sfills = "fill: " * cols * "; " * "fill" * opacs * "; "
-    sstrokes = "stroke: " * cols * "; " * "stroke" * opacs * "; "
-    styles = sfills * sstrokes
-    news = olds * styles
-    replace(step1, olds  => news)
-end
\ No newline at end of file
diff --git a/src/scale.jl b/src/scale.jl
index 0818be6..c6b0664 100644
--- a/src/scale.jl
+++ b/src/scale.jl
@@ -1,4 +1,3 @@
-
 """
     set_scale_sketch(s::T, pixels::Int) where {T <: Length}
     set_scale_sketch(s::T, pixels::Int) where {T <: Force}
diff --git a/src/streamline_convolution.jl b/src/streamline_convolution.jl
index d91eeae..81e3cc2 100644
--- a/src/streamline_convolution.jl
+++ b/src/streamline_convolution.jl
@@ -1,12 +1,16 @@
 """
-    draw_streamlines(center, xs, ys, f_xy, h)
+    draw_streamlines(pt, f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true, 
+        h = 0.1s, nsteps = 10, probability = 0.0001)
+
 where
-    center is local origo, a Point, typically equal to global O.
-    xs, ys are quantity iterators, e.g. (-5.0:0.007942811755361398:5.0)m
+    pt         point on screen corresponding to local origo
+    centered   if true, physwidth and physheight extend symmetrically from origo
+               if false, physwidth and physheight extend on the positive x and y axes.
+                The streamlines are drawn to the right (+x) and above (+y) pt.
+    physwidth, physheight Together with 'scale_to_pt', determines input to f_xy which correspond to pixels.
     f_xy is a function of quantities x and y where the range is also a quantity
-    h is a step-size quantity. If the integration variable is time,
+    h is a step-size quantity, default time. If the integration variable is time,
         it defines the length of time betweeen points. It can be negative.
-Optional keyword arguments:
     nsteps is the number of steps to trace, forwards and backwards from each
         randomly selected point
     probability is the number of streamlines divided by length(xs)*length(ys)
@@ -14,9 +18,13 @@ Optional keyword arguments:
 Draws a few random streamlines, starting at a random number of points picked from (xs, ys).
 Streamlines are found by starting at the randomly picked centre of their extent.
 They may extend outside of the rectangle defined by xs and ys.
+A better visualization might be based on where particles enter.
 """
-function draw_streamlines(center, xs, ys, f_xy, h; nsteps = 10, probability = 0.0001)
+function draw_streamlines(pt, f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true, 
+                          h = 0.1s, nsteps = 10, probability = 0.0001)
     gsave()
+    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight, centered)
+
     cx = xs[1]::Quantity{Float64}
     cy = ys[1]
     @assert cx isa Quantity{Float64}
@@ -32,10 +40,11 @@ function draw_streamlines(center, xs, ys, f_xy, h; nsteps = 10, probability = 0.
             # Find the streamline path
             rk4_steps!(f_xy, vxf, vyf, h)
             rk4_steps!(f_xy, vxb, vyb, -h)
-            # Reorder to most recent position to oldest position
+            # Reorder to most recent position -> oldest position
             vx = vcat(reverse(vxf), vxb)
             vy = vcat(reverse(vyf), vyb)
-            trace_diminishing(center, vx, vy)
+            # Draw it
+            trace_diminishing(pt, vx, vy)
         end
     end
     grestore()
@@ -62,8 +71,8 @@ f_xy    Function of x and y for visualization
 n_xy    Noise function of x and y. The spectrum should be adapted for f_xy.
 x_mid   Coordinate for which to evaluate. Streamlines forward and back are found for convolution.
 y_mid   Coordinate for which to evaluate
-f_s Sampling frequency
-f_0 Frequency of interest
+f_s     Sampling frequency
+f_0     Frequency of interest
 
 
 The resulting complex number ŝ yields:
@@ -288,7 +297,7 @@ end
 
 
 """
-    streamlines_add!(v_xy::Extrapolation, streamlinepixels::BitMatrix; centered = true, targetdensity = 0.42)
+    streamlines_matrix!(streamlinepixels::BitMatrix, v_xy; centered = true, targetdensity = 0.42)
 
 Input:
     v_xy: (Q, Q)  → (Q, Q)    Function taking a tuple of quantities, outputs tuple of quantities
@@ -299,7 +308,7 @@ Input:
 
 Input matrix streamlinepixels is modified in place.
 """
-function streamlines_add!(v_xy::Extrapolation, streamlinepixels::BitMatrix; centered = true, targetdensity = 0.42)
+function streamlines_matrix!(streamlinepixels::BitMatrix, v_xy; centered = true, targetdensity = 0.42)
     ny, nx = size(streamlinepixels)
     physwidth = nx * scale_pt_to_unit(m)
     physheight = ny * scale_pt_to_unit(m)
@@ -329,7 +338,30 @@ function streamlines_add!(v_xy::Extrapolation, streamlinepixels::BitMatrix; cent
         # Check against a reasoable upper limit, 2500 * 4000 pixels
         count > 10^7 && break
     end
-    nothing
+    streamlinepixels
+end
+
+"""
+    streamlines_matrix(v_xy; physwidth = 10.0m, physheight = 4.0m, centered = true, targetdensity = 0.42)
+
+Input:
+    v_xy: (Q, Q)  → (Q, Q)      Function taking a tuple of quantities, outputs tuple of quantities
+    centered          Bool      Keyword argument. If true, values at the centre of the matrix corresponds to (x,y) = (0, 0)m
+    targetdensity     Real      Keyword argument. Fraction of pixels covered by streamlines.
+    physwidth         Quantity  Together with 'scale_to_pt', determines input to v_xy which correspond to pixels.
+    physheight        Quantity
+    centered          Bool
+
+Output:
+    streamlinepixels: A matrix of booleans where false indicates unoccupied pixel. The matrix
+        follows the image manipulation convention: xs correspond to output rows.
+"""
+function streamlines_matrix(v_xy; physwidth = 10.0m, physheight = 4.0m, centered = true, targetdensity = 0.42)
+    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight, centered)
+    velocity_matrix = [v_xy(x, y) for y in ys, x in xs]
+    streamlinepixels = falses(size(velocity_matrix))
+    streamlines_matrix!(streamlinepixels, v_xy)
+    streamlinepixels
 end
 
 """
@@ -535,7 +567,7 @@ function convolution_matrix(matrix::Array{Quantity{Complex{Float64}, D, U}, 2};
                            centered = true, targetdensity = 0.42) where {D, U}
     v_xy = matrix_to_function(matrix)
     streamlinepixels = falses(size(matrix))
-    streamlines_add!(v_xy::Extrapolation, streamlinepixels::BitMatrix; centered, targetdensity)
+    streamlines_matrix!(streamlinepixels::BitMatrix, v_xy; centered, targetdensity)
     convolution_matrix(v_xy, streamlinepixels; centered)
 end
 
@@ -602,7 +634,7 @@ function convolution_matrix(f_xy, streamlinepixels::BitMatrix; centered = true)
         ny, nx = size(streamlinepixels)
         physwidth = nx * scale_pt_to_unit(m)
         physheight = ny * scale_pt_to_unit(m)
-        function_to_interpolated_function(f_xy; physwidth = physwidth, physheight = physheight)
+        function_to_interpolated_function(f_xy; physwidth, physheight)
     end
     # Prepare a noise function:
     # Simplex based continuous noise function
@@ -632,7 +664,8 @@ where
 
 Output phase and amplitude information for 'streamlinepixels'. Used for rendering line-integral-convolution showing a static vector field.
 """
-function convolution_matrix(f_xy::Extrapolation, n_xy::Extrapolation, streamlinepixels::BitMatrix; centered = true)
+function convolution_matrix(f_xy, n_xy::Extrapolation, streamlinepixels::BitMatrix; centered = true)
+    # used to be f_xy::Extrapolation
     ny, nx = size(streamlinepixels)
     physwidth = nx * scale_pt_to_unit(m)
     physheight = ny * scale_pt_to_unit(m)
diff --git a/test/Manifest.toml b/test/Manifest.toml
index ad079e8..b6e2b3a 100644
--- a/test/Manifest.toml
+++ b/test/Manifest.toml
@@ -2,15 +2,14 @@
 
 [[AbstractFFTs]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "8ed9de2f1b1a9b1dee48582ad477c6e67b83eb2c"
+git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
 uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-version = "1.0.0"
+version = "1.0.1"
 
 [[AbstractTrees]]
-deps = ["Markdown"]
-git-tree-sha1 = "33e450545eaf7699da1a6e755f9ea65f14077a45"
+git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
 uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
-version = "0.3.3"
+version = "0.3.4"
 
 [[Adapt]]
 deps = ["LinearAlgebra"]
@@ -76,9 +75,9 @@ version = "3.10.2"
 
 [[ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
-git-tree-sha1 = "4bffea7ed1a9f0f3d1a131bbcd4b925548d75288"
+git-tree-sha1 = "5e9769a17f17b587c951d57ba4319782b40c3513"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.10.9"
+version = "0.10.10"
 
 [[ColorVectorSpace]]
 deps = ["ColorTypes", "Colors", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "StatsBase"]
@@ -156,9 +155,9 @@ version = "0.8.3"
 
 [[Documenter]]
 deps = ["Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "21fb992ef1b28ff8f315354d3808ebf4a8fa6e45"
+git-tree-sha1 = "3ebb967819b284dc1e3c0422229b58a40a255649"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.26.2"
+version = "0.26.3"
 
 [[DocumenterTools]]
 deps = ["Base64", "DocStringExtensions", "Documenter", "FileWatching", "LibGit2", "Sass"]
@@ -218,10 +217,10 @@ uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
 version = "3.3.9+7"
 
 [[FileIO]]
-deps = ["Pkg"]
-git-tree-sha1 = "fee8955b9dfa7bec67117ef48085fb2b559b9c22"
+deps = ["Pkg", "Requires", "UUIDs"]
+git-tree-sha1 = "8800ec70aee7292931a3d3c10a3be3445b9c6141"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.4.5"
+version = "1.6.2"
 
 [[FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
@@ -288,9 +287,9 @@ version = "0.53.0+0"
 
 [[GeometryBasics]]
 deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
-git-tree-sha1 = "4d4f72691933d5b6ee1ff20e27a102c3ae99d123"
+git-tree-sha1 = "82853ebc70db4f5a3084853738c68fd497b22c7c"
 uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-version = "0.3.9"
+version = "0.3.10"
 
 [[Gettext_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
@@ -352,16 +351,16 @@ uuid = "a09fc81d-aa75-5fe9-8630-4744c3626534"
 version = "0.8.20"
 
 [[ImageIO]]
-deps = ["FileIO", "Netpbm", "PNGFiles"]
-git-tree-sha1 = "0d6d09c28d67611c68e25af0c2df7269c82b73c7"
+deps = ["FileIO", "Netpbm", "PNGFiles", "TiffImages", "UUIDs"]
+git-tree-sha1 = "015adb094b2d81167821df8a24c180a8eee22f2c"
 uuid = "82e4d734-157c-48bb-816b-45c225c6df19"
-version = "0.4.1"
+version = "0.5.3"
 
 [[ImageMagick]]
 deps = ["FileIO", "ImageCore", "ImageMagick_jll", "InteractiveUtils", "Libdl", "Pkg", "Random"]
-git-tree-sha1 = "02558f83932fde6ebd3ab007dbff6bd8740a8247"
+git-tree-sha1 = "f35fb0080591e6f271b7e1cf4c8323a4a78d5150"
 uuid = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
-version = "1.1.6"
+version = "1.1.7"
 
 [[ImageMagick_jll]]
 deps = ["JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg", "Zlib_jll", "libpng_jll"]
@@ -448,15 +447,15 @@ uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
 version = "2.10.0+3"
 
 [[LaTeXStrings]]
-git-tree-sha1 = "c7aebfecb1a60d59c0fe023a68ec947a208b1e6b"
+git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
-version = "1.2.0"
+version = "1.2.1"
 
 [[Latexify]]
 deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
-git-tree-sha1 = "3a0084cec7bf157edcb45a67fac0647f88fe5eaf"
+git-tree-sha1 = "537df0b32bc0a99620872cc6d06278e454da4533"
 uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
-version = "0.14.7"
+version = "0.14.10"
 
 [[LazyArtifacts]]
 deps = ["Artifacts", "Pkg"]
@@ -562,15 +561,15 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
 [[LoweredCodeUtils]]
 deps = ["JuliaInterpreter"]
-git-tree-sha1 = "f008f15264cc11de6de8cbdda3d4712dd152f0c3"
+git-tree-sha1 = "8c96709706ce27471655247ad9a931447d16dd62"
 uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
-version = "1.2.7"
+version = "1.2.9"
 
 [[Luxor]]
 deps = ["Base64", "Cairo", "Colors", "Dates", "FileIO", "ImageMagick", "Juno", "QuartzImageIO", "Random", "Rsvg"]
-git-tree-sha1 = "c018cde3ab7b0867d9f9bd49d290cf15dc9fb926"
+git-tree-sha1 = "8796dfbdde81ec4a811a2d972a6d5f8a7c4113d0"
 uuid = "ae8d54c2-7ccd-5906-9d76-62fc9837b5bc"
-version = "2.9.0"
+version = "2.10.0"
 
 [[MKL_jll]]
 deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
@@ -655,9 +654,9 @@ uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 
 [[OffsetArrays]]
 deps = ["Adapt"]
-git-tree-sha1 = "986e7c0a0ef1863be969d191957ac32cb17d0d79"
+git-tree-sha1 = "b3dfef5f2be7d7eb0e782ba9146a5271ee426e90"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.6.0"
+version = "1.6.2"
 
 [[Ogg_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -696,9 +695,9 @@ version = "8.42.0+4"
 
 [[PNGFiles]]
 deps = ["Base64", "CEnum", "ImageCore", "IndirectArrays", "OffsetArrays", "libpng_jll"]
-git-tree-sha1 = "aa6e87a2361c2fe5a63b1a6a4b567f13aa108991"
+git-tree-sha1 = "28f9667d6ca9168b7198a481847f0879e6c72ef1"
 uuid = "f57f5aa1-a3ce-4bc8-8ab9-96f992907883"
-version = "0.3.5"
+version = "0.3.6"
 
 [[PaddedViews]]
 deps = ["OffsetArrays"]
@@ -714,9 +713,9 @@ version = "1.42.4+10"
 
 [[Parsers]]
 deps = ["Dates"]
-git-tree-sha1 = "50c9a9ed8c714945e01cd53a21007ed3865ed714"
+git-tree-sha1 = "223a825cccef2228f3fdbf2ecc7ca93363059073"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "1.0.15"
+version = "1.0.16"
 
 [[Pixman_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -728,6 +727,12 @@ version = "0.40.0+0"
 deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
+[[PkgVersion]]
+deps = ["Pkg"]
+git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
+uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
+version = "0.1.1"
+
 [[PlotThemes]]
 deps = ["PlotUtils", "Requires", "Statistics"]
 git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
@@ -765,6 +770,12 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 deps = ["Printf"]
 uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 
+[[ProgressMeter]]
+deps = ["Distributed", "Printf"]
+git-tree-sha1 = "6e9c89cba09f6ef134b00e10625590746ba1e036"
+uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
+version = "1.5.0"
+
 [[Qt_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
 git-tree-sha1 = "7760cfea90bec61814e31dfb204fa4b81bba7b57"
@@ -803,15 +814,15 @@ version = "1.0.0"
 
 [[Requires]]
 deps = ["UUIDs"]
-git-tree-sha1 = "cfbac6c1ed70c002ec6361e7fd334f02820d6419"
+git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
 uuid = "ae029012-a4dd-5104-9daa-d747884805df"
-version = "1.1.2"
+version = "1.1.3"
 
 [[Revise]]
 deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
-git-tree-sha1 = "4697c220e6448fb5e9c55404ae1a13271d533b52"
+git-tree-sha1 = "b72fa706920b1421d581525de9f4e442b95ba254"
 uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
-version = "3.1.12"
+version = "3.1.14"
 
 [[Rsvg]]
 deps = ["Cairo", "Glib_jll", "Librsvg_jll"]
@@ -849,9 +860,9 @@ version = "0.3.2"
 
 [[SnoopCompile]]
 deps = ["Cthulhu", "FlameGraphs", "InteractiveUtils", "OrderedCollections", "Pkg", "Printf", "Profile", "Requires", "Serialization", "SnoopCompileCore", "YAML"]
-git-tree-sha1 = "7d5529fa0592dde81f89850696db734345ca6807"
+git-tree-sha1 = "ad31dfca1bbe91d2adaaf319e218ffee7a8eec9f"
 uuid = "aa65fe97-06da-5843-b5b1-d5d13cad87d2"
-version = "2.5.2"
+version = "2.6.0"
 
 [[SnoopCompileCore]]
 deps = ["Serialization"]
@@ -924,6 +935,12 @@ uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[[TiffImages]]
+deps = ["ColorTypes", "DocStringExtensions", "FileIO", "FixedPointNumbers", "IndirectArrays", "OffsetArrays", "OrderedCollections", "PkgVersion", "ProgressMeter"]
+git-tree-sha1 = "94d7c4d760d552a90a98a7150c4a32c7be25e4ee"
+uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69"
+version = "0.3.0"
+
 [[TimeZones]]
 deps = ["Dates", "EzXML", "Mocking", "Pkg", "Printf", "RecipesBase", "Serialization", "Unicode"]
 git-tree-sha1 = "4ba8a9579a243400db412b50300cd61d7447e583"
diff --git a/test/test_23.jl b/test/test_23.jl
index 4ec765c..cf00572 100644
--- a/test/test_23.jl
+++ b/test/test_23.jl
@@ -75,9 +75,7 @@ settext(str, O + (-WI/2 + EM, 1.5EM), markup = true)
 
 # Plot the complex-valued function, aka velocity vectors.
 B = begin
-        unclamped = ∇_rectangle(ϕ,
-            physwidth = physwidth,
-            physheight = physheight);
+        unclamped = ∇_rectangle(ϕ; physwidth, physheight)
         map(unclamped) do u
             hypot(u) > 0.5m/s ? NaN∙u : u
         end
diff --git a/test/test_26.jl b/test/test_26.jl
index c6fe004..c0e02a4 100644
--- a/test/test_26.jl
+++ b/test/test_26.jl
@@ -68,7 +68,7 @@ draw_legend(lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23), legend)
 @layer begin
     sethue(PALETTE[1])
     setline(8)
-    draw_streamlines(OB_26, xs, ys, fxy, h)
+    draw_streamlines(OB_26, fxy; h)
     setline(1)
 end
 
diff --git a/test/test_26.png b/test/test_26.png
index 77fd91b..c82610e 100644
Binary files a/test/test_26.png and b/test/test_26.png differ
diff --git a/test/test_27.jl b/test/test_27.jl
index 495a8df..dc277a5 100644
--- a/test/test_27.jl
+++ b/test/test_27.jl
@@ -79,7 +79,7 @@ draw_legend(lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23), botlegend)
 @layer begin
     sethue(PALETTE[1])
     setline(8)
-    draw_streamlines(OB, xs, ys, fxy, h)
+    draw_streamlines(OB, fxy)
     setline(1)
 end
 
@@ -94,4 +94,5 @@ settext(str, O + (-WI/2 + EM, 2 * EM) , markup = true)
 
 finish()
 set_scale_sketch()
+nothing
 end
diff --git a/test/test_27.png b/test/test_27.png
index 48916eb..a0a5b2e 100644
Binary files a/test/test_27.png and b/test/test_27.png differ
diff --git a/test/test_3.jl b/test/test_3.jl
index 986a36c..b1c6353 100644
--- a/test/test_3.jl
+++ b/test/test_3.jl
@@ -3,25 +3,29 @@ import MechanicalSketch: Point, rect, arrow, line
 import MechanicalSketch: EM, TXTOFFV, O, PT, Point, WI, HE
 import MechanicalSketch: RGB, HSL, color_with_lumin, background, sethue, text, fontsize
 import MechanicalSketch: settext, setfont, setopacity, FS, color_with_lumin, color_with_lumin2
-import MechanicalSketch: empty_figure, PALETTE
+import MechanicalSketch: empty_figure, PALETTE, finish, paint, ∙
+import MechanicalSketch: transform, @layer, do_action, getmatrix, grestore, gsave, scale
+import MechanicalSketch: box_fill_outline, label, row_height, pixelwidth, @import_expand
+import MechanicalSketch: scale_pt_to_unit, lumin, label_boxed
+if !@isdefined °
+    @import_expand ~m # Will error if m² already is in the namespace
+    @import_expand s
+    @import_expand °
+    @import_expand rad
+end
+#let
+    typcol = PALETTE[9];
+    bckcol = color_with_lumin(typcol, 0.5 * lumin(typcol))
+    empty_figure(joinpath(@__DIR__, "test_3.png");
+        backgroundcolor = bckcol, hue = "white")
 
-let
-    this_fig = empty_figure(joinpath(@__DIR__, "test_3.png"))
-    typcol = MechanicalSketch.HSL(get(PALETTE, 0.9))
-    bckcol = MechanicalSketch.RGB(MechanicalSketch.HSL(typcol.h, typcol.s, typcol.l * 0.5))
-    background(bckcol)
-    sethue("white")
-    arrow(O, Point(0, -5EM) )
-    arrow(O, Point(10EM, -2EM), arrowheadlength=EM, arrowheadangle = pi/12, linewidth = PT)
-    arrow(O, Point(8EM, 0.5EM), arrowheadlength=EM, arrowheadangle = pi/12, linewidth = PT)
+    # Straight arrow width
+    pt = O + ( -6EM, 2EM)
+    arrow(pt, pt + (0, -5EM) )
+    arrow(pt, pt + (10EM, -2EM), arrowheadlength=EM, arrowheadangle = pi/12, linewidth = PT)
+    arrow(pt, pt + (8EM, 0.5EM), arrowheadlength=EM, arrowheadangle = pi/12, linewidth = PT)
 
-    # center, radius, start (cw), end(cw)
-    setopacity(0.2)
-    for (i, endang) in enumerate(-π:π/4:π)
-    sethue(PALETTE[i])
-    arrow(O, 8EM, π, endang, arrowheadlength = 2EM,   arrowheadangle=pi/12, linewidth = i *PT)
-    end
-    setopacity(1)
+    # Palette
     ncol = length(PALETTE)
     dw =  WI / ncol
     for i = 1:ncol
@@ -37,14 +41,14 @@ let
                     ps + (0, EM - TXTOFFV))
         fontsize(FS/2)
         setfont("Calibri", FS / 2)
-
-        text("lumin", Point(xs, 0))
-        text("lumin2", Point(xs + dw /2, 0))
+        posy = HE / 2 - 22 * 0.5 * EM
+        text("lumin", Point(xs, posy))
+        text("lumin2", Point(xs + dw /2, posy))
         for (j, lum) in enumerate(2:-0.1:0)
             curcol = PALETTE[i]
             # Left rectancle, simple luminance
             sethue(color_with_lumin(curcol, lum * 100))
-            posy = j * 0.5 * EM
+            posy += 0.5 * EM
             ps = Point(xs, posy)
             rect( ps , dw / 2 , EM / 2 , :fill)
             str = string(round(lum, digits = 2))
@@ -63,5 +67,64 @@ let
                         ps + (dw/2 + dw/ 4, -TXTOFFV + 0.6 *EM) , markup=true)
         end
     end
-    MechanicalSketch.finish()
-end
\ No newline at end of file
+    # Arc arrows, left screen
+    sethue(PALETTE[3])
+    y = -HE / 4 + 0EM
+    x = - WI / 3
+    cpl = O + (x, y)
+    text("Unitless angles: Always clockwise", cpl + (0, EM), halign = :center)
+    startangle = 0.0
+    iter  = enumerate(range(-π , π, step = π / 4))
+    for (i, endangle) in iter
+        # center, radius, start (cw), end(cw)
+        radius = (3 + i / 4 )EM
+        arrow(cpl, radius, startangle, endangle;
+               arrowheadlength = EM, arrowheadangle = pi / 12, linewidth = 2PT)
+    end
+    for (i, endangle) in iter
+        radius = (3 + i / 4 )EM
+        str = "$(round(startangle / π; digits = 2)) to $(round(endangle / π; digits = 2))∙π"
+        endpoint = cpl + radius .* (cos(endangle), sin(endangle))
+        label_boxed(endpoint, str)
+    end
+
+    # Arc arrows with units, mid-screen
+    cpm = O + (0, y)
+    text("Uniful angles, specified clockwise", cpm + (0, EM), halign = :center)
+    startangle = 0.0∙rad
+    iter  = enumerate(range(-π, π, step = π / 4)∙rad)
+    rdeg(a) = round(°, a; digits = 0) |> typeof(1°)
+    for (i, endangle) in iter
+        radius = (3 + i / 4 )EM * scale_pt_to_unit(m)
+        arrow(cpm, radius, startangle, endangle;
+              arrowheadlength = EM, arrowheadangle = pi / 12, linewidth = 2PT, clockwise = true)
+    end
+    for (i, endangle) in iter
+        radius = (3 + i / 4 )EM
+        str = "$(rdeg(startangle)) to $(rdeg(endangle))"
+        endpoint = cpm + radius .* (cos(endangle), -sin(endangle))
+        label_boxed(endpoint, str)
+    end
+
+    # Arc arrows with units, right-screen
+    cpr = O + (-x, y)
+    text("Uniful angles, default counterclockwise", cpr + (0, EM), halign = :center)
+    startangle = 0.0∙rad
+    iter  = enumerate(range(-π, π, step = π / 4)∙rad)
+    rdeg(a) = round(°, a; digits = 0) |> typeof(1°)
+    for (i, endangle) in iter
+        radius = (3 + i / 4 )EM * scale_pt_to_unit(m)
+        arrow(cpr, radius, startangle, endangle;
+              arrowheadlength = EM, arrowheadangle = pi / 12, linewidth = 2PT)
+    end
+    for (i, endangle) in iter
+        radius = (3 + i / 4 )EM
+        str = "$(rdeg(startangle)) to $(rdeg(endangle))"
+        endpoint = cpr + radius .* (cos(endangle), -sin(endangle))
+        label_boxed(endpoint, str)
+    end
+
+
+    finish()
+
+#end
\ No newline at end of file
diff --git a/test/test_3.png b/test/test_3.png
index 76e1ed8..ec03b9f 100644
Binary files a/test/test_3.png and b/test/test_3.png differ
diff --git a/test/test_32.jl b/test/test_32.jl
index 90a2f58..6ec87e9 100644
--- a/test/test_32.jl
+++ b/test/test_32.jl
@@ -31,7 +31,7 @@ totwidth = totheight * WI / HE
 Δy = physheight * 1.1
 
 # Reused velocity field from earlier tests
-velocity_matrix = clamped_velocity_matrix(ϕ_32; physwidth = physwidth, physheight = physheight, cutoff = 0.5m/s);
+velocity_matrix = clamped_velocity_matrix(ϕ_32; physwidth, physheight, cutoff = 0.5m/s);
 # One complex matrix: Phase and amplitude for the visualization. This can generate cyclic movies
 complex_convolution_matrix = convolute_image_32(velocity_matrix) # 15.071 s (364518 allocations: 61.17 MiB)
 # The distribution of phase angles (complex argument) ought to be flat between -π and π. Let's check that:
@@ -43,16 +43,17 @@ histogrampoints = Point.(collect(phasetop * n_pixels / 2π), -EM .* relfreq);
 # We'll also prepare a uniform flow field, constant velocity 0.5 m/s.
 fxy_unif(x, y) = (0.5, 0.0)m∙s⁻¹
 # Phase and amplitude for the visualization. This can generate cyclic movies
-complex_convolution_matrix_uniform = convolute_image_32(fxy_unif, physwidth = physwidth, physheight = physheight)
+complex_convolution_matrix_uniform = convolute_image_32(fxy_unif; physwidth, physheight)
 
 # We'll also plot a vertically uniform flow field, horizontally increasing velocity from 0 to 0.5 m/s.
 fxy_lin(x, y) = (0.5 * (x / physwidth + 0.5), 0.0)m∙s⁻¹
 # Phase and amplitude for the visualization. This can generate cyclic movies
-complex_convolution_matrix_linear = convolute_image_32(fxy_lin, physwidth = physwidth, physheight = physheight)
+complex_convolution_matrix_linear = convolute_image_32(fxy_lin; physwidth, physheight)
 
 legend = BinLegend(;maxlegend = 0.2, minlegend = -1.0, noofbins = 256, 
                        colorscheme = reverse(Greys_9), 
                        nan_color = color_with_lumin(PALETTE[1], 80), name = Symbol("Value{Float64}"))
+                       
 # Define scene functions (parts of each image)
 
 # Rectangular flow field plot including a visual frame counter
diff --git a/test/test_33.jl b/test/test_33.jl
index 88761b3..acb4d0b 100644
--- a/test/test_33.jl
+++ b/test/test_33.jl
@@ -7,7 +7,7 @@ import MechanicalSketch: lenient_max
 import MechanicalSketch: place_image, @layer
 import MechanicalSketch: lic_matrix_current
 import MechanicalSketch: ColorSchemes, color_with_lumin, convolution_matrix, clamped_velocity_matrix
-import MechanicalSketch: BinLegend, BinLegendVector, streamlines_add!, draw_legend
+import MechanicalSketch: BinLegend, BinLegendVector, streamlines_matrix, draw_legend
 import MechanicalSketch: leonardo, fontsize
 import ColorSchemes:     Paired_6, Greys_9, isoluminant_cgo_80_c38_n256
 import Base.show
@@ -38,12 +38,11 @@ Oadj = O + (0, EM /3)
 
 # Reused velocity field from earlier tests
 max_velocity = 0.5m/s
-velocity_matrix = clamped_velocity_matrix(ϕ_33; physwidth = physwidth, physheight = physheight, cutoff = max_velocity);
+velocity_matrix = clamped_velocity_matrix(ϕ_33; physwidth, physheight, cutoff = max_velocity);
 
 v_xy = matrix_to_function(velocity_matrix)
 
-streamlinepixels = falses(size(velocity_matrix))
-streamlines_add!(v_xy, streamlinepixels)
+streamlinepixels = streamlines_matrix(v_xy; physwidth, physheight)
 
 cent1 = Oadj + (-Δx, Δy)
 legend1 = BinLegend(;maxlegend = 1.0, minlegend = 0.0, noofbins = 2, 
@@ -115,7 +114,6 @@ settext(str, ulp, markup = true)
 
 
 cent6 = Oadj + (Δx, -Δy)
-noofbins = 256
 legend6 = BinLegend(;maxlegend = max_velocity, binwidth,
         colorscheme = Paired_6, name = :Speed)
 mixmat = map(legend6.(speedmatrix), curmat) do col, lu
diff --git a/test/test_33.png b/test/test_33.png
index 6ae8d6b..0bc7639 100644
Binary files a/test/test_33.png and b/test/test_33.png differ
diff --git a/test/test_34.jl b/test/test_34.jl
index ec948bb..422ba09 100644
--- a/test/test_34.jl
+++ b/test/test_34.jl
@@ -42,7 +42,7 @@ cent5 = Oadj + (-Δx, -Δy)
 
 # Reused velocity field from earlier tests
 max_velocity = 0.5m/s
-velocity_matrix = clamped_velocity_matrix(ϕ_33; physwidth = physwidth, physheight = physheight, cutoff = max_velocity);
+velocity_matrix = clamped_velocity_matrix(ϕ_33; physwidth, physheight, cutoff = max_velocity);
 speedmatrix = hypot.(velocity_matrix)
 c1 = convolution_matrix(velocity_matrix)
 lum = 50
diff --git a/test/test_38.jl b/test/test_38.jl
new file mode 100644
index 0000000..221cae3
--- /dev/null
+++ b/test/test_38.jl
@@ -0,0 +1,72 @@
+import MechanicalSketch: @import_expand, empty_figure, WI, HE, EM, O, PT, finish
+import MechanicalSketch: settext, place_image, PALETTE, color_with_lumin
+import MechanicalSketch: circle, ∙, arrow_nofill, Point, @layer, sethue
+import MechanicalSketch: latexify, arrow, generate_complex_potential_vortex
+import MechanicalSketch: generate_complex_potential_source, ∇, color_matrix
+import MechanicalSketch: fontsize, draw_legend, draw_streamlines, ∇_rectangle
+import ColorSchemes:     isoluminant_cgo_80_c38_n256
+
+if !@isdefined N
+    @import_expand ~m # Will error if m² already is in the namespace
+    @import_expand s
+    @import_expand N
+    @import_expand °
+end
+include("test_functions_33.jl")
+ϕ = ϕ_33
+include("test_functions_39.jl")
+
+
+# Scaling and placement
+physwidth = 10.0m
+physheight = 4.0m
+framedheight = physheight * 1.2
+totheight = framedheight * 4
+totwidth = totheight * WI / HE
+Δx = totwidth / 4
+Δy = framedheight / 2
+Oadj = O + (-3EM, EM /3)
+cent1 = Oadj + ( Δx, 3Δy)
+cent2 = Oadj + ( Δx, Δy)
+cent3 = Oadj + ( Δx, -Δy)
+cent4 = Oadj + ( Δx, -3Δy)
+
+empty_figure(joinpath(@__DIR__, "test_38.png"); 
+      backgroundcolor = color_with_lumin(PALETTE[5], 90),
+      hue = "black", height = totheight);
+
+pt = O + (-WI / 2 + EM, -HE / 2 + 2EM)
+settext("Potential field, defined previously:", pt)
+expression = :(ϕ(𝐫));
+ptul, ptbr, latexscale = place_image(pt + (EM, 0), latexify(expression); height = 1.2EM, centered = false)
+pt += (0, 1.5EM + ptbr.y - ptul.y)
+
+settext("Steady flow velocity field at <b>r</b>:", pt, markup = true)
+expression = :( 𝐕(𝐫) = 𝐢∙u(𝐫) + 𝐣∙v(𝐫) = ∇ϕ(𝐫))
+ptul, ptbr = place_image(pt + (EM, 0), latexify(expression); scalefactor = latexscale, centered = false)
+
+𝐕 = ∇(ϕ)
+colmat, legend1 = color_matrix(𝐕; centered = true, maxlegend = 1.0m/s)
+ulp, _ = place_image(cent1, colmat)
+
+@layer begin
+    fontsize(25)
+    draw_legend(ulp + (0.5EM, 0.5EM) + (physwidth, 0.0m), legend1)
+end
+# To better indicate positive direction for this static image:
+@layer begin
+    sethue(PALETTE[8])
+    draw_streamlines(cent1, 𝐕, probability = 0.001)
+end
+
+
+pt += (0, 1.5EM + ptbr.y - ptul.y)
+settext("Local acceleration, steady flow:", pt, markup = true)
+# Strictly incorrect, use paranthesis or .
+expression = :(𝐚 = d𝐕 / dt = X) # 𝐕∙(∇𝐕) = 𝐕.∙∇∇ϕ = 𝐕.∙ ∇²ϕ) 
+ptul, ptbr, _ = place_image(pt + (EM, 0), latexify(expression); scalefactor = latexscale, centered = false)
+pt += (0, 1.5EM + ptbr.y - ptul.y)
+#𝐚 = ∇_rectangle(𝐕)
+
+
+finish()
diff --git a/test/test_38.png b/test/test_38.png
new file mode 100644
index 0000000..e11ac16
Binary files /dev/null and b/test/test_38.png differ
diff --git a/test/test_39.jl b/test/test_39.jl
new file mode 100644
index 0000000..da47bca
--- /dev/null
+++ b/test/test_39.jl
@@ -0,0 +1,35 @@
+import MechanicalSketch: @import_expand, empty_figure, WI, HE, EM, O, PT, finish
+import MechanicalSketch: settext, place_image, PALETTE, color_with_lumin
+import MechanicalSketch: circle, ∙, arrow_nofill, Point, @layer, sethue
+import MechanicalSketch: latexify, arrow
+if !@isdefined N
+    @import_expand ~m # Will error if m² already is in the namespace
+    @import_expand s
+    @import_expand N
+    @import_expand °
+end
+include("test_functions_39.jl")
+empty_figure(joinpath(@__DIR__, "test_39.png"); 
+      backgroundcolor = color_with_lumin(PALETTE[5], 90),
+      hue = "black", height = 10m);
+pt = O + (-WI / 2 + EM, -HE / 2 + 2EM)
+settext("<b>Flow around a cylinder</b>", pt, markup = true)
+pt += (0, EM)
+pt += (1m, -1m)
+r = 1.0m
+@layer begin
+    sethue(PALETTE[8])
+    circle(pt; r, action=:fill)
+end
+vr =  √2 / 2 ∙(r, r)
+arrow_nofill(pt + 0.5vr, pt + vr;  arrowheadlength = 0.5EM)
+settext("$r", pt + vr)
+arrow(pt, 2EM, 0°, 30°; arrowheadlength = 0.4EM) 
+pt += (0, 1.5EM)
+settext("θ", pt + vr)
+# 
+#expression = :( δ / δr ∙ (δψ / δθ)  + δ / δθ ∙ (-δψ / δr)  = 0)
+
+
+
+finish()
\ No newline at end of file
diff --git a/test/test_39.png b/test/test_39.png
new file mode 100644
index 0000000..54e3727
Binary files /dev/null and b/test/test_39.png differ
diff --git a/test/test_5.jl b/test/test_5.jl
index 3b92318..8c7d4ae 100644
--- a/test/test_5.jl
+++ b/test/test_5.jl
@@ -5,7 +5,7 @@ import MechanicalSketch: empty_figure, PALETTE, dimension_aligned
 
 let
 empty_figure(joinpath(@__DIR__, "test_5.png"))
-background(color_with_lumin(PALETTE[6], 30))
+background(color_with_lumin(PALETTE[6], 40))
 @layer begin
   sethue("black")
   setopacity(0.5)
diff --git a/test/test_5.png b/test/test_5.png
index e58f9bc..af03eb3 100644
Binary files a/test/test_5.png and b/test/test_5.png differ
diff --git a/test/test_functions_38.jl b/test/test_functions_38.jl
new file mode 100644
index 0000000..391f66b
--- /dev/null
+++ b/test/test_functions_38.jl
@@ -0,0 +1 @@
+println("Wow, such empty")
\ No newline at end of file
diff --git a/test/test_functions_39.jl b/test/test_functions_39.jl
new file mode 100644
index 0000000..391f66b
--- /dev/null
+++ b/test/test_functions_39.jl
@@ -0,0 +1 @@
+println("Wow, such empty")
\ No newline at end of file