NB. Chaotic Attractor Exhibiting Quasicrystalline Structure
NB. Can be used to duplicate figure qca00
NB. Cliff Reiter, October 12, 2001

NB. Running this script creates (about 5 min on 750MHz)
NB. a low resolution version of the image qca_test.bmp in
NB. in the J directory. See the notes at the end of this 
NB. regarding creation of high resolution images.

NB. this script requires additional scripts
NB.   fvj2\raster5.ijs and fvj2\dwin+.ijs
NB. 

require 'fvj2\raster5 fvj2\dwin+ misc files'

NB. find hyperfacets
hv=:((i.1 0)"0)`(#:@:i.@:(2&^))@.*                                 NB. hyper vertices
nhf=: 2&^@-~ * !                                  NB. number hyper facets
perms2=: ((i.@! A. i.)@{:@$ {"1"1 2 ])"2          NB. permutations of rows of matrices
hf=: [: ~.@:;@:(<@:~.@:(/:~"2)@perms2)  hv@[ ,"1"_ 1 hv@-~     NB. hyperfacets

Perms2=:3 : 0
0 Perms y.
:
n=.{:$y.
z=.i.0 1 1*$y.
i=.i.M=.!n-x.
for_m. i.*/n-i.x. do.
  z=.~.z, /:~"2 ,/((i+M*m) A. i.n) {"1"1 2"2 y.
end.
z
)

HF=:4 : 0
if. y.<7
  do. x. hf y.
  else. 3 Perms2 x.(hv@[ ,"1"_ 1 hv@-~ )y. end.
)

NB. vector utilities
vmag=:+/&.(*:"_)
unit=:% +/&.(*:"_)

mkU=: 3 : 'unit"1&.|: (,./ +._12 o. (>:i.<:>.-:y.)*/(2p1%y.)*i.y.),.,./2 o. ((2|y.)}.1p1 0)*/i.y.'

init_n=: 3 : 0    NB. initize on dimension
U=:mkU y.
UI=:%.U
P_2=:2&}. @(UI&(+/ . *))"1
P2=:2&{. @(UI&(+/ . *))"1
fc=:(y.-3) HF y.
pfc=:P_2 fc{_0.5 0.5
fs=:({~ 1: i.~  0: ~: ])@:*   NB. first sign
hpl=: (3 : '(*fs)(% vmag@:}:) (1,(#y.)# 0) %. 1,.~0, (**|)y.') :: 0:
hfc=:(**|)hpl"2 pfc
ht=:(}:"1 hfc) (}:@{.,>./@:({:"1))/.  hfc
hn=:}:"1 ht
hmax=:{:"1 ht
invq=: *./@:(<:&hmax)@:|@:(hn&(+/ . *))"1
$hn
)

sglat=: 3 : 0     NB. search for good lattice points
6 sglat y.
:
n=.#y.
init_n n
off=:(U +/ . * y.)&+"1
$tlat0=:,/^:([:{._2:+$@$) >{n#<_1 0 1
$ntlat=:pk tlat0
$glat=:tlat=:i.0 1
]bksz=:3^8
while. 0<#ntlat do.
  ]sz=:bksz <. # ntlat
  $m=:invq@:P_2@:off n upk ntlat0=:sz{.ntlat
  $glat=:glat,n upk nglat=: m # ntlat0
  $tlat=:tlat,ntlat0
  $ntlat=:sz}.ntlat
  $ntlat0=:,/(n upk nglat) +"1"1 _ tlat0
  $ntlat0=:pk ntlat0 #~ *./"1 x. >:| ntlat0
  $ntlat=:~.ntlat,ntlat0 -. tlat
  if. 0<+/m do. (1!:2)&2 ": ($glat),($tlat),($ntlat) end.
end.
$glat
)

pk=:_5&((5#33)&#.@(5&{.)@(16&+)\)"1
upk=: 1 : '_16&+@(m.&{.)@(,/"2)@:((5#33)&#:)'

pedges=:3 : 0
k=.0
o=.=i.#{.glat
z=.''
while. k<#y. do.
  z=.z,<v,:"1"1 v+"1 o#~(o+"1 v=.k{y.)e."1 _ y.
  k=.k+1
end.
P2@:off ; z
)
upedges=:3 : 0
k=.0
o=.=i.#{.glat
z=.''
while. k<#y. do.
  z=.z,<v,:"1"1 v+"1 o#~(o+"1 v=.k{y.)e."1 _ y.
  k=.k+1
end.
; z
)

NB. ------------------------------------------------
NB.* cln254 v Gives log-bias cummulative discretization  
cln254=:3 : 0
  o=./:~ ,y.
  m=.(}.~:}:) o
  p=.((#i.@#),#) m
  n=.p{o
  ($y.)$(n i. ,y.){ 0 , >:@<.@(253.99&*)@(+/\ % +/)@:^.@ (}.-}:) p
)

NB. translation vector
]tv=:0 0 0.696 0.879,-:%:5

projc=: P2@:off@:(invq@P_2@:off # ])"2

pr_proj_z=: 2 : 0     NB. argument y. has *.bmp extension 
800 800 1.8 u. pr_proj_z n. y.
:
VRAWH=:0 1{x.
mm=.2{x.
if. fexist (_3}.y.),'jba'
  do.
  'seed tVRA'=:3!:1^:_1 fread (_3}.y.),'jba'
  VRAWH=:|.$tVRA
  else.
  seed=: f^:(0{n.) 0.1*1+i.5
  tVRA=:(|.VRAWH)$0
  end.
(mm*_1 _1 1 1) vwin 0
(>.mm) sglat tv
255 vline"2 u. upedges glat
eVRA=:VRA
VRA=:tVRA
glaty=:~.,/glat +"1 _"_ 1-#:i.2^5
smoutput (#glat),#glaty
for_k. i. 1{n. do.
  x5=:f^:(1+i.2{n.) seed
  seed=:{:x5
  for_glatpt. glaty do.
    vfpixel (*./"1@:(<:&mm@|) # ]) u. glatpt+"1 x5
  end.
  cvra=.cln254 VRA
  (P256;|. cvra ) writebmp8 y.
  (P256;|. cvra >. eVRA) writebmp8 (_4}.y.),'e.bmp'
  (3!:1 seed;VRA) fwrite (_3}.y.),'jba'
  smoutput k
end.
)

f=: 1: | ] + }:@(+/)@(] +/ .*"1 2 (10 6 6$1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 _1 0 0 0 0 0 0 0 0 0 0 1 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 0 0 0 1 0 0 0 _1 0 0 0 0 0 0 _1 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 0 0 1 0 0 0 0 _1 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 0 1 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 1)"_)@:((, 1:)@(1: | (((5 5 5$0.00073634718147697607 _0.00021293450539273625 _0.00079023829548641883 0.00046496611876875723 0.00068554668192336954 _1.6930105180414578e_5 _0.00054428546088215474 0.00019425538759143589 0.00085028990344171397 0.00082239266468434849 _4.6498912429901668e_5 0.0004927713288858267 7.7129015770859878e_6 _0.00036927108035354978 _0.00033905243812040655 _0.00045433266223811797 3.094149383803973e_5 3.367886788554369e_5 4.0724462867457916e_5 0.00045603926894817438 0.0006519818172342907 _0.00014161067189798746 _5.0569307106145944e_5 8.1648036913618566e_5 0.00025854794239490108 _0.00058474201856622654 0.000240890708045001 0.00065012040142912466 0.00057357390557748623 5.6618725933710746e_5 _0.00040908150133591283 0.00056720242329612596 0.00097111607340268427 _0.00045216974751294573 0.00038304926304430628 _9.104683817461772e_5 _0.00022421631424149162 _0.00040359952804036397 0.00070272755874574888 0.00074206651443808024 _8.810547281811517e_5 _0.00078868879078765287 0.00050749157882061936 _0.0005890465598451407 _0.00010553453313702206 0.00028109456624937678 0.00035636492758650789 _0.00057467266877617011 _0.00052354744963678296 0.00073801022650488678 _0.00026213673433487931 0.00026790025948351275 0.0005996512170877784 0.00033799307545460131 0.00064960869462192864 _2.6682399238163999e_5 _0.00045109161310930891 0.00050325417658617175 0.00019293411932227588 0.00064373411391278196 _0.00076076033223491197 _9.890574412318488e_5 _0.00030884853034172997 _0.00081725486257789182 0.00039752322379269 _0.00082718865371047396 _0.00055970426421584956 _0.00094957212167752177 0.00054135057194284217 0.00047905058085830001 _0.00059689908995935606 _8.3008101287305449e_5 0.00088283448841686781 _0.00020076791308891765 _0.00030632154020138224 _0.00034613159331132808 0.00056630609982302976 _9.3392532405808977e_5 0.00035170076075168402 _0.00096532462534059272 _0.00021097837035635334 8.6523246100659841e_5 0.00019618871081265496 _0.00065634673336180556 0.0007804496994059996 _0.00098191601781352994 0.00093748846737434018 0.00036865599767179789 1.3421587658506064e_6 0.00055765454150517911 0.00049986688746127012 _0.00073723417682315618 _0.00069481192314514142 0.00029600531169772645 0.00096126356158523227 _4.3335786129660749e_5 _0.000344564968109657 0.00089657585220092521 0.00075033309846216176 0.00084837215565924744 0.00059080569966771246 _0.00032861813467798581 0.00091500521350186808 0.00049260833905321396 _0.00073165721421284335 _0.00096280137539170817 0.00019728350098923179 _0.00025620824401845069 _9.1963038865442467e_5 0.00037719868251792957 _0.00042175369935273531 _0.000414429546719817 0.00068260369792823078 0.00052033791180731525 _0.00068072815302827226 _0.00099807044490992139 _0.00056996761594895981 0.00055427538096371599 _0.00029368430682801113 4.7849614217737417e_5 0.00020845795721603365 _0.00044712252770282566 _0.00078832742831513977 0.00058091065155427041 _0.00063469169999909748)"_ +/ .* {:) +/ .* {.)@:((1: , 1&o. , 2&o. , 1&o.@+: , 2&o.@+:)"0)@(6.2831853071795862&*))@}:"1)@(] +/ .*"1 2 (10 6 6$1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 _1 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 0 1 0 0 0 _1 0 0 0 0 0 0 _1 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 0 0 1 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 0 0 0 1 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 _1 0 0 0 0 0 0 0 0 0 0 1 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 _1 0 0 0 0 0 0 1)"_)@(, 1:)

NB. makes a small picture qca_test.bmp
NB.                   and qca_teste.bmp (with rhombs)
NB. by running the script (takes about 5 min on 750 MHz)
NB. see below for higher resolution

600 600 1.5 projc pr_proj_z 200 1 10000 'qca_test.bmp'

NB. uncomment and run the last line below for a higher resolution;
NB. It will be 1600 x 1600 pixels with window into projection
NB. plane having radius 3.5; using 20000 iterations at a block
NB. and repeating the blocks 50 times. Several hours to complete
NB. but ctrl-break and rerunning the same command will continue
NB. the computation where it was last saved and you may stop
NB. whenever the quality is satisfactory.

NB. 1600 1600 3.5 projc pr_proj_z 2000 50 20000 'qca_high.bmp'