NB. 3x+1 utilities
NB. Jeffrey Dumont & Clifford Reiter (reiterc@lafayette.edu)
NB. December 2000
NB. functions giving complex generalizations of the 3x+1 function and 
NB. utilities for visualizing various stopping times and creating
NB. some images assoicated with "Visualizing Generalized 3x+1 Function Dynamics"
NB. images require the additional script raster4.ijs or raster4+.ijs
NB. last ten lines create images -- increase numbers for higher resolution

load 'fvj2\raster4+.ijs'

mod2a=: *:@(1&o.)@(1r2p1&*)   NB. sin^2

mod2b=: 2r1&|                 NB. Re & Im remainder
mod2f=: |@(2r1&|&.>:)&.+.     NB. sawtooth
mod2g=: mod2a&.+.             NB. sin^2 on Re & Im parts
mod2w=: -:@(1:-_12&o.@o.)     NB. winding

T1=: 1 : '-:@((]+[* 3: ^ ]) u.) f.'
TX=: 1 : '-:@((]+[* (m."_) ^ ]) mod2a) f.'
c3x=: + 0.25"_ * 1: - 2&o.@o. * >:@+:

SIG=:1 : 0     NB. fct arg is "T"
sigs=.0 1&+@}: , u.@{:
sigl=.sigs^:((1&{ < 255"_)*. {. <:&| {:)^:_
(((~:>:)@{:*1&{)@sigl@(],0:,]) f. :: 0:)"0
)

TOT=:1 : 0     NB. fct arg is "T"
sigs=.>:@{. , u.@{:
sigl=.sigs^:(({. < 255"_)*. 1: ~: {:)^:_
(((~:>:)@{:*{.)@sigl@(0:,]) f. :: 0:)"0
)

LAM1=: 1 : 0    NB. fct arg is "mod2"
lams=: ({.+u.@(u. T1)@{:),>:@(1&{),(u. T1)@{:
laml=: lams^:((1&{ < 255"_)*. >:/@:(3 2&^)@:|@}:)^:_
(((~:>:)@{:*1&{)@laml@(|@u.,1:,]) f. :: 0:)"0
)


ESC=:1 : 0     NB. fct arg is "T"
sigs=.(>:@{. , u.@{:) :: (>:@{. , _:)
sigl=.sigs^:(({. < 255"_)*. ( |@{: < 1e10"_)*.(~:>:)@{:)^:_
({.@sigl@(0:,]) f. )"0
)

$pal3x=:255,0,~254{.<.,/ ( >:-:(i.%-)22) */(hue (i.%])12)

zl_clur=: 4 : 0
w=.-~/9 o.y.
h=.-~/11 o.y.
(h*(i:%j.)<.0.5+x.*h%w) +/ ({.y.)+w*(i.%<:)1+x.
)

zl_cccr=: 4 : 0
w=.--/y.
({.y.)+w*((i:%j.) +/ (i:%])) <.-:x.
)

mk_strip=:1 : 0
2000 u. mk_strip y.
:
for_k. Stripe_starts do.
  b=. u. x. zl_clur k+Stripe_offset
  (pal3x;b) writebmp8 y.,(":k),'.bmp'
end.
if. 1<#Stripe_starts do.
  for_k. Stripe_starts do.
    'p b'=. readbmp8 y.,(":k),'.bmp'
    try. B=.B,.b catch. B=.b end.
  end.
(pal3x;B) writebmp8 y.,'.bmp'
end.
)

mk_zoom=:2 : 0
1000 u. mk_zoom n. y.
:
for_k. Zoom_lev do.
  b=. u. x. zl_cccr n.+0 1*1.5*10^-k
  (pal3x;b) writebmp8 y.,(":k),'.bmp'
end.
)

Zoom_lev=:i.2               NB. 2 zoom levels
Stripe_starts=: _6+12*i.1   NB. 1 stripe, starts at _6
Stripe_offset=: 0 12j2      NB. stripe offset/width/upper corner

NB. Note the directory referred to in the path needs to exist
path=:'\3x+1\'
300 mod2a T1 ESC  mk_strip path,'mod2a_esc'
300 mod2a T1 SIG  mk_strip path,'mod2a_sig'
300 mod2a T1 TOT  mk_strip path,'mod2a_tot'
300 mod2a  LAM1 mk_strip   path,'mod2a_lam1'

200 mod2a T1 SIG  mk_zoom 27 path,'mod2a_sig_z27z'

300 mod2w T1 ESC  mk_strip path,'mod2w_esc'
300 (2.5 TX) ESC  mk_strip path,'tx_5r2'