# sacks-spiral.tcl
 # Author:      Gerard Sookahet
 # Date:        18 Feb 2012
 # Description: Plot Sacks prime spiral and Sacks divisor spiral 

 package require Tk
 bind all <Escape> {exit}

 proc SpiralMain { N } {
  set w .sp
  catch {destroy $w}
  toplevel $w
  wm withdraw .
  wm title $w "Sacks spiral"
  wm geometry $w +100+10

  set dim [expr {int(sqrt($N) + 10)}]
  set mid [expr {$dim/2}]
  pack [canvas $w.c -width $dim -height $dim -bg black]

  set f1 [frame $w.f1 -relief sunken -borderwidth 2]
  pack $f1 -fill x
  button $f1.bu -text "Sacks spiral" -width 12 -bg blue -fg white \
        -command "PlotSacks $w $N $mid prime"
  button $f1.bd -text "Divisor spiral" -width 12 -bg blue -fg white \
        -command "PlotSacks $w $N $mid divisor"
  button $f1.bq -text Quit -width 5 -bg blue -fg white -command exit
  eval pack [winfo children $f1] -side left
 }

 proc PlotSacks {w N mid type} {
  $w.c delete all
  set pix [image create photo]
  $w.c create image 0 0 -anchor nw -image $pix
  set xo $mid
  set yo $mid

  set n 1
  set 2pi 6.28318531
  set M [expr {$N/4}]

  if {$type == "prime"} {
    set cmap1 #00FFFF
    set cmap2 #606060
    while {$n < $M} {
       set sqn [expr {sqrt($n)}] 
       set 2pisqn [expr {$2pi*$sqn}]
       set x [expr {int(-cos($2pi*$sqn)*$sqn) + $xo}]
       set y [expr {int(sin($2pi*$sqn)*$sqn) + $yo}]
       if [IsPrime $n] {$pix put $cmap1 -to $x $y} else {$pix put $cmap2 -to $x $y}
       incr n
       update idletasks
    }
  } else {
    set cmap #030303
    while {$n < $M} {
       set sqn [expr {sqrt($n)}] 
       set 2pisqn [expr {$2pi*$sqn}]
       set x [expr {int(-cos($2pi*$sqn)*$sqn) + $xo}]
       set y [expr {int(sin($2pi*$sqn)*$sqn) + $yo}]
       $pix put [colormap [NbDivisor $n]] -to $x $y
       incr n
       update idletasks
    }
  }
 }

 # Primality testing
 proc IsPrime { n } {
  if {$n==1} {return 0}
  set max [expr {int(sqrt($n))}]
  set d 2
  while {$d <= $max} {
       if {$n%$d == 0} {return 0}
       incr d
  }
  return 1
 }
 # Return the number of divisors of an integer
 proc NbDivisor { n } {
  set max [expr {int(sqrt($n))}]
  set nd 0
  for {set i 2} {$i <= $max} {incr i} {
     if {$n%$i == 0} {incr nd}
  }
  return $nd
 }
 # Arbitrary color table
 proc colormap { n } {
  set lcolor {#030303 #CD0000 #CD4F39 #EE4000 #EE6A50 #FF7F00 #EE9A00 \
              #FF8C69 #FFC125 #EEEE00 #EED5B7 #D2691E #BDB76B #00FFFF \
              #7FFFD4 #FFEFD5 #AB82FF #E066FF
  }
  return [lindex $lcolor $n]
 }
 # The maximum integer. The canvas is sized from its square root
 SpiralMain 100000