SSS

jednotky a rozmery

dc=UnitConvert[Quantity[1.,"Meter"],Quantity[1.,"AstronomicalUnit"]]
vc={0.000210805}(*vc=UnitConvert[Quantity[1.,"Meters/Second"],Quantity[1,"AstronomicalUnit/Year"]]*)
yrs=UnitConvert[Quantity[1.,"Year"]]
mc=UnitConvert[Quantity[1,"Kilogram"],Quantity[1,"SolarMass"]]
gc0=UnitConvert[Quantity[1,"GravitationalConstant"]]
gc=gc0[[1]]((dc[[1]])^3)((yrs[[1]])^2)/mc[[1]]
G= gc;(*grav.konst*)

počateční pozice a rychlosti jednotlivych planet

(* ? DATUM ? *)

sun={1.9891*^30×mc[[1]],dc[[1]]×-7.0299*^8,dc[[1]]×-7.5415*^8,dc[[1]]×2.38988*^7,vc[[1]]×14.1931,vc[[1]]×-6.9255,vc[[1]]×-0.31676,Yellow,""};
mercury={3.302*^23×mc[[1]],dc[[1]]×2.60517*^10,dc[[1]]×-6.1102*^10,dc[[1]]×-7.3616*^9,vc[[1]]×34796.,vc[[1]]×22185.2,vc[[1]]×-1379.78,Orange,"☿"};
venus={4.8685*^24×mc[[1]],dc[[1]]×7.2129*^10,dc[[1]]×7.9106*^10,dc[[1]]×-3.0885*^9,vc[[1]]×-25968.7,vc[[1]]×23441.6,vc[[1]]×1819.92,Green,"♀"};
earth={5.9736*^24×mc[[1]],dc[[1]]×-2.91204*^10,dc[[1]]×1.43576*^11,dc[[1]]×2.39614*^7,vc[[1]]×-29699.8,vc[[1]]×-5883.3,vc[[1]]×0.050215,Blue,"♁"};
mars={6.4185*^23×mc[[1]],dc[[1]]×-2.47064*^11,dc[[1]]×-1.03161*^10,dc[[1]]×5.8788*^9,vc[[1]]×1862.73,vc[[1]]×-22150.6,vc[[1]]×-509.6,Red,"♂"};
jupiter={1.8986*^27×mc[[1]],dc[[1]]×2.67553*^11,dc[[1]]×7.0482*^11,dc[[1]]×-8.911*^9,vc[[1]]×-12376.3,vc[[1]]×5259.2,vc[[1]]×-509.6,Brown,"♃"};
saturn={5.9846*^26×mc[[1]],dc[[1]]×6.999*^11,dc[[1]]×1.16781*^12,dc[[1]]×-4.817*^10,vc[[1]]×-8792.6,vc[[1]]×4944.9,vc[[1]]×263.754,Cyan,"♄"};
uranus={8.6832*^25×mc[[1]],dc[[1]]×2.2993*^12,dc[[1]]×-1.90411*^12,dc[[1]]×-3.6864*^10,vc[[1]]×4293.6,vc[[1]]×4928.1,vc[[1]]×-37.32,Pink,"♅"};
neptune={1.0243*^26×mc[[1]],dc[[1]]×2.65363*^12,dc[[1]]×-3.6396*^12,dc[[1]]×1.37957*^10,vc[[1]]×4356.6,vc[[1]]×3233.3,vc[[1]]×-166.986,Magenta,"♆"};
pluto={1.27*^22×mc[[1]],dc[[1]]×-1.31126*^12,dc[[1]]×-4.2646*^12,dc[[1]]×8.3563*^11,vc[[1]]×5316.6,vc[[1]]×-2484.6,vc[[1]]×-1271.99,Black,"♇"};

>>>>> definice fce slunecni soustava <<<<<

