Flow visualisation experiments, parse quantities from file, variable …

…discipline, table improvement. modified: Manifest.toml modified: Project.toml modified: README.md new file: resource/Reaction forces.txt modified: src/MechanicalSketch.jl modified: src/arrow.jl modified: src/autodiff_unitfu.jl modified: src/curves.jl modified: src/flow.jl modified: src/rope.jl modified: src/scale.jl new file: src/streamline_convolution.jl modified: src/table.jl new file: test/clean_util.jl modified: test/runtests.jl modified: test/test_1.jl modified: test/test_10.jl modified: test/test_10.png modified: test/test_11.jl modified: test/test_11.png modified: test/test_12.jl modified: test/test_12.png modified: test/test_13.jl modified: test/test_13.png modified: test/test_14.jl modified: test/test_14.png modified: test/test_15.jl modified: test/test_16.jl modified: test/test_17.jl modified: test/test_18.jl modified: test/test_19.jl modified: test/test_2.jl modified: test/test_20.jl modified: test/test_20.png modified: test/test_21.jl modified: test/test_21.png modified: test/test_22.jl modified: test/test_22.png modified: test/test_23.jl modified: test/test_23.png modified: test/test_24.jl new file: test/test_24.png new file: test/test_25.jl new file: test/test_25.png new file: test/test_26.jl new file: test/test_26.png new file: test/test_27.jl new file: test/test_27.png new file: test/test_28.jl new file: test/test_28.png modified: test/test_3.jl modified: test/test_4.jl modified: test/test_5.jl modified: test/test_5.png modified: test/test_6.jl modified: test/test_6.png modified: test/test_7.jl modified: test/test_7.png modified: test/test_8.jl modified: test/test_8.png modified: test/test_9.jl modified: test/test_9.png new file: test/test_functions_24.jl new file: test/test_functions_25.jl new file: test/test_functions_26.jl new file: test/test_functions_27.jl new file: test/test_functions_28.jl new file: test/test_functions_4.jl
hustf · Nov 15, 2020 · 19b0c0a · 19b0c0a
1 parent d43eeb2
commit 19b0c0a
Show file tree

Hide file tree

Showing 68 changed files with 2,459 additions and 455 deletions.
diff --git a/Manifest.toml b/Manifest.toml
diff --git a/Project.toml b/Project.toml
@@ -20,10 +20,10 @@ REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 Unitfu = "5ee08b94-2369-4f4a-b8c7-99333ba35fb0"
 
 [compat]
-MechanicalUnits = "0.3.1"
 julia = "1.3.0"
 
 [extras]

diff --git a/README.md b/README.md
@@ -3,12 +3,10 @@ This package is for making sketches using Luxor and MechanicalUnits.
 
 The sketch format is intended for rougly 1/3 of an A4 page, a typical figure for a technical report.
 
-The quantity dimensions (length, velocity, power, force, etc.) combines well with multiple dispatch. 
-
-An example is the 'arrow' function,  which represents 2d vectors.
+The quantity dimensions (length, velocity, power, force, etc.) combines well with multiple dispatch. An example is the 'arrow' function,  which shows 2d vectors with visual hints to the type of quantity.
 
 The package is developed progressively by writing scripts in the test folder, and putting the most useful pieces of script into functions.
 
-There's no intention to make this very general, and we would probably be better off reading all of Luxor's documentation.
+There's no intention to make this very general, rather to add functionality as the need arise.
 
-Sketches should look good and be in a consistent pallette, though. The test images are in .png format, but .svg is better for zooming.
+Sketches should look good and be in a consistent pallette, though. The test images are in .png format, but .svg is better for zooming. For post-editing figures, not all svg editors handle fonts. 'Inkscape' seems to work well. You might need to install additional fonts for some test images.
diff --git a/resource/Reaction forces.txt b/resource/Reaction forces.txt
diff --git a/src/MechanicalSketch.jl b/src/MechanicalSketch.jl
@@ -32,7 +32,8 @@ intersectboundingboxes, boundingboxesintersect, pointcrossesboundingbox,
 
 boxmap,
 
-circle, circlepath, ellipse, hypotrochoid, epitrochoid, squircle, center3pts, curve,
+circle, circlepath, currentdrawing, ellipse, hypotrochoid, epitrochoid,
+squircle, center3pts, curve,
 arc, carc, arc2r, carc2r, isarcclockwise, arc2sagitta, carc2sagitta,
 spiral, sector, intersection2circles,
 intersection_line_circle, intersectionlinecircle, intersectioncirclecircle, ispointonline,
@@ -118,15 +119,14 @@ Colors, parse
 
 using MechanicalUnits
 import MechanicalUnits.Unitfu: Power, oneunit, numtype
