# Butter side up?

I've written a piece of code that determines if a piece toast with butter lands on the butter-side or not, depending on its initial velocity and the table height.

The code is written in Mathematica and to be honest, it was my first time writing a piece of code in this language (which wasn't a simple calculation of an integral or something alike). It works and produces results I would expect but it incredibly slow and I don't even know where to start to correct it.

So what I'm looking for are some general hints on how to imporve when writing Mathematica code, style conventions (is there something like pep8 for Python?) and how to speed the code up.

ClearAll["Global*"]

MapMonitored[f_, args_List] := Module[{x = 0},
Monitor[MapIndexed[(x = #2[]; f[#1]) &, args],
ProgressIndicator[x/Length[args]]]]

xb:=trajectory[t]/.{b\[Rule]0.02,L\[Rule]0.1}+{-bb/2.0 \
Sin[ut],-bb/2.0 Cos[ut]};

xleftdown:=xb+{-LL/2.0 Cos[ut],LL/2.0 Sin[ut]};

xleftup:=xleftdown+{bb Sin[ut],bb Cos[ut]};

xrightup:=xleftup+{LL Cos[ut],-LL Sin[ut]};

xrightdown:=xrightup+{-bb Sin[ut],-bb Cos[ut]};

p=Block[{b=0.02,L=0.1},
Show[
ListLinePlot[
{Flatten[{xb,xleftdown,xleftup,xrightup,xrightdown,xb}]},
PlotStyle\[Rule]Blue
],
ListLinePlot[
{Flatten[xleftup],Flatten[xrightup]},
PlotStyle[Rule]Red
],
ListLinePlot[
{Flatten[xrightup],Flatten[xrightdown],Flatten[xb]},
PlotStyle\[Rule] Blue
],
ListPlot[
{0.0,trajectory[t]},
PlotStyle\[Rule]Red
],
AspectRatio\[Rule]1.0,
PlotRange\[Rule]{{-0.1,1.0},{0.05,-1.5}}
]
];

Return[p]];

uu[time_]=Flatten[u[time]/.sol];
uu0=Abs[uu[tst]-uu'[tst] tst+uu[tst]];
angle:=uu'[tst] t-uu[tst]+uu0;
xb:=wurfparabel[t,sol,tst,vtstar]+{-bb/2.0 Sin[angle],-bb/2.0 Cos[angle]};
xleftdown:=xb+{-LL/2.0 Cos[angle],LL/2.0 Sin[angle]};
xleftup:=xleftdown+{bb Sin[angle],bb Cos[angle]};
xrightup:=xleftup+{LL Cos[angle],-LL Sin[angle]};
xrightdown:=xrightup+{-bb Sin[angle],-bb Cos[angle]};

p=Block[{b=0.02,L=0.1},
Show[
ListLinePlot[
{Flatten[xb],Flatten[xleftdown],Flatten[xleftup],Flatten[xrightup],
Flatten[xrightdown],Flatten[xb]},
PlotRange\[Rule]{{-0.1,1.0},{0.05,-1.5}},
PlotStyle\[Rule]Blue
],
ListLinePlot[
{Flatten[xleftup],Flatten[xrightup]},
PlotRange\[Rule]{{-0.1,1.0},{0.05,-1.5}},
PlotStyle\[Rule]Red
],
ListLinePlot[
{Flatten[xrightup],Flatten[xrightdown],Flatten[xb]},
PlotRange\[Rule]{{-0.1,1.0},{0.05,-1.5}},
PlotStyle\[Rule]Blue
],
ListPlot[{Flatten[wurfparabel[t,sol,tst,vtstar]]},
PlotStyle\[Rule]Red],
AspectRatio\[Rule]1.0
]
];
Return[p]
];

];

Module[{bb = 0.02, LL = 0.1},

uu[time_] = Flatten[u[time] /. sol];
uu0 = Abs[uu[tst] - uu'[tst] tst + uu[tst]];
angle = uu'[tst] t - uu[tst] + uu0;
xb = wurfparabel[t, sol, tst,vtstar] + {-bb/2.0 Sin[angle], -bb/2.0 Cos[angle]};
xleftdown = xb + {-LL/2.0 Cos[angle], LL/2.0 Sin[angle]};
xleftup = xleftdown + {bb Sin[angle], bb Cos[angle]};
xrightup = xleftup + {LL Cos[angle], -LL Sin[angle]};
xrightdown = xrightup + {-bb Sin[angle], -bb Cos[angle]};

Return[
Abs[Min[xleftdown, xleftup, xrightup, xrightdown]]]
];

AngleAtGround[sol_, tst_, tt_] := Module[{},

uu[t_] = Flatten[u[t] /. sol];
uu0 = Abs[uu[tst] - uu'[tst] tst + uu[tst]];
angle = uu'[tst] tt - uu[tst] + uu0;

Return[Mod[angle[], 2.0 Pi]]];

If[Pi/2.0 <= phi, If[phi <= 3.0 Pi/2.0, 1.0, 0.0], 0.0];

wurfparabel[tt_, sol_, tst_, vtstar_] := Module[{g = 9.81},

beta = Abs[Flatten[u[tst] /. sol]];
xspace = vtstar[] tt Cos[beta];
xx0 = Abs[trajectory[tst][] - vtstar[] tst Cos[beta]];
yspace = -(Norm[v[tst]] tt Sin[beta] + g tt^2/2.0);
y0 = Abs[
trajectory[tst][[
2]] + (Norm[v[tst]] tst Sin[beta] + g tst^2/2.0)];
Return[{xspace + xx0, yspace + y0}]];

x[t_] = {{s[t] Cos[u[t]] + b/2 Sin[u[t]]}, {-s[t] Sin[u[t]] +
b/2 Cos[u[t]]}};

m = 0.1; l = 0.1; g = 9.81;
J = 1.0/12.0 m (l^2 + b^2);

T[t_] = Simplify[1.0/2.0 (m) Flatten[x'[t]].Flatten[x'[t]]] +
1.0/2.0 J u'[t]^2;
V[t_] = m g x[t][[2, 1]];
L[t_] = Simplify[T[t] - V[t]];

<< VariationalMethods

GetTheThingsINeed[StartingValuesList_, Lagrangian_: L,
PositionVector_: x] := Module[{g = 9.81},

b = 0.02;
eoms[t_] = Simplify[VariationalD[Lagrangian[t], s[t], t]];
eomu[t_] = Simplify[VariationalD[Lagrangian[t], u[t], t]];

TabelHeight = StartingValuesList[];

solution =
NDSolve[{eoms[t] == 0.0, eomu[t] == 0.0, u == 0.0, s == 0.01,
s' == StartingValuesList[], u' == 0.1}, {u, s}, {t,
0.0, 10.0}];

f[t_] = D[x[t][[1, 1]], {t, 2}] /. solution;

tstar = t /. FindRoot[f[t], {t, 0.1}];

trajectory[t_] := Flatten[x[t] /. solution];

v[t_] = Flatten[D[trajectory[t], t]];

tGroundTime =
t /. Flatten[
TabelHeight}, {t, 0.5}]];

AngleAtGround[solution, tstar, tGroundTime]]]];

velocities = Range[0.1, 1.8, 0.2];
TabelHeights = Range[0.3, 1.5, 0.8];
mesh = Tuples[{velocities, TabelHeights}];

plotmesh =
MapMonitored[GetTheThingsINeed, mesh] //
AbsoluteTiming; plotmesh[]

ListPlot[Lookup[
GroupBy[
Last -> First],
{1.0, 0.0}],
PlotStyle -> {Green, Red}, AxesLabel -> {"v", "h"}
]


Edit: I cleared the code up a bit.. Didn't impact performance at all.

Edit II: On a whim I just changed the line

tGroundTime =
t /. Flatten[
TabelHeight, t > 0.0}, t]];


to

tGroundTime =
t /. Flatten[

EDIT III: Cleared the code up again and moved as much as possible out of the GetTheThingsINeed out. Turns out this improves the performance again, now the code runs in about 0.8 sec compared to the 12 sec in the beginning.
• NSolve was searching for all possible positive solutions, whereas FindRoot uses something like Newton's method to find a single solution. – KraZug Apr 27 '19 at 19:15