Clear[fceSS,RG,solarsys];
fceSS[xroky_,xmetoda_]:=(

(*========================= URCENI SOUCASTI SOL.SYS. A SESTAVENI HAMILTONIANU ======================================*)
npart=Length[solarsys ];

mass =solarsys[[All,1]];(* pole s hmotnosti jednotlivych tels sol.sys. *)
gg=Table[Table[-G×mass[[ii]] mass[[jj]],{ii,1,npart}],{jj,1,npart}];(* matice hmotnosti (pro interakci) *)

(* sestaveni Hamiltonianu - interakace "kazdy s kazdym"*)
H :=(∑_(k=1)^npart)((((((px[k][t]))^2)+(((py[k][t]))^2)+(((pz[k][t]))^2))/(2*mass[[k]]))+(∑_(ii=1)^(k-1))gg[[ii,k]]/(√((((qx[k][t]-qx[ii][t]))^2)+(((qy[k][t]-qy[ii][t]))^2)+(((qz[k][t]-qz[ii][t]))^2))));

(* jednotlivé proměné souradnice a hybnosti vsech teles sol.sys. *)
varqx=Table[qx[jj][t],{jj,1,npart}];varqy=Table[qy[jj][t],{jj,1,npart}];varqz=Table[qz[jj][t],{jj,1,npart}];
varpx=Table[px[jj][t],{jj,1,npart}];varpy=Table[py[jj][t],{jj,1,npart}];varpz=Table[pz[jj][t],{jj,1,npart}];
vars=Join[varqx,varqy,varqz,varpx,varpy,varpz];(* souradnice a hybnosti dohromady *)

posvar=Flatten[Join[varqx,varqy,varqz]];(* pouze souradnice: nutne pro "SymplecticPartitionedRungeKutta" metodu*)
(*========================= URCENI SOUCASTI SOL.SYS. A SESTAVENI HAMILTONIANU ======================================*)

(*========================= SESTAVENI POHYBOVYCH ROVNIC ======================================*)
(* pohybové rovnice *)
eqqx=Table[qx[jj]'[t]==D[H,px[jj][t]],{jj,1,npart}];
eqqy=Table[qy[jj]'[t]==D[H,py[jj][t]],{jj,1,npart}];
eqqz=Table[qz[jj]'[t]==D[H,pz[jj][t]],{jj,1,npart}];

eqpx=Table[px[jj]'[t]==-D[H,qx[jj][t]],{jj,1,npart}];
eqpy=Table[py[jj]'[t]==-D[H,qy[jj][t]],{jj,1,npart}];
eqpz=Table[pz[jj]'[t]==-D[H,qz[jj][t]],{jj,1,npart}];

(* pohybové rovnice všechny dohromady*)
ODEs=Join[eqqx,eqqy,eqqz,eqpx,eqpy,eqpz];
(*========================= ESESTAVENI POHYBOVYCH ROVNIC ======================================*)

(*========================= POCATECNI PODMINKY ======================================*)
(* počáteční pozice všech těles *)
initqx=solarsys[[All,2]];icqx=Table[qx[jj][t]==initqx[[jj]]/.t->0,{jj,1,npart}];
initqy=solarsys[[All,3]];icqy=Table[qy[jj][t]==initqy[[jj]]/.t->0,{jj,1,npart}];
initqz=solarsys[[All,4]];icqz=Table[qz[jj][t]==initqz[[jj]]/.t->0,{jj,1,npart}];

(* počáteční hybnosti všech těles *)
initpx=Table[solarsys[[i,1]]*solarsys[[i,5]],{i,1,npart}];icpx=Table[px[jj][t]==initpx[[jj]]/.t->0,{jj,1,npart}];

initpy=Table[solarsys[[i,1]]*solarsys[[i,6]],{i,1,npart}];
icpy=Table[py[jj][t]==initpy[[jj]]/.t->0,{jj,1,npart}];

initpz=Table[solarsys[[i,1]]*solarsys[[i,7]],{i,1,npart}];
icpz=Table[pz[jj][t]==initpz[[jj]]/.t->0,{jj,1,npart}];

(* počáteční pozice a hybnosti všech těles dohromady *)
ics=Join[icqx,icqy,icqz,icpx,icpy,icpz];
(*========================= POCATECNI PODMINKY ======================================*)

(*========================= RESI SE NUMERICKY DIF. ROVNICE ======================================*)
tA[];
eqns = Join[ODEs,ics];
tstart=0.00;tint=1.00xroky;

If[xmetoda==1,
(* 1>>>> pomala ale velice presna metoda: *) startstep=1/12;(*integracni krok jeden den*)
Metoda={"SymplecticPartitionedRungeKutta","DifferenceOrder"->10,"PositionVariables"->posvar};,
(* 0>>>> rychle ale nepresna presna metoda: *)startstep=1/12;(*integracni krok jeden mesic*)Metoda=Automatic;
];
sol=NDSolve[eqns,vars,{t,tstart,tint},Method->Metoda,StartingStepSize->startstep,InterpolationOrder->All,MaxSteps->(10^7)];
draha=sol[[1]];
tB[];
(*========================= RESI SE NUMERICKY DIF. ROVNICE ======================================*)

);

solarsys = {sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,pluto};
fceSS[300,1]

nn=5;cas=2;epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.sol[[1]]/.t->0]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.sol[[1]],{n,1,nn}];
orbits=Table[ParametricPlot[orbitfunctions[[n]],{t,0,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
pic1=Show[Table[orbits[[jj]],{jj,1,nn}],PlotRange->{{-1.7,1.5},{-1.6,1.6}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},RotateLabel->True,Epilog->epilog]

nn=npart;cas=250;epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.sol[[1]]/.t->50]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.sol[[1]],{n,1,nn}];
orbits=Table[ParametricPlot[orbitfunctions[[n]],{t,0,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
pic2=Show[Table[orbits[[jj]],{jj,1,nn}],PlotRange->{{-32,45},{-32,45}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=250;epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qz[n][t]}/.sol[[1]]/.t->50]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t],qz[n][t]}/.sol[[1]],{n,1,nn}];
orbits=Table[ParametricPlot[orbitfunctions[[n]],{t,0,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
pic3=Show[Table[orbits[[jj]],{jj,1,nn}],PlotRange->{{-40,45},{-42.5,42.5}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["z [AU]",13]},ImageSize->300,Frame->True,AspectRatio->1,RotateLabel->True,BaseStyle->{15,Background->White},Epilog->epilog]

pic=GraphicsRow[{pic1,pic2},Spacings->-7]
Export["sss1.png",pic];

pic=GraphicsRow[{pic2,pic3},Spacings->-7]
Export["sss2.png",pic];

nn=5;cas=2;
epilog=Table[Text[solarsys[[n,9]],N[{t,(Sqrt[((qx[n][t])^2)+((qy[n][t])^2)+((qz[n][t])^2)])}/.sol[[1]]/.t->0.2]],{n,2,nn}];
orbitfunctions=Table[(Sqrt[((qx[n][t])^2)+((qy[n][t])^2)+((qz[n][t])^2)])/.sol[[1]],{n,1,nn}];
orbits=Table[Plot[orbitfunctions[[n]],{t,0,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
tra1=Show[Table[orbits[[jj]],{jj,1,nn}],PlotRange->All,FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=200;
epilog=Table[Text[solarsys[[n,9]],N[{t,Log[(Sqrt[((qx[n][t])^2)+((qy[n][t])^2)+((qz[n][t])^2)])]}/.sol[[1]]/.t->50]],{n,6,nn}];
orbitfunctions=Table[(Sqrt[((qx[n][t])^2)+((qy[n][t])^2)+((qz[n][t])^2)])/.sol[[1]],{n,1,nn}];
orbits=Table[LogPlot[orbitfunctions[[n]],{t,0,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
tra2=Show[Table[orbits[[jj]],{jj,1,nn}],PlotRange->{1.5,4},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

pic=GraphicsRow[{tra1,tra2},Spacings->-7]
Export["sss3.png",pic];

“milion let”

solarsys = {sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,pluto};
end=10000;fceSS[10000,1]

end=10000;

n=2;cas=2;
ParametricPlot[{qx[n][t],qy[n][t]}/.draha,{t,end-cas,end}]

nn=5;cas=2;epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]]/.t->0]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]],{n,1,nn}];
orbits=Table[ParametricPlot[orbitfunctions[[n]],{t,end-cas,end},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
pic1=Show[Table[orbits[[jj]],{jj,1,nn}],PlotRange->{{-1.7,1.5},{-1.6,1.6}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},RotateLabel->True,Epilog->epilog]

alien star 0

trozdel=147;

tint=300;impact=10000;starmass=1;vel=20;incl=30;Δt=1;Δtg=0.5;
miscdata={tint,impact,starmass,vel,incl,Δt,Δtg};
vc1=vc (10^3);
vxS=0;vyS=-vel×vc1[[1]]×Cos[Degree×incl];vzS=-vel×vc1[[1]]×Sin[Degree×incl];
tanϕ=(impact/(-vyS×tint/2));ϕ=ArcTan[tanϕ];radvel=vel×√((((Cos[incl×Degree]×Cos[ϕ]))^2)+((Sin[incl×Degree])^2));

distance=(Sqrt[(impact^2)+(((vyS×tint/2))^2)+(((vyS×tint/2))^2)])

star={starmass,-impact,-vyS×tint/2,-vzS×tint/2,vxS,vyS,vzS,Yellow,"X"};

Clear[sol,draha]
solarsys = {sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,pluto,star};
fceSS[300,1]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.sol[[1]]/.t->(cas-20)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p1=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-40,50},{-40,50}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]]/.t->(cas-40)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p2=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-40,50},{-40,50}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{t,Log[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])]}/.sol[[1]]/.t->50]],{n,6,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.sol[[1]],{n,1,nn}];
orbitsA=Table[LogPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[LogPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p3=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{1.5,9.1},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

trozdel=149;
nn=5;t1=140;t2=trozdel;t3=trozdel;t4=trozdel+10;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.sol[[1]]/.t->trozdel+3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p4=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-2,1.5},{-1.5,2}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},Axes->False,ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;t1=0;t2=3;t3=tint-t2;t4=tint;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]]/.t->t3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p5=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-1.9,1.5},{-1.5,1.9}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;cas=tint;
epilog=Table[Text[solarsys[[n,9]],N[{t,(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])}/.sol[[1]]/.t->20]],{n,2,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.sol[[1]],{n,1,nn}];
orbitsA=Table[Plot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,nn,1,-1}];
orbitsB=Table[Plot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,nn,1,-1}];
p6=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{0,2},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

