VC.06 Sources, Sinks, Swirls and Singularities
 BASICS

\(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};

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]; \)

B.1)  Using a 2D integral to measure flow across closed curves

B.1.a)

Explain this:
To calculate the net flow of a vector field 
       Field[x, y] = {m[x, y], n[x, y]} 
across the boundary C of a region R, you have your choice:
-> You can go through the labor of parameterizing C, and then calculate 
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy,
-> or, if the field has no singularities inside R, you can put 
       divField[x, y] = D[m[x, y], x] +  D[n[x, y], y] 
and calculate the 2D integral
       (∫_^)(∫_( R)) 
          divField[x, y] dx dy.

Back in the lesson on 2D integrals (Lesson 2.05), you met up with the Gauss-Green formula. The Gauss-Green formula says that if C is the boundary curve of a region R, then you are guaranteed that 
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
        = 
        (∫_^)(∫_( R)) 
            divField[x, y] dx dy.
Go with Field[x, y] = {m[x, y], n [x, y}, and put 
       divField[x, y] = D[m[x, y], x] +  D[n[x, y], y] 
and read off  
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
        = 
        (∫_^)(∫_( R)) 
            divField[x, y] dx dy.
Because (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy measures the net flow of Field[x, y] = {m[x, y], n[x, y]} across C, you are guaranteed that   
       (∫_^)(∫_( R)) 
          divField[x, y] dx dy
makes the same measurement.

B.1.b)

Here's the rectang1e R with corners at {-3, -1}, {3, -1}, {3, 1}, and {-3, 1}:

Rplot = 
      Show[Graphics[{Red, Thickness[0.01], 
            Line[{{-3, -1}, {3, -1}, {3, 1}, {-3, 1}, {
                  -3, -1}}]}], Axes -> True, 
        AxesLabel -> {"x", "y"}, AspectRatio -> (1/GoldenRatio)];

Use a 2D integral to measure the net flow of the vector field 
       Field[x, y] = {(x^3) +  y, x -  y} 
across the boundary curve C of this rectangle.

Enter the vector field:

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

Calculate the divergence, divField[x, y]:

\(Clear[divField]; 
divField[x_, y_] = D[m[x, y], x] + D[n[x, y], y]\)

Take another look at R:

Show[Rplot];

Calculate (∫_^)(∫_( R)) 
          divField[x, y] dx dy:

(∫_(-1)^1)
        (∫_(-3)^3) divField[x, y] \[DifferentialD]x 
          \[DifferentialD]y

Big time positive.  This means that the net flow of this vector field across the boundary of this rectangle is from inside to outside.
There must be a lot of sources of new fluid inside the rectangle.

B.1.c.i)

Take a look at divField[x, y] for the vector field
       Field[x, y] = {Sin[y] -  x, Cos[x] -  y}:

\(Clear[x, y, m, n, Field, divField]; 
    
{m[x_, y_], n[x_, y_]} = {Sin[y] - x, Cos[x] - y}; 

    \(Field[x_, y_] = {m[x, y], n[x, y]}; \)\)divField[x_, y_] = D[m[x, y], x] + D[n[x, y], y]

You look at this and note that 
       divField[x, y] < 0 
no matter what {x, y} is, and then you say:
"Good, this tells me that the flow of this vector field across any closed curve is from outside to inside." 
You are right.
Why are you right?

Take any closed curve C, and call R the region C encloses. You can calculate the net flow  
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
of Field[x, y] = {m[x, y], n[x, y]} across C by calculating the 2D integral
       (∫_^)(∫_( R)) 
          divField[x, y] dx dy.
Because 
       divField[x, y] < 0 
no matter what {x, y} you go with, you are guaranteed that
       (∫_^)(∫_( R)) 
            divField[x, y] dx dy < 0.
So:
       (∮_( C)) -n[x, y] dx  +  
          m[x, y] dy = 
        (∫_^)(∫_( R)) 
              divField[x, y] dx dy < 0
no matter what closed curve C you go with.
This tells you that the flow of this vector field across any closed curve is from outside to inside.
Good deal.

B.1.c.ii)

Take a look at divField[x, y] for the vector field
     Field[x, y] = {Sin[y] +  (x^5), Cos[x] +  (y^3)}

\(Clear[x, y, m, n, Field, divField]; 
    
{m[x_, y_], n[x_, y_]} = {Sin[y] + x\^5, Cos[x] + y\^3}; 

    \(Field[x_, y_] = {m[x, y], n[x, y]}; \)\)divField[x_, y_] = D[m[x, y], x] + D[n[x, y], y]

