0

Hi there I have a C++ code and I need to write it in C to do the simulation. Can anyone help? Thank you so much. and for random.hpp I will replace that with random.c (a random number generator in c)

```
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <string>
#include <vector>
using namespace std;
#include "random.hpp"
using namespace cpl;
Random rg; // random number generator object
const double pi = 4 * atan(1.0); // value of pi
double a = 1, // inner radius of shield
b = 1.5, // outer radius of shield
tau = 0.2, // mean free time for v=1
p_a = 0.05, // absorption probability
f = 0.95; // inelasticity parameter
double v_0 = 1; // initial speed
// define * to compute dot product
double operator * (const vector<double>& v1, const vector<double>& v2) {
double sum = 0;
for (int i = 0; i < v1.size(); i++)
sum += v1[i] * v2[i];
return sum;
}
// define C++ Neutron object
struct Neutron {
vector<double> r, v; // position and velocity
bool absorbed, escaped; // final fate
// constructor initalizes neutron
Neutron() {
r = vector<double>(2, 0);
double theta = 2 * pi * rg.rand();
v.push_back(v_0 * cos(theta));
v.push_back(v_0 * sin(theta));
absorbed = escaped = false;
}
// attempt a Monte Carlo transport step
bool step_succeeded() {
if (absorbed || escaped)
return false;
double r_mag = sqrt(r*r);
if (r_mag > b) { // outside the shield
escaped = true;
return false;
}
if (r_mag < a) // inside core region
move_to_a();
move_free(); // inside the shield
interact();
return true;
}
// move to inner wall of shield
void move_to_a() {
double r2 = r * r, v2 = v * v, rv = r * v;
double t = sqrt((a*a - r2)/v2 + rv*rv/v2/v2) - rv/v2;
r[0] += t * v[0];
r[1] += t * v[1];
}
// free motion in shield
void move_free() {
double r2 = r * r, v2 = v * v, rv = r * v;
double t_max;
if (rv < 0 && abs(rv) >= v2 * (r2 - a*a))
t_max = abs(rv)/v2 - sqrt(rv*rv/v2/v2 - (r2 - a*a)/v2);
else
t_max = sqrt((b*b - r2)/v2 + rv*rv/v2/v2) - rv/v2;
double t = - tau * log(rg.rand());
if (t > t_max)
t = t_max;
r[0] += t * v[0];
r[1] += t * v[1];
}
// interaction with nucleus in shield
void interact() {
if (rg.rand() < p_a) { // absorbed by nucleus
absorbed = true;
} else { // scatter from nucleus
double v_mag = f * sqrt(v*v);
double theta = 2 * pi * rg.rand();
v[0] = v_mag * cos(theta);
v[1] = v_mag * sin(theta);
}
}
};
int main() {
cout << " Neutron Reactor Shielding Monte Carlo\n"
<< " -------------------------------------\n"
<< " Enter number of particles to simulate: ";
int n;
cin >> n;
// set random number generator seed
rg.set_seed_time();
cout << " Using " << rg.get_algorithm()
<< " and seed " << rg.get_seed() << endl;
ofstream file("neutron.data");
int step_sum = 0, absorbed = 0, escaped = 0;
for (int i = 0; i < n; i++) {
Neutron neutron;
int steps = 0;
do {
file << neutron.r[0] << '\t' << neutron.r[1] << '\n';
++steps;
} while (neutron.step_succeeded());
file << '\n';
step_sum += steps;
if (neutron.absorbed)
++absorbed;
if (neutron.escaped)
++escaped;
}
file.close();
cout << " Number escaped = " << escaped << '\n'
<< " Number absorbed = " << absorbed << '\n'
<< " Averge number of steps = " << step_sum/double(n) << '\n'
<< " " << n << " trajectories in file neutron.data" << endl;
}
```

*Edited 3 Years Ago by mike_2000_17*: Fixed formatting