pic=GraphicsGrid[{{p1,p2,p3},{p4,p5,p6}},Spacings->-7]

Export["ssa0.png",pic];

alien star 1

tint=300;impact=1000;starmass=10;vel=20;incl=30;Δt=1;Δtg=0.5;
miscdata={tint,impact,starmass,vel,incl,Δt,Δtg};
vc1=vc (10^3);
vxS=0;vyS=-vel×vc1[[1]]×Cos[Degree×incl];vzS=-vel×vc1[[1]]×Sin[Degree×incl];
tanϕ=(impact/(-vyS×tint/2));ϕ=ArcTan[tanϕ];radvel=vel×√((((Cos[incl×Degree]×Cos[ϕ]))^2)+((Sin[incl×Degree])^2));

distance=(Sqrt[(impact^2)+(((vyS×tint/2))^2)+(((vyS×tint/2))^2)])

star={starmass,-impact,-vyS×tint/2,-vzS×tint/2,vxS,vyS,vzS,Yellow,"X"};

solarsys = {sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,pluto,star};
fceSS[300,1]

fceX[t_]:=(Sqrt[(((qx[1][t]-qx[11][t]))^2)+(((qy[1][t]-qy[11][t]))^2)+(((qz[1][t]-qz[11][t]))^2)])/.sol[[1]];
rozdelmin=Quiet[FindMinimum[fceX[t],{t,150}]]
trozdel=t/.rozdelmin[[2]];
Plot[fceX[t],{t,0,300},Epilog->{Line[{{trozdel,0},{trozdel,2000}}]}]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.sol[[1]]/.t->(cas-20)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p1=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-45,45},{-45,45}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]]/.t->(cas-40)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p2=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-40,50},{-40,50}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{t,Log[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])]}/.sol[[1]]/.t->50]],{n,6,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.sol[[1]],{n,1,nn}];
orbitsA=Table[LogPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[LogPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p3=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{1.5,7.1},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;t1=140;t2=trozdel;t3=trozdel;t4=trozdel+10;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.sol[[1]]/.t->trozdel+3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p4=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-5.7,-1.5},{-0.7,3.5}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},Axes->False,ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;t1=0;t2=3;t3=tint-t2;t4=tint;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]]/.t->t3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p5=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-1.9,1.5},{-1.5,1.9}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;cas=tint;
epilog=Table[Text[solarsys[[n,9]],N[{t,(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])}/.sol[[1]]/.t->20]],{n,2,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.sol[[1]],{n,1,nn}];
orbitsA=Table[Plot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,nn,1,-1}];
orbitsB=Table[Plot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,nn,1,-1}];
p6=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{0,2},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

