////////////////////////////////////////////////////////////
// DEFINITION OF THE SYSTEM

N := 100 ;;

SCALE := 10000 ;;

rose := |N>,|NE>,|E>,|SE>,<N|,<NE|,<E|,<SE|,():ring ;;

trans searchPos[v] = {
  x / (x==v) => return(^x)
} ;;

// Wind

type wind = [float]Moore ;;

fun new_wind(dir:posgbf, force:(`Strong|`Weak|`None)) = (

  let (d,md,orth,diag,mdiag) =
    if (force == `Weak) then (0.9,1.1,1.0,1.0,1.04)
    elif (force == `Strong) then (0.5,1.8,1.0,0.7,1.3)
    else (1.0,1.0,1.0,1.0,1.0) fi
  in
  let p = searchPos[v=dir](rose) in
    {| d     @ rose.(p),
       md    @ rose.(p+4),
       orth  @ rose.(p+2),
       orth  @ rose.(p-2),
       diag  @ rose.(p+1),
       diag  @ rose.(p-1),
       mdiag @ rose.(p+3),
       mdiag @ rose.(p-3)
    |}:Moore

) ;;

// State

type       system = [cell]Moore
and record cell   = {
  wind:wind,
  burnt:float,
  height:float,
  rate:float
}
and record burnable = cell + { burnt ~= 1.0 } + { rate ~= 0.0 }
and record burnt    = cell + { burnt  = 1.0 }
;;

// Specification of a configuration constructor
fun new_grid(f) = (
  reverse(fold((\i,acc.(
    reverse(fold((\j,acc.(  f(i-1,j-1) :: acc  )),seq:(),N)) :: acc
  )),seq:(),N)) following |E>, |N>
) ;;

fun export_grid[ofile="/tmp/toto"](color,g) = (
  stdout << "iteration: " << 'iteration << "\n";
  ofile << " { " ;
  fold((\i,acc.(
    fold((\j,acc.(  let c = color(g.((i-1)*|E> + (j-1)*|N>)) in
      if c <> WHITE then
        ofile << " { " << (i-1) << " " << (j-1) << " " << c.r << " " << c.g << " " << c.b << " } "
      fi
    )),<undef>,N)
  )),<undef>,N) ;
  ofile << " } " ;
  g
) ;;

fun color(v) =
  if   burnt(v)       then WHITE
  elif v.burnt == 0.0 then WHITE
  else mult_color(WHITE,GRAY(1.0-v.burnt))
  fi
;;

// Definition of an initial state
fun init(i,j) = (
  let w = new_wind(|E>,`Strong) in
  let r = if (i > N/5) && (i < 2*N/5) && (j > N/5) && (j < 2*N/5) then 0.01 else 0.1 fi in
  let b = if (i,j) == (N/2, N/2) then 0.1 else 0.0 fi in
  let h =
    let i' = (8.*i)/N in
    let (c1, c2, c3, c4, c5) = map((\c.((N/2.-abs(N/2.-j))*c/5.)),iota(5,seq:())) in
    if   (0 <= i') && (i' < 1) then c1
    elif (1 <= i') && (i' < 2) then (2-i')*c1 + (i'-1)*c2
    elif (2 <= i') && (i' < 3) then c2
    elif (3 <= i') && (i' < 4) then (4-i')*c2 + (i'-3)*c3
    elif (4 <= i') && (i' < 5) then c3
    elif (5 <= i') && (i' < 6) then (6-i')*c3 + (i'-5)*c4
    elif (6 <= i') && (i' < 7) then c4
    //elif (7 <= i') && (i' < 8) then (8-i')*c4 + (i'-7)*c5
    else                              (8-i')*c4 + (i'-7)*c5 fi
    in
    !! h >= 0. ;
    stdout << h << "\n" ;
    { wind = w, burnt = b, height = h, rate = r }
) ;;

g0 := new_grid(init) ;;

system(g0) ;;

////////////////////////////////////////////////////////////
// SPECIFICATION OF THE DYNAMICS

fun moore_coeff(d) = (
  if (d == <N|) || (d == |N>) || (d == <E|) || (d == |E>) then 1.0 else 1.0/sqrt(2.0) fi
) ;;

trans rules = {

  x / ((x.burnt != 1.0) && (neighborsexists((\y.(y.burnt != 0.0)),x))) => (
    let b = neighborsfold_indexed((\p,y,acc.(
      if y != <undef> then
        x.rate *
        moore_coeff(p-^x) *
        x.wind.(p-^x) *
        max(0.0,(1.0 + (x.height - y.height))) *
        y.burnt + acc
      else acc fi
    )),x.burnt,x) in
      x + { burnt = min(1.0,b) }
  );

} ;;

rules[prelude = export_grid(color),
      interlude = export_grid(color),
      fixpoint
     ](g0) ;;