length = 1.0e-10 
level = 10
dx = length/2^level
dx_half = dx*0.5
nLength = math.ceil(length/dx)
print(nLength)
sph_pos = {0.0,0.0,0.0}
nRad = 10
radius = nRad*dx
Dia = radius*2.0
nDia = nRad*2

dt =dx*dx
ref_pot = 0.0 
permit = 8.854187e-12 
moleDens0 = 1e-3
valence_sqr = 1
charge = 1.60217657e-19
k_b = 1.3805e-23
temp = 273
N_A = 6.02e23 

spr_pos = {-length/2.0,-length/2.0,0.0}
spr_pos_1 = {-length/8.0,0.0,0.0}
spr_pos_1 = {0.0,0.0,0.0}
spr_pos_2 = {length/8.0,0.0,0.0}
Q = 1
sig = radius/3.0
-- volume of charged sphere
volume = 4.0*math.pi*radius^3/3.0

function chargedens_sphere(x,y,z,t)
  tmp = math.exp(-r^2/(2*sig^2))/(sig^3*(2.0*math.pi)^(3.0/2.0))
  return Q*charge * tmp
  --r = math.sqrt((x-spr_pos[1])^2+(y-spr_pos[2])^2+(z-spr_pos[3])^2)
  --if (r<=radius) then
  --  return Q * charge / volume 
  --else
  --  return 0.0
  --end 
end

function analy_pot(x,y,z,t)  
  r = math.sqrt((x-spr_pos[1])^2+(y-spr_pos[2])^2+(z-spr_pos[3])^2)
  if (r <= radius ) then
    return 0
  else
    return fac/r
  end
end  

function chargedens_sphere_dipole(x,y,z,t)
  r1 = math.sqrt((x-spr_pos_1[1])^2+(y-spr_pos_1[2])^2+(z-spr_pos_1[3])^2)
  r2 = math.sqrt((x-spr_pos_2[1])^2+(y-spr_pos_2[2])^2+(z-spr_pos_2[3])^2)
  if (r1<=radius) then
    return Q * charge / volume 
  --elseif (r2<=radius) then
  --  return -Q * charge / volume 
  else
    return 0.0
  end 
  --tmp = math.exp(-r^2/(2*sig^2))/(sig^3*(2.0*math.pi)^(3.0/2.0))
  --return Q*charge * tmp
end

-- Parameter for analytical solution
fac = Q*charge/(4.0*math.pi*permit)
function analy_pot_dipole(x,y,z,t)  
  r1 = math.sqrt((x-spr_pos_1[1])^2+(y-spr_pos_1[2])^2+(z-spr_pos_1[3])^2)
  r2 = math.sqrt((x-spr_pos_2[1])^2+(y-spr_pos_2[2])^2+(z-spr_pos_2[3])^2)
  if (r1 ~= 0 and r2 ~= 0) then
    --return fac/r1 - fac/r2
    return fac/r1
  else
    return 0
  end  
end


