#
# field:
#
# A field represents a discretization of a scalar field (i.e., a
# real-valued function) on a grid.
#

printf("field procedures:  \
create_field, \
initialize_field, \
Gaussian_func, \
advance_field, \
get_profile_of_field, \
get_profile_of_field_at_time");


create_field := proc(grid)

#
# Create a field on a grid.
#
# The return value from create_field(), say u, is the new field.
# The value of u at grid[j] is u[j], j = 0, ..., n_cells.
#

  local n_cells;

  n_cells := get_n_cells_of_grid(grid);
  array(0..n_cells)
end:


initialize_field := proc(field, grid, func)

#
# Initialize a field on a grid with values determined by a function.
#
# The return value from create_field(), say u, is the new field.
# The value of u at grid[j] is func(grid[j]), j = 0, ..., n_cells.
#

  local n_cells, j;

  n_cells := get_n_cells_of_grid(grid);
  for j from 0 to n_cells do
    field[j] := func(grid[j]);
  od
end:


Gaussian_func := proc(center, width, amplitude)

#
# Return a bump function (i.e., a Gaussian) characterized by its
# center x_0, width dx, and amplitude A:
#
#		x -> A*exp(- [(x - x_0)/dx]^2).
#

  x -> evalf(amplitude*exp( - ((x - center)/width)^2 ))
end:

bump_func := Gaussian_func:


advance_field := proc(field, grid, stencil, l_bc, r_bc, dt)

#
# Advance the field one time step using the stencil and boundary conditions.
#

  local n_cells, u_l, u_m, dx, j;

  n_cells := get_n_cells_of_grid(grid);
  dx := get_spacing_of_grid(grid);

  u_l := l_bc(field[0], field[1], dx, dt);
  u_m := field[0];
  field[0] := stencil(u_l, u_m, field[1], dx, dt);

  for j from 1 to n_cells-1 do
    u_l := u_m;
    u_m := field[j];
    field[j] := stencil(u_l, u_m, field[j+1], dx, dt);
  od;

  u_l := u_m;
  u_m := field[n_cells];
  field[n_cells] := r_bc(u_l, u_m, dx, dt);
end:


get_profile_of_field := proc(field, grid)

#
# Return the field values along the grid as a list of pairs
# suitable for 2d plotting or writing to a file.
#

  local j;

  [seq([grid[j], field[j]], j=0..get_n_cells_of_grid(grid))]
end:


get_profile_of_field_at_time := proc(field, grid, time)

#
# Return the field values along the grid as a list of triples
# suitable for 3d plotting or writing to a file.
#

  local j;

  [seq([grid[j], time, field[j]], j=0..get_n_cells_of_grid(grid))]
end: