<< "Calculus`VectorAnalysis`"

<< "Graphics`PlotField`"

phi[r_, ri_] := (q/(Sqrt[DotProduct]))[r - ri, r - ri]

"A charge q=6 is at the point x=5; a sphere of radius r=1 centered at \
the origin is at potential phi=0, find the image charge and its location. \
Plot phi, and electric field."
pot1[x_, y_, z_] := 
      phi[{x, y, z}, {5, 0, 0}] /. {q -> 6}

"Image charge q'=-a/y q = -1/5 6, at y'= 1/6";
pot2[x_, y_, z_] := 
      phi[{x, y, z}, {(1/6), 0, 
            0}] /. {q -> -(6/5)}

tpot[x_, y_, z_] := 
      pot1[x, y, z] + pot2[x, y, z] /; x^2 + y^2 + z^2 > 1
tpot[x_, y_, z_] := 0 /; x^2 + y^2 + z^2 < 1

ContourPlot[pot1[x, y, 0], {x, -2, 8}, {y, -4, 6}, 
      ContourShading -> False, Contours -> 20, 
      PlotPoints -> 30, AspectRatio -> 1]

ContourPlot[pot1[x,y,0]+pot2[x,y,0],{x,-2,8},{y,-4,6},ContourShading->False,
  Contours->20,PlotPoints->50,AspectRatio->1]

Plot3D[tpot[x,y,0], {x,-2,2}, {y,-2,2},PlotPoints->30 ]

" Construct electric field ";
"lets work in three dimensions";
e1[Xx_, Yy_, 
        Zz_] = -Grad[pot1[Xx, Yy, Zz] + pot2[Xx, Yy, Zz]]
" express each component as a vector ";

efx[Xx_, Yy_, 
        Zz_] := e1[Xx, Yy, Zz][[1]] /; 
        Xx^2 + Yy^2 + Zz^2 ≥ 1
efx[Xx_, Yy_, Zz_] := 0 /; (Xx^2) + (Xy^2) + Zz^2 < 1
efy[Xx_, Yy_, 
        Zz_] := e1[Xx, Yy, Zz][[2]] /; 
        (Xx^2) + (Yy^2) + Zz^2 ≥ 1
efy[Xx_, Yy_, Zz_] := 0

PlotVectorField[{efx[Xx, Yy, 0], efy[Xx, Yy, 0]}, {Xx, -2.5, 
        3.0}, {Yy, -2.5, 3.0}, ScaleFactor -> 0.4]

" Charge density "

electricfield = -Grad[
          pot1[Xx, Yy, Zz] + pot2[Xx, Yy, Zz]]

rad = (Sqrt[(Xx^2) + (Yy^2) + (Zz^2)]);
normalvector = {(Xx/rad), (Yy/rad), (Zz/rad)}

sigma = -4 π electricfield . normalvector

" Convert this into spherical coordinates ";
sigma /. {Xx -> Sin[theta] Cos[phi], 
          Yy -> Sin[theta] Sin[phi], 
          Zz -> Cos[theta]};

output = Simplify[%]

Plot3D[output, {theta, 0, π}, {phi, 0, 2 π}]

NIntegrate[Sin[theta] output, {theta, 0, π}, {phi, 0, 2 π}, 
      WorkingPrecision -> 2]

N[6/5 (4 Pi)^2]