restart: with(plots): with(Statistics): # draw a stellate pentagon at location (x, y), with size s and transparency tr star := (x, y, s, tr) -> PLOT( POLYGONS( evalf([seq([x+s*cos((2*i)*2*Pi/5+Pi/10), y+s*sin((2*i)*2*Pi/5+Pi/10)], i=1..5)]), COLOR(RGB, rand()/10^12, rand()/10^12, rand()/10^12), TRANSPARENCY(tr) ) ): # draw trajectories (x(t), y(t)) between the fired time and the explosion time Xballs := proc(theta) plots[display]( seq( piecewise( theta > Tstart||n+Tend||n and theta < Tstart||n+Tend||n+1, star(op(eval(rhs~(sol||n), t=Tstart||n+Tend||n)), 10*(Tstart||n+Tend||n+1-theta), theta-(Tstart||n+Tend||n)), star(op(eval(rhs~(sol||n), t=Tstart||n+Tend||n)), 0, 1) ), n=1..N ), seq( piecewise( theta > Tstart||n+1 and theta < Tstart||n+Tend||n, plot([op(eval(rhs~(sol||n), t=tau)), tau=theta-1..theta], colorscheme=CS||n, thickness=5), plot([op(eval(rhs~(sol||n), t=tau)), tau=Tstart||n..max(theta, Tstart||n)], color=black, thickness=1, transparency=1) ), n=1..N ), axes=none ) end proc: randomize(): # some parameters will have values randomly choosen LaunchingAngle := RandomVariable(Uniform(Pi/2-Pi/10, Pi/2+Pi/10)): InitialVelocity := RandomVariable(Uniform(40, 70)): StartTime := RandomVariable(Uniform(0, 0.25)): ExplosionTime := RandomVariable(Normal(4, 1)): # these color schemes are use in procedure Xballs ColorSchemes := [ ["Black", "Blue", "Cyan", "White"], ["Black", "Purple", "Red", "Magenta"], ["Black", "Red", "Magenta", "Pink"], ["Black", "Red", "Gold", "Yellow", "White"], ["Black", "Blue", "Green", "Yellow"] ]: roll := rand(1..numelems(ColorSchemes)): # system to solve and formal solution StartFrom := [0\$2]: g := 9.81: de1 := diff(x(t), t\$2) = 0: de2 := diff(y(t), t\$2) = -g: ic := x(T0)=StartFrom, y(T0)=StartFrom, D(x)(T0)=V0*cos(alpha), D(y)(T0)=V0*sin(alpha): sol := dsolve([de1, de2, ic], [x(t), y(t)]): # N fireworks for different initial conditions Tstart||0 := 0: N := 50: for n from 1 to N do p := n-1: data := [ alpha = Sample(LaunchingAngle, 1), V0 = Sample(InitialVelocity, 1), T0 = Tstart||p + Sample(StartTime, 1) ]; Tend||n := Sample(ExplosionTime, 1); sol||n := eval(sol, data); Tstart||n := eval(T0, data); CS||n := ColorSchemes[roll()]; end do: # to large a file to be loaded by the green arrow, set "false" to "true" if false then animate( Xballs, [t], t=Tstart||0..max(seq(Tstart||n+Tend||n, n=1..20))+5, background=display(PLOT(POLYGONS([[-100,0], [100,0], [100,250], [-100,250]]))), view=[-100..100, 0..250], size=[600, 750], scaling=constrained, frames=100 ); end if: