########################################################################
# RNPL source for axisymmetric wave-equation in cylindrical coordinates
# (rho,z).  Second 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 seems to be stable
#
# Second order form: 
#    phi_tt = phi_zz + phi_rhorho + phi_rho / rho
# 
# 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  z0 := 0.5
parameter float  rho0 := 1
parameter float  epsdis := 0

########################################################################
# parameter 'idstr' defines the initial data type.  Currently implemented:
#
# idstr := 'ts'         Time-symmetric gaussian
# idstr := 'in'         Ingoing gaussian
# idstr := 'out'        Outgoing gaussian
# idstr := 'gts'        Distorted time symmetric gaussian: add'l params:
#   z0                  z-scale
#   rho0                rho-scale 
#
########################################################################
parameter string idstr  := "gts"

########################################################################
# Explicitly define the 'out_gf' attribute and give a default value 
# which enables(1)/disables(0) output of the corresponding grid function 
# at times specified by the vector parameter 'output'.  Unfortunately,
# you still have to keep track of which elements of 'out_gf' to set to 1
# to get the desired output.  The file '.rnpl.attributes', automatically
# generated by the RNPL compiler, can help, and can also be modified to 
# set 'out_gf' arbitrarily without re-compilation.
########################################################################
attribute int out_gf encodeone := [1 0 0]

rec coordinates t,rho,z

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

float Phi          on g1 at -1,0,1
float r            on g1 at 0

########################################################################
# Difference operators 
########################################################################

operator D_FW(f,t)   := (<1>f[0][0] -  <0>f[0][0])/dt
operator D_BW(f,t)   := (<0>f[0][0] - <-1>f[0][0])/dt
operator D_LF(f,t,t) := (<1>f[0][0] - 2*<0>f[0][0] + <-1>f[0][0])/(dt^2)
operator D_LF(f,rho,rho) := (<0>f[1][0] - 2*<0>f[0][0] + <0>f[-1][0])/(drho^2)
operator D_LF(f,z,z) := (<0>f[0][1] - 2*<0>f[0][0] + <0>f[0][-1])/(dz^2)
operator D_LF(f,t)   := (<1>f[0][0] - <-1>f[0][0])/(2*dt)
operator D_LF(f,rho) := (<0>f[1][0] - <0>f[-1][0])/(2*drho)
operator D_LF(f,z)   := (<0>f[0][1] - <0>f[0][-1])/(2*dz)

operator D_BW2(f,t)   := ( 3*<1>f[0][0] - 4*<0>f[0][0]  + <-1>f[0][0])/(2*dt)
operator D_FW2(f,rho) := (-3*<1>f[0][0] + 4*<1>f[1][0]  - <1>f[2][0])/(2*drho)
operator D_BW2(f,rho) := ( 3*<1>f[0][0] - 4*<1>f[-1][0] + <1>f[-2][0])/(2*drho)
operator D_CADV(f,rho) := ( <1>f[1][0] - <1>f[-1][0])/(2*drho)
operator D_FW2(f,z  ) := (-3*<1>f[0][0] + 4*<1>f[0][1]  - <1>f[0][2])/(2*dz)
operator D_BW2(f,z  ) := ( 3*<1>f[0][0] - 4*<1>f[0][-1] + <1>f[0][-2])/(2*dz)
operator D_CADV(f,z  ) := ( <1>f[0][1] - <1>f[0][-1])/(2*dz)

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

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

evaluate residual Phi {
   [1:1]      [1:Nz]    :=  QFIT(Phi,rho);
   [Nrho:Nrho][2:Nz-1]  :=   
       <0>r[0][0] * D_BW2(Phi,t) + rho *  D_BW2(Phi,rho) + z * D_CADV(Phi,z);
   [2:Nrho-1][2:Nz-1]   := 
       D_LF(Phi,t,t) = D_LF(Phi,z,z) + D_LF(Phi,rho,rho) + D_LF(Phi,rho) / rho;
   [2:Nrho-1][ 1: 1]      := 
       <0>r[0][0] * D_BW2(Phi,t) + rho *  D_CADV(Phi,rho) + z * D_FW2(Phi,z);
   [2:Nrho-1][Nz:Nz]      :=  
       <0>r[0][0] * D_BW2(Phi,t) + rho *  D_CADV(Phi,rho) + z * D_BW2(Phi,z);
}

########################################################################
# Initializations and update structure 
# 
# Note: RNPL generated initialization routine will generate only
# time symmetric data.  
########################################################################

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) }

looper iterative

auto update Phi