2
$\begingroup$

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.

The contour plot looks something like this...what are the areas?!

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]]
$\endgroup$
4
  • $\begingroup$ Off topic question, why there are x =... in s definition? It causes errors and problems with integration. $\endgroup$
    – Kuba
    Apr 22, 2014 at 6:32
  • $\begingroup$ @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. $\endgroup$
    – pirtle
    Apr 22, 2014 at 18:34
  • 1
    $\begingroup$ I know, piecewise is ok, but compare it to the part of code in my answer. x = is redundant. $\endgroup$
    – Kuba
    Apr 22, 2014 at 18:39
  • $\begingroup$ Ah, yes. I see that now. Duly noted. $\endgroup$
    – pirtle
    Apr 23, 2014 at 1:37

1 Answer 1

5
$\begingroup$
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

referrence link

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

enter image description here

poly = Cases[Normal@plot, Polygon[n_] :> n, ∞]
Graphics`Mesh`MeshInit[];
PolygonArea /@ poly // Total  
0.117933

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

$\endgroup$
1
  • $\begingroup$ That is so sleek! Thanks so much! $\endgroup$
    – pirtle
    Apr 22, 2014 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.