Model Roche’a można stosować dla ciał, w których materia jest bardzo skoncentrowana w centrum, tak, aby przybliżenie masy punktowej dla potencjału grawitacyjnego miało sens. Otoczka wtedy jest traktowana jako praktycznie nieważka. Model ma sens głównie dla gwiazd, dla których kontrast gęstości jest duży, czyli powiedzmy n>=3.

Dla obiektów o niskim kontraście gęstości, w szczególności stałej gęstości, powyższy model nie ma sensu. Natomiast niemal natychmiast wpadamy w pułapkę problemu różniczkowego o nieznanym (swobodnym) brzegu. Chcielibyśmy zadać warunki brzegowe dla równania różniczkowego cząstkowego, lub podać granice całkowania w wersji całkowej, ale nie znamy ich a priori. Wynikną one dopiero z rozwiązania problemu. W przypadku stałej gęstości okazuje się, że mamy szczęście: kształt powierzchni jest elipsoidą. Dlatego często niżej zaprezentowane przykłady są częścią gałęzi astrofizyki określanej jako teoria elipsoidalnych figur równowagi.  Przez jej historię przewija się śmietanka nauki światowej: Maclaurin, Jacobi, Riemann, Dedekind, Poincare, Chandrasekhar.

(*G=ρ=M=1*)

