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;

}
``````

## All 9 Replies

Beside the utterly ugly format you have presented us in this code, where's the problem in converting c++ to c? Post your conversion so far (of this code), so that we know you'd tried (and failed).

*Basically c is structed language and c++ is object oriented language that depend upon their objects and class.first i suggest you place your conversion code there so i can check then i resolve your problem and further proceed .okay
*

well now it can run but the result is not right. because now no matter how many neotron I type in, they all show escaped. So anyone can check the code for me to point out where I made wrong? Thanks.

It won't compile with Microsoft Visual Studio because M\$ has not implemented many things in the C99 standards, such as declaring variables in the middle of functions like you can in c++.

Just be told the allocation memory part is not quite right. Anyone got idea?

The dynamic allocation is unnecessary (and leaking), you could simply do this:

``````int i = 0;
for (i; i < n; i++) {

Neutron neutron;
initNeutron(&neutron);

int steps = 0;
do {
fprintf(file, "%f\t%f\n", neutron.r[0], neutron.r[1]);
++steps;
} while (step_succeeded(&neutron));

fprintf(file,"\n");
step_sum += steps;

if (neutron.absorbed)
++absorbed;
if (neutron.escaped)
++escaped;
}
``````

Another problem is that you translated this line of C++ code:

``````r = vector<double>(2, 0);
``````

into these lines of C code:

``````this->r[0] = 2;
this->r[1] = 0;
``````

You misinterpreted what the C++ code meant. The vector was being constructed with the second version of the constructors, see here, which takes the length of the vector as first argument (2) and then, the value with which to fill the vector (0). This means, the C code should be:

``````this->r[0] = 0;
this->r[1] = 0;
``````

Then, another big issue is with the random number generator. In the C++ code, you have this line:

``````double theta = 2 * pi * rg.rand();
``````

which takes a random number in the range `[0,1]` (produced by `rg.rand()`) and multiplies it with 2-pi in order to get a random angle between 0 and 2-pi. In your C code, you have:

``````double theta = 2 * pi * rand();
``````

which might look like the exact same thing to the untrained eye, but they are entirely different. The C `rand()` function generates a random integer number between 0 and `RAND_MAX` (which is usually a very large number, usually the max value of a 32bit integer). This would be the correct C code:

``````double theta = (2.0 * pi * rand()) / RAND_MAX;
``````

Another issue is within the `move_free` function:

``````if (rv < 0 && abs(rv) >= v2 * (r2 - a*a))
t_max = abs(rv)/v2 - sqrt(rv*rv/v2/v2 - (r2 - a*a)/v2);
``````

In C, the `abs()` function is only for `int` values. For `double` values, you need to use the `fabs()` function. This is confusing because in C++, with overloading, all this crazy naming business is obsolete, like in this case the C functions: `abs` / `labs` / `llabs` / `fabs` / `fabsf` / `fabsl` / `cabs` / `cabsf` / `cabsl`, for each respective type, `int` / `long int` / `long long int` / `double` / `float` / `long double` / `complex double` / `complex float` / `complex long double`.

Hi would you please show me how to use rand_MAX for entire code? and I changed

``````this->r[0] = 2;
this->r[1] = 0;
``````

to

``````this->r[0] = 0;
this->r[1] = 0;
``````

But the result just show now all neutron are absorbed.

and here is my code for now.

would you please show me how to use rand_MAX for entire code?

You could just create a global functions like this:

``````double drand() {
return ((double) rand()) / RAND_MAX;
};
``````

That function creates uniformly-distributed random floating-point variable between 0 and 1.

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, learning, and sharing knowledge.