// Animate a heat flow equation for 100 time steps and produce an output
// lattice each 5 time steps. The equation iterated at each timestep is
//		T += 4 h^2 (1-k) Laplacian( T );
// where h is the grid spacing and k is a relaxation parameter. 
// T is the temperature.
//
// No input ports
// Port First_Out: Temperatures on square grid every 5 iterations
// After 100 iterations the temperature range is approximately [-80 ... 80].
   
   int stepPerIter = 5;

   n := 50;
   
   float next[n, n];  		// space for next iteration
   float temps[n, n]; 		// current temperatures starting a 0

// Alias a heat source and sink in the grid
   source := subshape([20, 30], [29, 44], temps);
   sink   := subshape([30, 20], [39, 24], temps);
   
// Alias neighboring cells in the four directions plus center
// and the current temperatures.
   M := subshape([1, 1], [n-2, n-2], temps);  	// all but edge of temps
   
   N := subshape([0, 1], [n-3, n-2], next);  	// north
   S := subshape([2, 1], [n-1, n-2], next);  	// south
   E := subshape([1, 0], [n-2, n-3], next);  	// east
   W := subshape([1, 2], [n-2, n-1], next);  	// west
   C := subshape([1, 1], [n-2, n-2], next);  	// center
   
// The iteration is T += 4 h^2 (1-k) Laplacian( T );
   k := (float) 0.3;  // weight of C
   int i;
   for (i = 0; i <  100; i+=1) {
     source += (float) 10.;
     sink   -= (float) 10.;
     next = 0;  
     N += M; S += M; E += M; W += M;
     next *= (float) 0.25 * (1 - k);
     C += M * k;
     temps = next;
     if (i % stepPerIter == 0) {

// And produce a scalar lattice for rendering.
       First_Out := scalar_lattice_out(temps);

       write(First_Out, oport_First_Out);
       flush_output_port(oport_First_Out);
     }
   }