I need a function that tells me, what transmission a gyrating particle has through a rectangular aperture if you consider all phase angles. So for a specific guiding center position and gyration radius, only part of the circle can go through the opening. I have to evaluate this function a lot of times (in the scope of a multidimensional numerical integration), so I want it to be as effective as possible.

TLDR: What is the best (fastest, acceptable accuracy) approach, to calculate a fraction of a circle through a rectangular region?

The main option I have investigated so far:

Define an ImplicitRegion for a rectangle + a circle and then do NIntegrate (or RegionMeasure - which should effectively be the same but is slower) on the region.

The ImplicitRegion includes parameters for the rectangular opening (xA and yA), offsets of the coordinate system (xOff, yOff), guiding center positions (xGC and yGC) as well as a radius r.

CircleInRectangle[xA_, yA_, xOff_, yOff_, xGC_, yGC_, r_] =
  -xA/2 + xOff <= x <= xA/2 + xOff && 
  -yA/2 + yOff <= y <= yA/2 + yOff &&
  {x, y} \[Element] Circle[{xGC, yGC}, r]
  {x, y}]

A visualization for an examplatory parameter set looks like the following:

Central guiding center

I can then either do RegionMeasure of the region or a NIntegrate[1,{x,y} \[Element] region] divided by $2\pi r$. With those methods, I can also control the accuracy of my result (with PrecisionGoal), which is good!

A mediocre second option I investigated:

Create a list of equidistant (with regards to an angle phi of the circle) points on the circle and check for all of them, how many are inside the rectangle (with Count). The accuracy of course depends on the number of points I use.


For both options, I also wrap around a Piecewise function to directly return 0 if the circle is entirely out of range or if the gyration radius is very small, then I just return 1 or 0.


For some reason, the evaluation time of Option one doesn't significantly increase if I increase PrecisionGoal of the integration. For a set of multidimensional parameter points of N=9801, the calculation time is between 90 and 92 seconds (for PrecisionGoal 1 to 4)!

For option two, the time increases of course linearly with the number of steps. At approximately 500 steps, the time is equivalent to Option one. In general, I would expect the accuracy to be a little less than Option one. One could even give an analytic upper bound for the accuracy of option two. First estimations suggest, that I would need more than the 500 points to get a relative precision of $10^{-3}$.

Possible 3rd option

I could in principle make a lot of Which cases and then calculate the fraction of the circular arc analytically, but I think the number of cases is very high and complex, which is why I haven't tried that yet.

My question: What is the fastest (with still acceptable accuracy of maybe $10^{-3}$ relative) way of calculating the transmission. Is it a fine tuning of one of the above options, or something completely different?

A last note: Compile

Would compiling the functions help?


I think your question doesn't involve any integration. You can describe any union of circle arcs with Interval of angles (union of angle intervals).

IntersectionAngles[x_] := Piecewise[{
  { {0, 0},      x<=-1},
  { {0, 2\[Pi]}, x>= 1},
  { {ArcCos[x], 2\[Pi]-ArcCos[x]},  True}}];

IntervalMod[{min_,max_}] := Interval@@Piecewise[{
  { {{0,max-2\[Pi]}, {min,2\[Pi]}}, min<=2\[Pi]<=max},
  { {Mod[{min,max},2\[Pi]]},  True}

    arr = {xOff-xGC+xA/2,
    int = IntervalIntersection@@Table[
  r RegionMeasure[int]

In[641]:= Timing[CircleInRectangle[2, 2, 0, 0, 0, 0, 1.1]]

Out[641]= {0.000399, 3.13015}

Since you are working with just arithmetics, you can control accuracy and precision as much as you like. If you enter precise values, you will get precise answer:

In[681]:= Timing[CircleInRectangle[15/8, 2, -1/3, -1, -1/16, -Sqrt[23/123], 1]]

Out[681]= {0.001364, (3 \[Pi])/2 - ArcCos[2/3] - ArcCos[Sqrt[23/123]]}
  • $\begingroup$ Thank you for the quick answer. There is a minor issue (piecewise needs another {}), but other than that, a very interesting solution. I have to look into it at test timing results. $\endgroup$ – Synchrotronus Jun 24 '19 at 14:56
  • $\begingroup$ I had some time to investigate your function. For very specific parameter sets, I got a RegionMeasure = 1 instead of 0. I had to provide RegionMeasure with the dimenions = 1, otherwise it would have counted the number of points for some reason (dim = 0). I am accepting your answer as it is very much faster than my solutions! Thank you $\endgroup$ – Synchrotronus Jun 25 '19 at 11:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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