// CE0 / Simple Cellular Automaton on a Torus '08.02
// by Tsutomu HIGO  URL:http://www.asahi-net.or.jp/~nj2t-hg/

#version 3.6;
global_settings {
  assumed_gamma 2.2
  max_trace_level 5
  }

light_source {<-300, 1000, -2000> color <1, 1, 1>}

camera {
  location <0, -140, -400>
  angle 84
  look_at <6, 76, -10>
  }

#declare Nmay = 1800;
#declare Nmax = 200;
#declare Rs=Nmax/pi;
#declare Rl=Nmay*1.732/2/2/pi;
#declare Cp = array[2*Nmax][2]
  #declare Nx = 0; #while (Nx < 2*Nmax)
    #if (Nx=0) #declare Cp[Nx][0] = 1; #else #declare Cp[Nx][0] = 0; #end
  #declare Nx = Nx+1; #end

#macro Obje (Rss, M, Mv)
  sphere {<Rss, 0, 0>, 0.4
    pigment {color rgb <Mv/M, 0, 1-Mv/M>}
    finish {
      reflection 0
      phong 1
      brilliance 0.5
      ambient 1
      }
    }
#end

#declare Cmod = 12;
union {
  object {Obje (Rs, Cmod, 1) rotate y*180/Nmax*0 translate x*Rl rotate z*360/Nmay*0}
  #declare Ny = 1; #while (Ny < Nmay*0.84)
    #declare Nx = 0; #while (Nx < 2*Nmax)
      #declare Cp[Nx][1] = mod(Cp[Nx][0]+Cp[mod(Nx+1,2*Nmax)][0], Cmod);
      #if (Cp[Nx][1] > 0) object {Obje (Rs, Cmod, Cp[Nx][1]) rotate y*180/Nmax*(Nx+0.5*Ny) translate x*Rl rotate z*360/Nmay*Ny} #end
    #declare Nx = Nx+1; #end
    #declare Nx = 0; #while (Nx < 2*Nmax)
      #declare Cp[Nx][0] = Cp[Nx][1];
    #declare Nx = Nx+1; #end
  #declare Ny = Ny+1; #end
  no_shadow
  rotate x*90  rotate y*10  rotate z*340
  }