NB. Apollonian circle packing algorithms
NB. J script by Cliff Reiter
NB. tested with J601
NB. September 2006

coinsert 'jgl2'
require 'trig gl2'

WIN_WH=:1000 1000    NB. default winx width and height

NB. In order to be able to save image files, 
NB. load image3.ijs from the image addon
NB. load '~addons/image3/image3.ijs'

NB. Hue y
NB. gives a pure hue running through red-yellow-green-cyan-blue-magenta-red
NB. as y runs from 0 to 1
Hue=:<.@(255.9999&*)@((-.,])@(1&|)+/ . *{&(#:7|3^10-i.8)@(0 1&+)@<.)@(6&*)"0

NB. a default palette for the circle
$circle_pal=:48>.}:Hue (i.%<:)16

NB. alternate white palette
NB. circle_pal=:,:3#255

NB. function that defines the next window name
nx_WIN=:3 : 0
try. WIN_num=: WIN_num + 1 catch. WIN_num=: 0 end.
WIN_nam=:'p',": WIN_num
)

NB. $b=:dwin_image ''  NB. save the current window in image3 form
dwin_image=:3 : 0      NB. capture window as raster array
wd 'psel ',WIN_nam,';'
glsel 'g'
z=.256 256 256#:(|.WIN_WH)$glqpixels 0 0,WIN_WH
)

NB. Opens a drawing window & define SC to convert given coords to pixels locations
NB. (xy_low_left,xy_up_right) dwin 'name'
dwinx=: 3 : 0           NB. pixel sized window drawing
0 0 1 1 dwinx y
:
'a c b d'=.x        NB. note non-alphabetic order
SC=:WIN_WH&*@((-&(a,d))"1 %((b-a),c-d)"1)
sz=.":WIN_WH
PN=:y,';   Window bounds are x: ',(":a),' to ',(":b),'   y: ',(":c),' to ',":d
nx_WIN ''
z=.'pc ',WIN_nam,' closeok;pn "',PN,'";cc g isigraph;setxywhx g 2 2 ',sz,';pas 2 2'
wd z,';pshow;'
)

NB. draw circles given as curvature 
NB. and curvature time center
drcircles=:3 : 0
wd 'psel ',WIN_nam
glrgb 0 0 0
glpen 1 0
for_yy. y do.
  yy=.x:^:_1 yy
  cen=. SC tcen=.(}. yy)%{.yy
  r=.cen-~x:^:_1 SC (%{.yy)+tcen
  glrgb circle_pal{~(#circle_pal)|round {.yy
  glbrush ''
  glellipse (cen-r),+:r
  glpaint ''
end.
)

NB. draw circles given as curvature 
NB. and curvature time center
NB. draw text giving curvatures for 1-3 digit #'s
drtcircles=:3 : 0
wd 'psel ',WIN_nam
glrgb 0 0 0
glpen 1 0
for_yy. y do.
  yy=.x:^:_1 yy
  cen=. SC tcen=.(}. yy)%{.yy
  r=.cen-~x:^:_1 SC (%{.yy)+tcen
  glrgb circle_pal{~(#circle_pal)|round {.yy
  glbrush ''
  glellipse (cen-r),+:r
  nd=.<. 10 ^. {.yy
  if. nd e. 0 1 2 do.
   text=.":round {. yy
   select. nd
    case. 0 do.
    glfont '"Square721 BT" ',":1>.{.r
    case. 1 do.
    glfont '"Square721 BT" ',":1>.0.75*{.r
    case. 2 do.
    glfont '"Square721 BT" ',":1>.<.0.5*{.r
    end.
    ext=.glqextent text
    gltextxy round cen--:ext
    gltext ":round {. yy
  end.
  glpaint ''
end.
)

NB. rounding utility
round=:<.@(0.5&+)

NB. takes 4 radius (first negative) and 
NB. attempts to create complex curvature-center 
NB. form for the Descartes configurations
to_cf=:3 : 0
'R A B C'=.y
r=.-%R
a=.%A
b=.%B
c=.%C
ct=.((*:r-a)+(*:a+b)-*:r-b)%2*(r-a)*a+b
st=.%:1-*:ct
t=.(r-a)j. 0
cf=.(R,0),(A,t),:B,t+(a+b)*(-ct) j. st
ct=.((*:r-a)+(*:a+c)-*:r-c)%2*(r-a)*a+c
st=.%:1-*:ct
cf=.cf,C,t+(a+c)*-ct j. st
)

NB. convert complex-center form to curvature rational form 
NB. round-off remains problematic
to_rcf=:3 : 0
7 to_rcf y
:
0 2 3{"1 ,"2 (*>&(1r10^x)@|)x:+.*/\"1 to_cf y
)

NB. uses the matrix products to build more Descartes
NB. configurations from those given and avoids those already
NB. done; input boxed list of "done" and "todo" configurations
NB. initially applied to empty;singleton Descartes configuration
more_desc=:3 : 0
'done todo'=.y
smoutput #todo
k=.2000>L todo
$m0=.(/:~)"2,/SS&( +/ . *)"2 k#todo
done=.~.done,todo
todo=.~.m0-.done
done;todo
)

NB. version for the 0011 special configuration
more_desc_0011=:3 : 0
'done todo'=.y
smoutput #todo
k=.2000>L todo
k=.k*. */"1] 1.2>|({:"1 todo)%{."1 todo
$m0=.(/:~)"2,/SS&( +/ . *)"2 k#todo
done=.~.done,todo
todo=.~.m0-.done
done;todo
)

NB. the matrices performing the Descartes swaps
]S1=:_1 2 2 2 (0)}=i.4
]S2=:2 _1 2 2 (1)}=i.4
]S3=:2 2 _1 2 (2)}=i.4
]S4=:2 2 2 _1 (3)}=i.4
SS=:S1,S2,S3,:S4

NB. measures the "size" of the configurations
L=:+/@:{.@:|:"2

NB. gives curvatures from Diophantine solutions to 
NB. ((n^2)+m^2)=d1 *d2
nmd1d2_to_abcd=:3 : 0
'n0 m0 d1 d2'=.y
z=.n0,(d1-n0),(d2-n0),d1+d2-n0+2*m0
)

ic=:(i.0 4 3);,:to_rcf _2 3 6 7
$&.> md=:more_desc^:_ ic
$d=:/:~~.,/;md
#d=:d#~2000>{."1 d
$circle_pal=:48>.}:Hue (i.%<:)16
(_1 _1 1 1*0.51) dwinx 'ac'
drtcircles d

NB. to save the image with the image addon
NB. z=:dwin_image ''
NB. 9 z write_image 'ac2367.png'

other_things_to_try4=: 0 :0

NB. _1 2 2 3
ic=:(i.0 4 3);,:to_rcf _1 2 2 3
$&.> md=:more_desc^:_ ic
$d=:/:~~.,/;md
#d=:d#~2000>{."1 d
$circle_pal=:48>.}:Hue (i.%<:)16
(_1 _1 1 1*1.01) dwinx 'ac'
drtcircles d

NB. 0 0 1 1 
]r0011=: 0 1,0 _1 0,1 0 1,:1 0 _1x
]ic=:(i.0 4 3);,:r0011
$&.> md=:more_desc_0011^:_ ic
$d=:/:~~.,/;md
#d=:d#~2000>{."1 d
$circle_pal=:48>.}:Hue (i.%<:)16
(_1 _1 1 1*1.0) dwinx 'ac'
drtcircles d
)