Motion in a Central Field

Power-law potentials with screening

(Go down and play with controls; Requires Mathematica 6 )

U[r_,n_,α_,b_]:=-α Sign[n]((((b^2)+(r^2)))^(n/2)); (* Screened central potential: Attractive for α<0, repulsive for α>0 *)

CentralFieldOrbit[n_,α_,b_,x0y0_,vx0y0_,tMax_]:= (* This subroutine solves Newtonean equations of motion and plots trajectories, depending on input parameters *)
Quiet[Module[{nds,x,y,localx0y0},
localx0y0=If[x0y0==={0,0},{.0001,.0001},x0y0]; (* Avoid singularity at r=0, if any *)nds =NDSolve[Join[Thread[{x''[t],y''[t]}==-{D[U[Sqrt[((x[t])^2) + ((y[t])^2)],n,α,b], x[t]],D[U[Sqrt[((x[t])^2) + ((y[t])^2)],n,α,b],y[t]]}],
{Thread[{x[0],y[0]}==localx0y0],Thread[{x'[0],y'[0]}==vx0y0]}],{x,y},{t,0,tMax},PrecisionGoal->6];
ParametricPlot[{x[t],y[t]}/.nds[[1]],{t,0,nds[[1,1,2,1,1,2]]}, (* Plot until the actual stop of NDSolve at a singularity rather than until tMax *)
Frame->True,Axes -> False, Prolog -> {GrayLevel[.5], PointSize[0.02], Point[{0,0}]},ImageSize->{700,700},PlotRange->{{-2,2},{-2,2}},
FrameTicks -> None,
PlotPoints -> Ceiling[8 tMax+4],MaxRecursion->ControlActive[2, Automatic],
PlotStyle->{Thickness[.005]}]
]]

Manipulate[
CentralFieldOrbit[n,α,b,x0y0,vx0y0,tMax],
{{tMax,20,"maximal time"}, 0.1, 100},
{{n,-1,"force exponent n"},-4,4, Appearance -> "Labeled"},
{{α, -1,"force constant α"},-10,10, Appearance -> "Labeled"},
{{b, 0,"screening b"},0,10, Appearance -> "Labeled"},
{{x0y0,{1, 0.84}, "initial\nposition"},{-2,-2},{2,2},ImageSize->Medium},
{{vx0y0,{-0.2,-0.66},  "initial\nvelocity"},{-2,-2},{2,2},ImageSize->Medium},
ControlPlacement->{Top,Top,Top,Top,Left,Left},SaveDefinitions->True,AutorunSequencing->{4,5},TrackedSymbols:>{tMax,n,α,b,x0y0,vx0y0}]