VC.10 3D Surface Measurements
TUTORIALS

\(Off[General::spell]; 
Off[General::spell1]; 

    \(Off[ParametricPlot3D::ppcom]; \)\)Off[ParametricPlot::ppcom]; Off[Plot::plnr]; SetOptions[Limit, Analytic -> True];

CMView = {2.7, 1.6, 1.2};

$Post  =  Chop;

Vector Drawer

BeginPackage["Vector3D`", "Graphics`Colors`"];

\(Axes3D::"usage" = 
      "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8]."; 

    \(Perpend::"usage" = 
      "Perpend[a] returns a unit vector which is perpendicular to the vector a."; \)\)Arrow::"usage" = 
      "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. \
Arrow[a, Tail -> tail] gives a vector from tail to a."; Vector::"usage" = 
      "Vector[a, b] gives a 2D or 3D vector from a to b."; ArrowHead::"usage" = 
      "ArrowHead[a, vec] produces an arrowhead with its tip placed at the \
point a, pointing in the direction of the vector vec."; Tail::"usage" = 
      "Tail -> point puts the tail of the vector at point."; Norm::"usage" = 
      "Norm[vec] returns the Euclidean length of vec."; Normalize::"usage" = 
      "Normalize[vec] returns a unit vector in the direction of vec, \
provided vec is nonzero."; Aperture::"usage" = 
      "Aperture is the ratio of the radius of the base to the length of the \
head of a vector."; TipSize::"usage" = 
      "TipSize specifies an absolute size for the inner tip length of the \
head of a vector."; TipRatio::"usage" = 
      "TipRatio is the ratio of the inner length of the arrowhead to the \
length of the shaft of a vector."; HeadRatio::"usage" = 
      "HeadRatio is the ratio of the outer length of the arrowhead to the \
length of the shaft of a vector."; HeadSize::"usage" = 
      "HeadSize specifies an absolute size for the head of a vector."; 
    ScaleFactor::"usage" = 
      "ScaleFactor specifies the amount to scale a vector in length and may \
be either a number or a pure function."; ZeroVectorPointSize::"usage" = 
      "ZeroVectorPointSize is the size of the point used to represent a \
zero vector."; Options[Vector] = 
      Options[Arrow] = 
        Options[ArrowHead] = {HeadRatio -> 0.18, 
            TipRatio -> 0.14, Aperture -> 0.3, 
            ScaleFactor -> 1, HeadType -> Polygon, 
            ShaftQ -> True, EdgesQ -> True, 
            VectorColor -> RGBColor[0, 0, 1], ShaftWidth -> 0.005, 
            ZeroVectorPointSize -> 0.01}; Begin["`Private`"]; SetOptions[ParametricPlot, AspectRatio -> Automatic]; SetOptions[Plot, AspectRatio -> Automatic]; SetOptions[Graphics, AspectRatio -> Automatic]; Axes3D[u_, v_] := 
      Graphics3D[{{Blue, Line[{{-(u/3), 0, 0}, {u, 0, 0}}]}, 
          Text["x", {u + v, 0, 0}], {Blue, 
            Line[{{0, -(u/3), 0}, {0, u, 0}}]}, 
          Text["y", {0, u + v, 0}], {Blue, 
            Line[{{0, 0, -(u/3)}, {0, 0, u}}]}, 
          Text["z", {0, 0, u + v}]}]; Axes3D[u_] := Axes3D[u, (u/8)]; Norm[a_?VectorQ] := N[(Sqrt[a . a])]; Normalize[a_?VectorQ] := If[a . a == 0, a, (a/Norm)[a]]; Perpend[a_?VectorQ] := 
      Normalize[{a[[2]], 
            -a[[1]]}] /; 
        Length[a] == 2Perpend[a_?VectorQ] := 
      Normalize[
          If[a[[1]] == 0, {1, 0, 0}, {0, 
              a[[3]], 
              -a[[2]]}]] /; 
        Length[a] == 3Base3D = 
      Table[{Cos[((2 k π)/8.) ], Sin[((2 k π)/8.) ]}, {k, 0, 
          8}]; Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = 
            {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, 
                  TipRatio, Aperture, ScaleFactor, VectorColor, 
                  ZeroVectorPointSize} /. {opts} /. Options[Vector]; 
          If[from == to, 
            Return[Graphics[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]; tr = 0.75 hr]; 
          If[NumberQ[ts], tr = (ts/Norm)[diff]]; 
          Graphics[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              With[{trans = ap hr Norm[diff] Perpend[diff]}, {vc, 
                  Thickness[sw], 
                  ht[{tip, tip - hr diff + trans, tip - tr diff, 
                      tip - hr diff - trans, tip}]}]}]] /; 
        Length[from] == Length[to] == 2Vector[from_?VectorQ, to_?VectorQ, opts___] := 
      Module[{sq, eq, sw, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, 
            tip}, {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = 
            {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, 
                  ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts}
                   /. Options[Vector]; 
          If[from == to, 
            Return[Graphics3D[{PointSize[zvps], vc, Point[from]}]]]; 
          diff = If[NumberQ[N[sf]], sf diff, sf[diff]]; tip = from + diff; 
          If[NumberQ[hs], hr = (hs/Norm)[diff]]; tr = (hr/2); 
          Graphics3D[{
              If[sq, {vc, Thickness[sw], Line[{from, tip - tr diff}]}, {}], 
              If[eq, {}, EdgeForm[]], {SurfaceColor[vc], 
                With[{temp = Perpend[diff]}, 
                  (Polygon[Append[#1, tip]]&)/@
                    Partition[
                      (tip - hr diff + #1&)/@
                        ((ap hr Norm[diff] Base3D) . {temp, 
                              Normalize[temp\[Cross]diff]}), 2, 1]]}}]] /; 
        Length[from] == Length[to] == 3Arrow[a_, Tail -> b_, opts___] := Vector[b, a + b, opts]Arrow[a_, opts___] := Vector[Table[0, {Length[a]}], a, opts]\(ArrowHead[a_, b_, opts___] := 
      Vector[a - b, a, ShaftQ -> False, opts]; 

End[]; 

    \(EndPackage[]; \)\)

Graphics`Colors`GosiaGreen = RGBColor[0,  0.392187,  0];

Spherical Setup

\(Clear[polar]; 

    polar[{x_, y_, z_}] := {\@\(x\^2 + y\^2 + z\^2\), ArcTan[y\/x], 
        ArcCos[z\/\@\(x\^2 + y\^2 + z\^2\)]}; 
Clear[X, length]; 

    length[X_] := \@\(X . X\); 
\(rr = 3; \)\)R = 5; scale = 0.8; red = RGBColor[1, 0, 0]; blue = RGBColor[0, 0, 1]; color1 = RGBColor[0.8, 0.2, 0.8]; color2 = RGBColor[0.2, 0.8, 0.2]; black = RGBColor[0, 0, 0]; \(Clear[axes]; 

    axes := Graphics3D[{Line[{X, {1.5\ rr, 0, 0}}], 
          Line[{X, {0, 1.5\ rr, 0}}], Line[{X, {0, 0, 2\ rr}}]}]; 

    Clear[xyplane]; 

    xyplane := 
      Graphics3D[
        Polygon[{{\(-2\)\ rr, \(-2\)\ rr, 0}, {\(-2\)\ rr, 2\ rr, 0}, {2\ rr, 
              2\ rr, 0}, {2\ rr, \(-2\)\ rr, 0}}]]; 
\(X = {0, 0, 0}; \)\)origin = Graphics3D[{black, Point[X]}]; Z = {2, 3, 5}; \(Clear[point]; 

    point := Graphics3D[{red, PointSize[0.03], Point[scale\ Z]}]; 

    \(Z0 = {0, 0, length[0.8\ Z]}; \)\)\(Clear[zpoint]; 
zpoint = Graphics3D[{black, Point[Z0]}]; 

    \(XY = scale\ Z\ {1, 1, 0}; \)\)\(Clear[xypoint]; 
xypoint := Graphics3D[{black, Point[XY]}]; 

    \(X0 = {length[scale\ Z\ {1, 1, 0}], 0, 0}; \)\)xpoint = Graphics3D[{black, Point[X0]}]; Y0 = {0, length[scale Z {1, 1, 0}], 0}; ypoint = Graphics3D[{black, Point[Y0]}]; \(Clear[ray]; 
ray := Graphics3D[{blue, Line[{X, scale\ Z}]}]; 

    Clear[ray1]; 
ray1 := Graphics3D[{color1, Line[{scale\ Z, Z}]}]; 

    Clear[xyray]; 
xyray := Graphics3D[{red, Line[{X, Z\ {1, 1, 0}}]}]; 

    \(rays := {ray, xyray}; \)\)xyaux = 
      (length[0.8 Z] (0.8 Z {1, 1, 0}))/
        length[0.8 Z {1, 1, 0}]; \(Clear[zarc]; 

    zarc = Graphics3D[{blue, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 0, 
                \(polar[0.8\ Z]\)[[3
                  ]],  .1}]]}]; 
Clear[zarc1]; 

    zarc1 = Graphics3D[{color1, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 
                \(polar[0.8\ Z]\)[[3]], 
                1,  .1}]]}]; 
Clear[zarc2]; 

    zarc2 = Graphics3D[{color1, 
          Line[Table[
              \(Cos[t]\ \) {0, 0, length[0.8\ Z]} + Sin[t]\ xyaux, {t, 
                \(-0.8\), 0,  .1}]]}]; 
Clear[xyarc]; 

    xyarc = Graphics3D[{blue, 
          Line[Table[
              Cos[t]\ X0 + Sin[t]\ Y0, {t, 0, 
                \(polar[XY]\)[[2
                  ]],  .1}]]}]; 
Clear[xyarc1]; 

    xyarc1 = Graphics3D[{color1, 
          Line[Table[Cos[t]\ X0 + Sin[t]\ Y0, {t, \(- .5\), 0,  .1}]]}]; 

    Clear[xyarc2]; 

    xyarc2 = Graphics3D[{color1, 
          Line[Table[
              Cos[t]\ X0 + Sin[t]\ Y0, {t, 
                \(polar[XY]\)[[2]], 
                \(polar[XY]\)[[2]] + 
                  0.3,  .1}]]}]; 
Clear[stick]; 

    stick = Graphics3D[{red, Line[{scale\ Z, XY}]}]; 

    \(xaxislabel = Text["x", {0.363, 0.26}]; \)\)yaxislabel = Text["y", {0.73, 0.3}]; zaxislabel = Text["z", { .5, 0.67}]; stickplot = 
      Show[xyplane, axes, ray, point, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; rho = Text["r", {0.53, 0.45}]; rlabel = 
      Show[xyplane, axes, ray, point, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {rho, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; phi = Text["s", {0.55, 0.57}]; slabel = 
      Show[xyplane, zpoint, axes, ray, point, zarc, zarc1, zarc2, 
        Boxed -> False, ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {phi, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; theta = Text["t", { .5, 0.3}]; tlabel = 
      Show[xyplane, axes, xypoint, xpoint, ray, stick, point, xyarc, xyarc1, 
        xyarc2, Boxed -> False, ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {theta, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity]; all = 
      Show[xyplane, axes, ray, stick, point, zpoint, xpoint, xypoint, xyarc, 
        xyarc1, xyarc2, zarc, zarc1, zarc2, Boxed -> False, 
        ViewPoint -> {2.747, 1.594, 1.168}, 
        Epilog -> {theta, phi, rho, xaxislabel, yaxislabel, zaxislabel}, 
        DisplayFunction -> Identity];

Crisscross

\(Clear[x]; 

    crisscross = 
      Plot[{Cos[x], Sin[x]}, {x, π\/8, \(3\ π\)\/8}, 
        PlotStyle -> {{RGBColor[1, 0, 0], Thickness[0.02]}}, 
        Axes -> None, AspectRatio -> 1\/3, 
        DisplayFunction -> Identity]; \)

T.1) Measuring flow across surfaces: 
Gauss's 3D formula versus calculation by surface integrals

T.1.a.i) The situation suggests a calculation by Gauss's 3D formula

Here's a surface:

\(Clear[x, y, z, r, s, t]; 
    
{x[r_, s_, t_], y[r_, s_, t_], z[r_, s_, t_]} = {
        0.6\ r\ Sin[s]\ Cos[t], r\ Sin[s]\ Sin[t], 2\ r\ Cos[s]}; 

    \({{slow, shigh}, {tlow, thigh}} = {{0, π}, {0, 2\ π}}; \)\)Clear[surfaceplotter]; surfaceplotter[s_, 
        t_] = {0, 1, 0} + {x[3, s, t], y[3, s, t], z[3, s, t]}; surface = 
      ParametricPlot3D[
        Evaluate[surfaceplotter[s, t]], {s, slow, shigh}, {t, tlow, thigh}, 
        Boxed -> False, BoxRatios -> Automatic, 
        ViewPoint -> CMView, 
        AxesLabel -> {"x", "y", "z"}];

The surface is the skin of a solid egg.
Measure the flow of 
       Field[x, y, z] = {-4 x, (y^2), 2  z} 
across this surface. 
Determine whether the net flow is from inside to outside or from outside to inside.

This is a natural for Gauss's 3D formula because the surface is all of the outside skin of a solid region.
Enter the vector field:

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {\(-4\)\ x, y\^2, 
        2\ z}; 
Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

Calculate divField[x, y, z]:

\(Clear[divField]; 

    divField[x_, y_, z_] = 
      D[m[x, y, z], x] + D[n[x, y, z], y] + D[p[x, y, z], z]\)

Call the surface you see above C and call the solid consisting of all points inside and on the surface R.  Gauss's formula tells you that flow of Field[x, y, z] across C is 
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            divField[x, y, z] dx dy 
              dz\)\)\).
To calculate this integral, take advantage of the parameterization of the surface.  Checking the plotting instructions above you see that R is described by

\(Clear[x, y, z, r, s, t]; 
    
{x[r_, s_, t_], y[r_, s_, t_], z[r_, s_, t_]} = {
        0.6\ r\ Sin[s]\ Cos[t], r\ Sin[s]\ Sin[t], 2\ r\ Cos[s]}\)

With:

{{rlow, rhigh}, {slow, shigh}, {tlow, thigh}} = {{0, 3}, {0, π}, {
          0, 2 π}}

Yyou know that the net flow of Field[x, y, z] across C is measured by 
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            divField[x, y, z] dx dy 
              dz\)\)\)
        = 
        (∫_rlow^rhigh)
            (∫_slow^shigh)
              (∫_tlow
                ^thigh) divField[x[r, s, t], y[r, s, t], z[r, s, t]] 
                (Subscript[V, xyz])[r, s, t] dt ds 
                dr
Unleash Mathematica:

\(Clear[gradx, grady, gradz, Vxyz]; 

    gradx[r_, s_, t_] = {D[x[r, s, t], r], D[x[r, s, t], s], 
        D[x[r, s, t], t]}; 

    grady[r_, s_, t_] = {D[y[r, s, t], r], D[y[r, s, t], s], 
        D[y[r, s, t], t]}; 

    gradz[r_, s_, t_] = {D[z[r, s, t], r], D[z[r, s, t], s], 
        D[z[r, s, t], t]}; 

    Vxyz[r_, s_, t_] = 
      TrigExpand[Det[{gradx[r, s, t], grady[r, s, t], gradz[r, s, t]}]]\)

Good, this never goes negative for the r, s, and t used here.
Here comes the measurement of the net flow of Field[x, y, z] across C  
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            divField[x, y, z] dx dy 
              dz\)\)\)
        = 
        (∫_rlow^rhigh)
            (∫_slow^shigh)
              (∫_tlow
                ^thigh) divField[x[r, s, t], y[r, s, t], z[r, s, t]] 
                (Subscript[V, xyz])[r, s, t] dt ds 
                dr

(∫_rlow^rhigh)
        (∫_slow^shigh)
          (∫_tlow^thigh)
             divField[x[r, s, t], y[r, s, t], z[r, s, t]] Vxyz[r, s, t] 
            \[DifferentialD]t \[DifferentialD]s \[DifferentialD]r

Negative.
Because Gauss's 3D formula says
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            divField[x, y, z] dx dy 
              dz\)\)\)
        = 
        (∫_^)(∫_( R)^) 
            Field . outerunitnormal dA
        = net flow of Field[x, y, z] across C, 
this negative result tells you that the net flow of this vector field is against the outernormals. 
The bottom line: 
The net flow of this vector field across this surface is from outside to inside. There must be lots of sinks inside the surface.

T.1.a.ii)

Go with the same surface as in part i) above, but this time measure the flow of the 3D vector field
       Field[x, y, z] = {x +  Sin[y], y +  4  (e^(-x)), -2 z +  (x^3)}
          
across this surface. 
Determine whether the net flow across the surface is from inside to outside or from outside to inside.

Enter the vector field:

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {x + Sin[y], 
        y + 4\ E\^\(-z\), \(-2\)\ z + x\^3}; 

    Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

Calculate divField[x, y, z]:

\(Clear[divField]; 

    divField[x_, y_, z_] = 
      D[m[x, y, z], x] + D[n[x, y, z], y] + D[p[x, y, z], z]\)

No sources or sinks.
The net flow of this vector field across this surface is 0.  In fact the net flow of this vector field across any surface that is the whole skin of a solid region is 0.  The outside-to-inside flow is exactly balanced by the inside-to-outside flow.
You could've done this one by hand.

T.1.a.iii)

Go with the same surface as in part i) above, but this time measure the flow of the 3D vector field
       Field[x, y, z] = {x +  Sin[y], y +  4  (e^(-x)), -5 z +  (x^3)}
          
across this surface. 
Determine whether the net flow is from inside to outside or from outside to inside.

Enter the vector field:

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {x + Sin[y], 
        y + 4\ E\^\(-z\), \(-5\)\ z + x\^3}; 

    Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

Calculate divField[x, y, z]:

\(Clear[divField]; 

    divField[x_, y_, z_] = 
      D[m[x, y, z], x] + D[n[x, y, z], y] + D[p[x, y, z], z]\)

Negative; every point is a sink.
The net flow of this vector field across this surface is from outside to inside.  In fact, the net flow of this vector field across any surface that is the whole skin of a solid region is from outside to inside.
You could've done this one by hand.

T.1.b) The situation dictates a direct calculation by a surface integral

Here's a surface:

\(Clear[x, y, z, u, v]; 
x[u_, v_] = 3\ v; 
\(y[u_, v_] = \(-2\)\ u\ v; 
    \)\)z[u_, v_] = 3 u - v; {ulow, uhigh} = {0, 2}; {vlow, vhigh} = {0, 1.3}; surface = 
      ParametricPlot3D[
        Evaluate[{x[u, v], y[u, v], z[u, v]}], {u, ulow, uhigh}, {v, vlow, 
          vhigh}, ViewPoint -> CMView, 
        AxesLabel -> {"x", " y", "z"}];

Go with 
       f[x, y, z] = ((((x +  4) (y +  3) z))^2), 
and measure the flow of the gradient field of f[x, y, z] across this surface.

You can't use Gauss's formula because this surface is not the skin of a solid region, so you've got to go with a direct calculation via the surface integral.
Enter the function and its gradient field:

\(Clear[f, x, y, z, gradf]; 

    f[x_, y_, z_] = \((\((x + 4)\)\ \((y + 3)\)\ z)\)\^2; 

    gradf[x_, y_, z_] = {D[f[x, y, z], x], D[f[x, y, z], y], D[f[x, y, z], z]}
        \)

Enter the vector field:

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = gradf[x, y, z]; 

    Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

Reenter the parameterization of the surface and calculate normal[u, v]:

\(Clear[x, y, z, u, v]; 
x[u_, v_] = 3\ v; 
\(y[u_, v_] = \(-2\)\ u\ v; 
    \)\)z[u_, v_] = 3 u - v; {ulow, uhigh} = {0, 2}; {vlow, vhigh} = {0, 1.3}; \(Clear[normal]; 

    normal[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]\)

Calculate  
       (∫_ulow^uhigh)
          (∫_vlow
            ^vhigh) Field[x[u, v], y[u, v], z[u, v]] . normal[u, v] 
            dv du:

flowacross = 
      (∫_ulow^uhigh)
          (∫_vlow
            ^vhigh) Field[x[u, v], y[u, v], z[u, v]] . normal[u, v] 
            \[DifferentialD]v \[DifferentialD]u; N[flowacross]

Strongly positive.
To interpret this, look at:

{u, ulow, uhigh}

{v, vlow, vhigh}

scalefactor = 0.3; \(normals = 
      Table[Arrow[normal[u, v], Tail -> {x[u, v], y[u, v], z[u, v]}, 
          ScaleFactor -> scalefactor], {u, 0.5, 1.5, 1}, {v, 0.2, 1.0, 
          0.8}]; 


    Show[surface, normals, ViewPoint -> CMView, Boxed -> False, 
      PlotRange -> All, AxesLabel -> {"x", "y", "z"}]; \)

The normals are pointing out from the side of the surface you are sitting on.
Because the flow measurement is hugely positive, this tells you that the net flow of gradf[x, y, z] across this surface is strongly from the side you can't see to the side you can see.

T.2)  Using Gauss's formula to avoid a calculational nightmare: 
Calculating flow across an oddball surface by calculating the flow across a substitute surface

(Subscript[C, 1]) is the top half of a sphere of radius 2 centered at the origin:

\(Clear[x1, y1, z1, s, t]; 
x1[s_, t_] = 2\ Sin[s]\ Cos[t]; 

    \(y1[s_, t_] = 2\ Sin[s]\ Sin[t]; \)\)z1[s_, t_] = 2 Cos[s]; {{slow, shigh}, {tlow, thigh}} = {{0, (π/2)}, {0, 2 π}}; 
    C1plot = 
      ParametricPlot3D[{x1[s, t], y1[s, t], z1[s, t]}, {s, slow, shigh}, {t, 
          tlow, thigh}, Boxed -> False, ViewPoint -> CMView, 
        AxesLabel -> {"x", "y", "z"}, 
        PlotRange -> All];

(Subscript[C, 2]) is an oddball surface that fits under (Subscript[C, 1]) but agrees with (Subscript[C, 1]) on their common boundary curve.

\(Clear[x2, y2, z2, r, t]; 
x2[r_, t_] = r\ Cos[t]; 

    \(y2[r_, t_] = r\ Sin[t]; \)\)z2[r_, t_] = 0.2 r Sin[(t/2)]^2 ((E^((Sqrt[4 - (r^2)]))) - 1); 
    {{rlow, rhigh}, {tlow, thigh}} = {{0, 2}, {0, 2 π}}; C2plot = 
      ParametricPlot3D[{x2[r, t], y2[r, t], z2[r, t]}, {r, rlow, rhigh}, {t, 
          tlow, thigh}, Boxed -> False, ViewPoint -> CMView, 
        AxesLabel -> {"x", "y", "z"}, 
        PlotRange -> All];

Take a look at the surfaces together from underneath:

Show[C1plot, C2plot, ViewPoint -> {1, 1, -1}];

Together (Subscript[C, 1]) and (Subscript[C, 2]) make the skin of a solid region.

T.2.a.i)

How do you know that when you go with a 3D vector field Field[x, y, z] with divField[x, y, z] = 0 throughout the solid region whose top skin is (Subscript[C, 1]) and whose bottom skin is (Subscript[C, 2]), then
       (∫_^)(∫_( (Subscript[C, 1]))^) 
          Field . topunitnormal dA
        = 
        (∫_^)(∫_( (Subscript[C, 2]))^) 
            Field . topunitnormal dA,
so that the flow of Field[x, y, z] across both surfaces is the same?

Remember that (Subscript[C, 1]) and (Subscript[C, 2]) share a common boundary curve. 
Make a solid region R whose top skin is (Subscript[C, 1]) and whose bottom skin is (Subscript[C, 2]).  Agree that C stands for the skin of R.
You are armed with the fact that divField[x, y, z] = 0 throughout R. Use this fact and Gauss's formula to see that
       \(0 = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
              0 dx dy 
                dz\)\)\) 
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
              divField[x, y, z] dx dy 
                dz\)\)\)\)
        = 
        (∫_^)(∫_( C)^) 
            Field . outerunitnormal dA
        = 
        (∫_^)(∫_( (Subscript[C, 1]))^) 
              Field . topunitnormal dA - 
          (∫_^)(∫_( (Subscript[C, 2]))^) 
              Field . topunitnormal dA,
because along (Subscript[C, 1]), the outer unit normal of C agrees with the top unit normal of (Subscript[C, 1]), but along (Subscript[C, 2]) the outer unit normal of C agrees with the negative of the top unit normal of (Subscript[C, 2]).
The upshot:
       (∫_^)(∫_( (Subscript[C, 1]))^) 
            Field . topunitnormal dA = 
        (∫_^)(∫_( (Subscript[C, 2]))^) 
            Field . topunitnormal dA.
Accordingly,
       flow of Field[x, y, z] across (Subscript[C, 1]) 
        = 
        (∫_^)(∫_( (Subscript[C, 1]))^) 
            Field . topunitnormal dA
        = 
        (∫_^)(∫_( (Subscript[C, 2]))^) 
            Field . topunitnormal dA
        = flow of Field[x, y, z] across (Subscript[C, 2]).

T.2.a.ii)

What calculational nightmare does this help you to avoid?

Because the parameterization of C involves screwy functions like
       z2[r, t]=0.2r ((Sin[(t/2)])^( 2))((E^((Sqrt[4 -  (r^2)]))) -  1), 
calculating 
       (∫_^)(∫_( (Subscript[C, 1]))^) 
          Field . topunitnormal dA
will probably be a nightmare.  Its normals make you want to vomit.  Take a look:

\(Clear[normal2]; 

    normal2[r_, t_] = 
      D[{x2[r, t], y2[r, t], z2[r, t]}, r]
          \[Cross]D[{x2[r, t], y2[r, t], z2[r, t]}, t]\)

These normals will make you choke, and they will probably make Mathematica choke too. 
But you know that
       flow of Field[x, y, z] across (Subscript[C, 1]) 
        = 
        (∫_^)(∫_( (Subscript[C, 1]))^) 
            Field . topunitnormal dA
        = 
        (∫_^)(∫_( (Subscript[C, 2]))^) 
            Field . topunitnormal dA
        = flow of Field[x, y, z] across (Subscript[C, 2]),
so instead of confronting the nightmarish calculation of 
       (∫_^)(∫_( (Subscript[C, 2]))^) 
          Field . topunitnormal dA,
you can opt for the simpler calculation of 
       (∫_^)(∫_( (Subscript[C, 1]))^) 
          Field . topunitnormal dA.

T.2.a.iii)

Illustrate this idea by calculating the flow of the 3D vector field
       Field[x, y, z] = {x -  y, -2 y -  x, z +  3} 
across (Subscript[C, 2]).

Enter the vector field:

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {x - y, \(-2\)\ y - x, 
        z + 3}; 

    Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

Check divField[x, y, z]:

\(Clear[divField]; 

    divField[x_, y_, z_] = 
      D[m[x, y, z], x] + D[n[x, y, z], y] + D[p[x, y, z], z]\)

Good. 
Now you know you can measure the flow of this vector field across the squirrely surface (Subscript[C, 2]) by measuring the flow of this vector field across the clean surface (Subscript[C, 1]).  
Here you go:
       (∫_^)(∫_( (Subscript[C, 1]))^) 
          Field . topunitnormal dA
        = 
        (∫_0
            ^(2 π))(∫_0
              ^((π/2)))Field[x1[s, t], y1[s, t], z1[s, t]] . 
                normal1[s, t] ds dt
where normal1[s, t] is:

\(Clear[normal1]; 

    normal1[s_, t_] = 
      TrigExpand[
        D[{x1[s, t], y1[s, t], z1[s, t]}, s]
            \[Cross]D[{x1[s, t], y1[s, t], z1[s, t]}, t]]\)

You can see that normal1[s, t] is an top normal by looking at the third slot, which is 2  Sin[2  s].  In the parameterization of (Subscript[C, 1]), s runs from 0 to (π/2) and 2  Sin[2  s] remains nonnegative for 0 ≤ s ≤ (π/2):

Plot[2 Sin[2 s], {s, 0, (π/2)}, 
      AspectRatio -> (1/(2 π))];

Now you're ready to calculate:
       flow of Field[x, y, z] across (Subscript[C, 2])
        = flow of Field[x, y, z] across (Subscript[C, 1])
        = 
        (∫_^)(∫_( (Subscript[C, 1]))^) 
            Field . topunitnormal dA
        = 
        (∫_0
            ^(2 π))(∫_0
              ^((π/2)))Field[x1[s, t], y1[s, t], z1[s, t]] . 
                normal1[s, t] ds dt

(∫_0
        ^(2 π))(∫_0
          ^((π/2)))Field[x1[s, t], y1[s, t], z1[s, t]] . normal1[s, t] 
          \[DifferentialD]s \[DifferentialD]t

Strong flow from low to high across (Subscript[C, 2]) (as well as (Subscript[C, 1])).
The direct calculation of the flow of this field across (Subscript[C, 2]) is not practical.
Gauss was a master at avoiding calculational nightmares.

\(Clear[k]; 
∑\+\(k = 1\)\%1000 k\)

T.3) Using Gauss's formula to take advantage of singularities: Calculating flow across the skin of a solid region by calculating the flow across a substitute sphere

Here's an elliptical cylinder with top and bottom:

radial1 = {1, -2, -1}; radial2 = {1, 1, -1}; core = {0.5, 2, 2}; \(Clear[sideplotter, topplotter, bottomplottter, r, u, t]; 

    sideplotter[u_, t_] = Cos[t]\ radial1 + Sin[t]\ radial2 + u\ core; 

    \({{ulow, uhigh}, {tlow, thigh}} = {{\(-1\), 2}, {0, 2\ π}}; \)\)topplotter[r_, t_] = 
      r Cos[t] radial1 + r Sin[t] radial2 + uhigh core; bottomplotter[r_, t_] = 
      r Cos[t] radial1 + r Sin[t] radial2 + ulow core; {rlow, rhigh} = {0, 1}; sides = 
      ParametricPlot3D[
        Evaluate[sideplotter[u, t]], {u, ulow, uhigh}, {t, tlow, thigh}, 
        PlotPoints -> {2, Automatic}, DisplayFunction -> Identity]; 
    top = 
      ParametricPlot3D[
        Evaluate[topplotter[r, t]], {r, rlow, rhigh}, {t, tlow, thigh}, 
        PlotPoints -> {2, Automatic}, DisplayFunction -> Identity]; 
    \(bottom = 
      ParametricPlot3D[
        Evaluate[bottomplotter[r, t]], {r, rlow, rhigh}, {t, tlow, thigh}, 
        PlotPoints -> {2, Automatic}, DisplayFunction -> Identity]; 
    

Show[sides, top, bottom, Boxed -> False, 
      ViewPoint -> CMView, 
      AxesLabel -> {"x", "y", "z"}, PlotRange -> All, 
      DisplayFunction -> $DisplayFunction]; \)

