implemented u_com and tried making a ghost particle for each collisio…

…n and failed (core dumped)

include/openmc/particle.h

@@ -54,7 +54,7 @@ class Particle : public ParticleData {
   //! simply as a secondary particle.
   //! \param src Source site data
   void from_source(const SourceSite* src);
-
+  void initilze_ghost_particle(Particle& p,Direction u, double E);
   // Coarse-grained particle events
   void event_calculate_xs();
   void event_advance();

prints

@@ -0,0 +1,36 @@
+  fmt::print("post energy = {0} pre energy = {1}\n",p.E(),p.E_last());
+  fmt::print("event mt = {}\n",p.event_mt());
+  fmt::print("post weight = {0} pre weight = {1}\n",p.wgt(),p.wgt_last());
+  fmt::print("pos = {0} , {1} , {2}\n",p.r().x,p.r().y,p.r().z);
+  fmt::print("u = {0} , {1} , {2}\n",p.u().x,p.u().y,p.u().z);
+  fmt::print("u_last = {0} , {1} , {2}\n",p.u_last().x,p.u_last().y,p.u_last().z);
+  fmt::print("distance = {}\n",distance);
+  fmt::print("xs = {}\n",p.macro_xs().total);
+  fmt::print("number of nuclear = {}\n",n);
+  fmt::print("name of 1# nuclicde = {}\n",nuc->name_);
+  fmt::print("name of 2# nuclicde = {}\n",nuc2->name_);
+  fmt::print("reaction MT = {}\n",rx->mt_);
+  fmt::print("d type = {}\n",typeid(d.get()).name());
+  fmt::print("energy distribution at 0 {}\n",d_->angle().get_energy(0));
+  fmt::print("energy distribution at 1 {}\n",d_->angle().get_energy(1));
+  fmt::print("post energy = {}\n",p.E);
+  fmt::print(p.E);
+
+  // Neutron velocity in LAB
+  Direction v_n = vel * p.u_last();
+
+ // // Sample velocity of target nucleus
+ // Direction v_t {};
+ // if (!p.neutron_xs(i_nuclide).use_ptable) {
+ //   v_t = sample_target_velocity(*nuc, p.E_last(), p.u_last(), v_n,
+ //     p.neutron_xs(i_nuclide).elastic, kT, p.current_seed());
+ // }
+
+  // Velocity of center-of-mass
+  Direction v_cm = (v_n ) / (awr + 1.0);  //(v_n + awr * v_t) / (awr + 1.0);
+
+  // Transform to CM frame
+  v_n -= v_cm;
+
+  // Find speed of neutron in CM
+  vel = v_n.norm();

src/particle.cpp

@@ -60,6 +60,22 @@ double Particle::speed() const
   return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma);
 }
 
