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

Question

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]]
Graphics`Mesh`MeshInit[];
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!!!

Answer 1

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},
  Method -> "AdaptiveMonteCarlo"
  ]

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

Answer 2

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}} *)