Help needed convert C++ to C

[hxc210]

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;

}

[Lucaci Andrew]

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).

[asifalizaman]

*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
*

[hxc210]

here is the code I done so far. Can be compiled but the result isn't right.

http://codepad.org/oqau7J83

[hxc210]

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.

http://codepad.org/6kiNnL5O

[Ancient Dragon]

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++.

[hxc210]

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

[mike_2000_17]

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.

[hxc210]

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.

http://codepad.org/vBruw3v4

[mike_2000_17]

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.