-@import_expand kW # Don't use ~ since that would also import W, which we use elsewhere.
+@import_expand kW # Don't use ~ since that would also import WI, which we use elsewhere.
 
-import LinearAlgebra: norm
-# TODO: using, not import.
-import REPL.TerminalMenus
+# TODO: Cleanup
+import REPL.TerminalMenus ## Where used?
 import ColorSchemes
 import ColorSchemes: getinverse, get
 import ColorSchemes: HSL, HSLA, RGB, RGBA, HSV, LCHuv, LCHuvA
-import Base: -, +, *, /, abs
+import Base: -, +, *, /, hypot
 import Colors: @colorant_str
 import FileIO: @format_str, File, save
 import ForwardDiff
@@ -168,23 +168,28 @@ For line thickness, we are used to points.
 """
 const PT = Int(round(EM / 12))
 """
-We want to make high quality figures for A4 with 5 cm total margin. Width and height are
+WI is the pixel widt of the figure. See HE and origo, O.
+
+We want to make quality figures for A4 with 5 cm total margin. Width and height are
     a4_w_72dpi, _ = Luxor.paper_sizes["A4"]
     marginfrac = 50 / 210
     w_300f = a4_w_72dpi * (1 - marginfrac) * 300 / 72
-    W = Int(round(w_300f))
-    H = Int(round(w_300f * 2 / 3))
-"""
-const W = 1889
-const H = 1259
-global SCALEDIST = 20m / H
-global SCALEVELOCITY = 70m/s / H
-global SCALEFORCE = 20kN / H
+    WI = Int(round(w_300f))
+    HE = Int(round(w_300f * 2 / 3))
 
-
-# Note, this could be extended by modulus, such that points to the left of the screen
-# would appear at right. But that could make some strange bugs too.
-# CONSIDER generalize to quantity, let scale(x) do the transformation?
+"""
+const WI = 1889
+"HE is the pixel height of the figure. See WI and orgo, O"
+const HE = 1259
+"Used for scaling quantities to pixels. Use setscale_dist to change"
+global SCALEDIST = 20m / HE
+"Used for scaling quantities to pixels. Use setscale_velocity to change"
+global SCALEVELOCITY = 70m/s / HE
+"Used for scaling quantities to pixels. Use setscale_force to change"
+global SCALEFORCE = 20kN / HE
+
+
+# CONSIDER generalize to quantity, let scale_sketch(x) do the transformation?
 Point(x::T, y::T) where T<:Length = Point(x / SCALEDIST, y / scaledisty())
 Point(x::T, y::T) where T<:Velocity = Point(x / SCALEVELOCITY , y / scalevelocityy())
 Point(x::T, y::T) where T<:Force = Point(x / SCALEFORCE , y / scaleforcey())
@@ -221,7 +226,9 @@ const ForceTuple = NTuple{2, <: Force}
 -(p1::Point, shift::Tuple{Real, Quantity}) = p1 - Point(shift[1], shift[2])
 *(p1::Point, shift::Tuple{Real, Quantity}) = p1 * Point(shift[1], shift[2])
 /(p1::Point, shift::Tuple{Real, Quantity}) = p1 / Point(shift[1], shift[2])
-
+# extend function in base with Point
+hypot(p::Point) = hypot(p.x, p.y)
+hypot(p::QuantityTuple) = hypot(p[1], p[2])
 """
 Rotations given with units occur around the positive z axis, when y is up and x is to the right.
 We don't dispatch on the dimension of angles, because the dimension is NonDims
@@ -231,34 +238,36 @@ polyrotate!(f, ang::Angle) = polyrotate!(f, - ustrip( ang |> rad))
 
 
 """
-    empty_figure(filename = "HiThere.png")
+    empty_figure(filename = "HiThere.png";
+    backgroundcolor = color_with_luminance(PALETTE[8], 0.1),
+    hue = PALETTE[8] )
 
 Establish a drawing sized for A4 300 dpi figures,
 black on white figure, line width 3 pt  default.
 """
-function empty_figure(filename = "HiThere.png")
-    fig = Drawing(W, H, filename)
+function empty_figure(filename = "HiThere.png";
+        backgroundcolor = color_with_luminance(PALETTE[8], 0.1),
+        hue = PALETTE[8] )
+    fig = Drawing(WI, HE, filename)
     # Font for the 'toy' text interface
     # (use e.g. JuliaMono for unicode symbols like ∈. There are no nice fonts for text AND math.
     fontface("Calibri")
     # 1 pt font size = 12/72 inch - by the book.
     # Letter spacing works differently here than in Word, so we adjust a little.
     fontsize(FS)
     setfont("Calibri", FS)
-    background("white")
-    sethue("black")
+    background(backgroundcolor)
+    sethue(hue)
     setline(0.5PT)
     setdash("solid")
     # Origo at centre
     origin()
     # Scale and rotation - x right, y up ('z' is in.). And rotations may
     # clockwise. Just deal with it. Or stick to using dimensions.
-
-    setmatrix([1, 0, 0, 1, W / 2, H / 2])
+    setmatrix([1, 0, 0, 1, WI / 2, HE / 2])
     fig
 end
 
-
 """
     text(t, pt::Point, angle::T; kwargs) where {T <: Angle}
 For angles with unit, use rotation around z axis.