You look at this and note that 
       divField[x, y] > 0 
unless {x, y} = {0, 0}, and then you say:
"Good, this tells me that the flow of this vector field across any closed curve is from inside to outside." 
You are right.
Why are you right?

Take any closed curve C, and call R the region C encloses.  You can calculate the flow  
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
of Field[x, y] = {m[x, y], n[x, y]} across C by calculating the 2D integral
       (∫_^)(∫_( R)) 
          divField[x, y] dx dy.
Because
       divField[x, y] > 0 
except at one point, you are guaranteed that
       (∫_^)(∫_( R)) 
            divField[x, y] dx dy > 0.
So:
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
        = 
        (∫_^)(∫_( R)) 
              divField[x, y] dx dy > 0
no matter what closed curve C you go with.
This tells you that the flow of this vector field across any closed curve is from inside to outside.
DivField[x, y] is a really handy tool.

B.2)  Sources, sinks, and the divergence of a vector field

B.2.a)  The meaning of the sign of
           divField[x, y] = D[m[x, y], x] +  D[n[x, y], y]

Given a vector field 
       Field[x, y] = {m[x, y], n[x, y]},
you calculate 
       divField[x, y] = D[m[x, y], x] +  D[n[x, y], y].
How does the sign of   
       divField[x, y] = D[m[x, y], x] + D[n[x, y], y] 
tell you whether {x, y} is a source of new fluid or a sink (drain) for old fluid?

If divField[(Subscript[x, 0]), (Subscript[y, 0])] > 0, then the point {(Subscript[x, 0]), (Subscript[y, 0])} is a source of new fluid.
If divField[(Subscript[x, 0]), (Subscript[y, 0])] < 0, then the point {(Subscript[x, 0]), (Subscript[y, 0])} is a sink for old fluid.
Here's why: 
Take a small circle C centered at {(Subscript[x, 0]), (Subscript[y, 0])}.  Calculate the flow of Field[x, y] across C by calculating 
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
        = 
        (∫_^)(∫_( R)) 
            divField[x, y] dx dy.
Here's the kicker:  
If 
       divField[(Subscript[x, 0]), (Subscript[y, 0])] > 0, 
then it is positive for all {x, y}'s close to {(Subscript[x, 0]), (Subscript[y, 0])}. So if C is so small that
       divField[(Subscript[x, 0]), (Subscript[y, 0])] > 0
at all {x, y}'s inside C, then you see that
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
        = 
        (∫_^)(∫_( R)) 
              divField[x, y] dx dy > 0.
This means that if divField[(Subscript[x, 0]), (Subscript[y, 0])] > 0, then the net flow of Field[x, y] across small circles centered at {(Subscript[x, 0]), (Subscript[y, 0])} is from inside to outside.  
The upshot: 
If divField[(Subscript[x, 0]), (Subscript[y, 0])] > 0, then the point {(Subscript[x, 0]), (Subscript[y, 0])} is a source of new fluid.
Similarly, if divField[(Subscript[x, 0]), (Subscript[y, 0])] < 0, then the net flow of Field[x, y] across small circles centered at {(Subscript[x, 0]), (Subscript[y, 0])} is from outside to inside.  
So:
If divField[(Subscript[x, 0]), (Subscript[y, 0])] < 0, then the point {(Subscript[x, 0]), (Subscript[y, 0])} is a sink for old fluid.
Check this out:
Go with Field[x, y] = {3  x y, y} and look at divField[x, y]:

\(Clear[x, y, Field, m, n, divField]; 
    
{m[x_, y_], n[x_, y_]} = {3\ x\ y, y}; 

    \(Field[x_, y_] = {m[x, y], n[x, y]}; \)\)divField[x_, y_] = D[m[x, y], x] + D[n[x, y], y]

See whether the point {0, 2} is a source or a sink:

divField[0, 2]

Positive.  This tells you that the point {0, 2} is a source.
Take a look at this vector field on a small circle centered at {0, 2} to see whether this calculation agrees with reality:

