8 #ifndef WARPX_TWO_PRODUCT_FUSION_UTIL_H
9 #define WARPX_TWO_PRODUCT_FUSION_UTIL_H
48 void TwoProductFusionComputeProductMomenta (
49 const amrex::ParticleReal& u1x_in,
50 const amrex::ParticleReal& u1y_in,
51 const amrex::ParticleReal& u1z_in,
52 const amrex::ParticleReal& m1_in,
53 const amrex::ParticleReal& u2x_in,
54 const amrex::ParticleReal& u2y_in,
55 const amrex::ParticleReal& u2z_in,
56 const amrex::ParticleReal& m2_in,
57 amrex::ParticleReal& u1x_out,
58 amrex::ParticleReal& u1y_out,
59 amrex::ParticleReal& u1z_out,
60 const amrex::ParticleReal& m1_out,
61 amrex::ParticleReal& u2x_out,
62 amrex::ParticleReal& u2y_out,
63 amrex::ParticleReal& u2z_out,
64 const amrex::ParticleReal& m2_out,
65 const amrex::ParticleReal& E_fusion,
68 using namespace amrex::literals;
72 constexpr amrex::ParticleReal inv_csq = 1._prt / ( c_sq );
74 const amrex::ParticleReal E_rest_in = (m1_in + m2_in)*c_sq;
76 const amrex::ParticleReal E_rest_out = (m1_out + m2_out)*c_sq;
79 const amrex::ParticleReal g1_in = std::sqrt( 1._prt
80 + (u1x_in*u1x_in + u1y_in*u1y_in + u1z_in*u1z_in)*inv_csq );
81 const amrex::ParticleReal g2_in = std::sqrt( 1._prt
82 + (u2x_in*u2x_in + u2y_in*u2y_in + u2z_in*u2z_in)*inv_csq );
85 const amrex::ParticleReal p1x_in = u1x_in * m1_in;
86 const amrex::ParticleReal p1y_in = u1y_in * m1_in;
87 const amrex::ParticleReal p1z_in = u1z_in * m1_in;
88 const amrex::ParticleReal p2x_in = u2x_in * m2_in;
89 const amrex::ParticleReal p2y_in = u2y_in * m2_in;
90 const amrex::ParticleReal p2z_in = u2z_in * m2_in;
92 const amrex::ParticleReal p_total_sq = powi<2>(p1x_in+p2x_in) +
93 powi<2>(p1y_in+p2y_in) +
94 powi<2>(p1z_in+p2z_in);
97 const amrex::ParticleReal E_lab = (m1_in * g1_in + m2_in * g2_in) * c_sq;
100 const amrex::ParticleReal E_star_sq = E_lab*E_lab - c_sq*p_total_sq;
107 const amrex::ParticleReal E_star_f_sq = powi<2>(std::sqrt(E_star_sq)
108 - E_rest_in + E_rest_out + E_fusion);
114 const amrex::ParticleReal E_ratio = std::sqrt(E_star_f_sq)/((m1_out + m2_out)*c_sq);
115 const amrex::ParticleReal p_star_f_sq = m1_out*m2_out*c_sq * ( powi<2>(E_ratio) - 1._prt )
116 + powi<2>(m1_out - m2_out)*c_sq*0.25_prt * powi<2>( E_ratio - 1._prt/E_ratio );
120 amrex::ParticleReal px_star, py_star, pz_star;
125 amrex::ParticleReal p1x_out, p1y_out, p1z_out;
127 const amrex::ParticleReal mass_g = m1_in * g1_in + m2_in * g2_in;
128 const amrex::ParticleReal vcx = (p1x_in+p2x_in) / mass_g;
129 const amrex::ParticleReal vcy = (p1y_in+p2y_in) / mass_g;
130 const amrex::ParticleReal vcz = (p1z_in+p2z_in) / mass_g;
131 const amrex::ParticleReal vc_sq = vcx*vcx + vcy*vcy + vcz*vcz;
135 if ( vc_sq > std::numeric_limits<amrex::ParticleReal>::min() )
137 const amrex::ParticleReal gc = 1._prt / std::sqrt( 1._prt - vc_sq*inv_csq );
138 const amrex::ParticleReal g_star = std::sqrt(1._prt + p_star_f_sq / (m1_out*m1_out*c_sq));
139 const amrex::ParticleReal vcDps = vcx*px_star + vcy*py_star + vcz*pz_star;
140 const amrex::ParticleReal factor0 = (gc-1._prt)/vc_sq;
141 const amrex::ParticleReal factor = factor0*vcDps + m1_out*g_star*gc;
142 p1x_out = px_star + vcx * factor;
143 p1y_out = py_star + vcy * factor;
144 p1z_out = pz_star + vcz * factor;
154 const amrex::ParticleReal p2x_out = p1x_in + p2x_in - p1x_out;
155 const amrex::ParticleReal p2y_out = p1y_in + p2y_in - p1y_out;
156 const amrex::ParticleReal p2z_out = p1z_in + p2z_in - p1z_out;
159 u1x_out = p1x_out/m1_out;
160 u1y_out = p1y_out/m1_out;
161 u1z_out = p1z_out/m1_out;
162 u2x_out = p2x_out/m2_out;
163 u2y_out = p2y_out/m2_out;
164 u2z_out = p2z_out/m2_out;
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_INLINE void RandomizeVelocity(amrex::ParticleReal &ux, amrex::ParticleReal &uy, amrex::ParticleReal &uz, const amrex::ParticleReal vp, amrex::RandomEngine const &engine)
Function to perform scattering of a particle that results in a random velocity vector with given magn...
Definition: ParticleUtils.H:162
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44