Quasicrystalline structure in Just Intonations

just.typ

#set page(height: auto, width: auto)
#set text(font: "New Computer Modern", weight: 450)
#set par(leading: 0pt)

#import "@preview/cetz:0.3.4": canvas, draw
#import calc: *
#let dot(v1, v2) = v1.zip(v2, exact: true).map(((x, y)) => x*y).sum()
#let scale(s, v) = v.map(x => x*s)
#let add(v1, v2) = v1.zip(v2, exact: true).map(((x, y)) => x+y)

#let colors = (0,1,2).map(i => color.oklch(90%, 30%, i * 120deg + 54deg, 40%))
/*
We are on a plane P:
  x ln 2 + y ln 3 + z ln 5 = c.
Let's find a basis for that, e.g. (-ln 3, ln 2, 0), (0, -ln 5, ln 3) will do.
But a normalized basis is better.
*/
#let v1 = (-ln(3), ln(2), 0)
#let v2 = (0, -ln(5), ln(3))

#let n1 = sqrt(dot(v1, v1))
#let v1 = scale(1/n1, v1) // normalized v1
#let d12 = dot(v1, v2)
#let v2 = add(v2, scale(-d12, v1)) // orthogonalized v2
#let n2 = sqrt(dot(v2, v2))
#let v2 = scale(1/n2, v2) // normalized v2

// Convert between RR^3 and the plane P
#let toR3(s, t) = add(scale(s, v1), scale(t, v2))
#let toR2(..vec) = (dot(v1, vec.pos()), dot(v2, vec.pos()))

/*
Now we wish to intersect planes of integral coordinates with P,
which will become three sets of parallel lines.
Every intersection of the lines signify a parallelogram.
Using s and t as coordinates on P to prevent confusion, we have
  (s v1 + t v2)[i] = n
where [i] fetches the i-th coordinate. Let's compute the intersection.
*/
#let intersect(dir1, lvl1, dir2, lvl2) = {
  if dir1 == dir2 {
    panic("Identical direction!")
  }
  // Regular 2x2 matrix
  // a11 a12 | b1
  // a21 a22 | b2
  let a11 = v1.at(dir1)
  let a12 = v2.at(dir1)
  let a21 = v1.at(dir2)
  let a22 = v2.at(dir2)
  let det = a11 * a22 - a12 * a21 // TODO precompute
  return ((lvl1 * a22 - lvl2 * a12)/det, (-lvl1 * a21 + lvl2 * a11)/det)
}

#let letters = "FCGDAEB"
#let octaves = (-1,0,0,1,1,2,2)
#let display(oct, fif, thd) = {
  let fif = fif + thd * 4
  let oct = oct - thd * 2
  let sharps = div-euclid(fif, 7)
  let sharp-octave = sharps * 4
  let letter = rem-euclid(fif, 7)
  let letter-octave = octaves.at(letter)

  let sharps-content = if sharps > 0 {
    $sharp$ * sharps
  } else {
    $flat$ * -sharps
  }
  let downs-content = if thd > 0 {
    $arrow.b$ * thd
  } else {
    $arrow.t$ * -thd
  }
  set align(center)
  [
    #letters.at(letter)#super(sharps-content + downs-content + $zwj$) \
    #set text(size: 0.8em)
    $#(sharp-octave + letter-octave + oct)$ octs
  ]
}

#let Radius = 6
#figure(canvas(length: 4cm, {
  import draw: *
  let allPoints = (:)

  for dir1 in range(3) {
  for dir2 in range(dir1) {
    let dir3 = rem(- dir1 - dir2, 3) // funny trick specific to mod 3
    for lvl1 in range(-5, 6) {
    for lvl2 in range(-5, 6) {
      let posR2 = intersect(dir1, lvl1, dir2, lvl2)
      if dot(posR2, posR2) > Radius*Radius {
        continue
      }
      let posR3 = toR3(..posR2)
      let lvl3 = ceil(posR3.at(dir3))
      let points = ()
      for (eps1, eps2) in ((0,0), (0,1), (1,1), (1,0)) {
          // The 3D lattice point is at (lvl1 + eps1, lvl2 + eps2, lvl3), reordered
          let pt = (none, none, none)
          pt.at(dir1) = lvl1 + eps1
          pt.at(dir2) = lvl2 + eps2
          pt.at(dir3) = lvl3

          let pt2 = toR2(..pt)
          // By the way, see if allPoints has it
          if repr(pt) not in allPoints {
            allPoints.insert(repr(pt), pt2)
          }

          // Project down to P
          points.push(pt2)
      }
      line(..points, close: true, fill: colors.at(dir3), stroke: 2pt)
    }
    }
  }
  }
  // Draw lattice points
  for (pt3, pt2) in allPoints.pairs() {
    circle(pt2, radius: 2em, stroke: 2pt, fill: white)
    content(pt2, display(..eval(pt3)))
  }


  // circle((0,0), radius: Radius, stroke: 2cm + white)
}))