(* 
Master di II livello in Scienza e Tecnologia Spaziale
     METODI NUMERICI PER L' ASTRONOMIA

Il software Mathematica per il calcolo simbolico, numerico e grafico.

Tor Vergata, 3 febbraio 2015

	Corrado Falcolini
	Università degli Studi Roma Tre 
*)

Problema 2: visualizzazione delle superfici a velocità zero definite tramite l'integrale di Jacobi.

(* Nel problema ristretto dei 3 corpi (una particella di massa trascurabile che si muova 
sotto l'azione gravitazionale di due masse (Subscript[m, 1]) ed (Subscript[m, 2]) in orbita circolare intorno al 
centro di massa con uguale velocità angolare) fissiamo le unità di misura riscrivendo 
 le masse (Subscript[μ, 1]) e (Subscript[μ, 2]) tali che (Subscript[μ, 1])+(Subscript[μ, 2])=1 ed il moto medio n=1 in un sistema di riferimento Oxyz
tale che (Subscript[m, 1]) ed (Subscript[m, 2]) siano fisse in posizione (-(Subscript[μ, 2]),0,0) e ((Subscript[μ, 1]),0,0). 
Tenendo conto allora del potenziale centrifugo e gravitazionale si può calcolare 
l'integrale di Jacobi: (Subscript[C, j]) = 2 U - (v^2) = (n^2)((x^2)+(y^2))+2((Subscript[μ, 1])/(Subscript[r, 1])+(Subscript[μ, 2])/(Subscript[r, 2]))−((OverDot[x])^2)-((OverDot[y])^2)−((OverDot[z])^2) 
che permette di definire le superfici a velocità zero che delimitano regioni dove
il moto è accessibile.
Vogliamo visualizzare tali superfici e loro sezioni piane.  *)

n=1;

mu2=1-mu1;

r1[x_,y_,z_]:=Sqrt[(x+mu2)^2+y^2+z^2];

r2[x_,y_,z_]:=Sqrt[(x-mu1)^2+y^2+z^2];

v0[x_,y_,z_]:=n^2*(x^2+y^2)+2*(mu1/r1[x,y,z]+mu2/r2[x,y,z])-cj;

Plot3D[(v0[x,y,0]+cj)/.mu1->0.8,{x,-2,2}, {y,-2,2}]

FindMinimum[v0[x,y,0]+cj /. mu1->0.8,{{x,2},{y,2}}]

FullForm[a==0]

ContourPlot[v0[x,y,0]==0/.{mu1->0.8,cj->3.4}//Evaluate,{x,-2,2},{y,-2,2}]

Table[k,{k,1,10}]

Table[k^2,{k,1,10}]

Table[k^2,{k,1,100,5}]

Table[ContourPlot[v0[x,y,0]/.{mu1->0.8,cj->2.84+k }//Evaluate,{x,-2,2},{y,-2,2}],{k,0,2,.1}]

Table[ContourPlot[v0[x,y,0]==0/.{mu1->0.8,cj->2.84+k }//Evaluate,{x,-2,2},{y,-2,2}],{k,0,2,.1}]

Clear[x]

Clear[y]

Clear[z]

Manipulate[ContourPlot[v0[x,y,0]==0/.{mu1->0.8,cj->2.84+k}//Evaluate,{x,-5,5},{y,-5,5},PlotPoints->80],{{k,1},0,5}]

ContourPlot3D[v0[x,y,z]==0/.{mu1->0.8,cj->3}//Evaluate,{x,-3,3}, {y,-3,3},{z,-3,3},ViewPoint->{0,0,2}]

ContourPlot3D[(v0[x,y,z]==0)/.{mu1->0.8,cj->3.7}//Evaluate,{x,-2,2}, {y,-2,2},{z,-2,2}]

Options[ContourPlot3D]