(* Content-type: application/mathematica *) (*** Wolfram Notebook File ***) (* http://www.wolfram.com/nb *) (* CreatedBy='Mathematica 7.0' *) (*CacheID: 234*) (* Internal cache information: NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 145, 7] NotebookDataLength[ 400926, 8753] NotebookOptionsPosition[ 394305, 8534] NotebookOutlinePosition[ 394830, 8552] CellTagsIndexPosition[ 394787, 8549] WindowFrame->Normal*) (* Beginning of Notebook Content *) Notebook[{ Cell[CellGroupData[{ Cell[TextData[StyleBox["You can use a single program for everything: \ Mathematica 7", FontSlant->"Italic"]], "Title", CellChangeTimes->{{3.4419825218125*^9, 3.441982528640625*^9}, { 3.441982575859375*^9, 3.441982597703125*^9}}], Cell["Roman Shcherbakov, 30 Jan 2009, group meeting", "Subtitle", CellChangeTimes->{{3.44198261159375*^9, 3.4419826255*^9}}], Cell[CellGroupData[{ Cell["Manipulate: N-body simulation of a merger (single processor)", "Section", CellChangeTimes->{{3.441985587921875*^9, 3.441985610453125*^9}}], Cell[BoxData[{ RowBox[{ RowBox[{"M1", "=", "5"}], ";"}], "\n", RowBox[{ RowBox[{ RowBox[{"MinRad", "=", "5"}], ";"}], "\n", RowBox[{"(*", " ", RowBox[{ "return", " ", "the", " ", "positions", " ", "and", " ", "velocities", " ", "of", " ", "the", " ", "particles"}], " ", "*)"}]}], "\n", RowBox[{ RowBox[{ RowBox[{"posandvel", "[", RowBox[{ RowBox[{"{", RowBox[{"a_", ",", "b_", ",", "c_"}], "}"}], ",", RowBox[{"{", RowBox[{"r_", ",", "dr_"}], "}"}], ",", "m_"}], "]"}], ":=", "\[IndentingNewLine]", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", RowBox[{ "xp", ",", "yp", ",", "pos", ",", "unitvel", ",", "vel", ",", "xandy"}], "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{"xp", ",", "yp"}], "}"}], "=", RowBox[{"norms", "[", RowBox[{"{", RowBox[{"a", ",", "b", ",", "c"}], "}"}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"xandy", "=", RowBox[{"Select", "[", RowBox[{ RowBox[{"Flatten", "[", RowBox[{ RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{"x", ",", "y"}], "}"}], ",", RowBox[{"{", RowBox[{"x", ",", RowBox[{"-", "r"}], ",", "r", ",", "dr"}], "}"}], ",", RowBox[{"{", RowBox[{"y", ",", RowBox[{"-", "r"}], ",", "r", ",", "dr"}], "}"}]}], "]"}], ",", "1"}], "]"}], ",", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"MinRad", "^", "2"}], "<", RowBox[{"(", RowBox[{"#", ".", "#"}], ")"}], "<", RowBox[{"r", "^", "2"}]}], ")"}], "&"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"pos", "=", RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"xp", " ", RowBox[{"#", "[", RowBox[{"[", "1", "]"}], "]"}]}], "+", RowBox[{"yp", " ", RowBox[{"#", "[", RowBox[{"[", "2", "]"}], "]"}]}]}], ")"}], "&"}], ",", "xandy"}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"unitvel", "=", RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"-", RowBox[{"#", "[", RowBox[{"[", "1", "]"}], "]"}]}], " ", "yp"}], "+", RowBox[{ RowBox[{"#", "[", RowBox[{"[", "2", "]"}], "]"}], " ", "xp"}]}], ")"}], "/", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"#", "[", RowBox[{"[", "1", "]"}], "]"}], "^", "2"}], "+", RowBox[{ RowBox[{"#", "[", RowBox[{"[", "2", "]"}], "]"}], "^", "2"}]}], "]"}]}], ")"}], "&"}], ",", "xandy"}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"vel", "=", RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{ RowBox[{"#", "[", RowBox[{"[", "1", "]"}], "]"}], "*", RowBox[{"Sqrt", "[", RowBox[{"m", "/", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"#", "[", RowBox[{"[", "2", "]"}], "]"}], ".", RowBox[{"#", "[", RowBox[{"[", "2", "]"}], "]"}]}], "]"}]}], "]"}]}], "&"}], ",", RowBox[{"Transpose", "[", RowBox[{"{", RowBox[{"unitvel", ",", "pos"}], "}"}], "]"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"{", RowBox[{"pos", ",", "vel"}], "}"}]}]}], "]"}]}], "\n", RowBox[{"(*", " ", RowBox[{ "calculate", " ", "the", " ", "normal", " ", "vectors", " ", "to", " ", "the", " ", "velocities"}], " ", "*)"}]}], "\n", RowBox[{ RowBox[{"norms", "[", RowBox[{"{", RowBox[{"x_", ",", RowBox[{"0", "|", "0."}], ",", RowBox[{"0", "|", "0."}]}], "}"}], "]"}], ":=", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{"0", ",", "1", ",", "0"}], "}"}], ",", RowBox[{"{", RowBox[{"0", ",", "0", ",", "1"}], "}"}]}], "}"}]}], "\n", RowBox[{ RowBox[{"norms", "[", "v_List", "]"}], ":=", "\[IndentingNewLine]", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", RowBox[{ "w", ",", "u", ",", "wy", ",", "wz", ",", "ux", ",", "uy", ",", "uz"}], "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"w", "=", RowBox[{"{", RowBox[{"0", ",", "wy", ",", "wz"}], "}"}]}], ";", "\[IndentingNewLine]", RowBox[{"u", "=", RowBox[{"{", RowBox[{"ux", ",", "uy", ",", "uz"}], "}"}]}], ";", "\[IndentingNewLine]", RowBox[{"w", "=", RowBox[{"w", "/.", RowBox[{"First", "[", RowBox[{"NSolve", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"w", ".", "v"}], "\[Equal]", "0"}], ",", RowBox[{ RowBox[{"w", ".", "w"}], "\[Equal]", "1"}]}], "}"}], ",", RowBox[{"{", RowBox[{"wy", ",", "wz"}], "}"}]}], "]"}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"u", "=", RowBox[{"u", "/.", RowBox[{"First", "[", RowBox[{"NSolve", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"u", ".", "v"}], "\[Equal]", "0"}], ",", RowBox[{ RowBox[{"u", ".", "w"}], "\[Equal]", "0"}], ",", RowBox[{ RowBox[{"u", ".", "u"}], "\[Equal]", "1"}]}], "}"}], ",", "u"}], "]"}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"{", RowBox[{"w", ",", "u"}], "}"}]}]}], "]"}]}], "\n", RowBox[{ RowBox[{"CenterOfMass", "[", RowBox[{ RowBox[{"{", RowBox[{"X1_", ",", "Y1_", ",", "Z1_"}], "}"}], ",", RowBox[{"{", RowBox[{"X2_", ",", "Y2_", ",", "Z2_"}], "}"}], ",", RowBox[{"{", RowBox[{"M1_", ",", "M2_"}], "}"}]}], "]"}], ":=", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"M1", " ", "X1"}], "+", RowBox[{"M2", " ", "X2"}]}], ")"}], "/", RowBox[{"(", RowBox[{"M1", "+", "M2"}], ")"}]}], ",", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"M1", " ", "Y1"}], "+", RowBox[{"M2", " ", "Y2"}]}], ")"}], "/", RowBox[{"(", RowBox[{"M1", "+", "M2"}], ")"}]}], ",", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"M1", " ", "Z1"}], "+", RowBox[{"M2", " ", "Z2"}]}], ")"}], "/", RowBox[{"(", RowBox[{"M1", "+", "M2"}], ")"}]}]}], "}"}]}], "\n", RowBox[{ RowBox[{"Initialize", "[", RowBox[{ "Dist_", ",", "MF_", ",", "XX1_", ",", "XX2_", ",", "pos_", ",", "pos2_", ",", "vel_", ",", "vel2_", ",", "V1_", ",", "V2_", ",", "ss_", ",", "tcol_", ",", "icol_", ",", "vp_", ",", "TIME_", ",", "ir_", ",", "tr_", ",", "sep_"}], "]"}], ":=", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", RowBox[{ "Rad1", ",", "Rad2", ",", "dr", ",", "intpl", ",", "Distlocal", ",", "MFlocal", ",", "M2local", ",", "poslocal", ",", "pos2local", ",", "vellocal", ",", "vel2local", ",", "XX1local", ",", "XX2local", ",", "V1local", ",", "V2local", ",", "sslocal", ",", "tcollocal", ",", "icollocal", ",", "vplocal", ",", "TIMElocal"}], "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{ "Distlocal", ",", "Rad1", ",", "Rad2", ",", "dr", ",", "MFlocal", ",", "XX2local", ",", "V2local", ",", "intpl", ",", "sslocal", ",", "tcollocal", ",", "icollocal", ",", "vplocal", ",", "TIMElocal"}], "}"}], "=", RowBox[{"{", RowBox[{ "100", ",", "tr", ",", "ir", ",", "sep", ",", "MF", ",", "XX2", ",", "V2", ",", "Automatic", ",", ".001", ",", "Green", ",", "Cyan", ",", RowBox[{"{", RowBox[{"1.3", ",", RowBox[{"-", "2.4"}], ",", "2."}], "}"}], ",", "0"}], "}"}]}], ";", "\[IndentingNewLine]", RowBox[{"intpl", "=", "V2local"}], ";", "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"If", "[", RowBox[{ RowBox[{"intpl", "===", "Automatic"}], ",", RowBox[{"intpl", "=", "V2local"}]}], "]"}], ";"}], "*)"}], "\[IndentingNewLine]", RowBox[{"XX1local", "=", RowBox[{"V1local", "=", RowBox[{"{", RowBox[{"0.", ",", "0.", ",", "0."}], "}"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"M2local", "=", RowBox[{"MFlocal", " ", "M1"}]}], ";", "\[IndentingNewLine]", RowBox[{ RowBox[{"{", RowBox[{"poslocal", ",", "vellocal"}], "}"}], "=", RowBox[{"posandvel", "[", RowBox[{ RowBox[{"{", RowBox[{"0", ",", "0", ",", "1"}], "}"}], ",", RowBox[{"{", RowBox[{"Rad1", ",", "dr"}], "}"}], ",", "M1"}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{ RowBox[{"{", RowBox[{"pos2local", ",", "vel2local"}], "}"}], "=", RowBox[{"posandvel", "[", RowBox[{"intpl", ",", RowBox[{"{", RowBox[{"Rad2", ",", "dr"}], "}"}], ",", "M2local"}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"pos2local", "=", RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"#", "+", "XX2"}], ")"}], "&"}], ",", "pos2local"}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"vel2local", "=", RowBox[{"Map", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"#", "+", "V2"}], ")"}], "&"}], ",", "vel2local"}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"{", RowBox[{ "Distlocal", ",", "MFlocal", ",", "V1local", ",", "V2local", ",", "sslocal", ",", "tcollocal", ",", "icollocal", ",", "vplocal", ",", "TIMElocal", ",", RowBox[{"Sequence", "@@", RowBox[{"Realign", "[", RowBox[{ "XX1local", ",", "XX2local", ",", "M2local", ",", "poslocal", ",", "pos2local"}], "]"}]}], ",", "vellocal", ",", "vel2local"}], "}"}]}]}], "\[IndentingNewLine]", "]"}]}], "\n", RowBox[{ RowBox[{ RowBox[{"Draw", "[", RowBox[{ "XX1_", ",", "XX2_", ",", "MF_", ",", "pos_", ",", "pos2_", ",", "Dist_", ",", "TIME_", ",", "vp_", ",", "ss_", ",", "tcol_", ",", "icol_"}], "]"}], ":=", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", RowBox[{"cm", "=", RowBox[{"CenterOfMass", "[", RowBox[{"XX1", ",", "XX2", ",", RowBox[{"{", RowBox[{"M1", ",", RowBox[{"M1", "*", "MF"}]}], "}"}]}], "]"}]}], "}"}], ",", RowBox[{"Graphics3D", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{"tcol", ",", RowBox[{"PointSize", "[", "ss", "]"}], ",", RowBox[{"Point", "[", "pos", "]"}], ",", RowBox[{"PointSize", "[", ".025", "]"}], ",", RowBox[{"Point", "[", "XX1", "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"icol", ",", RowBox[{"PointSize", "[", "ss", "]"}], ",", RowBox[{"Point", "[", "pos2", "]"}], ",", RowBox[{"PointSize", "[", ".025", "]"}], ",", RowBox[{"Point", "[", "XX2", "]"}]}], "}"}]}], "}"}], ",", RowBox[{"BaseStyle", "\[Rule]", RowBox[{"{", RowBox[{"GrayLevel", "[", "1", "]"}], "}"}]}], ",", RowBox[{"AxesStyle", "\[Rule]", RowBox[{"GrayLevel", "[", ".5", "]"}]}], ",", RowBox[{"BoxStyle", "\[Rule]", RowBox[{"GrayLevel", "[", ".5", "]"}]}], ",", RowBox[{"PlotRange", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"-", "Dist"}], "+", RowBox[{"cm", "[", RowBox[{"[", "1", "]"}], "]"}]}], ",", RowBox[{"Dist", "+", RowBox[{"cm", "[", RowBox[{"[", "1", "]"}], "]"}]}]}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"-", "Dist"}], "+", RowBox[{"cm", "[", RowBox[{"[", "2", "]"}], "]"}]}], ",", RowBox[{"Dist", "+", RowBox[{"cm", "[", RowBox[{"[", "2", "]"}], "]"}]}]}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{"-", "Dist"}], "+", RowBox[{"cm", "[", RowBox[{"[", "3", "]"}], "]"}]}], ",", RowBox[{"Dist", "+", RowBox[{"cm", "[", RowBox[{"[", "3", "]"}], "]"}]}]}], "}"}]}], "}"}]}], ",", RowBox[{"Axes", "\[Rule]", "True"}], ",", RowBox[{"ViewPoint", "\[Rule]", "vp"}], ",", RowBox[{"Background", "\[Rule]", RowBox[{"GrayLevel", "[", "0", "]"}]}], ",", RowBox[{"PlotLabel", "\[Rule]", RowBox[{"\"\