point = {0, 2}; radius = 0.1; pointplot = Graphics[{PointSize[0.04], Point[point]}]; \(Clear[x, y, t]; 
{x[t_], y[t_]} = point + radius\ {Cos[t], Sin[t]}; 

    \(tlow = 0; \)\)thigh = 2 π; smallcircle = 
      ParametricPlot[{x[t], y[t]}, {t, tlow, thigh}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        DisplayFunction -> Identity]; scalefactor = 0.05; \(fieldplot = 
      Table[Arrow[Field[x[t], y[t]], Tail -> {x[t], y[t]}, \ 
          ScaleFactor -> scalefactor], {t, tlow, thigh, 
          \(thigh - tlow\)\/16}]; 

    \(Show[pointplot, smallcircle, fieldplot, Axes -> True, 
      AxesLabel -> {"x", "y"}, AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)\)

Confirm by looking at the normal components of the field vectors on the curve:

Clear[normal];
normal[t_]={y'[t],-x'[t]};
Clear[normalcomponent];
normalcomponent[t_] = 
        Field[x[t], y[t]] . normal[t]/normal[t] . normal[t] normal[t];
jump = (π/5);
actualflowacross = 
        Table[Arrow[normalcomponent[t], Tail -> {x[t], y[t]},  
            ScaleFactor -> scalefactor], {t, 0, 2 π - jump, jump}]
          ;

Show[pointplot, smallcircle, actualflowacross, Axes -> True, 
        AxesLabel -> {"x", "y"}, AspectRatio -> Automatic, 
        DisplayFunction -> $DisplayFunction];

Yessiree, Bob.  The plot shows a lot more flow from inside to outside than from outside to inside.  
Just as you would expect on a small circle centered at a source.

B.2.b) Sources and sinks

Here's a vector field:

\(Clear[Field, m, n, x, y]; 
    
{m[x_, y_], n[x_, y_]} = {Sin[x]\ Cos[y], Sin[y]\ Cos[x]}; 

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

Give a sample plot of some of the sources and sinks in this vector field.

Here's divField[x, y]:

\(Clear[divField]; 
divField[x_, y_] = D[m[x, y], x] + D[n[x, y], y]\)

A point {x, y} is a source if divField[x, y] > 0, and {x, y} is a sink if divField[x, y] < 0.
Here comes the plot:

sourcesandsinks = 
      Show[Table[
          If[N[divField[x, y]] > 0, 
            Graphics[{PointSize[0.015], Red, Point[{x, y}]}], 
            Graphics[{PointSize[0.025], GrayLevel[0.5], Point[{x, y}]}]], {x, 
            -5, 5, 0.25}, {y, -4, 4, 0.25}], Axes -> True, 
        AxesLabel -> {"x", "y"}];

Alternate squares of sources and sinks.
Think of the sources as little individual springs feeding the flow.  
Think of the sinks as tiny little holes through which fluid seeps out as the flow goes by.

B.2.c.i) All sources inside C

If every point inside a closed curve (like a deformed circle) C is a source of a given vector field, and if the vector field has no singularities inside C, then how do you know that the net flow of the given vector field across C is automatically from inside to outside?

Answer:

If every point inside a closed curve C is a source of a given vector field, then 
-> new fluid is oozing out of each point inside C, and 
-> there there is no place within C to absorb excess outside-to-inside flow.
The result:
If every point inside a closed curve C is a source of a given vector field, then the flow of this vector field across C is automatically from inside to outside.
For example, look at:

\(Clear[Field, m, n, x, y, divField]; 
    
{m[x_, y_], n[x_, y_]} = {x\^3 - y, y\^3 + x}; 

    \(Field[x_, y_] = {m[x, y], n[x, y]}; \)\)divField[x_, y_] = D[m[x, y], x] + D[n[x, y], y]

Unless {x, y} = {0, 0}, divField[x, y] > 0. 
This tells you that all points {x, y} except {0, 0} are sources for Field[x, y] and the lone exception is not a sink.
This vector field has no singularities. 
So:
On the basis of this information you can say with confidence and authority that the flow of this vector field across any closed curve is from inside to outside.

B.2.c.ii) All sinks inside C

If every point inside a closed curve (like a deformed circle) C is a sink of a given vector field, and if the vector field has no singularities inside C, then how do you know that the net flow of the given vector field across C is automatically from outside to inside?

Answer:

If every point inside a closed curve C is a sink of a given vector field, then 
-> old fluid is soaking into each point inside C, and 
-> there there is no place within C to generate excess inside-to-outside flow.
The result:
If every point inside a closed curve C is a sink of a given vector field, then the flow of this vector field across C is automatically from outside to inside.
For example, look at:

\(Clear[Field, m, n, x, y, divField]; 
    
{m[x_, y_], n[x_, y_]} = {y - x\^3, x - y\^7}; 

    \(Field[x_, y_] = {m[x, y], n[x, y]}; \)\)divField[x_, y_] = D[m[x, y], x] + D[n[x, y], y]

