2.05 2D Integrals and the Gauss-Green Formula
TUTORIALS

Off[General::spell]; Off[General::spell1]; Off[Plot::plnr]; Off[ParametricPlot3D::ppcom]; SetOptions[Limit, Analytic -> True]; \(Needs["Graphics`FilledPlot`"]; 


    If[MemberQ[{"Global`", \ "Miscellaneous`SIUnits`", \ 
          "Graphics`Colors`"}, 
\ \ \ \ Context[Gray]], \ 
      Remove["*`Gray"]]; 
<< "Miscellaneous`Units`"; 

    Remove["*`Gray"]; 
<< "Graphics`Colors`"; 

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

    
$Pre\  = \ Chop; 
CMView\  = \ {2.7, \ 1.6, \ 1.2}; \)

T.1)  Using a 2D integral to measure area

T.1.a) Area measurements by double integrals

Here's a little region R in the xy-plane:

\(Clear[high, low, x]; 

    high[x_] = 2\ E\^\(\(-0.2\)\ x\)\ Sin[\(π\ x\)\/4]; 

    low[x_] = 0.05\ x\^3\ \((x - 4)\); 


    regionplot = 
      Plot[{high[x], low[x]}, {x, 0, 4}, 
        PlotStyle -> {{Blue, Thickness[0.01]}}, 
        AspectRatio -> Automatic, AxesLabel -> {"x", "y"}, 
        Epilog -> {Text["y = high[x]", {3, 1.1}], 
            Text["y = low[x]", {3, \(-1.2\)}], Text["R", {2, 1}]}]; 
    \)

Use a single integral to measure the area of the region R between the two curves, and then use a 2 D integral to make the same measurement.

Without 2 D integrals, you just calculate
       (∫_0^4) high[x] - low[x] dx:

NIntegrate[high[x] - low[x], {x, 0, 4}]

With 2 D integrals you set f[x, y] = 1 for all x's and y's and calculate
       ∫(∫_R) f[x, y] dx 
              dy = 
        (∫_0^4)
            (∫_(low[x])^(high[x]))f[x, y] dy 
              dx:

\(Clear[y]; 
NIntegrate[1, {x, 0, 4}, {y, low[x], high[x]}]\)

The reason this works is that when you go with
       f[x, y] = 1 
for all x's and y's, then ∫(∫_R) f[x, y] dx 
            dy measures the volume of the flat solid one unit high whose vertical sides have the same shape as R:

\(Clear[r]; 

    top = ParametricPlot3D[{x, low[x] + r\ \((high[x] - low[x])\), 1}, {x, 0, 
          4}, {r, 0, 1}, PlotPoints -> {Automatic, 2}, 
        DisplayFunction -> Identity]; 


    bottom = ParametricPlot3D[{x, low[x] + r\ \((high[x] - low[x])\), 0}, {x, 
          0, 4}, {r, 0, 1}, PlotPoints -> {Automatic, 2}, 
        DisplayFunction -> Identity]; 


    side1 = ParametricPlot3D[{x, low[x], r}, {x, 0, 4}, {r, 0, 1}, 
        PlotPoints -> {Automatic, 2}, DisplayFunction -> Identity]; 
    

side2 = 
      ParametricPlot3D[{x, high[x], r}, {x, 0, 4}, {r, 0, 1}, 
        PlotPoints -> {Automatic, 2}, DisplayFunction -> Identity]; 
    

CMView = {2.7, 1.6, 1.2}; 

    \(Show[top, bottom, side1, side2, ViewPoint -> CMView, 
      DisplayFunction -> $DisplayFunction]; \)\)

Because this solid is one unit high, the volume measurement of the solid also gives you the area measurement of its base. 
The upshot:
When you go with f[x, y] = 1 for all x's and y's, then
       ∫(∫_R) f[x, y] dx 
            dy
measures the area of R.

T.2)  Volume measurements with 2D integrals

