# ########################################################################
# 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.  
#
# 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 + 2*(d/d rho^2)( rho*Phirho )
#    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
# 
# Copyright 1997, Modifications by Steve Liebling, 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 := 1
parameter float rho_center := 10
parameter float z_center   := 10
parameter float a := 10
parameter float b := 10

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 r      on g1 at 0

attribute int out_gf encodeone :=  [ 1 0 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])

#added by sll   ---> ( d / d rho^2 )
operator D_rhosq(f,rho)   := ( rho[1]*<0>f[1][0] - rho[-1]*<0>f[-1][0] ) 
															    / 
									  ( rho[1]^2          - rho[-1]^2           )
operator D_rhosqCN(f,rho) := MU(D_rhosq(<0>f[0][0],rho),t)

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

evaluate residual pi { 
#sll---added MU(Pi,t) term for proper outgoing radiation condition
   [2:Nrho-1][2:Nz-1]  := D_CN(pi,t)  - D_CN0(Phiz,z) - 2.0*D_rhosqCN(Phirho,rho);

   [2:Nrho-1][1:1]     := r * D_CN(pi,t) + rho * D_CN0(pi,rho)  + z * D_CNFW(pi,z) + MU(Pi,t);
   [Nrho:Nrho][1:1]    := r * D_CN(pi,t) + rho * D_CNBW(pi,rho) + z * D_CNFW(pi,z) + MU(Pi,t);
   [2:Nrho-1][Nz:Nz]   := r * D_CN(pi,t) + rho * D_CN0(pi,rho)  + z * D_CNBW(pi,z) + MU(Pi,t);
   [Nrho:Nrho][Nz:Nz]  := r * D_CN(pi,t) + rho * D_CNBW(pi,rho) + z * D_CNBW(pi,z) + MU(Pi,t);
   [Nrho:Nrho][2:Nz-1] := r * D_CN(pi,t) + rho * D_CNBW(pi,rho) + z * D_CN0(pi,z)  + MU(Pi,t);
   [1:1][Nz:Nz]        := r * D_CN(pi,t) + rho * D_CNFW(pi,rho) + z * D_CNBW(pi,z) + MU(Pi,t);
   [1:1][1:1]          := r * D_CN(pi,t) + rho * D_CNFW(pi,rho) + z * D_CNFW(pi,z) + MU(Pi,t);
   [1:1][2:Nz-1]       := QFIT(pi,rho);
}
evaluate residual Phirho { 
   [1:1][1:Nz]       := <1>Phirho[0][0];
   [2:Nrho-1][1:Nz]  := D_CN(Phirho,t) - D_CN0(Pi,rho);
   [Nrho:Nrho][1:Nz] := D_CN(Phirho,t) - D_CNBW(Pi,rho);
}
evaluate residual Phi {
#  [1:Nrho][1:Nz]    := D_CN(Phi,t) = Pi;
   [1:Nrho][1:Nz]    := D_CN(Phi,t) = MU(Pi,t);
}

evaluate residual Phiz { 
# sll--eliminated QFIT at rho=1 and extended domain of below from rho=2 to rho=1
   [1:Nrho][2:Nz-1]  := D_CN(Phiz,t) - D_CN0(Pi,z);
   [1:Nrho][1:1]     := D_CN(Phiz,t) - D_CNFW(Pi,z);
   [1:Nrho][Nz:Nz]   := D_CN(Phiz,t) - D_CNBW(Pi,z);
}

########################################################################
# 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-rho_center)/a)^2 + ((z-z_center)/b)^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, Phiz