# Area between Contours in ContourPlot

I feel slightly foolish for asking this because I am so close, but I'm having trouble, so I will anyway.

I asked this question two days ago regarding finding the lengths of contours. Now, I'd like to find the areas between them. I have read and played around with the information here, here and here, but am having difficulties.

Having extracted the individual contours of f[x,y] via the method given by chris, I would like to find the areas bounded by the contours and the plot region. Is there a simple way to do this?

For reference, I am using the contour plot of the following function:

q[r_] := Piecewise[{{25/(0.1*1), r < 0.1}, {25/r, r >= 0.1}}]
phi[r_, t_] := (Pi/2) + q[r]*t
v[r_, t_] := q[r]*r*Cos[phi[r, t]]
s[x_] := Piecewise[{{x = -1, x < 0}, {x = 1, x >= 0}}]

f[x_,y_] := s[x]*v[Sqrt[x^2 + y^2],ArcTan[y/x]/q[Sqrt[x^2 + y^2]]]

ContourPlot[f[x, y], {x, -1, 1}, {y, -1, 1}, RegionFunction -> Function[{x, y}, x^2 + y^2 <= 1],
PlotPoints -> 100, Contours -> Range[-25, 25, 1]]

-
Off topic question, why there are x =... in s definition? It causes errors and problems with integration. – Kuba Apr 22 '14 at 6:32
@Kuba I simply wanted a function that was -1 when less than zero and 1 when greater than or equal to it. There's probably a more elegant method, but I'm new to Mathematica and a piecewise function was quick and dirty. – pirtle Apr 22 '14 at 18:34
I know, piecewise is ok, but compare it to the part of code in my answer. x =  is redundant. – Kuba Apr 22 '14 at 18:39
Ah, yes. I see that now. Duly noted. – pirtle Apr 23 '14 at 1:37

q[r_] := Piecewise[{{25/(0.1*1), r < 0.1}, {25/r, r >= 0.1}}]
phi[r_, t_] := (Pi/2) + q[r]*t
v[r_, t_] := q[r]*r*Cos[phi[r, t]]
s[x_] := Piecewise[{{-1, x < 0}, {1, x >= 0}}]

f[x_, y_] := s[x]*v[Sqrt[x^2 + y^2], ArcTan[y/x]/q[Sqrt[x^2 + y^2]]]


Here are two ways to go:

## 1

NIntegrate[Boole[(23 <= f[x, y] <= 24 && x^2 + y^2 <= 1)],
{x, -1, 1}, {y, -1, 1}]

0.118004


## 2

plot = RegionPlot[23 <= f[x, y] <= 24 && x^2 + y^2 <= 1, {x, -1, 1}, {y, -1, .1},
PlotPoints -> 100]


poly = Cases[Normal@plot, Polygon[n_] :> n, ∞]
GraphicsMeshMeshInit[];
PolygonArea /@ poly // Total

0.117933


You can work with precission by adjusting specific options for plot or integration.

-
That is so sleek! Thanks so much! – pirtle Apr 22 '14 at 18:31