Laboratorio di fisica computazionale
(Prova d'esame del giugno 2005)

Un corpo di massa m  viene abbandonato in quiete su un piano lungo s  inclinato di un angolo β rispetto all'orizzontale e dotato di coefficiente di attrito μ. 
   Alla fine del piano il corpo,dopo aver percorso un tratto orizzontale d in presenza dello stesso coefficiente d'attrito μ, urta centralmente ed elasticamente un secondo corpo di uguali dimensioni e massa sospeso ad un filo inestensibile e di massa trascurabile lungo l≥ s Sin β, sospeso ad un punto fisso O posto sulla verticale del punto di contatto tra le due sfere
- Scrivere una procedura che calcoli il moto dei due corpi
- Analizzare i possibili tipi di moto al variare dei parametri 
- Riportare in un grafico le traiettorie dei due corpi e animare il moto.

Generazione della figura

Clear[a, b, h]
Fig=
Module[{O = {1, 1}, a = {0, 0}, b = {1, 0}, c = {0, 1}, 
              dr = {0.01, 0.01},  d = {0.2, 0}, φ},
O = O + d;
φ = 
                ArcTan[c[[2]]/b[[1]]];
Show[Graphics[{
	Line[{a, b, c, a}],Circle[b, 
                          0.14, {π - φ, π}],
Line[{O, b + d}],Line[{b, b + d}],
PointSize[0.01],Point[O],PointSize[0.03],Point[c + dr],Point[b + d + {0, 0.01}],
Text[StyleForm["l",FontFamily -> "Times",FontWeight -> "Bold",FontSize -> 14],((b + d + O)/2),{-2, -2}],Text[StyleForm["s",FontFamily -> "Times",FontWeight -> "Bold",FontSize -> 14],((b + c)/2),{-1, -2}],Text[
                          StyleForm["β", 
                            FontFamily -> "Times", 
                            FontWeight -> "Bold", 
                            FontSize -> 14], b, {3, -1}],Text[
                          StyleForm["O", FontFamily -> "Times", 
                            FontWeight -> "Bold", 
                            FontSize -> 14], O, {0, -1}],Text[StyleForm["d",FontFamily -> "Times",FontWeight -> "Bold",FontSize -> 14],b + (d/2),{0, -1}]
}],PlotRange -> All,AspectRatio -> Automatic]];

Off[General::"spell1"]; 
Unprotect[g];
Clear[Subscript, T, Olds, β, d, μ, l, Moto, 
      g]
g = 9.81;
Protect[g];
Start[k_, s_, β_, 
        d_, μ_] := -d Sec[β] + (s + 
              d Sec[β]) (((-1 + ((2 Sin[β])/(μ Cos[β] + Sin[β]))))^k)
TOld[α_] := 
        Re[(2/π) 
            Csc[(α/2)] EllipticF[(α/2), 
              Csc[(α/2)]^2]];
T[α_] := 
        Module[{θ, γ = 
              0.05}, 
If[α < γ, 
            1 + (1/(2^2)) Sin[(α/2)]^2 + (1/(2^2)) ((3^2)/(4^2)) 
                Sin[(α/2)]^4, ((Sqrt[2])/π) 
              NIntegrate[
                (1/(Sqrt[Cos[θ] - Cos[α]])), {θ, 
                  0, α}, SingularityDepth -> 100, 
                MaxRecursion -> 100]]];
\(Moto[s_, β_, d_, μ_, 
        l_] := 
Module[{a, τ, v, 
          x, θ, θm, Δ}, 
a = 
          g {Sin[β] - μ\ Cos[β], \(-μ\)\ , 0, μ, 
              Sin[β] + μ\ Cos[β]\ }; 
x[1, 
            t_] = 1\/2\ a[\([1]\)]\ t\^2 - s; 
τ[
            1] = \@\(\(2\ s\)\/a[\([1]\)]\); 
v[1] = 
          a[\([1]\)] τ[1]; 
x[2, 
            t_] = \(1\/2\) a[\([2]\)]\ \((t - τ[1])\)\^2 + 
            v[1] \((t - τ[1])\); 
Δ = 
          v[1]\^2 + 
            2\ a[\([2]\)]\ d; 
If[Δ ≤ 
            0, 
τ[2] = τ[1] - 
              v[1]/a[\([2]\)]; 
{{Map[Evaluate, 
                Which[0 ≤ t < τ[1], 
                  x[1, t] Cos[β], τ[1] ≤ 
                    t ≤ τ[2], 
                  x[2, t], τ[2] ≤ t, 
                  x[2, τ[2]]]], 
Map[Evaluate, 
                Which[0 ≤ t < τ[1], \(-x[1, t]\) 
                    Sin[β], τ[1] ≤ t, 0]]}, {d, 
              0}, τ[2]}, 
τ[2] = τ[1] + 
              If[a[\([2]\)] == 0, 
                d\/v[1], \(\(-v[1]\)\  + \@Δ\)\/a[\([2]\)]]; 
v[2] = 
            v[1]\  + \ 
              a[\([2\ ]\)] \((τ[2] - τ[
                      1])\); 
θm = 
            ArcCos[1 - v[2]\^2\/\(2  g\ l\)]; 
τ[
              3] = τ[
                2] + π\ \@\(l\/g\)\ T[θm]; 

θ = θ /. \(NDSolve[
{\(θ''\)[t] == \(-\(g\/l\)\) 
                        Sin[θ[t]], \(θ'\)[τ[2]] == 
                      v[2]\/l, θ[τ[2]] == 
                      0}, θ, {t, τ[2], τ[
                      3]}]\)[\([1]\)]; 
x[3, t_] = 
            d; 
v[3] = \(-v[2]\); 
x[4, 
              t_] = \(1\/2\) a[\([4]\)]\ \((t - τ[3])\)\^2 + 
              v[3] \((t - τ[3])\) + 
              d; 
Δ = 
            v[3]\^2 - 
              2\ a[\([4]\)]\ d; 
If[Δ ≤ 
              0, 
τ[4] = τ[3] - 
                v[3]/a[\([4]\)]; 
Expand[{{Map[Evaluate, 
                    Which[0 ≤ t < τ[1], 
                      x[1, t] Cos[β], τ[1] ≤ 
                        t ≤ τ[2], 
                      x[2, t], τ[2] ≤ 
                        t ≤ τ[3], 
                      x[3, t], τ[3] ≤ 
                        t ≤ τ[4], 
                      x[4, t], τ[4] ≤ t, 
                      x[4, τ[4]]]], 
Map[Evaluate, 
                    Which[0 ≤ t < τ[1], \(-x[1, t]\) 
                        Sin[β], τ[1] ≤ t, 
                      0]]}, 
{Map[Evaluate, 
                    Which[0 ≤ t < τ[2], 
                      d, τ[2] ≤ t ≤ τ[3], 
                      d + l\ Sin[θ[t]], τ[3] ≤ t, 
                      d]], 
Map[Evaluate, 
                    Which[0 ≤ t < τ[2], 
                      0, τ[2] ≤ t ≤ τ[3], 
                      l \((1 - \ Cos[θ[t]])\), τ[3] ≤ 
                        t, 0]]}, τ[4]}], 
τ[
                4] = τ[3] + 
                If[a[\([4]\)] == 
                    0, \(-\(d\/v[
                          3]\)\), \(\(-v[3]\)\  - \@Δ\)\/a[\([4]\)]]; 
v[4] = 
              v[3] + a[\([4]\)]\ \((τ[4] - τ[
                        3])\); 
x[5, 
                t_] = \(1\/2\) a[\([5]\)] \((t - τ[4])\)\^2 + 
                v[4]\ \((t - τ[4])\); 
τ[
                5] = τ[4] - 
                v[4]\/a[\([5]\)]; 
Expand[{{Map[Evaluate, 
                    Which[0 ≤ t ≤ τ[1], 
                      x[1, t] Cos[β], τ[1] < 
                        t ≤ τ[2], 
                      x[2, t], τ[2] < t ≤ τ[3], 
                      x[3, t], τ[3] < t ≤ τ[4], 
                      x[4, t], τ[4] < t ≤ τ[5], 
                      x[5, t] Cos[β]]], 
Map[
                    Evaluate, 
                    Which[0 ≤ 
                        t ≤ τ[1], \(-x[1, t]\) 
                        Sin[β], τ[1] < t ≤ τ[4], 
                      0, τ[4] < t ≤ τ[5], \(-x[5, t]\) 
                        Sin[β]]]}, 
{Map[Evaluate, 
                    Which[0 ≤ t ≤ τ[2], 
                      d, τ[2] < t ≤ τ[3], 
                      d + l\ Sin[θ[t]], τ[3] < 
                        t ≤ τ[5], d]], 
Map[
                    Evaluate, 
                    Which[0 ≤ t < τ[2], 
                      0, τ[2] ≤ t ≤ τ[3], 
                      l \((1 - \ Cos[θ[t]])\), τ[3] ≤ 
                        t ≤ τ[5], 0]]}, τ[
                  5]}]
]
]
]\)

Clear[r, rp]
μ = 0.2;
β = (π/4);
s0 = 1;
d = 0.2;
l = s0  Sin[β];
{r[t_], rp[t_], (Subscript[T, 0])} = Moto[s0, β, d, μ, l];

ParametricPlot[r[t], {t, 0, (Subscript[T, 0])}, 
        PlotRange -> {{-s0 Cos[β], d + l}, {0, 
              s0 Sin[β]}}, 
        DisplayFunction -> Identity];
ParametricPlot[rp[t], {t, 0, (Subscript[T, 0])}, 
        PlotRange -> {{-s0 Cos[β], d + l}, {0, 
              s0 Sin[β]}}, 
        DisplayFunction -> Identity];
Show[%, %%, AspectRatio -> Automatic, 
        DisplayFunction -> $DisplayFunction];
Plot[Evaluate[r[t]], {t, 0, (Subscript[T, 0])}, 
        PlotRange -> {{0, (Subscript[T, 0])}, {-s0 Sin[β], 
              d + l}}];
Plot[Evaluate[rp[t]], {t, 0, (Subscript[T, 0])}, 
        PlotRange -> {{0, (Subscript[T, 0])}, {-s0 Sin[β], 
              d + l}}];

μ = 0.1;
β = (π/4);
s0 = 1;
d = 0.2;
l = s0  Sin[β];
PSize = 0.02;
x1 = -1.05 s0 Cos[β];
x2 = d + l;
ϵ = 
        PSize (x2 - x1)/(2 - PSize);
x2 = x2 + ϵ;
kmax = 
        If[μ == 0, 0, k = 1; 
          While[Start[k, s0, β, d, μ] > 0, ++k]; 
          k - 1];
P1 = -ϵ {Tan[(β/2)], 
            1};
P2 = 
        s0 {- Cos[β], Sin[β]} - ϵ {Sin[β], 
              Cos[β]};
P3 = {-s0 Cos[β] -  ϵ Sin[β], -ϵ};
P4 = {d + 
            l + ϵ, -ϵ};
Sfondo = 
        Graphics[Line[{P1, P2, P3, P1, P4}]];
{r[t_], rp[t_], (Subscript[T, 0])} = 
        Moto[s0, β, d, μ, l];
np = 100;
Δt = Moto[s0, β, d, μ, l][[3]]/np;
Do[
Clear[r, rp]; 
s = 
        Start[k, s0, β, d, μ]; 
{r[t_], rp[t_], 
          (Subscript[T, 0])} = Moto[s, β, d, μ, l]; 
fg[k] = 
        Table[
Show[
            Sfondo, 
Graphics[{PointSize[PSize], 
                Point[r[t]], Point[rp[t]], 
                Line[{{d, l}, 
                    rp[t]}]}], 
PlotRange -> {{x1, 
                  x2}, {-0.1 l, Max[l, 1.05  s0 Sin[β]]}}, 
            AspectRatio -> Automatic], {t, 0, 
            (Subscript[T, 0]), Δt}], 
{k, 0, 
        kmax}]