########################################################################
# RNPL source for axisymmetric wave-equation in cylindrical coordinates
# (rho,z).  First order form and single uniform mesh (including (0,0)
# used.  Radiation conditions imposed on outer-rho and both z-boundaries.
# Regularity imposed at rho=0.  Scheme unstable, instability probably 
# due to bad boundary/regularity conditions, particularly at (0,zmin),
# (0,zmax)
#
# Second order form: 
#    phi_tt = phi_zz + phi_rhorho + phi_rho / rho
# 
# ########################################################################
# RNPL source for axisymmetric wave-equation in cylindrical coordinates
# (rho,z).  First order form and single uniform mesh (including (0,0)
# used.  Radiation conditions imposed on outer-rho and both z-boundaries.
# Regularity imposed at rho=0.  Scheme unstable, instability probably 
# due to bad boundary/regularity conditions, particularly at (0,zmin),
# (0,zmax)
#
# Second order form: 
#    phi_tt = phi_zz + phi_rhorho + phi_rho / rho
# 
# First  order form:  (this code)
#    Pi := phi_t
#    Phirho := phi_rho
#    Phiz := phi_z
#
#    Pi_t := Phiz_z + Phirho_rho + Phirho / rho
#    Phirho_t := Pi_rho
#    Phiz_t := Pi_z
#
# Ingoing/outgoing/time symmetric initial data: signum = +1/-1/0
#
#    Pi  := signum * (Phi + rho * Phi_rho + z * Phi_z)/r;
# 
# Crank-Nicholson (time) differencing
#
########################################################################
#
# Copyright 1996, Matthew W Choptuik, The University of Texas at Austin
########################################################################


########################################################################
# Parameters, grid and grid functions
########################################################################

# This is how to set the memory size
system parameter int memsiz := 2000000

parameter float rhomin 
parameter float rhomax
parameter float zmin 
parameter float zmax
parameter float amp
parameter float delta
parameter float r0
parameter float epsdis := 0
parameter float signum := 0

rec coordinates t,rho,z

uniform rec grid g1 [1:Nrho][1:Nz] {rhomin:rhomax} {zmin:zmax}

float Phi    on g1 at 0,1
float pi     on g1 at 0,1
float Phirho on g1 at 0,1
float Phiz   on g1 at 0,1
float Phirhobyrho  on g1 at  0,1
float r      on g1 at 0

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

########################################################################
# Difference operators 
########################################################################
operator D_CN(f,t) := (<1>f[0][0] - <0>f[0][0]) / dt
operator MU(f,t)   := (<1>f[0][0] + <0>f[0][0]) / 2
operator D_0(f,rho) := (<0>f[1][0] - <0>f[-1][0])/(2 * drho)
operator D_0(f,z) := (<0>f[0][1] - <0>f[0][-1])/(2 * dz)
operator D_ADV0(f,z) := (<1>f[0][1] - <1>f[0][-1])/(2 * dz)
operator D_FW(f,rho) := (-3 * <0>f[0][0] + 4 * <0>f[1][0]  - <0>f[2][0])/(2 * drho)
operator D_BW(f,rho) := ( 3 * <0>f[0][0] - 4 * <0>f[-1][0] + <0>f[-2][0])/(2 * drho)
operator D_FW(f,z  ) := (-3 * <0>f[0][0] + 4 * <0>f[0][1]  - <0>f[0][2])/(2 * dz)
operator D_BW(f,z  ) := ( 3 * <0>f[0][0] - 4 * <0>f[0][-1] + <0>f[0][-2])/(2 * dz)
operator D_CN0(f,rho) := MU(D_0(<0>f[0][0],rho),t)
operator D_CN0(f,z) := MU(D_0(<0>f[0][0],z),t)
operator D_CNFW(f,rho) := MU(D_FW(<0>f[0][0],rho),t)
operator D_CNFW(f,z) := MU(D_FW(<0>f[0][0],z),t)
operator D_CNBW(f,rho) := MU(D_BW(<0>f[0][0],rho),t)
operator D_CNBW(f,z) := MU(D_BW(<0>f[0][0],z),t)

