Binary Collisions

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

In this notebook we study some of the properties of nonrelativistic binary collisions.  Inelasticity and mass transfer are included.  Tools are developed to display the vector relationships and to produce both momentum and velocity diagrams.  Kinematic relationships are also plotted.

The reaction is described by the mass relationship (Subscript[m, 1]) + (Subscript[m, 2]) == (Subscript[m, 3]) + (Subscript[m, 4]), where (Subscript[m, 1]) is designated the projectile, (Subscript[m, 2]) the target, (Subscript[m, 3]) the ejectile, and (Subscript[m, 4]) the recoil particle.  If no mass is transferred between the projectile and the target, then (Subscript[m, 3]) -> (Subscript[m, 1]) and (Subscript[m, 4]) -> (Subscript[m, 2]).   The particle velocities are designated \(TraditionalForm\`{\(v\& ⇀ \)\_1, \(v\& ⇀ \)\_2, \(v\& ⇀ \)\_3, \(v\& ⇀ \)\_4}\).  For convenience we define the lab frame to be the rest frame of the target, such that (Subscript[v, 2]) -> 0 in the lab frame.  Other frames can be developed easily using Galilean transformations.  Conservation of momentum is then represented by the equation \(TraditionalForm\`\(p\& ⇀ \)\_1 + \(p\& ⇀ \)\_2 == \(p\& ⇀ \)\_3 + \(p\& ⇀ \)\_4\) and conservation of energy is represented by

\(TraditionalForm\`\(\(p\& ⇀ \)\_1·\(p\& ⇀ \)\_1\)\/\(2  m\_1\)\  + \ \(\(p\& ⇀ \)\_2·\(p\& ⇀ \)\_2\)\/\(2  m\_2\)\  + \ 
        Q\  = \ \(\(p\& ⇀ \)\_3·\(p\& ⇀ \)\_3\)\/\(2  m\_3\)\  + \ \(\(p\& ⇀ \)\_4·\(p\& ⇀ \)\_4\)\/\(2  m\_4\)\)

where the Q-value, Q, represents the internal energy liberated by the reaction.  Thus, the Q-value is positive for an exothermic reaction in which the final kinetic energy is greater than the initial kinetic energy , whereas  the Q-value is negative for an endothermic reaction in which some of the initial kinetic energy is converted into internal energy.  Elastic scattering requires the Q-value to vanish and for there to be no mass transfer.

Initialization

Options and packages

ClearAll["Global`*"];
Off[General::spell,  General::spell1];
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

Needs["Graphics`PlotField`"];
Needs["Graphics`Arrow`"]

Definitions

Establish symbolic names for the masses, momenta, and velocities that describe the  reaction.  These names can be used in replacement rules but should not be set to specific values.

masses = {m1, m2, m3, m4};
momenta = {p1, p2, p3, p4};
velocities = {v1, v2, v3, v4};

Implementation

Procedures for drawing vectors

describeVector[vector_?VectorQ, position_: {0, 0}, 
        style_: {HeadScaling -> Absolute, HeadLength -> 5.0, 
            HeadWidth -> 0.5}] := {vector, position, style}

drawVector[VectorStuff_] := 
      Module[{vector, position, style}, 
        vector = VectorStuff[[1]]; 
        position = VectorStuff[[2]]; 
        style = VectorStuff[[3]] /. 
            List -> Sequence; 
        Graphics[Arrow[position, position + vector, style]]]

drawVectors[VectorList_] := (drawVector[#1] &) /@ 
        VectorList

drawVectorSum[v1_?VectorQ, v2_?VectorQ] := 
        Module[{position1, position2, position3, style1, style2, style3, v3, 
            vector1, vector2, vector3}, position1 = position3 = {0, 0}; 
          position2 = v1; v3 = v1 + v2; 
          style1 = style2 = style3 = {HeadScaling -> Absolute, 
                  HeadLength -> 5.0, HeadWidth -> 0.5}; 
          vector1 = describeVector[v1, position1, style1]; 
          vector2 = describeVector[v2, position2, style2]; 
          vector3 = describeVector[v3, position3, style3]; 
          Show[{drawVectors[{vector1, vector2}], 
              Graphics[AbsoluteThickness[1]], drawVector[vector3]}, 
            AspectRatio -> Automatic]];

drawVectorDifference[v1_?VectorQ, v2_?VectorQ] := 
        Module[{position1, position2, position3, style1, style2, style3, v3, 
            vector1, vector2, vector3}, position1 = position2 = {0, 0}; 
          position3 = v1; v3 = v2 - v1; 
          style1 = style2 = style3 = {HeadScaling -> Absolute, 
                  HeadLength -> 5.0, HeadWidth -> 0.5}; 
          vector1 = describeVector[v1, position1, style1]; 
          vector2 = describeVector[v2, position2, style2]; 
          vector3 = describeVector[v3, position3, style3]; 
          Show[{drawVectors[{vector1, vector2}], 
              Graphics[AbsoluteThickness[1]], drawVector[vector3]}, 
            AspectRatio -> Automatic]];

drawCollision[p1i_?VectorQ, p2i_?VectorQ, p1f_?VectorQ, 
          p2f_?VectorQ] := 
        Module[{vector1i, vector2i, vector1f, vector2f}, 
          vector1i = describeVector[p1i, -p1i]; 
          vector2i = describeVector[p2i, -p2i]; 
          vector1f = describeVector[p1f]; vector2f = describeVector[p2f]; 
          Show[drawVectors[{vector1i, vector2i, vector1f, vector2f}], 
            AspectRatio -> Automatic]];

Procedures for calculating and drawing collisions

Establish equations expressing conservation of mass, momentum and energy

collisionEqs = {m1 + m2 == m3 + m4, p1 + p2 == p3 + p4, 
          p1 . (p1/(2 m1)) + p2 . (p2/(2 m2)) + qvalue == 
            p3 . (p3/(2 m3)) + p4 . (p4/(2 m4))};

momentumComponentRules = {p1 -> {p1x, p1y}, 
          p2 -> {p2x, p2y}, p3 -> {p3x, p3y}, 
          p4 -> {p4x, p4y}, p1x -> p1r Cos[theta1], 
          p1y -> p1r Sin[theta1], p2x -> p2r Cos[theta2], 
          p2y -> p2r Sin[theta2], p3x -> p3r Cos[theta3], 
          p3y -> p3r Sin[theta3]};

magnitude[p_?VectorQ] := (Sqrt[p . p]);
angle[p_?VectorQ] := 
      ArcTan[p[[2]], 
        p[[1]]]

cmVelocity[masses_, 
        momenta_] := ((momenta[[1]] + momenta[[2]])/(masses[[1]] + masses[[2]]))

LabToCM[masses_, momenta_] := 
        Module[{Vcm, velocities}, Vcm = cmVelocity[masses, momenta]; 
          velocities = (momenta/masses) - Vcm; 
          Simplify[{momenta -> masses velocities}]];

Define general solution and procedures for plotting momenta and angles

generalSolution[p3r, p4x, p4y] = 
        Solve[collisionEqs //. momentumComponentRules, {m4, p3r, p4x, 
            p4y}];

\(\(\(finalMomenta::usage\  = \ "finalMomenta[angle,sol,values] calculates the ejectile and recoil momenta given the scattering angle (in degrees) and the solution and initialization rules, sol and values.";\)\(
\)
    \)\)
finalMomenta[angle_, sol_, 
        values_] := {magnitude[p3], magnitude[p4]} //. 
              momentumComponentRules /. sol //. 
          values /. {theta3 -> angle Degree}

\(\(\(plotFinalMomenta::usage\  = \ "plotFinalMomenta[sol,values,thmax] plots the ejectile and recoil momenta for angles up to thmax (in degrees) given solution and initialization rules, sol and values.";\)\(
\)
    \)\)
\(\(plotFinalMomenta[sol_, values_, thmax_] := 
        Plot[Evaluate[finalMomenta[theta3, sol, values]], {theta3, 0, thmax}, 
          Frame -> True, GridLines -> Automatic, 
          FrameLabel -> {\*""(Subscript[θ, 3]) [deg]"", \*""(Subscript[p, 3])/(Subscript[p, 1]), (Subscript[p, 4])/(Subscript[p, 1])""}, 
          PlotLabel -> "scattered and recoil momenta"];\)\)

\(\(\(recoilAngle::usage\  = \ "recoilAngle[angle,sol,values] calculates the recoil angle (in degrees) given a scattering angle (in degrees) and the replacement rules that define a solution and the collision parameters.";\)\(
\)
    \)\)
recoilAngle[angle_, sol_, 
        values_] := ArcTan[p4x, -p4y]/Degree /. sol //. 
          values /. {theta3 -> angle Degree}

\(\(\(plotRecoilAngle::usage\  = \ "plotRecoilAngle[sol,values,thmax]  plots the recoil angle (in degrees) given solution and initialization rules, sol and values.  The maximum scattering angle plotted is thmax (in degrees).";\)\(
\)
    \)\)
\(\(plotRecoilAngle[sol_, values_, thmax_] := 
        Plot[Evaluate[recoilAngle[theta3, sol, values]], {theta3, 0.0001, 
            thmax}, Frame -> True, GridLines -> Automatic, 
          FrameLabel -> {\*""(Subscript[θ, 3]) [deg]"", \*""(Subscript[θ, 4]) [deg]""}, 
          PlotLabel -> "recoil angle"];\)\)

\(\(\(calcCollision::usage\  = \ "calcCollision[sol,values,LabAngle]  evaluates many quantities describing a collision given the solution replacement rules, sol, the list of rules defining the masses and initial conditions, values, the laboratory angle for particle 3, LabAngle.  The following quantities are calculated using this information: 
 LabMomenta = list of laboratory momenta 
 LabVelocities = list of laboratory velocities 
 cmMomenta = list of cm momenta 
 cmVelocities = list of cm velocities 
 mui,muf = initial, final reduced masses 
 vCM, pCM = velocity, momentum of CM 
 vi,vf = initial, final relative velocities 
 pi,pf = initial, final relative momenta in lab 
 ki,kf = initial, final relative momenta in cm ";\)\(
\)
    \)\)
calcCollision[sol_, values_, LabAngle_] := 
        Module[{}, 
          LabMomenta = momenta //. momentumComponentRules /. sol //. 
                values /. theta3 -> LabAngle; 
          LabVelocities = (LabMomenta/masses) /. sol //. values; 
          cmMomenta = momenta /. LabToCM[masses, momenta] //. 
                    momentumComponentRules /. sol //. values /. 
              theta3 -> LabAngle; 
          cmVelocities = (cmMomenta/masses) /. sol //. values; {mui, muf, 
              vi, vf, pi, pf, vCM, 
              pCM} = {((masses[[1]] masses[[2]])/(masses[[1]] + masses[[2]])), ((masses[[3]] masses[[4]])/(masses[[3]] + masses[[4]])), 
                  LabVelocities[[1]] - 
                    LabVelocities[[2]], 
                  LabVelocities[[3]] - 
                    LabVelocities[[4]], 
                  LabMomenta[[1]] - 
                    LabMomenta[[2]], 
                  LabMomenta[[3]] - 
                    LabMomenta[[4]], ((LabMomenta[[1]] + LabMomenta[[2]])/(masses[[1]] + masses[[2]])), 
                  LabMomenta[[1]] + 
                    LabMomenta[[2]]} /. 
                sol //. values; ki = mui vi; kf = muf vf;];

Procedures for drawing vector relationships

\(\(\(PlotCMmomenta::usage\  = \ "PlotCMmomenta[cmMomenta]  produces a plot which displays the initial and final momenta in the cm frame together with circles of constant momentum.";\)\(
\)
    \)\)
PlotCMmomenta[cmMomenta_] := 
        Module[{pi, pf, vectorList}, 
          pi = (1/2) (cmMomenta[[1]] - 
                  cmMomenta[[2]]); 
          pf = (1/2) (cmMomenta[[3]] - 
                  cmMomenta[[4]]); 
          vectorList = (describeVector[#1] &) /@ cmMomenta; {drawVectors[
              vectorList], 
            Graphics[{Circle[{0, 0}, magnitude[pi]], 
                Circle[{0, 0}, magnitude[pf]], AbsolutePointSize[3], 
                Point[{0, 0}]}]}];

\(\(\(VelocityPlot2::usage\  = \ "VelocityPlot2[LabVelocities,vCM,cmVelocities] produces two plots, named plot1 and plot2, which display velocity relationships in the initial and final states, respectively.  Each plot contains two lab velocity vectors and the cm velocity vector based at the origin, and the two velocity vectors in the cm frame based at the cm velocity.";\)\(
\)
    \)\)
VelocityPlot2[LabVelocities_, vCM_, cmVelocities_] := 
        Module[{initialVelocities, finalVelocities}, 
          initialVelocities = 
            Join[(describeVector[#1] &) /@ {LabVelocities[[1]], 
                  LabVelocities[[2]], 
                  vCM}, {describeVector[
                  cmVelocities[[1]], 
                  vCM], describeVector[
                  cmVelocities[[2]], 
                  vCM]}]; 
          finalVelocities = 
            Join[(describeVector[#1] &) /@ {LabVelocities[[3]], 
                  LabVelocities[[4]], 
                  vCM}, {describeVector[
                  cmVelocities[[3]], 
                  vCM], describeVector[
                  cmVelocities[[4]], 
                  vCM]}]; plot1 = drawVectors[initialVelocities]; 
          plot2 = drawVectors[finalVelocities];];

\(\(\(MomentumPlot2::usage\  = \ "MomentumPlot2[LabMomenta] produces two figures: plot1 shows the sum of initial momentum vectors and plot2 shows the sum of final momentum vectors.";\)\(
\)
    \)\)
MomentumPlot2[LabMomenta_] := 
        Module[{pCM}, 
          pCM = LabMomenta[[1]] + 
              LabMomenta[[2]]; 
          initialMomenta = {describeVector[
                LabMomenta[[1]]], 
              describeVector[
                LabMomenta[[2]], 
                LabMomenta[[1]]], 
              describeVector[pCM], 
              describeVector[
                cmMomenta[[1]], pCM], 
              describeVector[
                cmMomenta[[2]], pCM]}; 
          finalMomenta = {describeVector[
                LabMomenta[[3]]], 
              describeVector[
                LabMomenta[[4]], 
                LabMomenta[[3]]], 
              describeVector[pCM], 
              describeVector[
                cmMomenta[[3]], pCM], 
              describeVector[
                cmMomenta[[4]], pCM]}; 
          plot1 = drawVectors[initialMomenta]; 
          plot2 = drawVectors[finalMomenta];];

Simplification rules

MyAssumptions -> {m1 > 0, m2 > 0, m3 > 0, m4 > 0, mtotal > 0, 
          mui > 0, muf > 0, qvalue ∈ Reals, vi ∈ Reals, 
          vf ∈ Reals};
MySimplify = 
        Simplify[#, MyAssumptions] &;
MyFullSimplify = FullSimplify[#, MyAssumptions] &;

Rules for lab to cm transformations

Define reduced masses, relative and cm velocities

Clear[m1, m2, m3, m4, v1, v2, v2, v4, mui, muf, vi, vf, vCM, 
        pCM];

reducedMassRules = {mui -> ((m1 m2)/(m1 + m2)), 
          muf -> ((m3 m4)/(m3 + m4)), 
          mtotal -> m1 + m2};

vcmRules = {v1cm -> v1 - ((m1 v1 + m2 v2)/(m1 + m2)), 
            v2cm -> v2 - ((m1 v1 + m2 v2)/(m1 + m2)), 
            v3cm -> v3 - ((m3 v3 + m4 v4)/(m3 + m4)), 
            v4cm -> 
              v4 - ((m3 v3 + m4 v4)/(m3 + m4))} /. {v1 -> vi + v2, 
            v3 -> vf + v4}  // Simplify

vLabToCM = 
      Join[Simplify[Solve[{vi == v1 - v2, 
                vCM == ((m1 v1 + m2 v2)/(m1 + m2))}, {v1, 
                v2}][[1]]], 
        Simplify[Solve[{vf == v3 - v4, 
                vCM == ((m3 v3 + m4 v4)/(m3 + m4))}, {v3, 
                v4}][[1]]]]

masses velocities /. vLabToCM   // 
        Collect[#, {vi, vf, vcm}] &  // Simplify

massRules = Solve[{mtotal == m1 + m2 == m3 + m4, 
            mui == ((m1 m2)/(m1 + m2)), 
            muf == ((m3 m4)/(m3 + m4))}, {m1, m2, m3, 
            m4}][[1]]

Express cm momenta in terms of relative velocities and reduced masses

stuff = 
      Eliminate[{p1cm == m1 v1cm, p2cm == m2 v2cm, 
            p3cm == m3 v3cm, p4cm == m4 v4cm, 
            mtotal == m1 + m2 == m3 + m4, 
            mui == ((m1 m2)/mtotal), 
            muf == ((m3 m4)/mtotal)} /. vcmRules, {m1, m2, m3, m4, 
          v1, v2, v3, v4}]

pcmRules = Solve[
          stuff, {p1cm, p2cm, p3cm, 
            p4cm}][[1]]

Separate kinetic energy into internal and cm contributions

Solve[
        Eliminate[{Simplify[
              KiLab == (Collect[#1, {vi, 
                          vCM}] &)[((m1 (v1^2))/2) + ((m2 (v2^2))/2) /. 
                    vLabToCM]], mtotal == m1 + m2, 
            mui == ((m1 m2)/mtotal)}, {m1, m2}], 
        KiLab][[1]]

The first term is the kinetic energy associated with the motion of the center of mass,  whereas the second is the kinetic energy due to the internal motions relative to the center of mass.  The corresponding relationship for the final state is obviously obtained by the replacements mui -> muf, vi -> vf.  Therefore, we can define the following rules.

kineticEnergyRules = {KiLab -> KCM + KiCM, 
          KfLab -> KCM + KfCM, KCM -> ((mtotal (vCM^2))/2), 
          KiCM -> ((mui (vi^2))/2), 
          KfCM -> ((muf (vf^2))/2)};

Conservation of energy can now be analyzed using the relative velocities and the qvalue.  For example, we can solve for the relative velocity in the final state in terms of the initial relative velocity, the qvalue, and the reduced masses.  Furthermore, because the momenta of the outgoing particles must be equal and opposite in the center of mass or zero momentum frame, the cm velocities of each particle are inversely proportional to their masses.  Note that since these results are independent of scattering angle, we can develop simple relationships for the magnitudes of the cm vectors as follows.

vfRule = Solve[
            KfLab - KiLab - qvalue == 0 //. kineticEnergyRules, 
            vf][[2]] // 
        Simplify

{p3cm, p4cm} /. pcmRules /. vfRule

Simplify[(({v3cm, v4cm} (m3 + m4))/mtotal) /. vcmRules /. 
        vfRule]

Examples

Drawing vectors

Show[Graphics[AbsoluteThickness[1]], 
        drawVector[describeVector[{1, 3}]]] ;

drawVectorSum[{1, 2}, {3, -5}];

drawVectorDifference[{-10, 2}, {3, 5}];

drawCollision[{2, 0}, {-1, 0}, {1, 
          1}, {-2, -1}];

vectorList = 
      Join[(describeVector[#1] &) /@ {{1, 1}, {-1, 
              2}}, {describeVector[{2, 3}, {-2, -5}]}]

Show[drawVectors[vectorList]];

One dimensional collisions

General one-dimensional inelastic collision with mass transfer

It is easier to use scalar equations for one-dimensional collisions than to decide whether the scattering angle should be 0 or π.

collisionEqs1D = collisionEqs /. {Dot -> Times}

sol =  {p1, p2, p3, p4} /. 
          Solve[collisionEqs1D, {p3, p4, 
              m4}]//MySimplify

To select the appropriate solution, consider an elastic collision between equal masses.

sol /. {m3 -> m1, m4 -> m2, p2 -> 0, 
              qvalue -> 0}  // Simplify // PowerExpand // 
      Simplify

Recognizing that the first solution corresponds to no scattering (p3 -> p1) in this case, we choose the second solution.

momenta1D = 
        sol[[2]];

velocities1D  = (momenta1D/masses) /. {p1 -> m1 v1, 
                p2 -> m2 v2} // Simplify // PowerExpand // 
        Simplify

One-dimensional elastic collisions

momenta1D /. {m3 -> m1, m4 -> m2, qvalue -> 0} // 
          Simplify // PowerExpand // Simplify

LabVelocities = velocities1D /. {m3 -> m1, m4 -> m2, 
                  qvalue -> 0} // Simplify // PowerExpand // 
          Simplify // Collect[#, {v1, v2}] &

It is useful to examine a few limiting cases.

LabVelocities /. {m2 -> m1}

Thus, we find that when the masses are equal a one-dimensional elastic collision simply exchanges the velocities of the two particles.  For example, if (Subscript[m, 2]) is initially at rest, after the collision it acquires the original velocity of (Subscript[m, 1]) and (Subscript[m, 1]) is brought to rest.

Limit[LabVelocities, m2 -> ∞]

On the other hand, if the target is extremely massive its velocity is unaffected by the collision whereas the velocity of the projectile is reversed and is increased by twice the target velocity.  For example, suppose that (Subscript[m, 1]) represents a baseball and that (Subscript[m, 2]) represents the combined mass of the bat and batter and is much larger than (Subscript[m, 1]).  The speed of a line drive is then the speed of the pitch plus twice the bat speed.

Limit[LabVelocities, m1 -> ∞]

Alternatively, if the projectile is much more massive than the target, the recoiling target acquires twice the projectile velocity minus its original velocity.  For example, if (Subscript[m, 1]) is a large truck which strikes a small car, (Subscript[m, 2]), originally at rest at a stop light, the car recoils with twice the truck's velocity while the truck continues with its original velocity (assuming that the collision is approximately elastic).

One-dimensional inelastic collision without mass transfer

LabMomenta = momenta1D /. {m3 -> m1, m4 -> m2} // 
            Simplify // PowerExpand // Simplify

LabVelocities = velocities1D /. {m3 -> m1, m4 -> m2} // 
            Simplify // PowerExpand // Simplify

LabVelocities /. {v2 -> 0}  // Simplify

Elastic colinear collision for equal masses

values = {m1 -> 1, m2 -> m1, m3 -> m1, 
          qvalue -> 0, p1r -> 1, theta1 -> 0, p2r -> 0, 
          theta2 -> π};

sol = 
      Simplify[PowerExpand[generalSolution[p3r, p4x, 
                p4y][[1]] //. 
            values]]

calcCollision[sol, values, 60 Degree]

LabMomenta

LabMomenta[[3]] . 
      LabMomenta[[4]]

drawCollision @@ LabMomenta;

drawVectorSum[LabMomenta[[3]], 
        LabMomenta[[4]]];

Simplify[p3 . p4 //. momentumComponentRules /. sol //. 
        values]

VelocityPlot2[LabVelocities, vCM, cmVelocities];

Show[Graphics[RGBColor[0, 1, 0]], plot1, 
        Graphics[Circle[vCM, magnitude[vi]/2]], Graphics[RGBColor[1, 0, 0]], 
        plot2, Graphics[Circle[vCM, magnitude[vf]/2]], 
        AspectRatio -> Automatic];

It should be clear from this plot that the laboratory scattering angle is limited to 90  °, where the struck particle carries the full momentum and the projectile ends up at rest.

plotFinalMomenta[sol, values, 90. ];

plotRecoilAngle[sol, values, 90. ];

PlotCMmomenta[cmMomenta];

Show[%, AspectRatio -> Automatic];

Colinear inelastic collision between equal masses

Inelastic collisions can be classified as endothermic when kinetic energy is transferred into internal energy or as exothermic when internal energy is liberated and becomes kinetic energy.  Thus, endothermic reactions have negative qvalue and exothermic reactions positive qvalue.  For simplicity, we consider inelastic collisions between equal masses and examine the kinematical relationships between the laboratory frame, defined as the rest frame of particle 2, and the center-of-mass frame.

Endothermic reaction

values = {m1 -> 1, m2 -> 1, m3 -> 1, 
          qvalue -> -0.2, p1r -> 1, theta1 -> 0, 
          p2r -> 0, theta2 -> π};

sol = generalSolution[p3r, p4x, 
              p4y][[1]] //. 
          values;

calcCollision[sol, values, 20 Degree]

drawCollision @@ LabMomenta;

For an endothermic reaction, the lab angle between the two emerging particles (with equal masses) is now less than 90  °.  (Recall that this angle is 90  ° for an elastic collision between equal masses.)

PlotCMmomenta[cmMomenta];

Show[%, AspectRatio -> Automatic];

VelocityPlot2[LabVelocities, vCM, cmVelocities];

Show[Graphics[RGBColor[0, 1, 0]], plot1, 
        Graphics[Circle[vCM, magnitude[vi]/2]], Graphics[RGBColor[1, 0, 0]], 
        plot2, Graphics[Circle[vCM, magnitude[vf]/2]], 
        AspectRatio -> Automatic];

Because the outgoing velocities must lie on the inner circle, the velocity plot  demonstrates that there is are maximum scattering and recoil angles for this  endothermic reaction.  The maximum scattering angle is derived in a later section using conservation of energy and momentum, without assuming equal masses.

Exothermic reaction

values = {m1 -> 1, m2 -> 1, m3 -> 1, 
          qvalue -> +0.2, p1r -> 1, theta1 -> 0, 
          p2r -> 0, theta2 -> π};

sol = generalSolution[p3r, p4x, 
              p4y][[1]] //. 
          values;

calcCollision[sol, values, 20 Degree]

drawCollision @@ LabMomenta;

For an exothermic reaction, the lab angle between the two emerging particles (with equal masses) is now greater than 90  °.  (Recall that this angle is 90  ° for an elastic collision between equal masses.)

PlotCMmomenta[cmMomenta];

Show[%, AspectRatio -> Automatic];

VelocityPlot2[LabVelocities, vCM, cmVelocities];

Show[Graphics[RGBColor[0, 1, 0]], plot1, 
        Graphics[Circle[vCM, magnitude[vi]/2]], Graphics[RGBColor[1, 0, 0]], 
        plot2, Graphics[Circle[vCM, magnitude[vf]/2]], 
        AspectRatio -> Automatic];

Since the outgoing velocities may appear anywhere on the outer circle, the scattering and recoil angles for this exothermic reaction are not restricted.

Elastic colinear collision for unequal masses

We now consider an elastic collision between unequal masses, such that m2/m1 = mratio.  For simplicity we assume that the target is at rest in the lab.

values = {m1 -> 1, m2 -> mratio m1, m3 -> m1, 
          qvalue -> 0, p1r -> 1, theta1 -> 0, p2r -> 0, 
          theta2 -> π};

sol  =  generalSolution[p3r, p4x, 
                  p4y][[1]] //. 
              values  // Simplify // PowerExpand // Simplify

Projectile heavier than target

Suppose that the projectile is heavier than the target, such that mratio<1.  We then find that there is a maximum scattering angle, which is derived below.  Plots are made for the specific case mratio->0.5.

calcCollision[sol, Join[values, {mratio -> 0.5}], 
      20 Degree]

drawCollision @@ LabMomenta;

VelocityPlot2[LabVelocities, vCM, cmVelocities];

Show[Graphics[RGBColor[0, 1, 0]], plot1, Graphics[RGBColor[1, 0, 0]], 
        plot2, AspectRatio -> Automatic];

PlotCMmomenta[cmMomenta];

Show[%, AspectRatio -> Automatic];

Simplify[p3r //. momentumComponentRules /. sol //. 
        values]

Notice that if mratio<1, the range of theta3 must be restricted.

maxTheta[mratio_] := 
      (theta3/Degree) /. Solve[-1 + 2 (mratio^2) + Cos[2 theta3] == 0, 
            theta3][[2]]

thmax = N[maxTheta[0.5]]

plotFinalMomenta[sol, Join[values, {mratio -> 0.5}], 
        0.9999 thmax];

plotRecoilAngle[sol, Join[values, {mratio -> 0.5}], 
        0.9999 thmax];

Target heavier than projectile

Suppose that the projectile is heavier than the target, such that mratio>1.  The scattering angle is then unrestricted.  We examine below the specific case of mratio->2.  There is a maximum scattering angle, which is derived below.

plotFinalMomenta[sol, Join[values, {mratio -> 2.0}], 
        180. ];

plotRecoilAngle[sol, Join[values, {mratio -> 2.0}], 
        180. ];

Maximum scattering angle

We have encountered several situations in which there is a maximum possible scattering angle.  For simplicity we assume that the target is at rest in the lab.  Then for elastic collisions there is a maximum scattering angle when the projectile is heavier than the target.  Alternatively, for endothermic reactions there is a maximum scattering angle even when the projectile mass is less than or equal to the target mass.   To obtain a general symbolic relationship for the maximum scattering angle, we apply conservation of momentum and energy to the triangle formed by the incident momentum and the two exit momenta, as  shown below.

drawVectorSum[LabMomenta[[3]], 
        LabMomenta[[4]]];

Use conservation of energy and apply the law of cosines, where x = Cos[theta3].  We eliminate p4 and express p3 in units of p1 and the qvalue as a fraction of the initial kinetic energy, such that qvalue->-inelasticity (p1^2)/(2 m1).  Similarly, we define mratio->m2/m1 as the ratio of masses.

eq = 
      Simplify[Eliminate[{((p1^2)/(2 m1)) + qvalue == 
                    ((p3^2)/(2 m1)) + ((p4^2)/(2 m2)), 
                  (p4^2) == (p1^2) + (p3^2) - 2 p1 p3 x}, 
                p4][[1]] /. {qvalue -> -((inelasticity (p1^2))/(2 m1)), 
              m2 -> m1 mratio} /. {p1 -> 1, 
            m1 -> 1}]

xsol = 
      x /. Solve[eq, 
            x][[1]]

Now we maximize theta3 by minimizing x wrt p3.

p3sol = 
      Simplify[Solve[(∂_p3) xsol == 0, 
            p3][[1]]]

cos2 = (xsol^2) /. p3sol  // Simplify  // 
        Collect[#, inelasticity] &

Finally, recognizing that the expression above gives Cos[(Subscript[θ, 3])]^2, we construct a maximum angle function that recognizes that there are situations in which the scattering angle is not restricted.  This function depends upon the mass ratio and the relative inelasticity (negative is endothermic, positive is exothermic) and is displayed in two forms below.

Clear[thetaMax];
thetaMax[m_, x_] := 
      Module[{temp}, 
        temp = cos2 /. {mratio -> m, inelasticity -> x}; 
        If[temp > 0 && temp < 1, ArcCos[(Sqrt[temp])], π]]

\(\(plotThetaMax3D = 
        Plot3D[Evaluate[thetaMax[mratio, inelasticity]\/Degree], {mratio, 0, 
            1}, {inelasticity, 0, 0.5}, PlotPoints -> 25, 
          PlotRange -> {Automatic, Automatic, {0, 180}}, 
          AxesLabel -> {"mass ratio", "relative inelasticity", \*""Max[(Subscript[θ, 3])] (deg)""}, 
          PlotLabel -> "Maximum Scattering Angle", 
          DisplayFunction -> Identity];\)\)

Show[plotThetaMax3D, ViewPoint -> {3, 2, 1}, 
        BoxRatios -> {1, 1, 1}, PlotLabel -> {}, 
        DisplayFunction -> $DisplayFunction];

\(\(plotThetaMax = 
        Plot[Evaluate[
            Table[thetaMax[mratio, inelasticity]\/Degree, {mratio, 0, 2, 
                0.2}]], {inelasticity, \(-1\), 1}, 
          PlotRange -> {Automatic, {0, 180}}, PlotPoints -> 50, 
          Frame -> True, GridLines -> Automatic, 
          FrameLabel -> {"relative inelasticity", \*""Max[(Subscript[θ, 3])] (deg)""}, 
          PlotLabel -> "Maximum Scattering Angle"];\)\)

Oblique elastic collision between unequal masses

values = {m1 -> 2, m2 -> 1, m3 -> 2, qvalue -> 0, 
          p1r -> 1, theta1 -> 30 Degree, p2r -> 3, 
          theta2 -> -50 Degree};

sol = generalSolution[p3r, p4x, 
              p4y][[1]] //. 
          values;

calcCollision[sol, values, 20 Degree]

drawCollision @@ LabMomenta;

VelocityPlot2[LabVelocities, vCM, cmVelocities];

Show[Graphics[RGBColor[0, 1, 0]], plot1, Graphics[RGBColor[1, 0, 0]], 
        plot2, AspectRatio -> Automatic];

PlotCMmomenta[cmMomenta];

Show[%, AspectRatio -> Automatic];

MomentumPlot2[LabMomenta];

Show[Graphics[RGBColor[0, 1, 0]], plot1, Graphics[RGBColor[1, 0, 0]], 
        plot2, AspectRatio -> Automatic];