Finding area under a region of a Parametric Plot using mesh functions

Question

I have the following Mathematica code:

$Assumptions = {r \[Element] Reals, r >= 0, rh \[Element] Reals, 
   rh > 0, g \[Element] Reals, g > 0, l \[Element] Reals, l > 0};
M[rh_, g_] := (1 + rh^2/l^2) (rh^2 + g^2)^(3/2)/(2 rh^2)
f[r_, rh_, g_] := 1 - (2 M[rh, g] r^2)/(g^2 + r^2)^(3/2) + r^2/l^2
T[rh_, g_] := D[f[r, rh, g], r]/(4 \[Pi]) /. r -> rh
Veff[r_] := f[r, rh, g]*(L^2/r^2)
r0 = Solve[Veff'[r] == 0, r, Reals][[2]]
\[Lambda][rh_, 
   g_] = (Sqrt[-((r^2 f[r, rh, g])/L^2) Veff''[r]]) /. r0 // 
   Simplify;
l = 1;
(tp1 = {#[[1]], \[Lambda][rh1 = (rh /. #[[2]]), 0.02]} &@
      Maximize[{T[rh, 0.02], 0.01 < rh < 1}, rh] // FullSimplify) // N;
(tp2 = {#[[1]], \[Lambda][rh2 = (rh /. #[[2]]), 0.02]} &@
     Minimize[{T[rh, 0.02], rh1 < rh < 1}, rh] // FullSimplify) // N

pplt = ParametricPlot[{T[rh, 0.02], \[Lambda][rh, 0.02]}, {rh, 0.01, 
    8}, PlotRange -> {{0, 1}, {0, 10}}, AspectRatio -> 0.5, 
   PlotPoints -> 50, ImageSize -> 800, 
   PlotTheme -> {"Scientific", "BoldColor"}, 
   FrameLabel -> {"Temperature \[Rule]", "Lyapunov Exponent \[Rule]"},
    LabelStyle -> {Bold, Black, 22}, PlotStyle -> Thickness[0.004],
   ColorFunction -> 
    Function[{T, \[Lambda], rh}, 
     If[rh <= rh1, Blue, If[rh <= rh2, Red, Green]]], 
   ColorFunctionScaling -> False,
   
   Prolog -> {{Gray, Dashing[0.008], Line[{tp1, {tp1[[1]], 0}}], 
      Line[{tp2, {tp2[[1]], 0}}]}, {Darker[Black], 
      Line[{{((tp1 + tp2)/2)[[1]], 0.5}, {((tp1 + tp2)/2)[[1]], 
         9}}]}, {Gray, Thick, Dashing[0.001], Circle[tp1, {0.01, .2}],
       Circle[tp2, {0.01, .2}]}}];

P3 = Legended[pplt, 
  Placed[SwatchLegend[{Blue, Red, Green}, {"Small BH", 
     "Intermediate BH", "Large BH"}, LabelStyle -> {Bold, Black, 20}, 
    LegendMarkers -> "SphereBubble", LegendMarkerSize -> 20, 
    LegendFunction -> (Framed[#, RoundingRadius -> 10, 
        FrameStyle -> None, Background -> GrayLevel[0.95]] &), 
    LegendLabel -> 
     Placed["Black Holes", Left, 
      Rotate[Style[#, 16], 90 Degree] &]], {0.2, .6}]]

which gives me the following result:

How can I find the area under Region 1 and Region 2 using Mesh functions ? I also wish to find the position of the line (Basically the value of T) for which the Region 1 and Region 2 have equal areas. What approach should I use for that, without succumbing to trial-and-error methods?

Answer 1

Hope I understood your question correctly.

Using P3 as start, first take all the points

pts = Cases[Normal[P3], Line[p_] :> p, -1][[1]];

Now take three subsets defined by range tp1,tp2

blue = Select[pts, #[[1]] <= tp1[[1]] && #[[2]] >= tp1[[2]] &];
red = Select[pts,tp2[[1]] <= #[[1]] <= tp1[[1]] && tp2[[2]] <= #[[2]] <= tp1[[2]] &] ;
green = Select[pts ,  #[[1]] >=   tp2[[1]] &&   #[[2]] <= tp2[[2]] &] ;

Interpolate these points

ipblue = Interpolation[blue];
ipred = Interpolation[red];
ipgreen = Interpolation[green];

and integrate

{Ablue, Ared, Agreen} =NIntegrate[{ipblue[tp], ipred[tp], ipgreen[tp]}, {tp,tp2[[1]],tp1[[1]]}] // Quiet

The two areas you 're looking for follow to

{Ablue - Ared, Ared - Agreen}
(*{2.82274, 1.9862}*)

Hope it helps!

addendum

calculation of midpoint T

reg1[T_?NumericQ] :=NIntegrate[ ipblue[tp] - ipred[tp] , {tp, T, tp1[[1]]} ]  
reg2[T_?NumericQ] := NIntegrate[  ipred[tp] - ipgreen[tp] , {tp, tp2[[1]], T} ]  

line1 = Line[Map[{#, reg1[#]} &, Subdivide[tp2[[1]], tp1[[1]], 20]]];
line2 = Line[Map[{#, reg2[#]} &, Subdivide[tp2[[1]], tp1[[1]], 20]]];
pT=RegionIntersection[line1, line2](*Point[{{0.669968, 0.707078}}]*)

T=pT[[1, 1]][[1]] (*0.669968*)

More direct solution using FindRoot (NDSolve and NMinimize won't evaluate)

T =.
FindRoot[reg1[T] == reg2[T], {T, tp2[[1]], tp1[[1]]}]   
(*{T -> 0.669965} *)

Answer 2

After answering your question, here I'll add a supplement which determines tp1,tp2 using MeshFunctions

plot = ParametricPlot[{T[rh, 0.02], \[Lambda][rh, 0.02]}, {rh, 0.01, 
   8}, PlotRange -> {{0, 1}, {0, 10}}, AspectRatio -> 0.5, 
  PlotPoints -> 50, ImageSize -> 800, 
  MeshFunctions -> 
   Function[{x, y, rh}, 
    Evaluate[Derivative[1, 0][ T][rh, 1/50]/1(*Derivative[ 1,
     0][\[Lambda]][rh,1/50]*)]], Mesh -> {{0}}] 

get points

 pts = plot[[1]][[1, 1]]; (* all points*)
 index = Cases[plot, Point[ p_] :>  p, -1][[1]]  (* index meshpoints*)
 tp=pts[[index]] (*{tp1,tp2}*)
 (*{{0.274999, 1.5824}, {0.951599, 7.78498}}*) 

Show[{ListPlot[pts , PlotRange -> All],Graphics[{Red,PointSize[Large], Point[tp]}]}, PlotRange -> {{0, 1}, {0, 10}}]