kill(all);

load("myhelpers.mac");

/* redefine this to change dimensionality: */
n:[nx,ny]; 

dims:length(n);

assume(c>0);

if dims = 1 then n:[1];

esymmatrix(n, v, i,j):=ematrix(n,n,v,i,j)+ematrix(n,n,v,j,i);
wave_A:sum(n[i]*esymmatrix(dims+1, -c, 1+i,1),i,1,dims);

[wave_V, wave_D, wave_invV]:hypdiagonalize(wave_A);

wave_vm:makelist(concat(vm,i),i,1,dims);
wave_vp:makelist(concat(vp,i),i,1,dims);
wave_wm:append([um],wave_vm);
wave_wp:append([up],wave_vp);

wave_sm:makelist(concat(sm,i),i,1,length(wave_D));
wave_sp:makelist(concat(sp,i),i,1,length(wave_D));

wave_sminw:wave_invV.wave_wm;
wave_spinw:wave_invV.wave_wp;

wave_wmins:wave_V.wave_sm;
wave_wpins:wave_V.wave_sp;

wave_radbdryspinw:makelist( 
  if wave_D[i,i] >= 0 then wave_sminw[i,1] else 0,
  i, 1, length(wave_D));
wave_radbdrywp:fullhypsimp(wave_V.wave_radbdryspinw);

wave_dirbdryspinw:makelist( 
  if wave_D[i,i] >= 0 then wave_sminw[i,1] else ubc,
  i, 1, length(wave_D));
wave_dirbdrywp:fullhypsimp(wave_V.wave_dirbdryspinw);

print("Radiation boundary condition for the wave equation:");
print(wave_radbdrywp);

print("Dirichlet boundary condition for the wave equation:");
print(expand(wave_dirbdrywp));

wave_known_dirbdrywp:vstack([-n.wave_vm/2 + um/2 + ubc],
  -n*um/2 + n*ubc - 1/2*n*(n.wave_vm) + wave_vm);

assert(norm_2_squared(fullhypsimp(wave_known_dirbdrywp - wave_dirbdrywp))=0);

print("Homogeneous-dirichlet in characteristic:");
print(fullhypsimp(wave_invV.vstack(
    -columnvector(first(wave_wmins)), 
    rest(wave_wmins))
    ));

print("Homogeneous-Neumann in characteristic:");
print(fullhypsimp(wave_invV.vstack(
    columnvector(first(wave_wmins)), 
    -rest(wave_wmins))
    ));

/* ------------------------------------------------------------------------- */

wave_eigenvalues:makelist(wave_D[i,i], i, 1, length(wave_D));

if member(0, wave_eigenvalues) then
  wave_sflux:hyp_upwind_flux([-c,0,c], wave_D)
else
  wave_sflux:hyp_upwind_flux([-c,c], wave_D);
wave_wflux:ratsimp(wave_V.ev(wave_sflux, [sm=wave_sminw,sp=wave_spinw]));
wave_strongwflux:fullhypsimp(subst([c=cm], wave_A).wave_wm - wave_wflux);

/*
print("Wave equation flux in terms of characteristic variables:");
print(wave_sflux);


print("Weak flux divided by (-c), as implemented in StrongWaveOperator:");
print(hypsimp(ev(wave_wflux, cp=c, cm=c)/(-c)));
*/
print("Wave equation weak flux in terms of physical variables:");
print(wave_wflux);

print("Strong-form wave equation flux in terms of physical variables:");
/*print(wave_strongwflux);*/

print(fullhypsimp(ev(wave_strongwflux)));

/* Closed-form expression for upwind flux */
wave_knownstrongwflux:vstack(
  [(cp*n.wave_vp - cm*n.wave_vm + cp*up-cm*um)/2],
  (cp*n*(n.wave_vp)-cm*n*(n.wave_vm) + cp*n*up - cm*n*um)/2
  );

assert(norm_2_squared(fullhypsimp(wave_knownstrongwflux - wave_strongwflux))=0);