ofile := "heterocyst.jbv" ;;
K     := random(2.0) ;;
v     := 0.5 ;;
R     := 1.1 ;;
T     := 0.1 ;;
Smax  := 0.8 ;;
dt    := 0.001 ;;

record W = { J: float } ;;
record C = { a:`V|`H, c:float, s:float, p:`Right|`Left } ;;
record H = C + { a = `H } ;;
type anabaena = [W|C]seq ;;

S0 := ({ a = `H, c = 1.0, s = 1.0, p = `Left  }, { J = 0.0 },
       { a = `V, c = 1.0, s = 1.0, p = `Right }, { J = 0.0 },
       { a = `H, c = 1.0, s = 1.0, p = `Left  }) ;;

trans rules = {

  wall:W => { J = K * (left(wall).c - right(wall).c) } ;

  cell:H => cell; 

  cell:C =>
    let cell' = cell + { c = cell.c + dt * ((left(cell).J - right(cell).J) - v * cell.c),
                         s = cell.s * (R ** dt) }
    in
      if cell'.c < T
      then cell' + { a = `H, c = 1.0, s = 1.0 }
      elif cell'.s > Smax
      then if cell'.p == `Right
           then (cell' + { s = 3.0 * cell'.s / 4.0, p = `Left  }, { J = 0.0 },
                 cell' + { s = 1.0 * cell'.s / 4.0, p = `Right })
           else (cell' + { s = 1.0 * cell'.s / 4.0, p = `Left  }, { J = 0.0 },
                 cell' + { s = 3.0 * cell'.s / 4.0, p = `Right }) fi
      else cell' + { a = `V } fi

} ;;

fun output(S) = (
  if iteration % 1000 == 0 then
    fold_indexed((fun p, c, acc ->
                    if C(c)
                    then
                      ofile << " { " << (2.0 * acc + c.s) << " " << (iteration / 500);
                      if c.a == `H then ofile << " 1 0 0 } " else ofile << " 0 0 1 } " fi;
                      acc + c.s
                    else acc fi),
      0.0, S)
  fi;
  S
);;

ofile << " { ";;
rules[iter = 40000,
      prelude = output,
      interlude = output,
      postlude = output
     ](S0) ;;
ofile << " }\n" ;;