/* -*- mode: maxima -*- */
isolate_wrt_times:true$

/* Shifts in x and y directions. */
shift_x(expr, d) := op(expr)[args(expr)[1]+d, args(expr)[2]]$
shift_y(expr, d) := op(expr)[args(expr)[1], args(expr)[2]+d]$
shift[x] : shift_x; shift[y] : shift_y;

/* weights */
weight(var, direction, d) := shift[direction](subst(w, op(var), var), d)$

/* one-sided differences with weights */
d_px(var) := weight(var, x,  1/2) * (shift[x](var, 1) - var)$
d_mx(var) := weight(var, x, -1/2) * (var - shift[x](var, -1))$

d_py(var) := weight(var, y,  1/2) * (shift[y](var, 1) - var)$
d_my(var) := weight(var, y, -1/2) * (var - shift[y](var, -1))$

/* centered differences defined as sums of one-sided differences with weights */
D_x(foo) := d_px(foo) + d_mx(foo)$
D_y(foo) := d_py(foo) + d_my(foo)$

load("ssa.mac")$

denominator : (4 * dx**2 * dy**2)$

/* Clear the denominator */
e1 : ''lhs1 * denominator, expand;
e2 : ''lhs2 * denominator, expand;

/* define weights to be equal to 1 in the interior; give them PIK names otherwise */
a_ones : [w[i+1/2,j] = 1,
          w[i-1/2,j] = 1,
          w[i,j+1/2] = 1,
          w[i,j-1/2] = 1,
          N[i+1/2,j] = c_e,
          N[i-1/2,j] = c_w,
          N[i,j+1/2] = c_n,
          N[i,j-1/2] = c_s,
          w[i-1/2,j+1] = 1,
          w[i+1/2,j+1] = 1,
          w[i+1,j+1/2] = 1,
          w[i+1,j-1/2] = 1,
          w[i+1/2,j-1] = 1,
          w[i-1/2,j-1] = 1,
          w[i-1,j-1/2] = 1,
          w[i-1,j+1/2] = 1]$

a_names : [w[i+1/2,j] = aPP,
           w[i-1/2,j] = aMM,
           w[i,j+1/2] = bPP,
           w[i,j-1/2] = bMM,
           N[i+1/2,j] = c_e,
           N[i-1/2,j] = c_w,
           N[i,j+1/2] = c_n,
           N[i,j-1/2] = c_s,
           w[i-1/2,j+1] = aMn,
           w[i+1/2,j+1] = aPn,
           w[i+1,j+1/2] = bPe,
           w[i+1,j-1/2] = bMe,
           w[i+1/2,j-1] = aPs,
           w[i-1/2,j-1] = aMs,
           w[i-1,j-1/2] = bMw,
           w[i-1,j+1/2] = bPw]$

if interior then (e1 : at(e1, a_ones), e2 : at(e2, a_ones))
else (e1 : at(e1, a_names), e2 : at(e2, a_names))$

/* Finally compute coefficients: */
for m: -1 thru 1 do (for n: -1 thru 1 do (c1u[m,n] : combine(expand(coeff(e1, u[i+m,j+n]) / denominator))));
for m: -1 thru 1 do (for n: -1 thru 1 do (c1v[m,n] : combine(expand(coeff(e1, v[i+m,j+n]) / denominator))));

for m: -1 thru 1 do (for n: -1 thru 1 do (c2u[m,n] : combine(expand(coeff(e2, u[i+m,j+n]) / denominator))));
for m: -1 thru 1 do (for n: -1 thru 1 do (c2v[m,n] : combine(expand(coeff(e2, v[i+m,j+n]) / denominator))));