Go with the electric field 
       Field[x, y, z]={((2 x)/((((x^2)  +  (y^2)  +  (z^2)))^(3/2))),((2 y)/((((x^2)  +  (y^2)  +  (z^2)))^(3/2))),((2 z)/((((x^2)  +  (y^2)  +  (z^2)))^(3/2)))}
resulting from a charge of strength 2 placed at the point {0, 0, 0} which is inside this skin.

T.3.a)

Try to use Gauss's 3D formula to calculate the flow of this vector field across the skin of the surface plotted above.

Enter the vector field:

\(Clear[x, y, z, m, n, p, Field]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = {
        \(2\ x\)\/\((x\^2 + y\^2 + z\^2)\)\^\(3/2\), 
        \(2\ y\)\/\((x\^2 + y\^2 + z\^2)\)\^\(3/2\), 
        \(2\ z\)\/\((x\^2 + y\^2 + z\^2)\)\^\(3/2\)}; 

    Field[x_, y_, z_] = {m[x, y, z], n[x, y, z], p[x, y, z]}\)

Calculate divField[x, y, z]:

\(Clear[divField]; 

    divField[x_, y_, z_] = 
      Together[D[m[x, y, z], x] + D[n[x, y, z], y] + D[p[x, y, z], z]]\)

Good; Gauss's 3D formula tells you if you let R stand for the solid of which the plotted surface above is the outside skin, then you are guaranteed that the flow of this vector field across this surface is given by  
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            divField[x, y, z] dx dy 
              dz\)\)\)
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
                0 dx dy dz\)\) = 
          0\)\).
You happily report that the flow of Field[x, y, z] across the skin plotted above is 0 and go on to the next problem.

T.3.b)

Was the answer given in part a) correct?

No way.

T.3.c.i)

What went wrong?

Look at:

Field[0, 0, 0]

Field[x, y, z] goes nuts at {0, 0, 0}.
In short, {0, 0, 0} is a singularity of Field[x, y, z], and the damage comes from the fact that {0, 0, 0} is inside the plotted surface above.
Anytime you have a singularity on or in a solid region R, Gauss's formula
       \(∫\_\%\(∫\_\%\(∫\_\(\ R\)\%\ 
            divField[x, y, z] dx dy 
              dz\)\)\)
        = flow of Field[x, y, z] across C
has the possibility of failing. 
And it did fail in the "answer" given in part a).

T.3.c.ii)

Are there other singularities?

Look at:

Field[x, y, z]

A short examination of this formula for Field[x, y, z] shows that the only point at which the denominators are 0 is {0, 0, 0}.  There are no other singularities.

T.3.d)

One way to try to do the calculation of the flow of Field[x, y, z] across this surface is to let C stand for the whole surface and to try to do a direct calculation of
       (∫_^)(∫_( C)^) 
          Field[x, y, z] . outerunitnormal dA
by breaking it up in the form:
       (∫_^)(∫_( (Subscript[C, sides]))^) 
          Field[x, y, z] . outerunitnormal dA
       +
        (∫_^)(∫_( (Subscript[C, top]))^) 
            Field[x, y, z] . outerunitnormal dA
       +
        (∫_^)(∫_( 
              (Subscript[C, bottom]))^) 
            Field[x, y, z] . outerunitnormal dA.
This is a calculational nightmare because it involves setting up three different normal vectors and then hoping that Mathematica can do the resulting integrals.
Part of the art of mathematics is knowing how to avoid calculational nightmares.
How can you avoid this nightmare?

You center a little sphere called (Subscript[R, little]) at the lone singularity:

singularity = {0, 0, 0}; \(Clear[xlittle, ylittle, zlittle, s, t, littlesphereplotter]; 

    littleradius = 0.1; 

    \({xlittle[s_, t_], ylittle[s_, t_], zlittle[s_, t_]} = 
      singularity + {littleradius\ Sin[s]\ Cos[t], 
          littleradius\ Sin[s]\ Sin[t], littleradius\ Cos[s]}; \)\)littlesphereplot = 
      ParametricPlot3D[
        Evaluate[{xlittle[s, t], ylittle[s, t], zlittle[s, t]}], {s, 0, 
          π}, {t, 0, 2 π}, Boxed -> False, 
        ViewPoint -> CMView, PlotRange -> All];

The singularity at {0, 0, 0} is encapsulated inside the little sphere, and the little sphere is inside the original surface.
Now you make a new hollow solid, taking everything inside the original skin C, but rejecting everything inside the skin, (Subscript[C, little]), of the little sphere.  Let (Subscript[R, new]) stand for the new hollow solid and let (Subscript[C, new]) stand for its skin.  Note that there are no singularities inside (Subscript[R, new]) because the lone singularity lies inside the hollow part. Apply Gauss's formula to the new hollow solid (Subscript[R, new]) to see that
       (∫_^)(∫_( (Subscript[C, new]))^) 
          Field[x, y, z] . outerunitnormal dA
       \(\( = 
        ∫\_\%\(∫\_\%\(∫\_\(\ R\_new\)\%\ 
              divField[x, y, z] dx dy 
                dz\)\)\)\)
        = 0
because divField[x, y, z] = 0 throughout (Subscript[R, new]).
This tells you that 
       (∫_^)(∫_( (Subscript[C, new]))^) 
            Field[x, y, z] . outerunitnormal dA = 0.
But if you let (Subscript[C, little]) stand for the skin of the little sphere centered at the singularity, and you take the outer unit normal from the little sphere, then you get 
       0 = 
        (∫_^)(∫_( (Subscript[C, new]))^) 
            Field[x, y, z] . outerunitnormal dA
       \(\( = 
          ∫\_\%\(∫\_\(\ C\)\%\ 
              Field[x, y, z] . outerunitnormal dA\)\)
		\ \ 
        \ \ \ \(-
          \(∫\_\%\(∫\_\(\ 
                C\_little\)\%\ 
              Field[x, y, z] . outerunitnormal dA\)\)\)\)