pic=GraphicsGrid[{{p1,p2,p3},{p4,p5,p6}},Spacings->-7]

Export["ssa1.png",pic];

alien star 2

Clear[draha]
tint=300;impact=100;starmass=1;vel=30;incl=30;Δt=1;Δtg=0.5;
miscdata={tint,impact,starmass,vel,incl,Δt,Δtg};
vc1=vc (10^3);
vxS=0;vyS=-vel×vc1[[1]]×Cos[Degree×incl];vzS=-vel×vc1[[1]]×Sin[Degree×incl];
tanϕ=(impact/(-vyS×tint/2));ϕ=ArcTan[tanϕ];radvel=vel×√((((Cos[incl×Degree]×Cos[ϕ]))^2)+((Sin[incl×Degree])^2));

distance=(Sqrt[(impact^2)+(((vyS×tint/2))^2)+(((vyS×tint/2))^2)])

star={starmass,-impact,-vyS×tint/2,-vzS×tint/2,vxS,vyS,vzS,Yellow,"x"};

solarsys = {sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,pluto,star};
fceSS[300,1]

fceX[t_]:=(Sqrt[(((qx[1][t]-qx[11][t]))^2)+(((qy[1][t]-qy[11][t]))^2)+(((qz[1][t]-qz[11][t]))^2)])/.draha;
rozdelmin=Quiet[FindMinimum[fceX[t],{t,150}]]
trozdel=t/.rozdelmin[[2]];
Plot[fceX[t],{t,0,300},Epilog->{Line[{{trozdel,0},{trozdel,2000}}]}]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.draha/.t->(cas-20)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.draha,{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p1=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-110,50},{-110,50}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.draha/.t->(cas-40)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.draha,{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p2=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-40,50},{-40,50}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{t,Log[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])]}/.draha/.t->100]],{n,6,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.draha,{n,1,nn}];
orbitsA=Table[LogPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[LogPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p3=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{1.5,6.1},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;t1=140;t2=trozdel;t3=trozdel;t4=trozdel+10;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.draha/.t->trozdel+3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.draha,{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p4=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-3.5,0.7},{-0.5,3.7}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},Axes->False,ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;t1=0;t2=3;t3=tint-t2;t4=tint;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.draha/.t->t3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.draha,{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p5=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-1.9,1.5},{-1.5,1.9}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;cas=tint;
epilog=Table[Text[solarsys[[n,9]],N[{t,(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])}/.draha/.t->20]],{n,2,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.draha,{n,1,nn}];
orbitsA=Table[Plot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,nn,1,-1}];
orbitsB=Table[Plot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,nn,1,-1}];
p6=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{0,2},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

pic=GraphicsGrid[{{p1,p2,p3},{p4,p5,p6}},Spacings->-7]

Export["ssa2.png",pic];

alien star 3

Clear[draha]
tint=300;impact=3;starmass=0.1;vel=20;incl=30;Δt=1;Δtg=0.5;
miscdata={tint,impact,starmass,vel,incl,Δt,Δtg};
vc1=vc (10^3);
vxS=0;vyS=-vel×vc1[[1]]×Cos[Degree×incl];vzS=-vel×vc1[[1]]×Sin[Degree×incl];
tanϕ=(impact/(-vyS×tint/2));ϕ=ArcTan[tanϕ];radvel=vel×√((((Cos[incl×Degree]×Cos[ϕ]))^2)+((Sin[incl×Degree])^2));

distance=(Sqrt[(impact^2)+(((vyS×tint/2))^2)+(((vyS×tint/2))^2)])

star={starmass,-impact,-vyS×tint/2,-vzS×tint/2,vxS,vyS,vzS,Yellow,"x"};

solarsys = {sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,pluto,star};
fceSS[300,1]

Clear[fceX];
fceX[t_]:=(Sqrt[(((qx[1][t]-qx[11][t]))^2)+(((qy[1][t]-qy[11][t]))^2)+(((qz[1][t]-qz[11][t]))^2)])/.draha;
rozdelmin=Quiet[FindMinimum[{fceX[t],140≤t≤160},{t,150}]]
trozdel=t/.rozdelmin[[2]];
Plot[fceX[t],{t,140,150},Epilog->{Line[{{trozdel,0},{trozdel,2000}}]},PlotRange->{0,10}]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.sol[[1]]/.t->(cas-20)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p1=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-70,50},{-50,70}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]]/.t->(cas-40)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p2=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-40,70},{-40,70}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.draha,{n,1,nn}];
LogPlot[orbitfunctions[[11]],{t,140,150},PlotStyle->Lighter[solarsys[[11,8]],0.5],PlotRange->{1.5,4}]

