# How to output the level of the contour plot based on the desired area enclosed in the contour level?

I am a beginner using Mathematica. I am doing a physics experiment that involves this 3D surface; the function $$f(x, y)$$. I need a program that can give me the contour level in the $$z$$-axis which encloses an area that I specify.

So far I have managed to create an algorithm that can accept a contour level and then output the area that it encloses:

Here is the code: (credits to Kuba :Area between Contours in ContourPlot)

f[x_, y_] :=
1/0.005*2.86*10^(-5)*
Sin[Pi/2 -
ArcTan[(0.04367/35)*((Sin[2*Pi*35*y])^2)*2*(Pi/0.00763)*
Cos[2*(Pi/0.00763)*x]]]*(2*
Pi*35*(0.04367/35)*((Sin[2*Pi*35*y])^2)*Sin[2*(Pi/0.00763)*x]*
Cos[2*Pi*35*
y] + (2*9.81*(0.023 -
2*1.9*10^(-3) - (0.04367/35)*((Sin[2*Pi*35*y])^2)*
Sin[2*Pi/0.00763*x]))^0.5)

plot = RegionPlot[0.0036 <= f[x, y], {x, 0, 0.007632}, {y, 0, 1/35},
PlotPoints -> 100]

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


The value 0.0036, seen in RegionPlot, is the contour level in this case. And the program outputs an area value of 0.0000394488 when I run it. But how do I specify an area value and then get the contour level out of the program?

This is for my school project which is due soon; any help is immensely appreciated!!!

• The solution in the question was originally written by Kuba here. – C. E. Dec 8 '18 at 15:26
• Hi, thanks for the reply! But I referred to Kuba's code in writing my code above, but its essentially the same thing; they go from a specified contour level to the area, but for my experiment, I need to find the contour level that has a specified area enclosed in it. So its the reverse process. – Steve Gordan Dec 8 '18 at 15:31
• I am so sorry Sir! I had no idea. – Steve Gordan Dec 8 '18 at 15:43
• Write a function area[level] that returns the area of f[x, y] > level and use FindRoot[] to find the desired level – Michael E2 Dec 8 '18 at 16:06
• @MichaelE2 Thank you! However, I am terrible at coding. Please guide me further. This is what I have so far : f[x_, y_] := ..... area[level_] = RegionPlot[level <= f[x, y], {x, 0, 0.007632}, {y, 0, 1/35}, PlotPoints -> 100] poly = Cases[Normal@plot, Polygon[n_] :> n, [Infinity]] GraphicsMeshMeshInit[]; PolygonArea /@ poly // Total FindRoot[area[0.000039], {x = 0}] – Steve Gordan Dec 8 '18 at 16:45

Let me just add my five cents on this problem. The idea that I pursued was the same as Bob but I tried to make it faster by using one of the Monte Carlo methods of NIntegrate instead of RegionPlot.

This is my function for finding the area in terms of level threshold:

area[thresh_] := NIntegrate[
Boole[f[x, y] >= thresh],
{x, 0, xmax},
{y, 0, ymax},
]

area[0.0036] // AbsoluteTiming


{0.030328, 0.0000390467}

This is quite fast and the result is almost the same as the one given by OP's method. Unfortunately, I found that FindRoot, NMinimize etc. do not play nice with this function, I'm guessing because it is stochastic. The Jacobian is singular, FindRoot says, for example. Using a similarly stochastic method of NMinimize such as SimulatedAnnealing also didn't work, it is still unclear to me why.

So I had to resort to a workaround which rather put a blemish on the whole thing, here it is:

table = Table[{Quiet@area[x], x}, {x, 0.003, 0.004, 0.001 0.003}];
areaInTermsOfLevel = Interpolation[table];
areaInTermsOfLevel[0.00003930850798691388]


0.00360097

This works but yeah, it's not very pretty, and also evaluating area a thousand times is pretty inefficient. Another problem is finding the limits 0.003 and 0.004. The reason is that e.g. the largest level allowed is too large (0.005 say) then the area will be 0 for multiple values, and that will cause problems for Interpolation. (Bob's answer shows how to get the thresholds that should work, in theory. The smallest and largest values of f[x, y])

• THANK YOUU!!! What do I have to do achieve your coding prowess? Please tell me. The way you throw this jargon around is mesmerizing. I want to be like you. What's your secret? What did you study? – Steve Gordan Dec 8 '18 at 17:29
• @SteveGordan I mostly learned by reading and contributing here on Stack Exchange, a well as by studying physics at the university. – C. E. Dec 8 '18 at 18:19

As suggested by Micheal E2, use FindRoot

\$Version

(* "11.3.0 for Mac OS X x86 (64-bit) (March 7, 2018)" *)

Clear["Global*"]

f[x_, y_] :=
Evaluate@Simplify[
Rationalize[
1/0.005*2.86*10^(-5)*
Sin[Pi/2 -
ArcTan[(0.04367/35)*((Sin[2*Pi*35*y])^2)*2*(Pi/0.00763)*
Cos[2*(Pi/0.00763)*x]]]*(2*Pi*35*(0.04367/35)*((Sin[2*Pi*35*y])^2)*
Sin[2*(Pi/0.00763)*x]*
Cos[2*Pi*35*
y] + (2*9.81*(0.023 -
2*1.9*10^(-3) - (0.04367/35)*((Sin[2*Pi*35*y])^2)*
Sin[2*Pi/0.00763*x]))^0.5), 0]]

{xmin, xmax} = {0, 477/62500};
{ymin, ymax} = {0, 1/35};
{fmin, fmax} = #[{f[x, y], xmin <= x <= xmax, ymin <= y <= ymax}, {x,
y}] & /@ {NMinValue, NMaxValue}

(* {0.00244846, 0.00419155} *)

area[level_?NumericQ] := Module[{plot, polys},
plot = RegionPlot[level <= f[x, y], {x, xmin, xmax}, {y, ymin, ymax},
PlotPoints -> 100];
polys = Cases[Normal@plot, _Polygon, ∞];
RegionMeasure /@ polys // Total]


Checking with the level from your example

area1 = area[0.0036]

(* 0.0000394488 *)


Using FindRoot to reverse the process

threshold[rgnMeasure_?(0 < # < (xmax - xmin) (ymax - ymin) &)] :=
FindRoot[area[level] == rgnMeasure, {level, (fmax + fmin)/2}]


Unfortunately, this is quite slow

threshold[area1] // AbsoluteTiming

(* {110.345, {level -> 0.0036}} *)

• OHH YESSSS THANK YOU SOO MUCHHHHHHHHHH!!!!!!!!!!!! I am forever in your debt, Mr Bob Hanlon. – Steve Gordan Dec 8 '18 at 17:18
• Hi! Is there a way to decrease the precision or accuracy of the calculation to speed it up? Because I only need 4 significant figures – Steve Gordan Dec 10 '18 at 5:07
• @SteveGordan - In the RegionPlot replace the option PlotPoints -> 100 with MaxRecursion -> 3. That should cut the time to about 25% of the original time. If that is not enough, play with the options to FindRoot – Bob Hanlon Dec 10 '18 at 6:03
• Thanks! It's much faster now :))) – Steve Gordan Dec 10 '18 at 7:07