Unless {x, y} = {0, 0}, divField[x, y] < 0. 
This tells you that all points {x, y} except {0, 0} are sinks for Field[x, y] and the lone exception is not a source.
This vector field has no singularities.
So:
On the basis of this information you can say with confidence and authority that the flow of this vector field across any closed curve is from outside to inside.
The region inside C is like a big vacuum sucking up fluid.

B.2.c.iii) No sources or sinks inside C

If there are no sinks and there are no sources of a given vector field inside a closed curve (like a deformed circle) C, and if the vector field has no singularities inside C, then how do you know that the net flow of the given vector field across C is 0?

Answer:

If there are no sinks and there are no sources of a given vector field inside C, and there are no singularities inside C, then no new fluid is injected and no old fluid is sucked up inside C.
So:
-> What flows from inside to outside must be replaced by equal outside to inside flow.
-> What flows from outside to inside must be replaced by equal inside to outside flow.
The result:
If there are no sinks and there are no sources of a given vector field inside C, and there are no singularities inside C, then the net flow of the vector field across C is 0.
For example, look at:

\(Clear[Field, m, n, x, y, divField]; 
    
{m[x_, y_], n[x_, y_]} = {Cos[x]\ Cosh[y], Sin[x]\ Sinh[y]}; 

    \(Field[x_, y_] = {m[x, y], n[x, y]}; \)\)divField[x_, y_] = D[m[x, y], x] + D[n[x, y], y]

This tells you that this vector field has no sources or sinks.
This vector field has no singularities, so:
On the basis if this information you can say with confidence and authority that the net flow of this vector field across any closed curve is 0.

B.2.d) Divergence

Why do most folks call divField[x, y] the divergence of a vector field Field[x, y]?

The name fits.  
If divField[x, y] > 0, then new fluid is oozing out of the point 
{x, y} and diverges elsewhere.
If divField[x, y] < 0, then old fluid is sucked into the point 
{x, y} and converges onto this point.
If divField[x, y] = 0, then no new fluid is introduced, and no old fluid is sucked off as the flow passes by {x, y}.

B.3)  Flow-across-the-curve measurements in the presence of singularities

B.3.a) Singularities

Here are two closed curves, each parameterized correctly in the counterclockwise way:

\(Clear[x1, y1, x2, y2, t]; 
tlow = 0; 
\(thigh = 2\ π; \)\){x1[t_], y1[t_]} = {1, 0} + 4 {Cos[t], 0.6 Sin[t]}; \({x2[t_], y2[t_]} = {1, 0} + 2\ {Cos[t], Sin[t]}; 


    curves = ParametricPlot[{{x1[t], y1[t]}, {x2[t], y2[t]}}, {t, tlow, 
          thigh}, PlotStyle 
          -> {{Red, Thickness[0.01]}, {Red, Thickness[0.01]}}, 
        AxesLabel -> {"x", "y"}, 
        Epilog -> {Text["C1", {\(-0.6\), 0.7}], 
            Text["C2", {3.5, 1.6}]}]; \)

Go with 
       Field[x, y] = {m[x, y], n[x, y]}  
       ={((6 x)/((x^2) +  (y^2))),((6 y)/((x^2) +  (y^2)))},  
and look at the measurements 
              (∮_( (Subscript[C, 1]))) -n[x, y] dx  +  
        m[x, y] dy,
and 
              (∮_( (Subscript[C, 2]))) -n[x, y] dx  +  
        m[x, y] dy
of the flow of this vector field across each curve:

Clear[Field, m, n, x, y];
{m[x_, y_], n[x_, y_]} = {((6 x)/((x^2) + (y^2))), 
          ((6 y)/((x^2) + (y^2)))};
Field[x_, y_] = {m[x, y], n[x, y]};

flowacrossC1=NIntegrate[-n[x1[t], y1[t]] x1'[t]+m[x1[t], y1[t]] y1'[t],{t, tlow, thigh}]

flowacrossC2=NIntegrate[-n[x2[t], y2[t]] x2'[t]+m[x2[t], y2[t]] y2'[t],{t, tlow, thigh}]

The curves are different, but the flow of this vector field across the one curve is the same as the flow of this vector field across the other.
Was this an accident?

You got it right.

Take a look at the vector field.

Field[x, y]

Note the nasty singularity at {0, 0}.

Field[0, 0]

Look at divField[x, y]:

