Podajoci utezi, povezani z vzmetjo
Rok Zitko, 3/2011

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

ϵ=10^-8;

(* C0=2Sqrt[k l/(m1+m2)g *)
(* x=m1/m2 *)
(* calc[] racuna trajektorijo od zacetnega casa t0 do naslednjega dogodka (trka),
ob zacetnih pogojih ξinit (polozaj) in ξdinit (hitrosti). *)
calc[C0_,x_,t0_,ξ1init_,ξ2init_,ξ1dinit_,ξ2dinit_]:=Module[{tmax,ndsol,sol},
tmax=3;
tend=t0+tmax;
ClearAll[C1,C2];
eqs={
ξ1''[t]==-2+C1^2/2(ξ2[t]-ξ1[t]-1),
ξ2''[t]==-2-C2^2/2(ξ2[t]-ξ1[t]-1)
};
eqs=eqs/.{C1->C0 Sqrt[1+1/x],C2->C0 Sqrt[1+x]};
ndsol=NDSolve[
Join[eqs,
{ξ1[t0]==ξ1init,ξ2[t0]==ξ2init,ξ1'[t0]==ξ1dinit,ξ2'[t0]==ξ2dinit}],{ξ1[t],ξ2[t]},{t,t0,t0+tmax},
Method->{"EventLocator","Event"->(ξ1[t])(ξ1[t]-ξ2[t])==0,
"EventAction":>Throw[tend=t,"StopIntegration"]},
MaxStepSize->10^-3
];
sol=ndsol[[1]];
pl=Plot[{ξ1[t],ξ2[t]}/.sol,{t,t0,tend}];
ξ1end=(ξ1[t]/.sol)/.t->tend;
ξ2end=(ξ2[t]/.sol)/.t->tend;
ξ1dend=(D[ξ1[t]/.sol,t])/.t->tend;
ξ2dend=(D[ξ2[t]/.sol,t])/.t->tend;
(* Trk s podlago *)
If[Abs[ξ1end]<ϵ,ξ1dend=-ξ1dend];
(* Trk med kroglicama: elasticni odboj *)
If[Abs[ξ1end-ξ2end]<ϵ,
{ξ1dend,ξ2dend}={
((x-1)ξ1dend+2 ξ2dend)/(1+x),
((1-x)ξ2dend+2x ξ1dend)/(1+x)
}];
];

(* go[] racuna celotno trajektorijo za vec (nr) zaporednih dogodkov. *)
(* Parametri so C0 (brezdimenzijska moc vzmeti), x (razmerje mas) in v0 (zacetna hitrost v b.d. enotah). *)
go[C0_,x_,v0_,nr_]:=Module[{},
(* Zacetni pogoji: zacnemo v trenutku, po utez 1 prvic trci ob tla. *)
tend=0;
ξ1end=0;
ξ2end=1;
ξ1dend=v0;
ξ2dend=-v0;
plots=Table[calc[C0,x,tend,ξ1end,ξ2end,ξ1dend,ξ2dend];pl,{nr}];
Show[plots,PlotRange->{{0,6},All}]
];

Manipulate[go[cc,x,V0,20],
{{cc,1,"Konstanta vzmeti"},0,5},
{{x,1,"Masa 1/masa 2"},0,20},
{{V0,1,"Zacetna hitrost"},0,5}]