because the outer unit normal on the new solid points inside the little sphere.
This is big news because it tells you that
       (∫_^)(∫_( C)^) 
          Field[x, y, z] . outerunitnormal dA
        = 
        (∫_^)(∫_( (Subscript[C, little]))^) 
            Field[x, y, z] . outerunitnormal dA
As a result, you can avoid the calculation of the gruesome integral
       (∫_^)(∫_( C)^) 
          Field[x, y, z] . outerunitnormal dA
by calculating one single, easy integral:
       (∫_^)(∫_( (Subscript[C, little]))^) 
          Field[x, y, z] . outerunitnormal dA
        = 
        (∫_0^(2  π))(∫_0^π)  
            Field[(Subscript[x, little]), (Subscript[y, little]), (Subscript[z, little])] . normal[s, t] 
              ds dt
where normal[s, t] is given by:

\(Clear[normal]; 

    normal[s_, t_] = 
      TrigExpand[
        D[{xlittle[s, t], ylittle[s, t], zlittle[s, t]}, s]
            \[Cross]D[{xlittle[s, t], ylittle[s, t], zlittle[s, t]}, t]]\)

Make sure that these normals are outer normals:

scalefactor = 5; Show[littlesphereplot, 
      Table[Arrow[normal[s, t], 
          Tail -> {xlittle[s, t], ylittle[s, t], zlittle[s, t]}, 
          ScaleFactor -> scalefactor], {s, (π/4), ((3 π)/4), 
          (π/4)}, {t, 0, 2 π, (π/2)}], Boxed -> False, 
      ViewPoint -> CMView, PlotRange -> All];

They are outer normals.
Now you can correctly calculate the flow of Field[x, y, z] across the original cylindrical surface in one short, sweet calculation:

(∫_0
        ^π)(∫_0
          ^(2 π))Field[xlittle[s, t], ylittle[s, t], zlittle[s, t]] . 
            normal[s, t] \[DifferentialD]t \[DifferentialD]s

And you're out of here.

T.4) Combined electric fields and multiple singularities

If you understand B.5) and T.3), then you understand why it is that when you take a charge of strength q and place it at a point {a, b, c} inside the whole skin of a solid region R, then the flow of the resulting electric field across the whole skin is the same as its flow across any small sphere centered at {a, b, c}.
This calculates out to:

\(Clear[x, y, z, a, b, c, Field]; 

    Field[x_, y_, 
        z_] = {q\ \((x - a)\), q\ \((y - b)\), q\ \((z - c)\)}
          \/\((\((x - a)\)\^2 + \((y - b)\)\^2 + \((z - c)\)\^2)\)\^\(3/2\); 
    
Clear[x, y, z, r, s, t, q, normal]; 
center = {a, b, c}; 

    \({x[s_, t_], y[s_, t_], z[s_, t_]} = 
      center + {r\ Sin[s]\ Cos[t], r\ Sin[s]\ Sin[t], r\ Cos[s]}; \)\){{slow, shigh}, {tlow, thigh}} = {{0, π}, {0, 2 π}}; \(normal[s_, t_] = 
      TrigExpand[
        D[{x[s, t], y[s, t], z[s, t]}, s]
            \[Cross]D[{x[s, t], y[s, t], z[s, t]}, t]]; 

    ∫\_slow\%shigh
        \(∫\_tlow
          \%thigh Field[x[s, t], y[s, t], z[s, t]] . normal[s, t] 
          \[DifferentialD]t \[DifferentialD]s\)\)

This is the same as 4  π q.
The upshot:
If you place a charge of strength q inside the skin of a solid region, then flow of the electric field across this skin is 4  π q.
If q > 0, the flow is from inside to outside, and if q < 0, the flow is from outside to inside.

T.4.a)

Take any skin of any solid region you like.
Here is the electric field resulting from a charge of strength (Subscript[q, 1]) placed at a point {a, b, c} inside this skin.

\(Clear[a, b, c, q1, x, y, z, Field1]; 
point1 = {a, b, c}; 

    Field1[x_, y_, z_] = 
      \(q1\ \(({x, y, z} - point1)\)\)\/\((
            \((x - a)\)\^2 + \((y - b)\)\^2 + \((z - c)\)\^2)\)\^\(3/2\)\)

Here is the electric field resulting from a charge of strength (Subscript[q, 2]) placed at a different point {aa, bb, cc} inside the same skin.

\(Clear[aa, bb, cc, Field2]; 
point2 = {aa, bb, cc}; 

    Field2[x_, y_, z_] = 
      \(q2\ \(({x, y, z} - point2)\)\)\/\((
            \((x - aa)\)\^2 + \((y - bb)\)\^2 + \((z - cc)\)\^2)\)\^\(3/2
            \)\)

Measure the flow of the combined electric field 
       Field1[x, y, z] +  Field2[x, y, z] 
across this skin.

Again you try to practice the art of mathematics in an effort to avoid a calculational nightmare.
Check the divergence of Field1[x, y, z], Field2[x, y, z], and 
Field1[x, y, z] +  Field2[x, y, z]:

\(Clear[m1, n1, p1, divF1]; 
    
{m1[x_, y_, z_], n1[x_, y_, z_], p1[x_, y_, z_]} = Field1[x, y, z]; 

    divFie1d1[x_, y_, z_] = 
      Together[D[m1[x, y, z], x] + D[n1[x, y, z], y] + D[p1[x, y, z], z]]\)

\(Clear[m2, n2, p2, divField2]; 
    
{m2[x_, y_, z_], n2[x_, y_, z_], p2[x_, y_, z_]} = Field2[x, y, z]; 

    divField2[x_, y_, z_] = 
      Together[D[m2[x, y, z], x] + D[n2[x, y, z], y] + D[p2[x, y, z], z]]\)

\(Clear[m, n, p, divField1plusField2]; 
    
{m[x_, y_, z_], n[x_, y_, z_], p[x_, y_, z_]} = 
      Field1[x, y, z] + Field2[x, y, z]; 

    divField1plusField2[x_, y_, z_] = 
      Together[D[m[x, y, z], x] + D[n[x, y, z], y] + D[p[x, y, z], z]]\)

Good!
-> divField1[x, y, z] = 0 at all points inside the skin except at its singularity at {a, b, c};
-> divField2[x, y, z] = 0 at all points inside the skin except at its singularity at {aa, bb, cc}; 
-> The divergence of the combined electric field
       Field1[x, y, z] +  Field2[x, y, z] 
is 0 at all points inside the skin except at the two singularities at {a, b, c} and {aa, bb, cc}.
Now on to the calculation:
Center a little sphere (Subscript[C, 1]) at the singularity at {a, b, c} and center another little sphere (Subscript[C, 2]) at the singularity at {aa, bb, cc}.  If you understand what happened in T.3), then you'll understand that the flow of 
       Field1[x, y, z] +  Field2[x, y, z] 
across the skin is the same as the sum
       \(\(∫\_\%\(∫\_\(\ C\_1\)\%\ 
            \((Field1 + \ Field2)\) . outerunitnormal dA\)\)

        	\ \ \ \ \ \ 
        \(+\(∫\_\%\(∫\_\(\ C\_2\)\%\ 
              \((Field1 + \ Field2)\) . outerunitnormal 
                dA\)\)\)\)
       \(\( = 
          ∫\_\%\(∫\_\(\ C\_1\)\%\ 
              Field1 . outerunitnormal dA\)\)
	\ \ \ \ \ \ 
        \(+\(∫\_\%\(∫\_\(\ C\_1\)\%\ 
              \ Field2 . outerunitnormal dA\)\)\)
	\ \ \ \ \ 
        \ \(+
          \(∫\_\%\(∫\_\(\ C\_2\)\%\ 
              \ Field1 . outerunitnormal dA\)\)\)
	\ \ \ \ \ 
        \ \(+\(∫\_\%\(∫\_\(\ 
                C\_2\)\%\ \ 
              Field2 . outerunitnormal dA\)\)\)\)
Now, because divField1[x, y, z] = 0 inside (Subscript[C, 2]), and because Field1[x, y, z] has no singularity inside (Subscript[C, 2]), you know that
       (∫_^)(∫_( (Subscript[C, 2]))^)  
            Field1 . outerunitnormal dA = 0.
Similarly,
       (∫_^)(∫_( (Subscript[C, 1]))^)  
            Field2 . outerunitnormal dA = 0.
So the flow of the combined electric field
Field1[x, y, z] +  Field2[x, y, z] across the original skin reduces to
       \(\(∫\_\%\(∫\_\(\ C\_1\)\%\ 
            Field1 . outerunitnormal dA\)\)
	\ \ \ 
        \(+\(∫\_\%\(∫\_\(\ C\_2\)\%\ 
              Field2 . outerunitnormal dA\)\)\)\)
The calculation done at the very beginning tells you that this is just
       4  (Subscript[q, 1]) π +  4  (Subscript[q, 2]) π = 
        4  π (q1 +  q2).
By using your head instead of the machine, you side-stepped another gruesome calculation.
By the way, electrical folks and physicists really get off on this stuff.  
You should like it too because it is clean, useful mathematics.

T.5) Surface packaging: Parameterized, explicit and implicit

If you have a surface to work with, the best way it can come into your hands is in parameterized form 
       {x[u, v], y[u, v], z[u, v]}, 
where {u, v} varies in a region (Subscript[R, uv]).  In this form:
-> You can plot it with ParametricPlot3D. 
-> You can measure its surface area with the integral
       (∫_^)(∫_( (Subscript[R, uv]))^) 
          (Subscript[SA, xyz])[u, v] du dv
where (Subscript[SA, xyz])[u, v]  is given by:

\(Clear[x, y, z, u, v, cross, SAxyz]; 

    cross[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]; 
SAxyz[u_, v_] = \@\(cross[u, v] . cross[u, v]\)\)

-> You can integrate a function f[x, y, z] with respect to surface area on this surface by means of the integral
       (∫_^)(∫_( surface)^) 
          f[x, y, z] dA
        = 
        (∫_^)(∫_( (Subscript[R, uv]))^) 
            f[x[u, v], y[u, v], z[u, v]] (Subscript[SA, xyz])[u, v] du 
              dv.
-> You can measure the flow of a given vector field Field[x, y, z] across the surface by means of the integral
       (∫_^)(∫_( surface)^) 
          Field . unitnormal dA
        = 
        (∫_^)(∫_( (Subscript[R, uv]))^) 
            Field[x[u, v], y[u, v], z[u, v]] . normal du 
              dv.
where normal[u, v] is given by:

\(Clear[normal]; 

    normal[u_, v_] = 
      D[{x[u, v], y[u, v], z[u, v]}, u]\[Cross]D[{x[u, v], y[u, v], z[u, v]}, 
          v]\)

If you have the surface already parameterized, you can go right to work using these formulas to make measurements like a bandit.
But often, the surface is not given to you in parametric form. When this happens, you can be in the lucky or the (often) unlucky case.

T.5.a) The lucky case: An explicit formula

You are lucky if you get a surface in one of the forms 
       z = f[x, y], y = g[x, z], or x = h[y, z].

Why are you lucky when you get an explicit equation for a surface?

You are lucky because the explicit form carries with it a built-in parameterization.  
If the surface comes to you with the explicit formula 
       z = f[x, y], 
then you can parameterize it with the formulas
       x[u, v] = u,
       y[u, v] = v, and
       z[u, v] = f[u, v].
If the surface comes to you with the explicit formula 
       y = g[x, z], 
then you can parameterize it with the formulas
       x[u, v] = u,
       y[u, v] = g[u, v], and
       z[u, v] = v.
If the surface comes to you with the explicit formula 
       x = h[y, z], 
then you can parameterize it with the formulas
       x[u, v] = h[u, v],
       y[u, v] = u, and
       z[u, v] = v.

T.5.b) The (often) unlucky case: An implicit equation

You can be in deep doo-doo if you get a surface in the form  
       f[x, y, z] = constant.

Why can you be unlucky when you get an implicit equation for a surface?
What is the first thing you should try to do in order to try to work with surface specified with an implicit equation?

You can't do diddly-squat with an implicit equation until you come up with a parameterization or an explicit equation. 
So the first thing you should try to do is to try get a parameterization.  This is not always a pleasant task.  Sometimes it's impossible.