p3=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}],LogPlot[orbitfunctions[[11]],{t,140,155},PlotStyle->Lighter[solarsys[[11,8]],0.5]]}],PlotRange->{1.5,4},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{t,Log[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])]}/.draha/.t->50]],{n,6,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.draha,{n,1,nn}];
orbitsA=Table[LogPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[LogPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p3=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}],LogPlot[orbitfunctions[[11]],{t,0,300},PlotStyle->Lighter[solarsys[[11,8]],0.5]]}],PlotRange->{1.5,4},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

ParametricPlot[{qx[11][t],qy[11][t]}/.draha,{t,t1,t4},PlotStyle->solarsys[[11,8]],PlotRange->{{-3.5,1.5},{-2.5,2.5}}]

nn=5;t1=140;t2=trozdel;t3=trozdel;t4=trozdel+10;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.sol[[1]]/.t->trozdel+3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p4=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}],ParametricPlot[{qx[11][t],qy[11][t]}/.draha,{t,t1,t4},PlotStyle->solarsys[[11,8]]]}],PlotRange->{{-3.5,1.5},{-2.5,2.5}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},Axes->False,ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;t1=0;t2=3;t3=tint-t2;t4=tint;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]]/.t->t3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.sol[[1]],{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p5=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-1.9,1.5},{-1.5,1.9}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;cas=tint;
epilog=Table[Text[solarsys[[n,9]],N[{t,(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])}/.sol[[1]]/.t->20]],{n,2,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.sol[[1]],{n,1,nn}];
orbitsA=Table[Plot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,nn,1,-1}];
orbitsB=Table[Plot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,nn,1,-1}];
p6=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{0,2},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

pic=GraphicsGrid[{{p1,p2,p3},{p4,p5,p6}},Spacings->-7]

Export["ssa3.png",pic];

alien star 4

Clear[draha]
tint=300;impact=10;starmass=5;vel=30;incl=30;Δt=1;Δtg=0.5;
miscdata={tint,impact,starmass,vel,incl,Δt,Δtg};
vc1=vc (10^3);
vxS=0;vyS=-vel×vc1[[1]]×Cos[Degree×incl];vzS=-vel×vc1[[1]]×Sin[Degree×incl];
tanϕ=(impact/(-vyS×tint/2));ϕ=ArcTan[tanϕ];radvel=vel×√((((Cos[incl×Degree]×Cos[ϕ]))^2)+((Sin[incl×Degree])^2));

distance=(Sqrt[(impact^2)+(((vyS×tint/2))^2)+(((vyS×tint/2))^2)])

star={starmass,-impact,-vyS×tint/2,-vzS×tint/2,vxS,vyS,vzS,Yellow,"x"};

solarsys = {sun,mercury,venus,earth,mars,jupiter,saturn,uranus,neptune,pluto,star};
fceSS[300,1]

fceX[t_]:=(Sqrt[(((qx[1][t]-qx[11][t]))^2)+(((qy[1][t]-qy[11][t]))^2)+(((qz[1][t]-qz[11][t]))^2)])/.draha;
rozdelmin=Quiet[FindMinimum[fceX[t],{t,150}]]
trozdel=t/.rozdelmin[[2]];
Plot[fceX[t],{t,0,300},Epilog->{Line[{{trozdel,0},{trozdel,2000}}]}]

trozdel=140

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.draha/.t->(trozdel+12)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.draha,{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p1=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-110,50},{-110,50}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.draha/.t->(trozdel-5)]],{n,6,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.draha,{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p2=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-40,50},{-40,50}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,Axes->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=npart;cas=300;
epilog=Table[Text[solarsys[[n,9]],N[{t,Log[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])]}/.draha/.t->100]],{n,6,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.draha,{n,1,nn}];
orbitsA=Table[LogPlot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[LogPlot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p3=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{1.5,6.1},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

trozdel=146;

nn=5;t1=trozdel-30;t2=trozdel;t3=trozdel;t4=trozdel+10;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t],qy[n][t]}/.draha/.t->trozdel+1.5]],{n,1,npart}];
orbitfunctions=Table[{qx[n][t],qy[n][t]}/.draha,{n,1,npart}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,npart}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,npart}];
p4=Show[Join[{Table[orbitsA[[jj]],{jj,1,npart}],Table[orbitsB[[jj]],{jj,1,npart}]}],PlotRange->{{-17,5},{-5,17}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},Axes->False,ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

trozdel=146;
nn=5;t1=trozdel-3;t2=trozdel;t3=trozdel;t4=trozdel+3;
epilog=Table[Text[solarsys[[n,9]],N[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.draha/.t->t3]],{n,1,nn}];
orbitfunctions=Table[{qx[n][t]-qx[1][t],qy[n][t]-qy[1][t]}/.draha,{n,1,nn}];
orbitsA=Table[ParametricPlot[orbitfunctions[[n]],{t,t1,t2},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,1,nn}];
orbitsB=Table[ParametricPlot[orbitfunctions[[n]],{t,t3,t4},PlotStyle->solarsys[[n,8]]],{n,1,nn}];
p5=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{{-1.9,1.5},{-1.5,1.9}},FrameLabel->{Null,Null,Style["x [AU]",13],Style["y [AU]",13]},ImageSize->300,Frame->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

nn=5;trozdel=146;cas=tint;
epilog=Table[Text[solarsys[[n,9]],N[{t,(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])}/.draha/.t->20]],{n,2,nn}];
orbitfunctions=Table[(Sqrt[(((qx[1][t]-qx[n][t]))^2)+(((qy[1][t]-qy[n][t]))^2)+(((qz[1][t]-qz[n][t]))^2)])/.draha,{n,1,nn}];
orbitsA=Table[Plot[orbitfunctions[[n]],{t,0,trozdel},PlotStyle->Lighter[solarsys[[n,8]],0.5]],{n,nn,1,-1}];
orbitsB=Table[Plot[orbitfunctions[[n]],{t,trozdel,cas},PlotStyle->solarsys[[n,8]]],{n,nn,1,-1}];
p6=Show[Join[{Table[orbitsA[[jj]],{jj,1,nn}],Table[orbitsB[[jj]],{jj,1,nn}]}],PlotRange->{0,2},FrameLabel->{Null,Null,Style["t [rok]",13],Style["vzdalenost  od  Slunce  [AU]",13]},ImageSize->300,Frame->True,Axes->False,RotateLabel->True,AspectRatio->1,BaseStyle->{15,Background->White},Epilog->epilog]

pic=GraphicsGrid[{{p1,p2,p3},{p4,p5,p6}},Spacings->-7]

Export["ssa4.png",pic];