############################################################
# 1d wave equation 
# 
#    phi_tt = phi_xx
#
# solved in first order form 
#
#    pp_t  = pi_x
#    pi_t  = pp_x
# 
# where
#
#    pp := phi_x
#    pi := phi_t
# 
# using Crank-Nicholson with (implicit) Kreiss-Oliger 
# dissipation and (implicit) ingoing/outgoing radiation 
# conditions (O(h^2) forwards and backwards differences) at 
# the boundaries.  Initial data currently  
# Gaussian in pp with pi = idsignum * pp.
#
#   idsignum = -1, 0, 1 -> ingoing, time-symmetric, outgoing 
#                          initial data
#
############################################################

# set the memory size (only needed for Fortran)
system parameter int memsiz := 1000000

parameter float xmin       := 0.0
parameter float xmax       := 1.0
parameter float amp        := 1.0
parameter float xc         := 0.5
parameter float xwid       := 0.05
parameter float epsdis     := 0.5
parameter float idsignum   := 0.0

rect coordinates t, x

uniform rect grid g1 [1:Nx] {xmin:xmax}

float pp on g1 at 0,1
float pi on g1 at 0,1

attribute int out_gf encodeone := [1 1 0 0]

# Crank Nicholson time derivative (first forward difference)
operator D(f,t) := (<1>f[0] - <0>f[0]) / dt

# Forward averaging operator
operator MU(f,t) := (<1>f[0] + <0>f[0]) / 2

# Dissipation operator
operator DISS(f,t) := -epsdis / (16 * dt) *
        (6 *<0>f[0] + <0>f[2] + <0>f[-2] - 4 * (<0>f[1] + <0>f[-1]))

# O(h^2) centred, forward and backward first derivatives
operator D0(f,x) := (<0>f[1] - <0>f[-1]) / (2 * dx)

operator DB(f,x) := ( 3*<0>f[0] - 4*<0>f[-1] + <0>f[-2]) / (2 * dx)
operator DF(f,x) := (-3*<0>f[0] + 4*<0>f[1]  - <0>f[2])  / (2 * dx)

# Difference equations
evaluate residual pp { 
               [1:1] := D(pp,t)  =  MU(DF(pp,x),t);
               [2:2] := D(pp,t)  =  MU(D0(pi,x),t);
            [3:Nx-2] := D(pp,t)  =  MU(DISS(pp,t) + D0(pi,x),t);
         [Nx-1:Nx-1] := D(pp,t)  =  MU(D0(pi,x),t); 
             [Nx:Nx] := D(pp,t)  =  MU(-DB(pp,x),t);
                     }

evaluate residual pi { 
               [1:1] := D(pi,t) =  MU(DF(pi,x),t);
               [2:2] := D(pi,t) =  MU(D0(pp,x),t);
            [3:Nx-2] := D(pi,t) =  MU(DISS(pi,t) + D0(pp,x),t);
         [Nx-1:Nx-1] := D(pi,t) =  MU(D0(pp,x),t);
             [Nx:Nx] := D(pi,t) =  MU(-DB(pi,x),t);
                     }

# Intialize to ingoing/outgoing/time-symmetric gaussian pulse
initialize pp {[1:Nx] := amp*exp(-(x-xc)^2/xwid^2)}
initialize pi {[1:Nx] := idsignum * amp*exp(-(x-xc)^2/xwid^2)}

looper iterative

auto update pp, pi