operator QFIT(f,rho) := (<1>f[0][0] - 4 * <1>f[1][0]/3 + <1>f[2][0]/3)
operator D_ADVAVG(f,z  ) := (-<1>f[0][0] + (<1>f[0][1] + <1>f[0][-1])/2)
operator D_FEX(f,rho) := (-<0>f[0][0] + 2 * <0>f[1][0] - <0>f[2][0])

########################################################################
# Residual definitions (equations of motion)
########################################################################

evaluate residual pi { 
   [2:Nrho-1][1:1] := <0>r[0][0] * D_CN(pi,t) + rho * D_CN0(pi,rho) + z * D_CNFW(pi,z);
   [Nrho:Nrho][1:1] := <0>r[0][0] * D_CN(pi,t) + rho * D_CNBW(pi,rho) + z * D_CNFW(pi,z);

   [2:Nrho-1][Nz:Nz] := <0>r[0][0] * D_CN(pi,t) + rho * D_CN0(pi,rho) + z * D_CNBW(pi,z);
   [Nrho:Nrho][Nz:Nz] := <0>r[0][0] * D_CN(pi,t) + rho * D_CNBW(pi,rho) + z * D_CNBW(pi,z);

   [1:1][1:Nz] := QFIT(pi,rho);

#  [1:Nrho][2:2] := D_ADVAVG(pi,z);
#  [1:Nrho][Nz-1:Nz-1] := D_ADVAVG(pi,z);

   [Nrho:Nrho][1:1] := <0>r[0][0] * D_CN(pi,t) + rho * D_CNBW(pi,rho) + z * D_CNFW(pi,z);
   [Nrho:Nrho][2:Nz-1] := <0>r[0][0] * D_CN(pi,t) + rho * D_CNBW(pi,rho) + z * D_CN0(pi,z);
   [Nrho:Nrho][Nz:Nz] := <0>r[0][0] * D_CN(pi,t) + rho * D_CNBW(pi,rho) + z * D_CNBW(pi,z);

   [2:Nrho-1][2:Nz-1] := D_CN(pi,t) - D_CN0(Phirho,rho) - D_CN0(Phiz,z) - MU(Phirho,t) / rho;
}

evaluate residual Phirho { 
   [2:Nrho-1][1:Nz] := D_CN(Phirho,t) - D_CN0(Pi,rho);
   [1:1][1:Nz] := <1>Phirho[0][0];
   [Nrho:Nrho][1:Nz] := D_CN(Phirho,t) - D_CNBW(Pi,rho);
}

residual Phirhobyrho { 
   [2:Nrho][1:Nz] := <1>Phirhobyrho[0][0] = <1>Phirho[0][0]/rho;
   [1:1][1:Nz] := QFIT(Phirhobyrho,rho);
}

residual Phi {
   [1:Nrho][1:Nz] := D_CN(Phi,t) = Pi;
}

evaluate residual Phiz { 
   [2:Nrho][2:Nz-1] := D_CN(Phiz,t) - D_CN0(Pi,z);
   [2:Nrho][1:1] := D_CN(Phiz,t) - D_CNFW(Pi,z);
   [2:Nrho][Nz:Nz] := D_CN(Phiz,t) - D_CNBW(Pi,z);
   [1:1][1:Nz] := QFIT(Phiz,rho);
}

########################################################################
# Initializations and update structure 
########################################################################

initialize r      { [1:Nrho][1:Nz] := sqrt(rho^2 + z^2) }

initialize Phi    { [1:Nrho][1:Nz] := amp * exp(-((sqrt(rho^2 + z^2)-r0)/delta)^2) }

initialize Pi     { [2:Nrho-1][2:Nz-1] := signum * (Phi + 
                       rho * D_0(Phi,rho) + z * D_0(Phi,z))/r;
                    [1:1][1:Nz] := D_FEX(Pi,rho);
                    [Nrho:Nrho][1:Nz] := 0 
}

initialize Phirho { [2:Nrho-1][1:Nz] := D_0(Phi,rho);
                    [1:1][1:Nz] := 0;
                    [Nrho:Nrho][1:Nz] := 0
}

initialize Phiz   { [1:Nrho][2:Nz-1] := D_0(Phi,z);
                    [1:Nrho][1:1] := 0;
                    [1:Nrho][Nz:Nz] := 0
}

looper iterative

auto update pi, Phi, Phirho, Phirhobyrho, Phiz