\(Clear[divField]; 

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

Ah-ha! 
The vector field has no sources or sinks other than at the singularity at {0, 0}. 
Check out the position of the singularity relative to the curves:

singularity = {0, 0}; singularityplot = Graphics[{PointSize[0.03], Point[singularity]}]; 
    Show[curves, singularityplot];

The singularity is not between the curves.
And because there are no sources or sinks anywhere but at the singularity, you know that the same amount of fluid that flows across the inner curve also flows across the outer curve.
This is why when you measure the flow of the vector field
       Field[x, y] = {m[x, y], n[x, y]}  
       ={((6 x)/((x^2) +  (y^2))),((6 y)/((x^2) +  (y^2)))},
as above, you get   
       (∮_( (Subscript[C, 1]))) -n[x, y] dx  +  
        m[x, y] dy
        = 
        (∮_( (Subscript[C, 2]))) -n[x, y] dx  +  
          m[x, y] dy

B.3.b) Avoiding grisly parameterizations by encapsulating the singularities.

Look at the following curve C:

Cplot = 
      Graphics[{Thickness[0.01], Red, 
          Line[{{7, 0}, {4, 4}, {5, 5}, {-1, 6}, {-2, 3}, {-5, 
                2}, {-3, -4}, {1, -2}, {4, -3}, {7, 0}}]}]; 
    label = Graphics[Text["C", {2.5, 5}]]; Show[Cplot, label, Axes -> True, AxesLabel -> {"x", "y"}, 
      AspectRatio -> Automatic];

Calculate the flow of 
       Field[x, y]={((2 (x  -  2))/((((x  -  2))^2)  +  
                    (((y  -  1))^2))),((2 (y  -  1))/((((x  -  2))^2)  +  
                    (((y  -  1))^2)))}
across C without going to all the bother of parameterizing C.

This is definitely a case in which a little knowledge can save a lot of work. Only a bean-counting dweeb would find pleasure in parameterizing that silly polygonal curve C.

\(Clear[Field, m, n, x, y]; 

    m[x_, y_] = \(2\ \((x - 2)\)\)\/\(\((x - 2)\)\^2 + \((y - 1)\)\^2\); 

    \(n[x_, y_] = \(2\ \((y - 1)\)\)\/\(\((x - 2)\)\^2 + \((y - 1)\)\^2\); 
    \)\)Field[x_, y_] = {m[x, y], n[x, y]}

Note the singularity at {2, 1}.

Field[2, 1]

Look at divField[x, y]:

Together[D[m[x, y], x] + D[n[x, y], y]]

Good.
If there were no singularities inside C, then this information would be enough to tell you that the flow of this field across C is 0. But the nasty singularity (blow-up) at {2, 1} inside the curve C might be a sink or source. 
This won't hold you back because you can pull off a pretty neat stunt:
Encapsulate the singularity by centering a little circle (Subscript[C, 1]) around the singularity at {2, 1}, taking care that the little circle lies completely within C.  
Here's one:

singularity = {2, 1}; \(Clear[x1, y1, t]; 
    
{x1[t_], y1[t_]} = singularity + 0.5\ {Cos[t], Sin[t]}; 

    \(C1plot = 
      ParametricPlot[{x1[t], y1[t]}, {t, 0, 2\ π}, 
        PlotStyle -> {{Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; \)\)singularityplot = 
      Graphics[{Blue, PointSize[0.03], Point[singularity]}]; \(extralabel = Graphics[Text["C1", {1.3, 0.6}]]; 


    Show[Cplot, C1plot, singularityplot, label, extralabel, 
      PlotRange -> All, Axes -> True, 
      AxesLabel -> {"x", "y"}, AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)

Because
       divField[x, y] = D[m[x, y], x] +  D[n[x, y], y] = 0 
at all points between the circle (Subscript[C, 1]) and the ugly curve C, and because there are no singularities between the two curves, you can be certain that that the flow of this vector field across the little circle (Subscript[C, 1]) is the same as the flow of this vector field across C.  In other words,
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
        = 
        (∮_( (Subscript[C, 1]))) -n[x, y] dx  +  
          m[x, y] dy.  
The beauty of this is that you can calculate 
       (∮_( (Subscript[C, 1]))) -n[x, y] dx  +  
        m[x, y] dy  
with a couple flicks of your fingers:

flowacrossC1=NIntegrate[-n[x1[t], y1[t]] x1'[t]+m[x1[t], y1[t]] y1'[t],{t, 0, 2 π}]

The net flow across the little circle (Subscript[C, 1])  is from inside to outside. 
There was a gushing source at the singularity.
The flow of the given vector field across C is the same as the flow of of this field across (Subscript[C, 1]):

flowacrossC = flowacrossC1

And you did this without going to the trouble of parameterizing C.
A little knowledge saved a lot of work.
Does anyone other than a wuss really want to parameterize that silly curve?

B.3.c) Encapsulating two singularities

Stay with the same curve C as in part b):

Cplot = 
      Graphics[{Thickness[0.01], Red, 
          Line[{{7, 0}, {4, 4}, {5, 5}, {-1, 6}, {-2, 3}, {-5, 
                2}, {-3, -4}, {1, -2}, {4, -3}, {7, 0}}]}]; 
    \(label = Graphics[Text["C", {2.5, 5}]]; 


    Show[Cplot, label, Axes -> True, AxesLabel -> {"x", "y"}, 
      AspectRatio -> Automatic]; \)

Go with:

\(Clear[Field, m, n, x, y]; 

    m[x_, y_] = 
      \(x - 2\)\/\(\((x - 2)\)\^2 + \((y - 1)\)\^2\) - 
        \(2\ \((x + 3)\)\)\/\(\((x + 3)\)\^2 + \((y + 2)\)\^2\); 

    \(n[x_, y_] = 
      \(y - 1\)\/\(\((x - 2)\)\^2 + \((y - 1)\)\^2\) - 
        \(2\ \((y + 2)\)\)\/\(\((x + 3)\)\^2 + \((y + 2)\)\^2\); \)\)Field[x_, y_] = {m[x, y], n[x, y]}

Calculate the flow of this vector field across C without breaking into a heavy sweat.

Look at the vector field again.

Field[x, y]

Note the singularities at at {2, 1} and at {-3, -2}.

{Field[2, 1], Field[-3, -2]}

Now look at:

Together[D[m[x, y], x] + D[n[x, y], y]]

Good. 
The vector field has no sources or sinks aside from the singularities.
Encapsulate the singularities by centering a little circle (Subscript[C, 1]) at the singularity at {2, 1}, and centering another little circle (Subscript[C, 2]) at the singularity at {-3, -2} taking care that both little circles lie completely within C, and taking care that they do not invade each other's territory.  
Here are two very acceptable circles:

singularity1 = {2, 1}; singularity2 = {-3, -2}; \(Clear[x1, y1, x2, y2, t]; 
    
{x1[t_], y1[t_]} = singularity1 + 0.8\ {Cos[t], Sin[t]}; 

    \({x2[t_], y2[t_]} = singularity2 + 0.4\ {Cos[t], Sin[t]}; \)\)littlecircles = 
      ParametricPlot[{{x1[t], y1[t]}, {x2[t], y2[t]}}, {t, 0, 2 π}, 
        PlotStyle -> {{Thickness[0.01], Red}, {Thickness[0.01], Red}}, 
        DisplayFunction -> Identity]; singularityplot = {
        Graphics[{Blue, PointSize[0.02], Point[singularity1]}], 
        Graphics[{Blue, PointSize[0.02], Point[singularity2]}]}; \(extralabels = {Graphics[Text["C1", {1, 0.5}]], 
        Graphics[Text["C2", {\(-3\), \(-1.3\)}]]}; 


    Show[Cplot, littlecircles, singularityplot, label, extralabels, 
      PlotRange -> All, Axes -> True, 
      AxesLabel -> {"x", "y"}, AspectRatio -> Automatic, 
      DisplayFunction -> $DisplayFunction]; \)

Because 
       divField[x, y] = D[m[x, y], x] +  D[n[x, y], y] = 0 
at all points between the circles and the ugly curve C, and since there are no singularities between the circles and C, you can be certain that the net flow across C is given by
       (∮_( C)) -n[x, y] dx  +  
        m[x, y] dy
        = 
        (∮_( (Subscript[C, 1]))) -n[x, y] dx  +  
          m[x, y] dy +  
          (∮_( (Subscript[C, 2]))) -n[x, y] dx  + 
           m[x, y] dy.  
Here you go with just a couple flicks of your fingers:

flowacrossC1=NIntegrate[-n[x1[t], y1[t]] x1'[t]+m[x1[t], y1[t]] y1'[t],{t, 0, 2 π}];flowacrossC2=NIntegrate[-n[x2[t], y2[t]] x2'[t]+m[x2[t], y2[t]] y2'[t],{t, 0, 2 π}];

flowacrossC = flowacrossC1 + flowacrossC2

The net flow of this vector field across C is from outside to inside. At least one of the singularities must have been a big time sink.