(* Content-type: application/mathematica *) (*** Wolfram Notebook File ***) (* http://www.wolfram.com/nb *) (* CreatedBy='Mathematica 6.0' *) (*CacheID: 234*) (* Internal cache information: NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 145, 7] NotebookDataLength[ 51862, 1209] NotebookOptionsPosition[ 51485, 1193] NotebookOutlinePosition[ 51826, 1208] CellTagsIndexPosition[ 51783, 1205] WindowFrame->Normal ContainsDynamic->False*) (* Beginning of Notebook Content *) Notebook[{ Cell[BoxData[ RowBox[{"(*", " ", RowBox[{ RowBox[{ RowBox[{ "This", " ", "notebook", " ", "implements", " ", "a", " ", "simulation", " ", "of", " ", "the", " ", "time"}], "-", RowBox[{"dependent", " ", "behavior", " ", "of", " ", "a", " ", "star"}], "-", RowBox[{ "forming", " ", "disk", " ", "with", " ", "a", " ", "flat", " ", "rotation", " ", "curve"}]}], ",", " ", RowBox[{ RowBox[{ RowBox[{ RowBox[{ "following", " ", "the", " ", "procedure", " ", "of", " ", "Krumholz"}], " ", "&"}], " ", "Burkert", " ", "2010.", " ", "User"}], "-", RowBox[{ "settable", " ", "parameters", " ", "are", " ", "in", " ", "the", " ", "first", " ", "three", " ", RowBox[{"cells", "."}]}]}]}], " ", "*)"}]], "Input", CellChangeTimes->{{3.474996867725772*^9, 3.474996899417407*^9}, { 3.474997281755753*^9, 3.4749972884515963`*^9}}], Cell[BoxData[ RowBox[{ RowBox[{"(*", " ", RowBox[{"Parameters", " ", "to", " ", "control", " ", "the", " ", "run"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"nx", " ", "=", " ", "500"}], ";"}], " ", RowBox[{"(*", " ", RowBox[{"Disk", " ", "resolution", " ", "in", " ", "cells"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"xmin", " ", "=", " ", "0.1"}], ";"}], " ", RowBox[{"(*", " ", RowBox[{"Minimum", " ", "disk", " ", "radius", " ", RowBox[{"(", RowBox[{ RowBox[{"choose", " ", "a", " ", "small", " ", "number"}], " ", ">", " ", "0"}], " "}]}], "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"tmax", " ", "=", " ", "4"}], ";", " ", RowBox[{"(*", " ", RowBox[{ "Duration", " ", "of", " ", "run", " ", "in", " ", "dimensionless", " ", "time", " ", "units"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"tsave", " ", "=", " ", "0.05"}], ";"}], " ", RowBox[{"(*", " ", RowBox[{ "Time", " ", "interval", " ", "at", " ", "which", " ", "to", " ", "save", " ", "the", " ", "state"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Eta]", "=", RowBox[{"3", "/", "2"}]}], ";", " ", RowBox[{"(*", " ", RowBox[{"Turbulent", " ", "dissipation", " ", "rate"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"fg", "=", RowBox[{"1", "/", "2"}]}], ";"}], " ", RowBox[{"(*", " ", RowBox[{"Gas", " ", "fraction", " ", RowBox[{"(", "fixed", ")"}]}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"maxstep", " ", "=", "10000"}], ";", " ", RowBox[{"(*", " ", RowBox[{ "Maximum", " ", "number", " ", "of", " ", "steps", " ", "to", " ", "take"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"outbasename", "=", "\"\\""}], ";", " ", RowBox[{"(*", " ", RowBox[{"Base", " ", "name", " ", "of", " ", "output", " ", "file"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"kdiff", "=", "0.005"}], ";", " ", RowBox[{"(*", " ", RowBox[{"Magnitude", " ", "of", " ", "numerical", " ", "diffusion"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"outfreq", "=", "100"}], ";", " ", RowBox[{"(*", " ", RowBox[{"Frequency", " ", "of", " ", "status", " ", "update", " ", RowBox[{"(", RowBox[{ "set", " ", "to", " ", "0", " ", "for", " ", "no", " ", "updates"}], ")"}]}], " ", "*)"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Accretion", " ", "rate", " ", "onto", " ", RowBox[{"disk", "."}]}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"\[Chi]", "=", "0.1"}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Initial", " ", RowBox[{"conditions", ".", " ", "The"}], " ", "user", " ", "must", " ", "specify", " ", "a", " ", "function", " ", "sinit", " ", "that", " ", "give", " ", "the", " ", "initial", " ", "value", " ", "for", " ", RowBox[{ RowBox[{"s", "[", "x", "]"}], ".", " ", "An"}], " ", "example", " ", "is", " ", "provided", " ", "that", " ", "sets", " ", "the", " ", "initial", " ", "velocity", " ", "dispersion", " ", "to", " ", "twice", " ", "the", " ", "equilibirium", " ", RowBox[{"value", "."}]}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"sinit", "[", "x_", "]"}], ":=", RowBox[{"2", " ", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"\[Chi]", "/", RowBox[{"(", RowBox[{"\[Eta]", " ", "fg"}], ")"}]}], ")"}], "^", RowBox[{"(", RowBox[{"1", "/", "3"}], ")"}]}], "/", RowBox[{"Sqrt", "[", "2", "]"}]}]}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Parameters", " ", "to", " ", "control", " ", "the", " ", "ODE", " ", RowBox[{"solver", ".", " ", "Only"}], " ", "alter", " ", "these", " ", "if", " ", "you", " ", "know", " ", "what", " ", RowBox[{"you", "'"}], "re", " ", "doing"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"workprec", "=", "$MachinePrecision"}], ";", " ", RowBox[{"(*", " ", RowBox[{ RowBox[{"Working", " ", RowBox[{"precision", ".", " ", "Increasing"}], " ", "this", " ", "gives", " ", "more", " ", "accuracy"}], ",", " ", RowBox[{ "allowing", " ", "stiffer", " ", "systems", " ", "to", " ", "be", " ", "solved"}], ",", "\[IndentingNewLine]", RowBox[{ "but", " ", "slows", " ", "things", " ", "down", " ", "*", "a", " ", "lot", "*", RowBox[{".", " ", "Alter"}], " ", "only", " ", "if", " ", "you", " ", "start", " ", "getting", " ", "lots", " ", "of", " ", "ODE", " ", "boundary", " ", "condition", " ", RowBox[{"warnings", "."}]}]}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"maxODEstep", "=", "10000"}], ";", " ", RowBox[{"(*", " ", RowBox[{ "Maximum", " ", "number", " ", "of", " ", "steps", " ", "for", " ", "the", " ", "ODE", " ", "solver", " ", "to", " ", RowBox[{"take", ".", " ", "May"}], " ", "need", " ", "to", " ", "be", " ", "increased", " ", "for", " ", "very", " ", "stiff", " ", "systems"}], " ", "*)"}], "\[IndentingNewLine]"}]}]}]], "Input", CellChangeTimes->{{3.4749969021779947`*^9, 3.474997043816266*^9}, { 3.474997103944409*^9, 3.4749971401995173`*^9}, {3.474997234702406*^9, 3.474997272813674*^9}, {3.474997380749653*^9, 3.474997402547162*^9}, { 3.4749975162573757`*^9, 3.474997531585354*^9}, {3.474998625490718*^9, 3.474998670129097*^9}, 3.474999400015346*^9, {3.474999795240542*^9, 3.474999795871346*^9}, {3.474999836743084*^9, 3.474999838094569*^9}, { 3.475000142436895*^9, 3.475000142802535*^9}, 3.4750075091905823`*^9, { 3.475008216637669*^9, 3.47500825009105*^9}, {3.475008367827375*^9, 3.475008403041665*^9}, {3.475008561474432*^9, 3.475008561783391*^9}, { 3.47500869216505*^9, 3.4750087158911057`*^9}, {3.475010162986949*^9, 3.475010166296692*^9}, {3.475180383150744*^9, 3.475180385033599*^9}, { 3.475180439793973*^9, 3.475180440808362*^9}, 3.4751804848682947`*^9, 3.475180719741984*^9, {3.4751995742166862`*^9, 3.475199616695016*^9}, { 3.4751999726478767`*^9, 3.47519998945464*^9}, {3.475200052793407*^9, 3.475200058934082*^9}, {3.475200100325251*^9, 3.475200112519825*^9}, { 3.475203893084903*^9, 3.475203910449291*^9}, 3.4752076120547323`*^9, { 3.475209759423234*^9, 3.4752097595416937`*^9}, {3.475209846816435*^9, 3.4752098474538937`*^9}, 3.47521024383393*^9, {3.4752103664461193`*^9, 3.475210378243372*^9}, {3.4752111912578373`*^9, 3.475211224968378*^9}, { 3.475211325418968*^9, 3.47521133079608*^9}, 3.475213543643821*^9, { 3.475240393137054*^9, 3.475240410262514*^9}, {3.475240476357739*^9, 3.475240497060203*^9}, {3.47524060690458*^9, 3.475240609118251*^9}, 3.475240681844304*^9, {3.47524523904965*^9, 3.4752452409168663`*^9}, { 3.475250993794153*^9, 3.475250993953017*^9}, {3.475251072311255*^9, 3.475251072422171*^9}, {3.4752511241827183`*^9, 3.475251126133649*^9}, 3.475251260518099*^9, 3.47525130965331*^9, {3.475251547527857*^9, 3.475251548359254*^9}, 3.475252407115344*^9, 3.47525259814497*^9, 3.475252687800002*^9, {3.475252767964892*^9, 3.475252768188078*^9}, 3.475252936507454*^9, {3.475254473347845*^9, 3.475254531634942*^9}, { 3.475254730608322*^9, 3.47525474761478*^9}, 3.4752561576736794`*^9, { 3.475256771880575*^9, 3.4752567799340057`*^9}, {3.47525696322014*^9, 3.475256963755486*^9}, {3.4752581676576242`*^9, 3.4752581720646772`*^9}, 3.475258732809083*^9, {3.475258768816201*^9, 3.4752587783346167`*^9}, 3.475259316105179*^9, 3.47525948659663*^9, {3.475260628879097*^9, 3.4752606294902487`*^9}, {3.475260705897814*^9, 3.475260706097333*^9}, { 3.4752609846386538`*^9, 3.475260993444689*^9}, {3.475261407592259*^9, 3.475261425061809*^9}, 3.4752614596076183`*^9, 3.475261618012952*^9, { 3.4752616936742897`*^9, 3.475261693761615*^9}, {3.475261915760272*^9, 3.475261915958871*^9}, 3.475262092197974*^9, {3.475262284386971*^9, 3.4752622911692057`*^9}, {3.475264512633052*^9, 3.475264520151976*^9}, { 3.475264671230968*^9, 3.4752646720848722`*^9}, {3.475268075932087*^9, 3.4752680775229692`*^9}, {3.4752681514894342`*^9, 3.475268166000928*^9}, { 3.4752682086906013`*^9, 3.47526821855973*^9}, {3.475503970713833*^9, 3.4755039732238007`*^9}, 3.47550410000711*^9, 3.475504139333968*^9, { 3.475504173044417*^9, 3.475504186083543*^9}, 3.4755043125386667`*^9, 3.475504432905954*^9, {3.475504745021304*^9, 3.4755047451559258`*^9}, { 3.4755048285704527`*^9, 3.475504828674267*^9}, {3.47550524139965*^9, 3.475505290660439*^9}, 3.4755053249966*^9, {3.475505363116128*^9, 3.475505364618867*^9}, {3.4755054442114067`*^9, 3.475505444281752*^9}, { 3.4755055383614902`*^9, 3.47550553847998*^9}, {3.475505652199999*^9, 3.475505652478155*^9}, {3.4755056930714273`*^9, 3.475505732149128*^9}, { 3.475505794517658*^9, 3.4755057950994873`*^9}, {3.475505829923983*^9, 3.475505830451878*^9}, {3.475506048457738*^9, 3.475506049400505*^9}, { 3.4755061220734386`*^9, 3.475506145599457*^9}, {3.475506327382201*^9, 3.475506327460204*^9}, {3.475506803357852*^9, 3.475506823565814*^9}, { 3.4755068673185377`*^9, 3.475506867383564*^9}, {3.475506899349825*^9, 3.4755069288968287`*^9}, {3.475506960140777*^9, 3.4755069658278313`*^9}, { 3.475510499272355*^9, 3.4755104994404087`*^9}, {3.47552924411198*^9, 3.475529247302082*^9}, {3.475529578610878*^9, 3.475529579808934*^9}, { 3.475529620737352*^9, 3.4755296329448233`*^9}, {3.475529742858221*^9, 3.475529743167864*^9}, {3.47552985634262*^9, 3.475529862349218*^9}, { 3.475530098588661*^9, 3.475530098849536*^9}, {3.47553012965899*^9, 3.47553013001469*^9}, {3.475530710075275*^9, 3.475530711433372*^9}, { 3.475531126932589*^9, 3.475531127075162*^9}, {3.4755327049593477`*^9, 3.475532709197322*^9}, {3.475533354308632*^9, 3.475533354363552*^9}, { 3.4755339243426027`*^9, 3.4755339243960247`*^9}, {3.475590387199382*^9, 3.4755903947731237`*^9}, 3.4755919928156*^9, 3.475592272210308*^9, 3.4755927980998507`*^9, {3.475610851071475*^9, 3.475610854070095*^9}, { 3.47561562991635*^9, 3.475615632777451*^9}, {3.475861086955414*^9, 3.475861128361848*^9}, {3.4758612170248823`*^9, 3.475861269431497*^9}, { 3.4758613807343273`*^9, 3.475861436774602*^9}, {3.475861818619561*^9, 3.475861818728752*^9}, {3.475861941609796*^9, 3.475861941660232*^9}, { 3.4758619759032173`*^9, 3.475861985478351*^9}, {3.475862015992098*^9, 3.475862040012515*^9}, {3.4758620987343407`*^9, 3.475862113422594*^9}, { 3.4758621830124817`*^9, 3.475862183156332*^9}, {3.4758626451765842`*^9, 3.475862675028138*^9}, {3.4758627086548347`*^9, 3.475862708717319*^9}, { 3.475862826940114*^9, 3.475862827028674*^9}, {3.475862868181262*^9, 3.475862868259572*^9}, {3.475863800184835*^9, 3.4758638003989887`*^9}, 3.475863862303906*^9, {3.47586392591159*^9, 3.475863926029113*^9}, { 3.475864092517684*^9, 3.4758641135385933`*^9}, {3.475866114585766*^9, 3.475866123352128*^9}, {3.475868189125601*^9, 3.4758681969229183`*^9}, { 3.475877644678445*^9, 3.475877644925768*^9}, {3.475947060793048*^9, 3.4759470717171803`*^9}, {3.4759472864675703`*^9, 3.4759473650807543`*^9}, {3.475947425056711*^9, 3.475947484518352*^9}, { 3.4759475681125317`*^9, 3.475947568237443*^9}, {3.475947665357347*^9, 3.475947704283649*^9}, {3.4759478256974688`*^9, 3.4759479458161983`*^9}, { 3.475947988714519*^9, 3.4759481204229603`*^9}, {3.475948189016342*^9, 3.475948241644761*^9}, {3.475948513188047*^9, 3.475948527399261*^9}, { 3.475948748639737*^9, 3.475948748989387*^9}, {3.47595005281129*^9, 3.475950054863411*^9}, {3.475961631203006*^9, 3.475961669166939*^9}, { 3.475961849681617*^9, 3.475961851748334*^9}, {3.476023587956935*^9, 3.476023629603334*^9}, {3.476032166906905*^9, 3.476032171841279*^9}, { 3.4760419961806583`*^9, 3.476042008290184*^9}, {3.476049575474523*^9, 3.47604958398635*^9}, {3.476124905945891*^9, 3.4761249458775883`*^9}, { 3.476125101848288*^9, 3.476125101956749*^9}, {3.476125499502365*^9, 3.476125499606472*^9}, {3.4761258517399483`*^9, 3.476125851916665*^9}, { 3.4761270538925543`*^9, 3.476127054057143*^9}, {3.47613431966473*^9, 3.476134380030322*^9}}], Cell[BoxData[ RowBox[{"(*", " ", RowBox[{ "Below", " ", "is", " ", "the", " ", "computation", " ", "part", " ", "of", " ", "the", " ", RowBox[{"code", ".", " ", RowBox[{"Don", "'"}]}], "t", " ", "alter", " ", "this", " ", "unless", " ", "you", " ", "know", " ", "what", " ", RowBox[{"you", "'"}], "re", " ", RowBox[{"doing", "."}]}], " ", "*)"}]], "Input", CellChangeTimes->{{3.474997407683619*^9, 3.4749974248652353`*^9}, { 3.474999416308052*^9, 3.474999424008504*^9}}], Cell[BoxData[ RowBox[{ RowBox[{"(*", " ", RowBox[{ "Function", " ", "to", " ", "return", " ", "the", " ", "minmod", " ", "derivative", " ", "on", " ", "our", " ", "grid"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{ RowBox[{"minmodslope", "[", "vec_", "]"}], ":=", "\[IndentingNewLine]", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", RowBox[{"lslope", ",", " ", "rslope", ",", "prod", ",", "minmodsl"}], "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"lslope", "=", RowBox[{ RowBox[{"1", "/", RowBox[{"xgrid", "[", RowBox[{"[", RowBox[{"2", ";;", RowBox[{"nx", "-", "1"}]}], "]"}], "]"}]}], RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"vec", "[", RowBox[{"[", "i", "]"}], "]"}], "-", RowBox[{"vec", "[", RowBox[{"[", RowBox[{"i", "-", "1"}], "]"}], "]"}]}], ")"}], "/", RowBox[{"(", "dlnx", ")"}]}], ",", RowBox[{"{", RowBox[{"i", ",", "2", ",", RowBox[{"nx", "-", "1"}]}], "}"}]}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"rslope", "=", RowBox[{ RowBox[{"1", "/", RowBox[{"xgrid", "[", RowBox[{"[", RowBox[{"2", ";;", RowBox[{"nx", "-", "1"}]}], "]"}], "]"}]}], RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"vec", "[", RowBox[{"[", RowBox[{"i", "+", "1"}], "]"}], "]"}], "-", RowBox[{"vec", "[", RowBox[{"[", "i", "]"}], "]"}]}], ")"}], "/", RowBox[{"(", "dlnx", ")"}]}], ",", RowBox[{"{", RowBox[{"i", ",", "2", ",", RowBox[{"nx", "-", "1"}]}], "}"}]}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"prod", "=", RowBox[{"lslope", " ", "rslope"}]}], ";", "\[IndentingNewLine]", RowBox[{"minmodsl", "=", RowBox[{"Table", "[", "\[IndentingNewLine]", RowBox[{ RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"prod", "[", RowBox[{"[", "i", "]"}], "]"}], ">", "0"}], ",", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"Abs", "[", RowBox[{"lslope", "[", RowBox[{"[", "i", "]"}], "]"}], "]"}], "<", RowBox[{"Abs", "[", RowBox[{"rslope", "[", RowBox[{"[", "i", "]"}], "]"}], "]"}]}], ",", RowBox[{"lslope", "[", RowBox[{"[", "i", "]"}], "]"}], ",", RowBox[{"rslope", "[", RowBox[{"[", "i", "]"}], "]"}]}], "]"}], ",", "\[IndentingNewLine]", "0"}], "]"}], ",", " ", RowBox[{"{", RowBox[{"i", ",", RowBox[{"nx", "-", "2"}]}], "}"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"Flatten", "[", RowBox[{"{", RowBox[{ RowBox[{"minmodsl", "[", RowBox[{"[", "1", "]"}], "]"}], ",", "minmodsl", ",", RowBox[{"minmodsl", "[", RowBox[{"[", RowBox[{"nx", "-", "2"}], "]"}], "]"}]}], "}"}], "]"}]}]}], "\[IndentingNewLine]", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"censlope", "[", "vec_", "]"}], ":=", "\[IndentingNewLine]", RowBox[{"Module", "[", RowBox[{ RowBox[{"{", "cslope", "}"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"cslope", "=", RowBox[{ RowBox[{"1", "/", RowBox[{"xgrid", "[", RowBox[{"[", RowBox[{"2", ";;", RowBox[{"nx", "-", "1"}]}], "]"}], "]"}]}], " ", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"vec", "[", RowBox[{"[", RowBox[{"i", "+", "1"}], "]"}], "]"}], "-", RowBox[{"vec", "[", RowBox[{"[", RowBox[{"i", "-", "1"}], "]"}], "]"}]}], ")"}], "/", RowBox[{"(", RowBox[{"2", " ", "dlnx"}], ")"}]}], ",", RowBox[{"{", RowBox[{"i", ",", "2", ",", RowBox[{"nx", "-", "1"}]}], "}"}]}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"Flatten", "[", RowBox[{"{", RowBox[{ RowBox[{"cslope", "[", RowBox[{"[", "1", "]"}], "]"}], ",", "cslope", ",", RowBox[{"cslope", "[", RowBox[{"[", RowBox[{"nx", "-", "2"}], "]"}], "]"}]}], "}"}], "]"}]}]}], "\[IndentingNewLine]", "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Establish", " ", "grid", " ", "and", " ", "initial", " ", "s", " ", "and", " ", RowBox[{"s", "'"}], " ", "functions"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"dlnx", "=", RowBox[{ RowBox[{"-", RowBox[{"Log", "[", "xmin", "]"}]}], "/", RowBox[{"(", RowBox[{"nx", "-", "1"}], ")"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"xgrid", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"xmin", RowBox[{ RowBox[{"(", RowBox[{"1", "/", "xmin"}], ")"}], "^", RowBox[{"(", RowBox[{ RowBox[{"(", RowBox[{"i", "-", "1"}], ")"}], "/", RowBox[{"(", RowBox[{"nx", "-", "1"}], ")"}]}], ")"}]}]}], ",", RowBox[{"{", RowBox[{"i", ",", "nx"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"sgrid", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"sinit", "[", RowBox[{"xgrid", "[", RowBox[{"[", "i", "]"}], "]"}], "]"}], ",", RowBox[{"{", RowBox[{"i", ",", "nx"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"xsgrid", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", "i", "]"}], "]"}], ",", RowBox[{"sgrid", "[", RowBox[{"[", "i", "]"}], "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"i", ",", "nx"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"sfunc", "=", RowBox[{"Interpolation", "[", "xsgrid", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"spgrid", "=", RowBox[{"minmodslope", "[", "sgrid", "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"xspgrid", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", "i", "]"}], "]"}], ",", RowBox[{"spgrid", "[", RowBox[{"[", "i", "]"}], "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"i", ",", "nx"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"spfunc", "=", RowBox[{"Interpolation", "[", "xspgrid", "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Establish", " ", "grid", " ", "to", " ", "save", " ", "the", " ", "run", " ", "state"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"savectr", "=", "1"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"ssavegrid", "=", RowBox[{"Table", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"i", ",", "nx"}], "}"}], ",", RowBox[{"{", RowBox[{"j", ",", RowBox[{ RowBox[{"tmax", "/", "tsave"}], "+", "1"}]}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"ssavegrid", "[", RowBox[{"[", "savectr", "]"}], "]"}], "=", "sgrid"}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Define", " ", "h", " ", "coefficient", " ", "functions"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"h0", "[", "x_", "]"}], ":=", RowBox[{ RowBox[{ RowBox[{"-", "1"}], "/", RowBox[{"(", RowBox[{ RowBox[{"3", " ", "fg"}], "-", "1"}], ")"}]}], "/", RowBox[{"(", RowBox[{ RowBox[{"x", "^", "2"}], " ", RowBox[{ RowBox[{"sfunc", "[", "x", "]"}], "^", "2"}]}], ")"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"h1", "[", "x_", "]"}], ":=", RowBox[{ RowBox[{"-", "5"}], " ", "x", " ", RowBox[{ RowBox[{"spfunc", "[", "x", "]"}], "/", RowBox[{"(", RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"3", "fg"}], "-", "1"}], ")"}], " ", RowBox[{"sfunc", "[", "x", "]"}], " ", "x"}], ")"}]}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"H", "[", "x_", "]"}], ":=", RowBox[{ RowBox[{"\[Eta]", "/", "\[Chi]"}], " ", RowBox[{"(", RowBox[{"2", " ", "fg", " ", RowBox[{ RowBox[{"Sqrt", "[", "2", "]"}], "/", RowBox[{"(", RowBox[{ RowBox[{"3", " ", "fg"}], "-", "1"}], ")"}]}]}], ")"}], " ", RowBox[{ RowBox[{"sfunc", "[", "x", "]"}], "/", "x"}]}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Dump", " ", "initial", " ", "state"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Export", "[", RowBox[{ RowBox[{"outbasename", "<>", "\"\<.\>\"", "<>", RowBox[{"IntegerString", "[", RowBox[{ RowBox[{"savectr", "-", "1"}], ",", "10", ",", "5"}], "]"}]}], ",", RowBox[{"ssavegrid", "[", RowBox[{"[", "1", "]"}], "]"}], ",", "\"\\""}], "]"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"savectr", "=", RowBox[{"savectr", "+", "1"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Begin", " ", "main", " ", "loop"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"t", " ", "=", " ", "0.0"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"ctr", "=", "0"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"saveflag", "=", "0"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"kelast", "=", RowBox[{ RowBox[{"Total", "[", RowBox[{ RowBox[{"sgrid", "^", "3"}], " ", "xgrid"}], "]"}], " ", RowBox[{"Sqrt", "[", "2", "]"}], " ", RowBox[{"fg", "/", RowBox[{"(", RowBox[{"\[Pi]", " ", "\[Chi]"}], ")"}]}], " ", "dlnx"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"xshoot", "=", "0.3"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"dxshmax", "=", "0.02"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"dxshmin", "=", RowBox[{"10", "^", RowBox[{"-", "4"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"dxsh", "=", "dxshmax"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Tau]hlast", "=", RowBox[{"-", "xshoot"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Tau]phlast", "=", RowBox[{"-", "1"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"lastdir", "=", "0"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"While", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"t", "<", "tmax"}], ")"}], " ", "&&", " ", RowBox[{"(", RowBox[{"ctr", "<", "maxstep"}], ")"}]}], ",", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Solve", " ", "for", " ", "toque", " ", "at", " ", "this", " ", "time"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Tau]sol", "=", RowBox[{ RowBox[{"NDSolve", "[", RowBox[{ RowBox[{"{", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{ RowBox[{ RowBox[{"\[Tau]", "''"}], "[", "x", "]"}], "+", RowBox[{ RowBox[{"h1", "[", "x", "]"}], " ", RowBox[{ RowBox[{"\[Tau]", "'"}], "[", "x", "]"}]}], "+", RowBox[{ RowBox[{"h0", "[", "x", "]"}], RowBox[{"\[Tau]", "[", "x", "]"}]}]}], "\[Equal]", RowBox[{"H", "[", "x", "]"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"\[Tau]", "'"}], "[", "1", "]"}], "\[Equal]", RowBox[{"-", "1"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"\[Tau]", "[", "xmin", "]"}], "\[Equal]", RowBox[{"-", "xmin"}]}]}], "\[IndentingNewLine]", "}"}], ",", " ", "\[Tau]", ",", RowBox[{"{", RowBox[{"x", ",", "xmin", ",", "1"}], "}"}], ",", "\[IndentingNewLine]", RowBox[{"Method", "\[Rule]", RowBox[{"{", RowBox[{"\"\\"", ",", RowBox[{"\"\\"", "\[Rule]", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"\[Tau]", "'"}], "[", "xshoot", "]"}], "\[Equal]", "\[Tau]phlast"}], ",", RowBox[{ RowBox[{"\[Tau]", "[", "xshoot", "]"}], "\[Equal]", "\[Tau]hlast"}]}], "}"}]}]}], "}"}]}], ",", "\[IndentingNewLine]", RowBox[{"WorkingPrecision", "\[Rule]", "workprec"}], ",", " ", RowBox[{"MaxSteps", "\[Rule]", "maxODEstep"}]}], "\[IndentingNewLine]", "]"}], "[", RowBox[{"[", "1", "]"}], "]"}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Monitor", " ", "ODE", " ", "solution", " ", "error", " ", "and", " ", "adjust", " ", "fitting", " ", "point", " ", "if", " ", "necessary"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"errlo", "=", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"xmin", "+", RowBox[{"\[Tau]", "[", "xmin", "]"}]}], ")"}], "/", "xmin"}], "/.", "\[Tau]sol"}]}], ";", "\[IndentingNewLine]", RowBox[{"errhi", "=", RowBox[{ RowBox[{"1", "+", RowBox[{ RowBox[{"\[Tau]", "'"}], "[", "1", "]"}]}], "/.", "\[Tau]sol"}]}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Abs", "[", "errlo", "]"}], " ", "<", " ", RowBox[{"0.5", RowBox[{"Abs", "[", "errhi", "]"}]}]}], ")"}], " ", "&&", " ", RowBox[{"(", RowBox[{ RowBox[{"xshoot", "+", "dxsh"}], "<", "1"}], ")"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"xshoot", "=", RowBox[{"xshoot", "+", "dxsh"}]}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"lastdir", "\[Equal]", "1"}], ")"}], " ", "&&", " ", RowBox[{"(", RowBox[{"dxsh", ">", "dxshmin"}], ")"}]}], ",", RowBox[{"dxsh", "=", RowBox[{"dxsh", "/", "2"}]}]}], "]"}], ";", "\[IndentingNewLine]", RowBox[{"lastdir", "=", RowBox[{"-", "1"}]}], ";"}]}], "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Abs", "[", "errhi", "]"}], " ", "<", " ", RowBox[{"0.5", RowBox[{"Abs", "[", "errlo", "]"}]}]}], ")"}], " ", "&&", " ", RowBox[{"(", RowBox[{ RowBox[{"xshoot", "-", "dxsh"}], ">", "xmin"}], ")"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"xshoot", "=", RowBox[{"xshoot", "-", "dxsh"}]}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"lastdir", "\[Equal]", RowBox[{"-", "1"}]}], ")"}], " ", "&&", " ", RowBox[{"(", RowBox[{"dxsh", ">", "dxshmin"}], ")"}]}], ",", RowBox[{"dxsh", "=", RowBox[{"dxsh", "/", "2"}]}]}], "]"}], ";", "\[IndentingNewLine]", RowBox[{"lastdir", "=", "1"}], ";"}]}], "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Abs", "[", "errlo", "]"}], "<", RowBox[{"2", " ", RowBox[{"Abs", "[", "errhi", "]"}]}]}], ")"}], " ", "&&", " ", RowBox[{"(", RowBox[{ RowBox[{"Abs", "[", "errhi", "]"}], "<", RowBox[{"2", " ", RowBox[{"Abs", "[", "errlo", "]"}]}]}], ")"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"lastdir", "\[Equal]", "0"}], ")"}], " ", "&&", " ", RowBox[{"(", RowBox[{"dxsh", "<", "dxshmax"}], ")"}]}], ",", RowBox[{"dxsh", "=", RowBox[{"2", " ", "dxsh"}]}]}], "]"}], ";", "\[IndentingNewLine]", RowBox[{"lastdir", "=", "0"}], ";"}]}], "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", RowBox[{"\[Tau]hlast", "=", RowBox[{ RowBox[{"\[Tau]", "[", "xshoot", "]"}], "/.", "\[Tau]sol"}]}], ";", "\[IndentingNewLine]", RowBox[{"\[Tau]phlast", "=", RowBox[{ RowBox[{ RowBox[{"\[Tau]", "'"}], "[", "xshoot", "]"}], "/.", "\[Tau]sol"}]}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"Abs", "[", "errlo", "]"}], ">", "0.1"}], ")"}], " ", "||", " ", RowBox[{"(", RowBox[{ RowBox[{"Abs", "[", "errhi", "]"}], ">", "0.1"}], ")"}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"Print", "[", RowBox[{ "\"\\"", ",", " ", "errlo", ",", "\"\< \>\"", ",", " ", "errhi", ",", " ", "\"\< at step \>\"", ",", "ctr", ",", "\"\<, t = \>\"", ",", "t"}], "]"}], ";"}]}], "]"}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Compute", " ", "\[Tau]", " ", "and", " ", "its", " ", "derivatives", " ", "on", " ", "the", " ", "grid"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"\[Tau]grid", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{ RowBox[{"\[Tau]", "[", RowBox[{"xgrid", "[", RowBox[{"[", "i", "]"}], "]"}], "]"}], "/.", "\[Tau]sol"}], ",", RowBox[{"{", RowBox[{"i", ",", "nx"}], "}"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"\[Tau]pgrid", "=", RowBox[{"minmodslope", "[", "\[Tau]grid", "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"\[Tau]ppgrid", "=", RowBox[{"minmodslope", "[", RowBox[{"censlope", "[", "\[Tau]grid", "]"}], "]"}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Compute", " ", RowBox[{"ds", "/", "dT"}]}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"dsdtgrid", "=", RowBox[{ RowBox[{ RowBox[{"-", "2"}], " ", "\[Pi]", " ", "sgrid", " ", RowBox[{"\[Eta]", "/", RowBox[{"(", RowBox[{"3", "xgrid"}], ")"}]}]}], " ", "-", " ", RowBox[{"\[Pi]", " ", RowBox[{"\[Chi]", "/", RowBox[{"(", RowBox[{"3", " ", RowBox[{"Sqrt", "[", "2", "]"}], " ", "fg", " ", RowBox[{"sgrid", "^", "2"}], " ", RowBox[{"xgrid", "^", "2"}]}], ")"}]}], " ", "\[Tau]grid"}], " ", "-", " ", RowBox[{"5", " ", "\[Pi]", " ", "\[Chi]", " ", RowBox[{"spgrid", "/", RowBox[{"(", RowBox[{"3", " ", RowBox[{"Sqrt", "[", "2", "]"}], " ", "fg", " ", "sgrid"}], ")"}]}], " ", "\[Tau]pgrid"}], "-", RowBox[{"\[Pi]", " ", RowBox[{"\[Chi]", "/", RowBox[{"(", RowBox[{"3", " ", RowBox[{"Sqrt", "[", "2", "]"}], " ", "fg"}], ")"}]}], " ", "\[Tau]ppgrid"}]}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ RowBox[{"Compute", " ", "the", " ", "time", " ", "step"}], ",", " ", RowBox[{ "being", " ", "sure", " ", "not", " ", "to", " ", "overshoot"}]}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"dt", "=", RowBox[{"0.01", " ", RowBox[{"Min", "[", RowBox[{"Abs", "[", RowBox[{"sgrid", "/", "dsdtgrid"}], "]"}], "]"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"t", "+", "dt"}], ">", "tmax"}], ",", RowBox[{"dt", "=", RowBox[{"tmax", "-", "t"}]}]}], "]"}], ";", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{ RowBox[{"t", "+", "dt"}], ">", RowBox[{"tsave", " ", RowBox[{"(", RowBox[{"savectr", "-", "1"}], ")"}]}]}], ",", RowBox[{ RowBox[{"dt", "=", RowBox[{ RowBox[{"tsave", RowBox[{"(", RowBox[{"savectr", "-", "1"}], ")"}]}], " ", "-", "t"}]}], ";", " ", RowBox[{"saveflag", "=", "1"}]}]}], "]"}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Construct", " ", "s", " ", "and", " ", RowBox[{"s", "'"}], " ", "at", " ", "the", " ", "next", " ", "time", " ", "step"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"sgrid", "=", RowBox[{"sgrid", "+", RowBox[{"dt", " ", "dsdtgrid"}]}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Diffuse", " ", "kinetic", " ", "energy"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"kegrid", "=", RowBox[{ RowBox[{"sgrid", "^", "3"}], "/", "xgrid"}]}], ";", "\[IndentingNewLine]", RowBox[{"diffmat", "=", RowBox[{"SparseArray", "[", RowBox[{ RowBox[{"{", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"{", RowBox[{"1", ",", "1"}], "}"}], "\[Rule]", RowBox[{"1", "+", RowBox[{"kdiff", " ", RowBox[{"dt", "/", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", "1", "]"}], "]"}], "^", "2"}], " ", RowBox[{"dlnx", "^", "2"}]}], ")"}]}]}]}]}], " ", ",", RowBox[{ RowBox[{"Band", "[", RowBox[{ RowBox[{"{", RowBox[{"2", ",", "2"}], "}"}], ",", RowBox[{"{", RowBox[{ RowBox[{"nx", "-", "1"}], ",", RowBox[{"nx", "-", "1"}]}], "}"}]}], "]"}], "\[Rule]", RowBox[{"1", "+", RowBox[{"2", " ", "kdiff", " ", RowBox[{"dt", "/", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", RowBox[{"2", ";;", RowBox[{"nx", "-", "1"}]}], "]"}], "]"}], "^", "2"}], " ", RowBox[{"dlnx", "^", "2"}]}], ")"}]}]}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"{", RowBox[{"nx", ",", "nx"}], "}"}], "\[Rule]", RowBox[{"1", "+", RowBox[{"kdiff", " ", RowBox[{"dt", "/", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", "nx", "]"}], "]"}], "^", "2"}], " ", RowBox[{"dlnx", "^", "2"}]}], ")"}]}]}]}]}], ",", "\[IndentingNewLine]", " ", RowBox[{ RowBox[{"Band", "[", RowBox[{"{", RowBox[{"1", ",", "2"}], "}"}], "]"}], "\[Rule]", RowBox[{ RowBox[{"-", "kdiff"}], " ", RowBox[{"dt", "/", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", RowBox[{"1", ";;", RowBox[{"nx", "-", "1"}]}], "]"}], "]"}], "^", "2"}], " ", RowBox[{"dlnx", "^", "2"}]}], ")"}]}]}]}], ",", RowBox[{ RowBox[{"Band", "[", RowBox[{"{", RowBox[{"2", ",", "1"}], "}"}], "]"}], "\[Rule]", RowBox[{ RowBox[{"-", "kdiff"}], " ", RowBox[{"dt", "/", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", RowBox[{"2", ";;", "nx"}], "]"}], "]"}], "^", "2"}], " ", RowBox[{"dlnx", "^", "2"}]}], ")"}]}]}]}]}], "}"}], ",", RowBox[{"{", RowBox[{"nx", ",", "nx"}], "}"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"kenew", "=", RowBox[{"LinearSolve", "[", RowBox[{"diffmat", ",", "kegrid"}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"sgrid", "=", RowBox[{ RowBox[{"(", RowBox[{"kenew", " ", "xgrid"}], ")"}], "^", RowBox[{"(", RowBox[{"1", "/", "3"}], ")"}]}]}], ";", "\[IndentingNewLine]", RowBox[{"ketot", "=", RowBox[{ RowBox[{"Total", "[", RowBox[{ RowBox[{"sgrid", "^", "3"}], " ", "xgrid"}], "]"}], RowBox[{"Sqrt", "[", "2", "]"}], " ", RowBox[{"fg", "/", RowBox[{"(", RowBox[{"\[Pi]", " ", "\[Chi]"}], ")"}]}], " ", "dlnx"}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ "Generate", " ", "new", " ", "interpolating", " ", "functions", " ", "for", " ", "s", " ", "and", " ", RowBox[{"s", "'"}]}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"xsgrid", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", "i", "]"}], "]"}], ",", RowBox[{"sgrid", "[", RowBox[{"[", "i", "]"}], "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"i", ",", "nx"}], "}"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"sfunc", "=", RowBox[{"Interpolation", "[", "xsgrid", "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"spgrid", "=", RowBox[{"minmodslope", "[", "sgrid", "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"xspgrid", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"{", RowBox[{ RowBox[{"xgrid", "[", RowBox[{"[", "i", "]"}], "]"}], ",", RowBox[{"spgrid", "[", RowBox[{"[", "i", "]"}], "]"}]}], "}"}], ",", RowBox[{"{", RowBox[{"i", ",", "nx"}], "}"}]}], "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"spfunc", "=", RowBox[{"Interpolation", "[", RowBox[{"xspgrid", ",", RowBox[{"InterpolationOrder", "\[Rule]", "3"}]}], "]"}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Update", " ", "time", " ", "and", " ", "counter"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"t", "=", RowBox[{"t", "+", "dt"}]}], ";", "\[IndentingNewLine]", RowBox[{"ctr", "=", RowBox[{"ctr", "+", "1"}]}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", "Save", " ", "*)"}], "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{"t", "\[GreaterEqual]", RowBox[{"tsave", " ", RowBox[{"(", RowBox[{"savectr", "-", "1"}], ")"}]}]}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"ssavegrid", "[", RowBox[{"[", "savectr", "]"}], "]"}], "=", "sgrid"}], ";", "\[IndentingNewLine]", RowBox[{"Export", "[", RowBox[{ RowBox[{"outbasename", "<>", "\"\<.\>\"", "<>", RowBox[{"IntegerString", "[", RowBox[{ RowBox[{"savectr", "-", "1"}], ",", "10", ",", "5"}], "]"}]}], ",", RowBox[{"ssavegrid", "[", RowBox[{"[", "savectr", "]"}], "]"}], ",", "\"\\""}], "]"}], ";", "\[IndentingNewLine]", RowBox[{"savectr", "=", RowBox[{"savectr", "+", "1"}]}], ";"}]}], "\[IndentingNewLine]", "]"}], ";", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{"Write", " ", "status"}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{"outfreq", ">", "0"}], ",", "\[IndentingNewLine]", RowBox[{"If", "[", RowBox[{ RowBox[{"Divisible", "[", RowBox[{"ctr", ",", "outfreq"}], "]"}], ",", "\[IndentingNewLine]", RowBox[{ RowBox[{"stot", "=", RowBox[{"Total", "[", "sgrid", "]"}]}], ";", "\[IndentingNewLine]", RowBox[{"Print", "[", RowBox[{ "\"\\"", ",", "ctr", ",", "\"\<: t = \>\"", ",", "t", ",", "\"\<, dt = \>\"", ",", "dt", ",", " ", "\"\<, dke = \>\"", ",", " ", RowBox[{"ketot", "-", "kelast"}]}], "]"}]}]}], "]"}]}], "]"}], ";", "\[IndentingNewLine]", RowBox[{"kelast", "=", "ketot"}], ";"}]}], "\[IndentingNewLine]", "]"}], ";"}]}]}]], "Input", CellChangeTimes->{{3.474997417517845*^9, 3.474997464204397*^9}, { 3.4750083045192223`*^9, 3.475008305502061*^9}, {3.475199648483451*^9, 3.475199662771976*^9}, {3.475199718390427*^9, 3.475199951613265*^9}, { 3.475203652446167*^9, 3.475203657597926*^9}, {3.47520493490907*^9, 3.4752049365526323`*^9}, {3.4752079271900473`*^9, 3.4752079274219847`*^9}, {3.475208216251255*^9, 3.475208216307878*^9}, 3.4752097544447327`*^9, 3.475210075841832*^9, 3.4752134816366663`*^9, { 3.475240503755467*^9, 3.4752405038208313`*^9}, {3.4752516136778193`*^9, 3.475251631526338*^9}, {3.475252394616804*^9, 3.475252395821953*^9}, 3.4752524870150337`*^9, {3.475252760506626*^9, 3.4752527642093563`*^9}, { 3.4752528025696793`*^9, 3.4752528923117228`*^9}, 3.475252942710415*^9, { 3.475254453392529*^9, 3.475254453982201*^9}, {3.4752554114250383`*^9, 3.4752554311580257`*^9}, {3.475256198770466*^9, 3.4752562016909018`*^9}, { 3.475256232475298*^9, 3.475256257089217*^9}, {3.475259310129586*^9, 3.475259311801202*^9}, {3.475260610744933*^9, 3.475260621609857*^9}, { 3.47526078478251*^9, 3.4752607874124537`*^9}, {3.4752615224298677`*^9, 3.475261530827283*^9}, {3.475261837017815*^9, 3.475261845709169*^9}, { 3.475262033674716*^9, 3.475262036179125*^9}, {3.475262145155408*^9, 3.4752621529951677`*^9}, {3.475264817587757*^9, 3.47526481861684*^9}, { 3.475268124551646*^9, 3.475268124847033*^9}, {3.475503933974454*^9, 3.475503965765291*^9}, {3.475504178942194*^9, 3.475504180400734*^9}, { 3.4755047270931253`*^9, 3.475504731201713*^9}, {3.4755048185718937`*^9, 3.475504836531213*^9}, {3.475505124343555*^9, 3.475505147133534*^9}, { 3.4755052471707087`*^9, 3.475505261493829*^9}, 3.475505376031666*^9, { 3.475505450555887*^9, 3.4755055612347918`*^9}, {3.475505631968762*^9, 3.47550563395312*^9}, {3.4755056776113586`*^9, 3.475505677686504*^9}, { 3.475505868199876*^9, 3.4755059102632933`*^9}, {3.475505947389715*^9, 3.475505953916999*^9}, {3.4755061161734743`*^9, 3.475506117881803*^9}, { 3.47550615331826*^9, 3.475506154872521*^9}, {3.4755061985712423`*^9, 3.4755062278646393`*^9}, 3.475506323505481*^9, {3.475506545891642*^9, 3.4755065737815847`*^9}, {3.475506659428814*^9, 3.475506680553104*^9}, { 3.475506714345175*^9, 3.475506729955594*^9}, {3.475507649920438*^9, 3.475507654074243*^9}, {3.47552923314849*^9, 3.475529261546294*^9}, { 3.475529365941579*^9, 3.4755293691492243`*^9}, {3.475529585468748*^9, 3.475529638049492*^9}, {3.475529834122409*^9, 3.475529836255274*^9}, { 3.475529953167397*^9, 3.475530004120266*^9}, {3.475530068671484*^9, 3.475530080122314*^9}, {3.475530191920418*^9, 3.475530209488392*^9}, { 3.475530685320219*^9, 3.4755306884213123`*^9}, {3.4755308322540417`*^9, 3.475530891049656*^9}, {3.475530980706025*^9, 3.4755310093675213`*^9}, 3.475531905535264*^9, {3.475531953707556*^9, 3.4755319838664703`*^9}, 3.475533975967148*^9, {3.475590404729721*^9, 3.475590416273782*^9}, { 3.475591484774293*^9, 3.4755914937614317`*^9}, {3.475591651208124*^9, 3.475591903534731*^9}, {3.4755920344519167`*^9, 3.4755920349850483`*^9}, { 3.4755922022714663`*^9, 3.475592206450542*^9}, {3.4755922488812933`*^9, 3.47559229066127*^9}, {3.475592907733268*^9, 3.475592909356183*^9}, { 3.475593066740469*^9, 3.4755930888238487`*^9}, {3.475593268534079*^9, 3.47559327554251*^9}, {3.475593340672308*^9, 3.4755933444093113`*^9}, { 3.475593464631605*^9, 3.475593485629643*^9}, {3.47559358520163*^9, 3.475593622784289*^9}, {3.475593662177136*^9, 3.475593741124305*^9}, { 3.475610795933057*^9, 3.475610864739311*^9}, {3.475861164029917*^9, 3.4758611789925823`*^9}, {3.475861644897094*^9, 3.475861750801124*^9}, { 3.475861846959593*^9, 3.475861923038435*^9}, {3.475862061742535*^9, 3.475862067502318*^9}, {3.475862206290019*^9, 3.475862248670394*^9}, { 3.475862308274309*^9, 3.475862345226678*^9}, {3.475862520260133*^9, 3.4758625536903133`*^9}, {3.475862594696347*^9, 3.475862630304034*^9}, { 3.475862766281475*^9, 3.475862766456162*^9}, {3.4758629159125547`*^9, 3.475863142485787*^9}, {3.475863180428248*^9, 3.475863180500362*^9}, { 3.475863220331895*^9, 3.475863320321039*^9}, {3.475863361377619*^9, 3.4758634134818373`*^9}, {3.475863608800581*^9, 3.47586363405581*^9}, { 3.475863692341402*^9, 3.475863745582169*^9}, {3.475863982366568*^9, 3.475864082177923*^9}, {3.475879130897214*^9, 3.475879130976844*^9}, { 3.475947347298525*^9, 3.475947421040357*^9}, {3.475947603918736*^9, 3.4759476073466578`*^9}, {3.475947727848783*^9, 3.4759477435372066`*^9}, { 3.475947777048497*^9, 3.475947777114325*^9}, {3.4759478824490557`*^9, 3.475947883492934*^9}, {3.475948105762013*^9, 3.475948135175139*^9}, { 3.475961877433495*^9, 3.475961877640386*^9}, {3.476024560615917*^9, 3.476024619400283*^9}, {3.4761250412594852`*^9, 3.4761250876690903`*^9}, { 3.4761251944715242`*^9, 3.47612519453613*^9}, {3.4761253687490177`*^9, 3.4761253689615383`*^9}, 3.476125789267261*^9, {3.476134394627612*^9, 3.4761344886015244`*^9}}] }, WindowSize->{1216, 598}, WindowMargins->{{Automatic, -1984}, {421, Automatic}}, FrontEndVersion->"6.0 for Mac OS X x86 (32-bit) (June 19, 2007)", StyleDefinitions->"Default.nb" ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[568, 21, 916, 23, 43, "Input"], Cell[1487, 46, 12622, 221, 373, "Input"], Cell[14112, 269, 504, 11, 27, "Input"], Cell[14619, 282, 36862, 909, 2098, "Input"] } ] *) (* End of internal cache information *)