Centripetal acceleration in Kepler orbit in cartesian coordinate

I have some difficult in deriving the formula for centripetal acceleration of a planet in solar system moving in elliptical orbit with constant areal velocity $\dot S$.

The book report the formula in cartesian coordinate:

$$a_c=\frac{4\dot S^2(1+(\frac{cy}{b^2})^2)}{((x+c)^2+y^2)}\frac{a^4b}{(a^4+c^2x^2)^\frac 3 2}$$

With:
a = semi major axis
b = semi minor axis
$c = \sqrt{a^2-b^2}$ semi focal distance

I try using polar coordinate find $a_c$ and then return to cartesian but i get stuck.

I guess that the last part is the curvature and the first is $r^2\dot\theta$ but have no idea for the others.

In that figure, the Sun is at the origin and the origin is chosen to coincide with the right focal point of the ellipse. Now there is enough information to give a Cartesian description of the ellipse, $$\frac{(x+c)^2}{a^2} + \frac{y^2}{b^2} =1$$ where $(x,y)$ is the position of the planet at some instant in time. Use the formula (Stewart Calculus 7th edition, p891), $$\vec{a} = a_T \vec{T} + a_{N}\vec{N} = \frac{|\vec{v}(t)\times \vec{a}(t)|}{|\vec{v}(t)|}\vec{N} + \frac{|\vec{v}(t)\cdot \vec{a}(t)|}{|\vec{v}(t)|}\vec{T}$$ where $\vec{N}$ is a unit vector directed towards the center of curvature and hence normal to the velocity and $\vec{T}$ is a unit vector in the direction of velocity. The goal is to compute $a_N$. To avoid polar coordinates and what not we use $a_N = \kappa v^2$ (Stewart Calculus p891 7th ed). From page 881 of Stewart Calculus 7th ed. the curvature, $$\kappa(x) = \frac{|f''(x)|}{(1+f'(x)^2)^{\frac{3}{2}}}$$ which explains the power of $\frac{3}{2}$ you have. Theoretically speaking this approach is sketchy since $y$ is clearly not a function of $x$ whereas this formula works for functions only; however, there comes the beauty of symmetry, if we consider only the upper half of the ellipse, which is a graph of a function, the radius of curvature there, at any given abscissa $x$, will be the same for the lower half of the ellipse at the same abscissa.
Using implicit differentiation one gets, $$\frac{(x+c)}{a^2} + \frac{yy'}{b^2} = 0 \Rightarrow y' = -\frac{-b^2(x+c)}{ya^2}$$ Using implicit differentiation twice one gets, $$\frac{1}{a^2} + \frac{1}{b^2}(y'^2+yy'')=0 \Rightarrow |y''| = \frac{\frac{b^2}{a^2}+y'^2}{y}= \frac{\frac{b^2}{a^2}+\frac{b^4(x+c)^2}{y^2a^4}}{y} =\frac{b^4}{a^2 y^3}$$ Thus, $$y'^2+1= y|y''|-\frac{b^2}{a^2} +1= \frac{b^4}{a^2y^2}-\frac{b^2}{a^2}+1=\frac{b^4-b^2y^2+a^2y^2}{a^2y^2}=\frac{c^2y^2+b^4}{a^2y^2}$$ Therefore, $$\displaystyle \kappa(x) = \frac{b^4a}{(b^4+c^2y^2)^{\frac{3}{2}}}$$ Now, in polar coordinates $\vec{v} = \dot{\rho}\vec{e}_{\rho} + \rho\dot{\phi}\vec{e}_{\phi} \Rightarrow v^2 = \dot{\rho}^2 +(\rho\dot{\phi})^2$.
The area swept between the angles $\phi_0$ and $\phi$ is $S(\phi) = \frac{1}{2}\int_{\phi_0}^{\phi} (\rho(\phi))^2 d\phi$ hence $\dot{S} =\frac{1}{2} \rho^2 \dot{\phi}$ which is a constant due to conservation of angular momentum. From the figure I attached one sees that $\rho + \rho' = 2a =\rho + \sqrt{(\rho \sin{\phi})^2+(2c+\rho \cos{\phi})^2}$ hence $$\displaystyle \rho = \frac{a(1-e^2)}{1+e\cos{\phi}}$$ where $e=\frac{c}{a}$ is the eccentricity of the ellipse. Differentiating and using the formula obtained for $\dot{S}$ , $$\dot{\rho} = \frac{2e\dot{S}}{a(1-e^2)}\sin{\phi}=\frac{2e\dot{S}}{\rho a(1-e^2)}\rho\sin{\phi}=\frac{2e\dot{S}} {a(1-e^2)}\frac{y}{\sqrt{x^2+y^2}}$$ Now, to calculate $\rho \dot{\phi} = 2\frac{\dot{S}}{\rho}=2\frac{\dot{S}}{\sqrt{x^2+y^2}}$. Therefore, $$v^2 = \dot{\rho}^2+(\rho\dot{\phi})^2 = \frac{4\dot{S}^2}{x^2+y^2}+\frac{4e^2\dot{S}^2}{a^2(1-e^2)^2}\frac{y^2}{x^2+y^2}$$ after some careful algebraic changes, $$v^2 = \frac{4\dot{S}^2\left( 1+\left(\frac{cy}{b^2} \right)^2\right)}{x^2+y^2}$$ and finally $$a_c =\frac{4\dot{S}^2\left( 1+\left(\frac{cy}{b^2} \right)^2\right)}{x^2+y^2}\frac{b^4 a}{(b^4+c^2y^2)^{\frac{3}{2}}}$$