Przy wyprowadzeniu limitu Roche’a, wyznaczającego granicę przy której płynny satelita zostanie rozerwany przez siły pływowe pochodzące od ciała, które orbitruje, wykorzystamy poradę Feynmanna:


The same equations have the same solutions 

Okazuje się, że równania są niemal identyczne jak dla elipsoidy Maclaurina. 

Szkic innego wyprowadzenia: https://dlab.epfl.ch/wikispeedia/wpcd/wp/r/Roche_limit.htm

(*G=ρ=M=1*)

Clear[Ω,G,a,b,c,m,r]

R={x,y,z} (* promień wodzący *)

o={a,b,c} (* Osie elipsoidy *)

A={Ax,Ay,Az} (* Współczynniki do policzenia *)

b=c (* Zakładamy, że elipsoida jest obrotowa, ale inaczej niż dla elipsoidy Maclaurina ! *)

Block[{a=1,c=1/2},
ContourPlot3D[
x^2/a^2+y^2/b^2+z^2/c^2==1,{x,-2,2},{y,-2,2},{z,-2,2},
Axes->True,AxesOrigin->{0,0,0},Boxed->False,PlotPoints->64,AxesLabel->{"x,a", "y,b", "z,c"}
]
]

A = Table[a*b*c*Integrate[1/(Sqrt[(a^2 + u)*(b^2 + u)*(c^2 + u)]*(o[[i]]^2 + u)), {u, 0, Infinity}, 
     Assumptions -> a > 0 && b > 0 && c > 0&&a>c], {i, 1, 3}]

Φtidal=-3/2 Ω^2 x^2 (* Potencjał siły pływowej, (Ω^2)= GM/(d^3) *)

Φg=-Pi G ρ Sum[(o[[i]]^2-R[[i]]^2) A[[i]],{i,1,3}]//Simplify (* Potencjał WEWNĄTRZ elispoidy *)

const=(Φg+Φtidal/.{x->0,y->0,z->c}//Simplify) (* Stała C wyznaczona na północnym biegunie ,,księżyca'' *)

ellipsoid=Φg+Φtidal-const (* Równanie równowagi na powierzchni powinno być równoważne równaniu elipsoidy, czyli:

(x^2)/(a^2)+(y^2)/(b^2)+(z^2)/(c^2)==1

Uzyskana forma jest daleka od uczynienia powyższego oczywistym...
*)

(*CoefficientList[ellipsoid,{x,y,z}]//Simplify*)

α=D[ellipsoid,{x,2}]/2 (* Obliczamy wyraz stojący przy (x^2) *)

β=D[ellipsoid,{y,2}]/2(* Obliczamy wyraz stojący przy (y^2) *)

γ=D[ellipsoid,{z,2}]/2//Simplify (* Obliczamy wyraz stojący przy (z^2) *)

ζ=-(ellipsoid-α x^2- β y^2- γ z^2//Simplify) (* Obliczamy wyraz wolny *)

eqc=ζ/γ==c^2 //Simplify(* Spr. czy wyraz stojący przy (z^2) to, tak jak oczekujemy, (c^(2 ))*)

eqa=a^2==ζ/α//Simplify

eqb=b^2==ζ/β//Simplify

eqVol=4/3 Pi a b c ρ == m (* Zakładamy, że w wyniku deformacji pływowej masa/objętość pozostaje niezmienna *)

r=(m/(4/3 Pi ρ))^(1/3) (* Promień kuli w przypadku braku siły pływowej *)

{eqa,eqVol}/.c->ε a //Simplify(* ... i wprowadzamy spłaszczenie ε jak w modelu Roche'a  *)

Simplify[%,a>0]

Ω^2/.Solve[%[[1]],Ω][[1]]

FullSimplify[%,1>ε>0]

χ=%/(2 Pi G ρ)

Plot[χ,{ε,0,0.999}]

Simplify[D[χ,ε],1/4<ε<3/4]

ArcSech[0.5]

FindInstance[8-8 ε+(ε^2)-(ε^3)+(Sqrt[((1-ε)/(1+ε))]) (-4-6 (ε^2)+(ε^4)) ArcSech[ε]==0&&1/4<ε<3/4,ε]

N[%,32]

χ0=χ/.%[[1]]

N[χ0,32]

3/64 (* Ciekawostka: wynik jest bardzo bliski 3/64 *)

%-%%

CubeRoot[2/3/χ0]//N (* Słynny współczynnik w limicie Roche'a wyprowadzony ! *)

(* Jako mnemotechniczny sposób zapamiętania, graniczny kształt na wydłużenie 50%, co daje identyczny współczynnik *)

χ1=χ/.ε->1/2//FullSimplify//TrigToExp//FullSimplify

N[χ1/χ0-1]

CubeRoot[2/3/χ1]//N (* Słynny współczynnik w limicie Roche'a wyprowadzony ! *)

N[χ1,32]

0.090093/2 (* Chandrasekhar *)

CubeRoot[2/3/%]