Have you considered finding the intersections using an implicit form for the circles, $$\frac{x^2}{r^2} + \frac{y^2}{r^2} + ax + by + c = 0?$$ This representation doesn't have any coefficients that diverge as the circle approaches a straight line. To find intersections, you'll have to solve a quadratic equation whose leading coefficient could be zero or arbitrarily close to it, but the alternative form of the quadratic formula should be able to deal with that robustly.
You'll then have to do some jiggery-pokery to figure out whether the intersection points lie within the arcs. If the arc's bending angle is smaller than $\pi$, a projection onto the line joining the endpoints will suffice.
(Disclaimer: While all of this feels like it should work, I haven't analyzed it in any detail. Also, there could still be a problem when the circle is close to a line and you want the longer arc. But I can't imagine that's a case that would turn up in any practical application.)
Update: For a concrete example, here is the equation for a circular arc passing through the three points $(0,0)$, $(0.5, h)$, and $(1,0)$: $$\kappa^2 x^2 + \kappa^2 y^2 - \kappa^2 x - 2\eta y = 0,$$ where $$\begin{align}\kappa &= \frac{8h}{4h^2 + 1}, \\ \eta &= \frac{8h(4h^2-1)}{(4h^2+1)^2}.\end{align}$$ As you can see, the coefficients remain bounded as $h \to 0$.
Update 2: Wait, that equation becomes trivial if $h = 0$, which is bad. We really want something like $x^2/r + y^2/r + ax + by + c,$ i.e. multiply the previous expression through by $r$. Then for the same example, our equation becomes $$\kappa x^2 + \kappa y^2 - \kappa x - 2\eta' y = 0,$$ where $\eta' = (4h^2-1)/(4h^2+1)$. Here are some explicit values.
$h = 1/2$: $$2 x^2 + 2 y^2 - 2 x = 0,$$ $h = 0.01$: $$0.07997 x^2 + 0.07997 y^2 - 0.07997 x + 1.998 y = 0,$$ $h = 0$: $$2 y = 0.$$
By the way, in this format, the linear terms will always be simply $-2(x_0/r)x$ and $-2(y_0/r)y$, where the center of the circle is at $(x_0,y_0)$. As the center goes to infinity but the endpoints remain fixed, these coefficients remain bounded and nonzero (i.e. not both zero).
Hint:
When the surface is at height $h$ from the bottom, hence height $H-h$ from the top, the gravity center of the soda is at $\dfrac h2$ and that of air at $h+\dfrac{H-h}2=\dfrac{H+h}2$.
The height of the global gravity center is the average of these two heights, weighted by the weights of the two parts.
$$h_g=\frac{d_sh}{d_sh+d_a(H-h)}\frac h2+\frac{d_a(H-h)}{d_sh+d_a(H-h)}\frac{H+h}2=\frac12\frac{(d_s-d_a)h^2+d_aH^2}{(d_s-d_a)h+d_aH}\\
=\frac12\frac{\delta\left(\dfrac hH\right)^2+1}{\delta\dfrac hH+1}=\frac12\left(\delta\dfrac hH-1\right)+\frac1{\delta\dfrac hH+1}.$$
Find the root of the derivative of this expression.
Best Answer
Let me show the most general way to find the answer.
For a curve (or an infinitesimally thin wire with uniform density, i.e. uniform linear mass distribution), the center of mass is at the centroid of the curve.
In the general 2D case, the centroid of a parametric curve $\vec{s}(t) = \left ( x(t) , y(t) \right )$, $t_0 \le t \le t_1$ is at $( \hat{x} , \hat{y} )$, $$\begin{cases} \hat{x} = \frac{1}{L} \displaystyle \int_{t_0}^{t_1} x(t) \, \delta(t) \, dt \\ \hat{y} = \frac{1}{L} \displaystyle \int_{t_0}^{t_1} y(t) \, \delta(t) \, dt \end{cases} \tag{1}\label{NA1}$$ where $\delta(t) \, dt$ is the arc length parameter at $t$, $$\delta(t) \, dt = \sqrt{ \left( \frac{ d\, x(t) }{ d t } \right )^2 + \left( \frac{ d\, y(t) }{ d\, t} \right) ^2 } \, dt$$ and $L$ is the total length of the curve, $$L = \int_{t_0}^{t_1} \delta(t) \, dt$$
In this particular case, we have a circular arc, $$\begin{cases} x(\theta) = r \cos(\theta) \\ y(\theta) = r \sin(\theta) \end{cases}$$ and therefore $$\delta(\theta) \, d\theta = \sqrt{ \left(-r \sin(\theta)\right)^2 + \left(r \cos(\theta)\right)^2 } \, d\theta = \sqrt{ r^2 \left( (\sin\theta)^2 + (\cos\theta)^2 \right) } \, d\theta = \sqrt{ r^2 } \, d\theta = r \, d\theta$$ The arc distends one third of a full circle, or 120°. If we put the center of the circle at origin, and the midpoint of the arc on the positive $y$ axis, then $\theta$ ranges from $90°-120°/2 = 30°$ to $90°+120°/2 = 150°$, i.e. from $\theta = \pi/6$ radians to $\theta = 5 \pi/6$ radians.
The length $L$ of the circular arc we already know from geometry; it is one third of the perimeter of the circle of radius $r$, $$L = \frac{2 \pi r}{3}$$
Substituting these to $\eqref{NA1}$ we get $$\begin{cases} \hat{x} = \frac{3}{2 \pi r} \displaystyle\int_{\pi/6}^{5\pi/6} r \cos(\theta) \, r \, d\theta \\ \hat{y} = \frac{3}{2 \pi r} \displaystyle\int_{\pi/6}^{5\pi/6} r \sin(\theta) \, r \, d\theta \end{cases}$$ which simplify to $$\begin{cases} \hat{x} = \frac{3 r}{2 \pi} \displaystyle\int_{\pi/6}^{5\pi/6} \cos(\theta) \, d\theta = \frac{3 r}{2 \pi} \left(\Bigl[-\sin\theta \Bigr]_{\pi/6}^{5\pi/6} \right) \\ \hat{y} = \frac{3 r}{2 \pi} \displaystyle\int_{\pi/6}^{5\pi/6} \sin(\theta) \, d\theta = \frac{3 r}{2 \pi} \left(\Bigl[\cos\theta \Bigr]_{\pi/6}^{5\pi/6} \right) \end{cases}$$ Because $-\sin(\pi/6) - -sin(5\pi/6) = 0$, $\hat{x} = 0$. Which is completely expected, because we arranged the arc to be symmetric around the $y$ axis. Because $\cos(\pi / 6) - \cos(5\pi / 6) = \sqrt{3}/2 - -\sqrt{3}/2 = \sqrt{3}$, $$\hat{y} = \frac{3 r}{2 \pi} \sqrt{3} = \frac{3 \sqrt{3}}{2 \pi} r$$
In the case of $r = 3$, $$\hat{y} = \frac{9 \sqrt{3}}{2 \pi} \approx 2.48$$
This is in perfect agreement with King Tut's answer.