+void Particle::initilze_ghost_particle(Particle& p,Direction u_new, double E_new)
+{
+  type() = p.type();
+  wgt() = p.wgt_last();
+  wgt_last() = p.wgt_last();
+  r() = p.r();
+  r_last() = p.r();
+  u() = u_new;
+  u_last() = u_new;
+  //u_last() = p.u_last();
+  E() = E_new;  //settings::run_CE ? p.E_last() : p.g_last();
+  E_last() = E_new;
+  time() = p.time_last();
+  time_last() = p.time_last();
+}
+
 void Particle::create_secondary(
   double wgt, Direction u, double E, ParticleType type)
 {

src/tallies/tally_scoring.cpp

@@ -2367,44 +2367,69 @@ void score_collision_tally(Particle& p)
 
   double distance;
   distance = p.r().x*p.r().x+p.r().y*p.r().y+p.r().z*p.r().z;
-  //fmt::print("collison happened\n");
-  if(p.E_last()>2e5)
-  {
-  //fmt::print("post energy = {0} pre energy = {1}\n",p.E(),p.E_last());
-  //fmt::print("event mt = {}\n",p.event_mt());
-  //fmt::print("post weight = {0} pre weight = {1}\n",p.wgt(),p.wgt_last());
-  //fmt::print("pos = {0} , {1} , {2}\n",p.r().x,p.r().y,p.r().z);
-  //fmt::print("u = {0} , {1} , {2}\n",p.u().x,p.u().y,p.u().z);
-  //fmt::print("u_last = {0} , {1} , {2}\n",p.u_last().x,p.u_last().y,p.u_last().z);
-  //fmt::print("distance = {}\n",distance);
-  //fmt::print("xs = {}\n",p.macro_xs().total);
-
   const auto& mat {model::materials[p.material()]};
   int n = mat->nuclide_.size();
   const auto& nuc {data::nuclides[0]};
-  const auto& nuc2 {data::nuclides[1]};
-  //fmt::print("number of nuclear = {}\n",n);
-  //fmt::print("name of 1# nuclicde = {}\n",nuc->name_);
-  //fmt::print("name of 2# nuclicde = {}\n",nuc2->name_);
-  int i_nuclide=0;
+  double awr = nuc->awr_;
   const auto& rx {nuc->reactions_[0]};
-  //fmt::print("reaction MT = {}\n",rx->mt_);
+
   auto& d = rx->products_[0].distribution_[0];
   auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
-  //fmt::print("d type = {}\n",typeid(d.get()).name());
-  //fmt::print("energy distribution at 0 {}\n",d_->angle().get_energy(0));
-  //fmt::print("energy distribution at 1 {}\n",d_->angle().get_energy(1));
+
+  Direction u_lab {0-p.r().x,0-p.r().y,0-p.r().z};
+  u_lab = u_lab/u_lab.norm();
+  double cos_lab = u_lab.dot(p.u_last());
+  double theta_lab = std::acos(cos_lab);
+  double sin_lab = std::sin(theta_lab);
+  double mu_COM = (std::sqrt(awr*awr - sin_lab*sin_lab)*cos_lab - sin_lab*sin_lab)/awr;
+  double vel = std::sqrt(p.E_last());
+  double theta_pdf = d_->angle().get_pdf_value(p.E_last(),mu_COM,p.current_seed());
+  double E_ghost = p.E_last()*(1+awr*awr+2*awr*mu_COM)/(1+awr)/(1+awr);
+  Particle ghost_particle=Particle();
+  ghost_particle.initilze_ghost_particle(p,u_lab,E_ghost);
+  ghost_particle.event_calculate_xs();
+
+  //fmt::print("collison happened\n");
+  if(p.E_last()>2e5)
+  {
+  fmt::print("pos = {0} , {1} , {2}\n",p.r().x,p.r().y,p.r().z);
+  fmt::print("u_last = {0} , {1} , {2}\n",p.u_last().x,p.u_last().y,p.u_last().z);
+  fmt::print("u_lab = {0} , {1} , {2}\n",u_lab.x,u_lab.y,u_lab.z);
+
+
+  if(cos_lab > 1)
+  {
+  cos_lab = 1;
+  if(cos_lab>1.0001)
+  fmt::print("-----------------------error cos_lab = {}-----------------------\n",cos_lab);
   }
-  //fmt::print("post energy = {}\n",p.E);
-  //fmt::print(p.E);
+  fmt::print("cos_lab = {}\n",cos_lab);
+  fmt::print("awr = {}\n",awr);
+
+
+  fmt::print("mu_COM = {}\n",mu_COM);
+
+
+  fmt::print("theta pdf = {}\n",theta_pdf);
+
+
+  fmt::print("ghost_particle E = {}\n",ghost_particle.E());
+  fmt::print("ghost_particle u = {0} , {1} , {2}\n",ghost_particle.u().x,ghost_particle.u().y,ghost_particle.u().z);
+
+  fmt::print("ghost_particle XS = {}\n",ghost_particle.macro_xs().total);
+
+  }
+
+
   if (p.type() == ParticleType::neutron || p.type() == ParticleType::photon) {
     flux = 0;//p.wgt_last() / p.macro_xs().total;
     if(p.event_mt() == 2)
     {
-      flux = 0.5*p.wgt()*exp(-(distance-1)*p.macro_xs().total)/(2*3.14*distance*distance);
+      flux = 0.5*ghost_particle.wgt()*exp(-(distance-1)*ghost_particle.macro_xs().total)/(2*3.14*distance*distance);
     }
   }
 
+
   for (auto i_tally : model::active_collision_tallies) {
     const Tally& tally {*model::tallies[i_tally]};