T.2.a.i)

Here's a plot of the surface
       z = f[x, y] = 12 - x + 2  y 
above the rectangle R in the xy-plane consisting of all points {x, y, 0} with
      -1 ≤ x ≤ 3 and -1 ≤ y ≤ 2:

\(\(Clear[f, x, y]; 
f[x_, y_] = 12 - x + 2\ y; 
    
{{a, b}, {c, d}}\  = \ {{\(-1\), 3}, {\(-1\), 2}}; 


    surfaceplot\  = \ 
      Plot3D[f[x, y], {x, a, b}, {y, c, d}, 
DisplayFunction -> Identity]; 

    
Rplot\  = \ 
      Graphics3D[
\ \ \ \ \ \ \ \ 
        Polygon[{{a, c, 0}, {a, d, 0}, {b, d, 0}, {b, c, 0}}]]; 

    \ \ \ \ \ \ \ \ 
label\  = \ Graphics3D[Text["R", {2, 1.5, 0}]]; 

    \ \ \ 
spacer\  = \ 0.2; 

    threedims\  = \ 
\ \ \ \ \ \ 
      Graphics3D[{

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{a, 0, 0}, {b, 0, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["x", {b\  + \ spacer, 0, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, c, 0}, {0, d, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["y", {0, \ d\  + \ spacer, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, 0, 0}, {0, 0, 16}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["z", {0, 0, 16\  + \ 3\ spacer}]}]; 


    CMView\  = \ {2.7, 1.6, 1.2}; 


    Show[surfaceplot, Rplot, label, threedims, 
\ \ \ \ \ \ \ \ 
      ViewPoint -> CMView, PlotRange -> All, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ 
\)\)

Measure the volume of the solid consisting of all points between or on the rectangle and the surface plotted above.

That's duck soup. The top surface is just a plot of f[x, y] for points {x, y, 0} in the rectangle R in the xy-plane. The rectangle R consists of those points {x, y, 0} with a ≤ x ≤ b and c ≤ y ≤ d. The volume measurement is just
       (∫_c^d)
          (∫_a^b) f[x, y] dx 
            dy:

Integrate[f[x, y], {y, c, d}, {x, a, b}]

That's all there is to it.

T.2.a.ii)

Here's a plot of the surface
       z = f[x, y] = 5 - 0.3 ((x^2) + (y^2))
above region R in the xy-plane consisting of all points {x, y, 0} with {x, y} inside the ellipse
     ((((x/3)))^2)+((((((y - 1))/2)))^2)=1:

\(\(Clear[f, x, y, r, t]; 
    
{x[r_, t_], y[r_, t_]}\  = \ {0, 1}\  + \ r\ {3\ Cos[t], 2\ Sin[t]}; 

    f[x_, y_]\  = \ 5\  - \ 0.3\ \((x\^2 + y\^2)\); 


    surfaceplot\  = \ 
      ParametricPlot3D[{x[r, t], y[r, t], f[x[r, t], y[r, t]]}, 

        \ \ \ \ \ \ {t, 0, 2\ Pi}, {r, 0, 1}, 
\ \ \ \ \ \ 
        DisplayFunction -> Identity]; 
\ \ \ \ \ \ 

    Rplot\  = \ 
      ParametricPlot3D[{x[r, t], y[r, t], 0}, 

        \ \ \ \ \ \ {t, 0, 2\ Pi}, {r, 0, 1}, PlotPoints -> {Automatic, 2}, 

        \ \ \ \ \ \ DisplayFunction -> Identity]; 
\ \ \ \ \ \ 

    label\  = \ Graphics3D[Text["R", {\(-1\), 2, 0}]]; 
\ \ \ 

    spacer\  = \ 0.2; 

    threedims\  = \ 
\ \ \ \ \ \ 
      Graphics3D[{

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, 
            Line[{{\(-3\), 0, 0}, {3.5, 0, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["x", {3.5\  + \ spacer, 0, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, 
            Line[{{0, \(-1\), 0}, {0, 3.5, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["y", {0, \ 3.5\  + \ spacer, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Blue, Line[{{0, 0, 0}, {0, 0, 5.5}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["z", {0, 0, 5.5\  + \ spacer}]}]; 

    CMView\  = \ {2.7, 1.6, 1.2}; 


    Show[surfaceplot, Rplot, label, threedims, 
\ \ \ \ \ \ \ \ 
      ViewPoint -> CMView, PlotRange -> All, 
\ \ \ \ \ \ \ \ 
      DisplayFunction -> $DisplayFunction]; \ 
\)\)

Measure the volume of the solid consisting of all points between or on the ellipse and the surface plotted above.

The top surface is just a plot of f[x, y] for points {x, y, 0} inside the ellipse R in the xy-plane.  The volume measurement is just
       ∫(∫_R) f[x, y] dx 
            dy,
which you can easily calculate with the Gauss-Green formula once you remember the ellipse
       ((((x/3)))^2)+((((((y - 1))/2)))^2)=1
(which is the boundary of R) is parameterized in the counterclockwise way by:

\(Clear[x, y, t]; 
{x[t_], y[t_]} = {0, 1} + {3\ Cos[t], 2\ Sin[t]}\)

with t running from 0 to 2  π.  Now you are ready to turn the winning team of Mathematica and the Gauss-Green formula loose:
The Gauss-Green formula says
       ∫(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt
where a = 0 and b = 2  π. 
To calculate
       ∫(∫_R) f[x, y] dx 
            dy,
you just go with
       m[x, y] = 0 and
       n[x, y] = (∫_0^x) f[s, y] ds:

\(Clear[n, m, s]; 
m[x_, y_] = 0; 


    n[x_, y_] = ∫\_0\%x f[s, y] \[DifferentialD]s\)

This gives you f[x, y] = (∂_x) n[x, y] - (∂_y) m[x, y]:

Expand[{D[n[x, y], x]  -  D[m[x, y], x],  f[x, y]}]

Now you know that
       ∫(∫_R) f[x, y] dx 
            dy
        = (∫_a^b) m[x[t], y[t]] x'[t] + 
        n[x[t], y[t]] y'[t] dt 
where {x[t], y[t]} is the parameterization of the the boundary of R used above.

a = 0; \(b = 2\ π; 


    NIntegrate[m[x[t], y[t]]\ \(x'\)[t] + n[x[t], y[t]]\ \(y'\)[t], {t, a, b}]
      \)

The volume measures out to about 70.2 cubic units; that's all there is to it.

T.3)  Calculation strategies for minimum effort on your part

T.3.a) This situation signals integrating with respect to y first

Here's a region R plotted in the xy-plane.

\(Clear[x]; 

    Rplot = Plot[{x\ \((x - 2)\)\ E\^x, x\ \((2 - x)\)}, {x, 0, 2}, 
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        AxesLabel -> {"x", "y"}, 
        Epilog -> Text["R", {1, \(-1\)}]]; \)

How does this plot and its coding signal to you the minimum effort way of calculating  
       ∫(∫_R) f[x, y] dx 
            dy 
for a given function f[x, y]?

Things go best when you slice perpendicularly to the x-axis. In fact, when you look at the code above you can see that when you slice perpendicularly to the x-axis at a fixed x, then
       x (x - 2) (E^x) ≤ y ≤ x (2 - x):

\(Clear[yhigh, ylow, x, y]; 
yhigh[x_] = x\ \((2 - x)\); 

    ylow[x_] = x\ \((x - 2)\)\ E\^x; 

Clear[xslicer]; 

    xslicer[x_] = Graphics[Line[{{x, ylow[x]}, {x, yhigh[x]}}]]; 


    Show[Rplot, xslicer[0.8]]; \)

This tells you to integrate with respect to y first to get
       ∫(∫_R) f[x, y] dx 
            dy
        = 
        ∫(∫_R) f[x, y] dy 
              dx
        = 
        (∫_0^2)
            (∫_(ylow[x])^(yhigh[x]))f[x, y] dy 
              dx
because x runs from 0 to 2.
Here is the calculation in the case f[x, y] = 5  x + y:

\(Clear[f, x, y]; 
f[x_, y_] = 5\ x + y; 


    first = ∫\_\(ylow[x]\)\%\(yhigh[x]\)f[x, y] 
          \[DifferentialD]y\)

second = NIntegrate[first, {x, 0, 2}]

Or you can try to get it with one instruction:

N[(∫_0^2)
          (∫_(ylow[x])^(yhigh[x]))f[x, y] \[DifferentialD]y 
            \[DifferentialD]x]

Slicing perpendicularly to the y-axis would have resulted in a calculational nightmare; so you wouldn't want to try to integrate with respect to x first. 

Using the Gauss-Green formula would have been somewhat inconvenient because you would have to do quite a bit of thinking to come up with the required counterclockwise parameterization of the boundary.

T.3.b) This situation signals integrating with respect to x first

Here's another region R plotted in the xy-plane.

\(Clear[y]; 

    Rplot = 
	
      ParametricPlot[{{y\ \((y - 3)\), y}, {y\ \((3 - y)\)\ Sin[y]\^2, y}}, {
          y, 0, 3}, 
		PlotStyle -> {{Red, Thickness[0.01]}}, 
        AxesLabel -> {"x", "y"}, 
        Epilog -> Text["R", {0.3, 1}]]; \)

How does this plot and its coding signal to you the minimum effort way of calculating ∫(∫_R) f[x, y] dx 
            dy for a given function f[x, y]?

You get a calculationally clean answer when you slice perpendicularly to the y-axis. In fact, when you look at the code above you can see that when you slice perpendicularly to the y-axis at a fixed y, then
       y (y - 3) ≤ x ≤ 
        y (3 - y) Sin[y]^2.
Take a look:

\(Clear[xlow, xhigh, y]; 
xhigh[y_] = y\ \((3 - y)\)\ Sin[y]\^2; 

    xlow[y_] = y\ \((y - 3)\); 

Clear[yslicer]; 

    yslicer[y_] = Graphics[Line[{{xlow[y], y}, {xhigh[y], y}}]]; 


    Show[Rplot, yslicer[2.3]]; \)

This tells you to integrate with respect to x first to get
       ∫(∫_R) f[x, y] dx 
            dy
        = 
        (∫_0^3)
            (∫_(xlow[y])^(xhigh[y]))f[x, y] dx 
              dy
because y runs from 0 to 3.
Here is the calculation in the case f[x, y] = Sin[x + y]:

\(Clear[f, x, y]; 
f[x_, y_] = Sin[x + y]; 


    first = ∫\_\(xlow[y]\)\%\(xhigh[y]\)f[x, y] 
          \[DifferentialD]x\)

second = NIntegrate[first, {y, 0, 3}]

Or you can try a single instruction:

N[(∫_0^3)
          (∫_(xlow[y])^(xhigh[y]))f[x, y] \[DifferentialD]x 
            \[DifferentialD]y]

Slicing perpendicularly to the x-axis would have resulted in a calculational nightmare; so you wouldn't want to try to integrate with respect to y first. 
Using the Gauss-Green formula would have been inconvenient because you would have to do a bit of thinking to come up with the required counterclockwise parameterization of the boundary.

T.3.c) This situation signals using the Gauss-Green formula

Here's another region R plotted in the xy-plane.

\(Clear[x, y, t]; 
x[t_] = Cos[t]\ \((2 + Cos[t] - 1\/6\ Cos[2\ t])\); 

    
y[t_] = Sin[t]\ \((1.5 - Sin[t] + 1\/8\ Sin[3\ t])\); 


    Rplot = ParametricPlot[{x[t], y[t]}, {t, 0, 2\ π}, 
		
        PlotRange -> All, PlotStyle -> {{Red, Thickness[0.01]}}, 
        AxesLabel -> {"x", "y"}, 
        Epilog -> Text["R", {1, \(-1\)}]]; \)

How does this plot and its code signal to you the minimum effort way of calculating ∫(∫_R) f[x, y] dx 
            dy for a given function f[x, y]?

The boundary of R is given parametrically, and that's what Gauss-Green is set-up to use.  Before you use it, you have to check that the boundary curve is swept out exactly once in the counterclockwise direction.

\(\(ParametricPlot[{x[t], y[t]}, {t, 0, π\/2}, 
	
      PlotRange -> All, PlotStyle -> {{Red, Thickness[0.01]}}, 
      AxesLabel -> {"x", "y"}]; \)\)

\(\(ParametricPlot[{x[t], y[t]}, {t, π\/2, π}, 
	
      PlotRange -> All, PlotStyle -> {{Red, Thickness[0.01]}}, 
      AxesLabel -> {"x", "y"}]; \)\)

Yep, exactly once and counterclockwise.
Here's how to use the Gauss-Green formula to calculate 
       ∫(∫_R) f[x, y] dx 
            dy 
in the case that f[x, y] = x - 2  y: 
The Gauss-Green formula says
       ∫(∫_R)(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt
where a = 0 and b = 2  π. 
To calculate
       ∫(∫_R)(x - 2  y) dx 
            dy,
you just say
       f[x, y] = x - 2  y
       m[x, y] = 0 and
       n[x, y] = (∫_0^x) f[s, y] ds:

\(Clear[x, y, n, m, s]; 
f[x_, y_] = x - 2\ y; 
m[x_, y_] = 0; 

    n[x_, y_] = ∫\_0\%x f[s, y] \[DifferentialD]s\)

This gives you f[x, y] = (∂_x) n[x, y] - (∂_x) m[x, y]:

{(∂_x) n[x, y] - (∂_x) m[x, y], f[x, y]}

Now you know that
       ∫(∫_R) f[x, y] dx 
            dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt 
where {x[t], y[t]} is a parameterization of the boundary of R used above.

a = 0; b = 2 π; \(Clear[x, y, t]; 
x[t_] = Cos[t]\ \((2 + Cos[t] - 1\/6\ Cos[2\ t])\); 

    
y[t_] = Sin[t]\ \((1.5 - Sin[t] + 1\/8\ Sin[3\ t])\); 


    NIntegrate[m[x[t], y[t]]\ \(x'\)[t] + n[x[t], y[t]]\ \(y'\)[t], {t, a, b}]
      \)

That's it. 
If you don't go with the Gauss-Green formula, you've got to slice perpendicularly to the x-axis or the y-axis. Either way, you're going to have trouble because of the bulges and hollows on the boundary of R.

Show[Rplot];

T.3.d) This situation signals using the Gauss-Green formula

Given that the region R consists of everything inside the ellipse
       (((((x - 1)/5)))^2)+(((((y + 2)/3)))^2)=1 
in the xy-plane and
       f[x, y] = y Sin[x], 
how does the set up signal to you the minimum effort way of calculating
       ∫(∫_R) f[x, y] dx 
            dy?

This is a natural for the Gauss-Green formula because you can write down a quick parameterization of the ellipse
       (((((x - 1)/5)))^2)+(((((y + 2)/3)))^2)=1. 
Just go with the parameterization below and see the ellipse unfold right before your eager eyes:

\(Clear[x, y, t]; 
{x[t_], y[t_]} = {1 + 5\ Cos[t], \(-2\) + 3\ Sin[t]}; 
    

Rplot = 
      ParametricPlot[{x[t], y[t]}, {t, 0, 2\ π}, 
		
        PlotStyle -> {{Red, Thickness[0.01]}}, 
        AxesLabel -> {"x", "y"}]; \)

Here's how to use the Gauss-Green formula to calculate 
       ∫(∫_R) f[x, y] dx 
            dy 
in the case that f[x, y] = y Sin[x]: 
The Gauss-Green formula says
       ∫(∫_R)(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt
where a = 0 and b = 2  π. 
To calculate
       ∫(∫_R) y Sin[x] dx 
            dy
You just say
       f[x, y] = y Sin[x]
       m[x, y] = 0 and
       n[x, y] = (∫_0^x) f[s, y] ds:

\(Clear[x, y, n, m, t]; 
f[x_, y_] = y\ Sin[x]; 
\(m[x_, y_] = 0; 
    \)\)n[x_, y_] = (∫_0^x) f[t, y] \[DifferentialD]t

This gives you f[x, y] = (∂_x) n[x, y] - (∂_y) m[x, y]:

{(∂_x) n[x, y] - (∂_y) m[x, y], f[x, y]}

Now you know that
       ∫(∫_R) f[x, y] dx 
            dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt
where {x[t], y[t]} is the parameterization of the boundary of R used above.

a = 0; b = 2 π; \(Clear[x, y, t]; 
{x[t_], y[t_]} = {1 + 5\ Cos[t], \(-2\) + 3\ Sin[t]}; 
    

NIntegrate[
      m[x[t], y[t]]\ \(x'\)[t] + n[x[t], y[t]]\ \(y'\)[t], {t, a, b}, 
      AccuracyGoal -> 3]\)

Quick and sweet.

T.3.e) This situation signals integrating with respect to y first

Here's another region R plotted in the xy-plane.

\(Clear[high, low, x]; 
high[x_] = 3 + x + Cos[2\ x]; 

    low[x_] = x\ Sin[x]; 


    Rplot = Plot[{high[x], low[x]}, {x, 1, 4}, 
		PlotStyle -> Red, 
        AxesLabel -> {"x", "y"}, 
        Epilog -> {{Red, Line[{{1, low[1]}, {1, high[1]}}]}, {Red, 
              Line[{{4, low[4]}, {4, high[4]}}]}, Text["R", {3, 3}]}]; 
    \)

How does this plot and its coding signal to you the minimum effort way of calculating ∫(∫_R) f[x, y] dx 
            dy for a given function f[x, y]?

You get a calculationally clean answer when you slice perpendicularly to the x-axis. In fact, when you look at the code above you can see that when you slice perpendicularly to the y-axis at a fixed x, then
       low[x] ≤ x ≤ high[x].
Take a look:

\(Clear[xslicer]; 

    xslicer[x_] = Graphics[Line[{{x, low[x]}, {x, high[x]}}]]; 


    Show[Rplot, xslicer[3.6]]; \)

This tells you to integrate with respect to y first to get
       ∫(∫_R) f[x, y] dx 
            dy
        = 
        ∫(∫_R) f[x, y] dy 
              dx
        = 
        (∫_1^4)
            (∫_(low[x])^(high[x]))f[x, y] dy 
              dx
because x runs from 1 to 4.
Here is the calculation in the sample case f[x, y] = (E^(2  x  -  y)):

\(Clear[f]; 
f[x_, y_] = E\^\(2\ x - y\); 

    first = ∫\_\(low[x]\)\%\(high[x]\)f[x, y] \[DifferentialD]y\)

second = NIntegrate[first, {x, 1, 4}]

Or you can try one instruction:

\(Clear[f, x, y]; 
f[x_, y_] = E\^\(2\ x - y\); 

    N[∫\_1\%4
          \(∫\_\(low[x]\)\%\(high[x]\)f[x, y] \[DifferentialD]y 
            \[DifferentialD]x\)]\)

Slicing perpendicularly to the y-axis would have resulted in a calculational nightmare; so you wouldn't want to try to integrate with respect to x first.  
Using the Gauss-Green formula would have been inconvenient because you would have to do quite a bit of work to come up with the required counterclockwise parameterization of the boundary.

T.4)  Gauss-Green when you have a clockwise parameterization

T.4.a)

What happens when you try to use the Gauss-Green formula but you find that your parameterization moves in the clockwise way rather than the required counterclockwise way?

Go ahead and use the clockwise parameterization to make the calculation, but adjust your answer by multiplying it by -1.
The reason this works is that
       (∫_a^b) g[t] dt = 
        -(∫_b^a) g[t] dt.

T.5)  Gauss's normal law in 2D and using it to plan bombing runs

Look at this:

\(\(Clear[x, bell, y]; 

    bell[x_]\  = E\^\(-\(x\^2\/2\)\)\/\@\(2\ π\); 

    bellshaped\  = \ 
      Plot3D[bell[x]\ bell[y], \ {x, \(-3\), 3}, {y, \(-3\), 3}, 

        DisplayFunction -> Identity]; 

spacer\  = \ 0.2; 
h\  = \ 3.5; 

    threedims\  = \ 
\ \ \ \ \ \ 
      Graphics3D[{

          \ \ \ \ \ \ \ \ \ \ \ \ {Red, Line[{{\(-h\), 0, 0}, {h, 0, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["x", {h\  + \ spacer, 0, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Red, Line[{{0, \(-h\), 0}, {0, h, 0}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["y", {0, \ h\  + \ spacer, 0}], 

          \ \ \ \ \ \ \ \ \ \ \ \ {Red, Line[{{0, 0, 0}, {0, 0, 0.18}}]}, 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
          Text["z", {0, 0, 0.18\  + \ spacer/10}]}]; \ 

    \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 

    CMView\  = \ {2.7, \ 1.6, \ 1.2}; 
\ 

    Show[threedims, bellshaped, PlotRange -> All, 
\ \ \ \ 
      ViewPoint -> CMView, Boxed -> False, BoxRatios -> {1, 1, 1/2}, 

      \ \ \ \ DisplayFunction -> $DisplayFunction]; \ \ \ \ \ \ \ \ \ \ \ 
    \)\)

T.5.a.i)

What is this?

It's the bell-shaped surface! 
You get it by taking the bell-shaped curve function
       bell[x]=((1/(Sqrt[2  π])))(E^(((-(x^2))/2)))
and plotting the surface
       z = f[x, y] = bell[x] bell[y].

T.5.a.ii)

What's the version of Gauss's normal law corresponding to the bell-shaped surface?

When a numerical measurement is normally distributed, then the fraction of the time that this measurement comes out between given numbers a and b is estimated by
       (∫_a^b) normal[x, mean, dev] dx 
where the parameters mean and dev (= standard deviation) come from the data.
Here
       normal[x, mean, dev]=((bell[((x - mean)/dev)])/dev).

Some measurements come in pairs of numbers {x, y}. 
The version of Gauss's normal law corresponding to bell-shaped surfaces is this:
If you know that
→   the number in the x slot is normally distributed with mean = xmean and standard deviation = xdev and it's known that
→   the number in the y slot is normally distributed with mean = ymean and standard deviation = ydev, and
→   the measurements x and y have absolutely no influence on each other, 
then you can estimate the fraction of these measurements that come out in a given region R in the xy-plane by calculating the double integral
       ∫(∫_R) f[x, y] dy 
            dx
where
       f[x, y] = normal[x, xmean, xdev] normal[y, ymean, ydev].

T.5.b.i) Gauss's normal law and the bomb damage assessment

You've got a target sitting at {0, 0} and you want to knock it out with an iron bomb. The Northrup bombing sights are pretty good. In fact, if you set the sight to hit the target at {0, 0} and examine the point {x, y} at which an individual bomb lands, you will find that the resulting coordinates x and y
→   have no influence on each other;
→   are both normally distributed; and 
→   both average out to {0, 0}.
You will also find that the standard deviation of the x measurement is not the same as the standard deviation of the y measurement.
For a specific bomb and airplane and altitude, it turns out that 
       xdev = 100 feet and ydev = 300 ft.
And to destroy the target you are trying to knock out, the bomb must hit within 150 feet of the target. 
Estimate the fraction of these bombs dropped on this target that actaully destroy the target.

Load up Mathematica:

{xmean, ymean} = {0, 0}

{xdev, ydev} = {100, 300}

\(Clear[x, bell]; 
bell[x_] = E\^\(-\(x\^2\/2\)\)\/\@\(2\ π\)\)

\(Clear[normal, mean, dev]; 

    normal[x_, mean_, dev_] = bell[\(x - mean\)\/dev]\/dev\)

In order for a bomb to detroy the target at {0, 0}, the bomb has to land inside the region R inside the circle 
       (x^2) + (y^2) = (150^2).
To estimate the fraction of bombs that do this you calculate
       ∫(∫_R) f[x, y] dy 
            dx
where:

\(Clear[f]; 

    f[x_, y_] = normal[x, xmean, xdev]\ normal[y, ymean, ydev]\)

It is very easy to parameterize the circle (x^2) + (y^2) = (150^2) in accordance with the restrictions imposed by the Gauss-Green formula:

\(Clear[t]; 
{x[t_], y[t_]} = 150\ {Cos[t], Sin[t]}\)

with 0 ≤ t ≤ 2  π.
The Gauss-Green formula says
       ∫(∫_R)(∂_x) n[x, y] - 
        (∂_y) m[x, y] dx dy
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt
where a = 0 and b = 2  π. 
To calculate
       ∫(∫_R) f[x, y] dx 
            dy,
you just say
       m[x, y] = 0 and
       n[x, y] = (∫_0^x) f[s, y] ds:

\(Clear[n, m, s]; 
m[x_, y_] = 0; 

    n[x_, y_] = ∫\_0\%x f[s, y] \[DifferentialD]s\)

This gives you
       f[x, y] = (∂_x) n[x, y] - (∂_y) m[x, y]:

{(∂_x) n[x, y] - (∂_y) m[x, y], f[x, y]}

Now you know that
       ∫(∫_R) f[x, y] dx 
            dy
        = 
        ∫(∫_R)
              ((∂_x) n[x, y] - (∂_y) m[x, y]) 
              dx dy
       (because m[x, y] = 0)
        = 
        (∫_a^b) m[x[t], y[t]] x'[t] + 
          n[x[t], y[t]] y'[t] dt 
where {x[t], y[t]} is a parameterization of the the boundary of R used above.

a = 0; \(b = 2\ π; 


    onehit = NIntegrate[
        m[x[t], y[t]]\ \(x'\)[t] + n[x[t], y[t]]\ \(y'\)[t], {t, a, b}, 
        AccuracyGoal -> 3]\)

A little more than 1 out of 4 bombs dropped on the target will actually destroy the target.

T.5.b.ii)

In part i) above, what are your chances of destroying the target if you bomb it twice?
How many times do you have to bomb the target to get 95  % confidence of destroying it?

The fraction of the time that you expect to hit the target with one bomb is:

onehit

The fraction of the time you expect to miss with one bomb is:

onemiss = 1 - onehit

The fraction of the time you expect to miss with two bombs is:

(onemiss^2)

The fraction of the time you expect to miss with three bombs is:

(onemiss^3)

The fraction of the time you expect to destroy the target with k bombs is:

\(Clear[destroy, k]; 
destroy[k_] = 1 - onemiss\^k\)

To find out how many bombs you have to drop to get 95  % confidence of destroying the target, look at the following table:

Table[{k "bombs", destroy[k]}, {k, 1, 12}]

It will take 9 bombs to be 95  % sure that the target is destroyed.