@@ -277,5 +286,5 @@ include("table.jl")
 include("curves.jl")
 include("flow.jl")
 include("autodiff_unitfu.jl")
-
+include("streamline_convolution.jl")
 end # module
diff --git a/src/arrow.jl b/src/arrow.jl
@@ -4,22 +4,32 @@
 Draw a single force arrow with filled head.
 
     p        starting position.
-    v        scalar force.
+    f        scalar force quantity.
     α = 0°   is direction to the right, positive values rotate around z axis.
 
 """
 arrow(p::Point, f::F; kwargs...) where F <: Force = arrow(p, f, 0 * f; kwargs...)
 
 
 """
-    arrow(p::Point, F::ForceTuple; kwargs...)
+    arrow(p::Point, f::ForceTuple; kwargs...)
 
 Draw a single velocity arrow - no filled arrow head.
 
 p     starting position (at the normal between vanes tip and arrow spine)
-F     is a 2-tuple here, but you can optionally rotate the x-axis by keyword argument α.
+F     is a force tuple referred to an x-axis which can optionally be rotated around z by
+      keyword argument α.
+components = false: No components forces on the x and y axes are drawn.
 """
-arrow(p::Point, f::ForceTuple; kwargs...) = arrow(p, f[1], f[2]; kwargs...)
+function arrow(p::Point, f::ForceTuple; components = true, kwargs...)
+    if !components
+        αcalc = atan(f[2], f[1]) |> rad |> °
+        F = hypot(f)
+        arrow(p, F, 0.0 * f[2]; kwargs..., α = αcalc)
+    else
+        arrow(p, f[1], f[2]; kwargs...)
+    end
+end
 
 """
     arrow(p::Point, Fx::F, Fy::F;
@@ -52,29 +62,36 @@ function arrow(p::Point, Fx::F, Fy::F;
     avglumin = (luminback + luminfront) / 2
     curcol = get_current_RGB()
     avgcol = color_with_luminance(curcol, avglumin)
-    # Draw the components, both only if the resultant is not coincident with one of them.
-    if norm(ΔFx) > 0 && norm(ΔFy) > 0
+    # Draw the components, but only if the resultant is not coincident with one of them.
+    if hypot(ΔFx) > 0 && hypot(ΔFy) > 0
         gsave()
         setdash("longdashed")
         deltalumin = luminback - luminfront
         lesscontrast = color_with_luminance(get_current_RGB(), luminfront + deltalumin)
         sethue(lesscontrast)
-        shaftlength = norm(ΔFx)
+        shaftlength = hypot(ΔFx)
         arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3
         arrow_nofill(p, p + ΔFx, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT / 2)
-        shaftlength = norm(ΔFy)
+        shaftlength = hypot(ΔFy)
         arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3
         arrow_nofill(p, p + ΔFy, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT / 2)
         grestore()
     end
 
     # Draw the resultant arrow
     endpoint = p + ΔF
-    shaftlength = norm(ΔF)
+    shaftlength = hypot(ΔF)
     arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3
-    arrow(p, endpoint, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT)
+    if !isapprox(shaftlength, 0.0)
+        arrow(p, endpoint, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT)
+    else
+        # with shaft length approximately zero, draw two opposing arrows with minimum possible shaft lengths
+        tiny = oneunit(shaftlength) # almost zero
+        arrow(p + ( tiny, 0), endpoint, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT)
+        arrow(p + (-tiny, 0), endpoint, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT)
+    end
 
-    # add a label with the euclidean length of the resultant, at the tip of it.
+    # add an optional label with the euclidean length of the resultant, at the tip of it.
     if labellength
         sethue(avgcol)
         fontface("Calibri bold")
@@ -165,7 +182,7 @@ function setlabel(endpoint, vx, vy, α)
     Δvx = (vx * cos(α), vx *sin(α))
     Δvy = (-vy * sin(α), vy *cos(α))
     Δv = Δvx + Δvy
-    l = norm(Δv)
+    l = hypot(Δv)
     rounded = round(typeof(l), l, digits = 1)
     direction = atan(Δv[2], Δv[1]) |> °
     labdir = direction + 90°
@@ -194,7 +211,7 @@ function curved_point_arrow_with_vanes(p, Δv, endpoint, linewidth)
     arrowheadangle = pi/24
 
     isapprox(Δv, Point(0,0)) && throw(error("can't draw velocity arrow between two identical points"))
-    shaftlength = sqrt( Δv.x^2 + Δv.y^2)
+    shaftlength = hypot(Δv)
     arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3
 
     shaftangle = atan(p.y - endpoint.y, p.x - endpoint.x)
@@ -315,4 +332,3 @@ function arrow_nofill(startpoint::Point, endpoint::Point;
     end
     grestore()
 end
-
diff --git a/src/autodiff_unitfu.jl b/src/autodiff_unitfu.jl
@@ -99,6 +99,7 @@ end
 
 """
     ∇_rectangle(CQ_to_Q;
+        cutoff = NaN,
         physwidth = 20.0,
         width_relative_screen = 2.0 / 3,
         height_relative_width = 1.0 / 3)
@@ -111,22 +112,33 @@ where
     Q is a real valued Quantity
 
 Calculates the gradient given the extents of a rectangle and how they relate to screen pixels.
+
+Optional argument cutoff: In case given, the magnitude (vector length) is limited to that
+value without changing vector direction.
 """
 function ∇_rectangle(CQ_to_Q;
+    cutoff = NaN,
     physwidth = 20.0,
     width_relative_screen = 2.0 / 3,
     height_relative_width = 1.0 / 3)
 
     physheight = physwidth * height_relative_width
 
     # Discretize per pixel
-    nx = round(Int64, W * width_relative_screen)
+    nx = round(Int64, WI * width_relative_screen)
     ny = round(Int64, nx * height_relative_width)
     # 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)))
 
     # Matrix of 2-element static vectors, one per pixel
     xyq = [SA[ix * SCALEDIST, -iy * SCALEDIST] for iy = pixitery, ix = pixiterx]
-    ∇_matrix(CQ_to_Q, xyq)
+    if isnan(cutoff)
+        ∇_matrix(CQ_to_Q, xyq)
+    else
+        unclamped = ∇_matrix(CQ_to_Q, xyq)
+        map(unclamped) do u
+            hypot(u) > cutoff ? cutoff * u / hypot(u) : u
+        end
+    end
 end
diff --git a/src/curves.jl b/src/curves.jl
@@ -1,12 +1,19 @@
+"""
+    circle(p::Point, r::Length, action)
 
-circle(p::Point, r::Length, action) = circle(p, scale(r), action)
+    where
+
+action can be :nothing, :fill, :stroke, :fillstroke, :fillpreserve, :strokepreserve, :clip.
+    The default is :nothing.
+"""
+circle(p::Point, r::Length, action) = circle(p, scale_sketch(r), action)
 
 function ellipse(p::Point, w::Q, h::Q, action=:none)  where {Q<:Length}
-    ellipse(p, scale(w), scale(h), action)
+    ellipse(p, scale_sketch(w), scale_sketch(h), action)
 end
 function squircle(center::Point, hradius::Q, vradius::Q, action=:none;
     kwargs...) where {Q<:Length}
-    squircle(center, scale(hradius), scale(vradius),
+    squircle(center, scale_sketch(hradius), scale_sketch(vradius),
     action;
     kwargs...)
 end
@@ -16,4 +23,34 @@ end
 arc, carc, arc2r, carc2r, sector, pie, curve, circlepath,
 hypotrochoid, epitrochoid, spiral, intersection2circles,
 intersectioncirclecircle, circlepointtangent
-=#
+=#
+
+
+
+"""
+    diminishingtrace(center, vx, vy)
+
+center is local origo, a Point, typically equal to global O.
+vx, vy are vectors of corresponding coordinate quantities.
+
+Plots a polyline until NaN or end of vectors.
+The first segment is drawn with opacity 1, decreasing to one step before translucent(zero)
+at the last segment
+"""
+function diminishingtrace(center, vx, vy)
+    gsave()
+    n = length(vx)
+    opacity = 1.0
+    xo, yo  = vx[1], vy[1]
+    for (x, y) in zip(vx[2:end], vy[2:end])
+        isnan(x) && break
+        isnan(y) && break
+        setopacity(opacity)
+        move(center + (xo, yo))
+        line(center + (x, y))
+        do_action(:stroke)
+        xo, yo = x, y
+        opacity -= 1.0 / (n - 1)
+    end
+    grestore()
+end