Clear[Ω,G,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=a (* Zakładamy, że elipsoida jest obrotowa *)

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

Φc=-1/2 Ω^2 (x^2+y^2) (* Potencjał siły odśrodkowej *)

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

ΦgKula=Limit[Φg,c->a,Direction->"FromBelow",Assumptions->a>0] (* Dla c=a potencjał powinien przejść do tego dla jednorodnej kuli *)

ΦgKula/.{x->0,y->0,z->0}/.ρ->M/(4/3Pi a^3) (* Spr. w centrum *)

ΦgKula/.{x->a,y->0,z->0}/.ρ->M/(4/3Pi a^3) (* Spr. na powierzchni *)

const=(Φg+Φc/.{x->0,y->0,z->c}//Simplify) (* Stała C wyznaczona na biegunie *)

ellipsoid=Φg+Φc-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 *)

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

(* Pozostałe 2 warunki są identyczne, bo założyliśmy, że elipsoida jest obrotowa *)

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

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

eq1==eq2

eq3=4/3 Pi a b c ρ == M (* Zakładamy, że w wyniku obrotu masa/objętość pozostaje niezmienna *)

R0=(M/(4/3 Pi ρ))^(1/3) (* Promień kuli w przypadku braku rotacji *)

Solve[eq3,c][[1]]

eq1/.%//Simplify (* Eliminujemy oś biegunową za pomocą warunku na niezmienną objętość *)

Ω^2/.Solve[%,Ω][[1]] (* Wyliczamy kwadrat prędkości kątowej *)

X=FullSimplify[%,a>0] (* Upraszczamy .... *)

X/.M->4/3 Pi a^2 c ρ (* Zamieniamy masę z powrotem na gęstość... *)

%/.c->ε a (* ... i wprowadzamy spłaszczenie ε jak w modelu Roche'a  *)

Simplify[%,a>0]

(* Wynik końcowy odbiega od podręcznikowego, ale okazuje się równoważny *)

ArcSec[1/ε]-ArcCos[ε]//FullSimplify

Plot[(( ε (-3 ε (Sqrt[1-(ε^2)])+(1+2 (ε^2)) ArcSec[(1/ε)]))/(((1-(ε^2)))^(3/2)))//Sqrt,{ε,0,1},AxesLabel->{"ε=c/a","Ω/(Sqrt[2 π G ρ])"},BaseStyle->FontSize->18,PlotStyle->{Red,Thick},AxesStyle->Thick,Epilog->{Arrow[{{1,0.47},{0.5,0.47}}],Text["Obrót",{0.75,0.47},{0,1}]},AspectRatio->1,ImageSize->600]

(( ε (-3 ε (Sqrt[1-(ε^2)])+(1+2 (ε^2)) ArcSec[(1/ε)]))/(((1-(ε^2)))^(3/2)))==Ω^2/(2Pi G ρ)/.{Ω->Quantity[2Pi,"Days^-1"],G->Quantity[1,"GravitationalConstant"],ρ->Quantity[1,"EarthMass"]/(4/3Pi Quantity[1,"EarthEquatorialRadius"]^3)}

FindRoot[%,{ε,0.99}]

(1-ε)^-1/.%

(1-Quantity[1,"EarthPolarRadius"]/Quantity[1,"EarthEquatorialRadius"]//N)^-1

%/%%

Solve[(( ε (-3 ε (Sqrt[1-(ε^2)])+(1+2 (ε^2)) ArcSec[(1/ε)]))/(((1-(ε^2)))^(3/2)))==χ^2,ε,Reals]

Series[(( ε (-3 ε (Sqrt[1-(ε^2)])+(1+2 (ε^2)) ArcSec[(1/ε)]))/(((1-(ε^2)))^(3/2))),{ε,1,1}]//Normal

1-15/8Ω^2/(2Pi G ρ)/.{Ω->Quantity[2Pi,"Days^-1"],G->Quantity[1,"GravitationalConstant"],ρ->Quantity[1,"EarthMass"]/(4/3Pi Quantity[1,"EarthEquatorialRadius"]^3)}

(1-%)^-1

3/4/(8/15)

N[%]

Block[{M=1,ρ=1},Plot[{(((M^(1/3)) ((((3/π)))^(1/3)) (ε^(2/3)))/((2^(2/3)) (ρ^(1/3)))),(((M^(1/3)) ((((3/π)))^(1/3)))/((2^(2/3)) (ε^(1/3)) (ρ^(1/3)))),R0},{ε,0,1},AxesLabel->{"ε=c/a","c,a"},BaseStyle->FontSize->18,PlotStyle->{Red,Blue,{DotDashed,Thin,Black}},AxesStyle->Thick,Epilog->{Arrow[{{1,1},{0.5,1}}],Text["Obrót",{0.75,1},{0,-1}]},ImageSize->600,
PlotLegends->{"c","a","(Subscript[R, 0])"}]
]

(* Wzór na moment bezwładności elipsoidy obrotowej I = 2/5 M a^2 *)

Block[{M=1,ρ=1},Plot[((((M^(1/3)) ((((3/π)))^(1/3)))/((2^(2/3)) (ε^(1/3)) (ρ^(1/3)))))^2((( ε (-3 ε (Sqrt[1-(ε^2)])+(1+2 (ε^2)) ArcSec[(1/ε)]))/(((1-(ε^2)))^(3/2)))//Sqrt),{ε,0,1},AxesLabel->{"ε=c/a","J/(Subscript[J, 0])"},BaseStyle->FontSize->18,PlotStyle->{Red,Thick},AxesStyle->Thick,Epilog->{Arrow[{{1,0.47},{0.5,0.47}}],Text["Obrót",{0.75,0.47},{0,1}]},ImageSize->600]
]

Clear[Ω]

G=M=ρ=1

Ω=12/10

mclaurin=FindRoot[{eq1,eq3},{{a,R0},{c,1.1R0}},WorkingPrecision->64,MaxIterations->256]

mel =x^2/a^2+y^2/b^2+z^2/c^2==1/.mclaurin

ContourPlot3D[Evaluate@mel,{x,-1.3,1.3},{y,-1.3,1.3},{z,-1.3,1.3},PlotPoints->64]

mel2d=x^2/a^2+z^2/c^2==1/.mclaurin

Xs=(a/.mclaurin)/2

Zs=Z/.Solve[Xs^2/a^2+Z^2/c^2==1,Z][[2]]/.mclaurin

g={2Pi G ρ A[[1]]Xs,2Pi G ρ A[[3]]Zs}/.mclaurin

w=ContourPlot[{Evaluate@mel2d,x^2+z^2==R0^2},{x,-0.9,0.9},{z,-0.9,0.9},ContourStyle->{Red,{Black,Dashed}},Epilog->{
{Blue,Arrow[{{Xs,Zs},{Xs,Zs}+{Ω^2 Xs,0}/3-g/3}]},
Arrow[{{Xs,Zs},{Xs,Zs}+{Ω^2 Xs,0}/3}],Arrow[{{Xs,Zs},{Xs,Zs}-g/3}],Point[{0,0}]}]