I need help

[cash]

I have a project that deals with a nuclear wall containing neutrons, and I'm confused on how I should set up my code? Here is a link of the project
https://www.chegg.com/homework-help/questions-and-answers/assignment-wall-nuclear-reactor-problem-descritpion-byproduct-nuclear-reaction-taking-plac-q35388290?trackid=kuMqsP0j

  #include <stdlib.h>
#include <iostream>
#include <iomanip>
using namespace std;
unsigned int seed = (unsigned)time(0);
double random(unsigned int &seed);
int main (){
    
    int particles = 1000, direction = 1, oldd, collision = 0;
    double entry_angle, deflection_angle=0;
    const int max_collsion = 7;
 
    double pctEscaped=0;
    while (particles != 0){
        oldd=1;
    
        collision=0;
        
        while( collision < max_collsion  && pctEscaped<.02 ){
            
               direction = .6+(.7*random(seed)); //pick a number .6 through 1.3
               entry_angle = -70+(140*random(seed)); //pick a angle -70 through 70
               deflection_angle = 10+(160*random(seed));//pick a angle 10 through 170
            if (oldd!=direction)
            {//compare old direction to direction
                collision++;
                oldd = direction;
            }
            
            
        }
       
    }
    

    
    return 0;
    
}

double random(unsigned int &seed)
{
    const int MODULUS = 15749;
    const int MULTIPLIER = 69069;
    const int INCREMENT = 1;
    seed = ((MULTIPLIER*seed)+INCREMENT)%MODULUS;
    return double (seed)/MODULUS;
}

[dutch]

I interpret the problem like this.

#include <iostream>
#include <vector>
#include <random>
#include <algorithm>
#include <chrono>
#include <cmath>

const int NumberOfNeutrons = 1000000;
const int CollisionsUntilAbsorbed = 7;
const double ToRad = acos(0) / 90.0;

int main() {
    std::vector<double> dist;

    auto seed = std::chrono::high_resolution_clock::now().time_since_epoch().count();
    auto rnd = std::default_random_engine(seed);

    auto distrib_entry     = std::uniform_real_distribution<>( -70.0 * ToRad,  70.0 * ToRad);
    auto distrib_collision = std::uniform_real_distribution<>(-160.0 * ToRad, 160.0 * ToRad);
    auto distrib_distance  = std::uniform_real_distribution<>(   0.6,   1.3);

    for (int n = 0; n < NumberOfNeutrons; ++n) {

        double total_dist = 0.0;
        double angle = distrib_entry(rnd);

        for (int i = 0; i < CollisionsUntilAbsorbed; ++i) {

            double ad = distrib_collision(rnd);
            if (ad < 0) ad -= 10.0 * ToRad;
            else        ad += 10.0 * ToRad;

            angle += ad;

            total_dist += distrib_distance(rnd) * std::cos(angle);
        }

        dist.push_back(total_dist);
    }

    std::sort(dist.begin(), dist.end());

    std::cout << dist[NumberOfNeutrons * 98 / 100] << '\n';
}

[Duthomhas]

You have the right idea: simulate some large number of particles (you have particles = 1000, good job!) working its way through a wall. The trick is to find the boundary between too-thick and too-thin.

You will need two loops (good job!). The OUTER loop ends when enough particles have hit the maximum number of collisions. (Good job!)

The INNER loop simulates a single particle. For each particle you must determine a random entry vector. Do this before you enter the inner loop.

Next you must track the travel vector as it bounces through the wall. At each bounce, calculate a new vector and travel distance.

ALSO track the cumulative ORTHOGONAL thickness of the wall as you bounce through it. You’ll need to remember your trigonometry (SOH-CAH-TOA). In other words, sum the X-components of the distance traveled by the particle at each bounce.

When you hit your inner loop exit condition (seven collisions!), then you know that the cumulative orthogonal distance (wall thickness) was sufficient for that particular particle.

Now, the nasty trick. You want to permit up to 2% of the particles to escape. You cannot do this with a loop condition.

Instead, consider the following:

    distance (wall thickness when particle hits 7th collision) -->
    | ~~~~~> | ~~~~~> |
    ^entry   ^min     ^max
    0        n        N

What percentage of the particles escape with a wall distance of n?
What percentage of the particles escape with a wall distance of N?

So the question is, where between n and N do only 2% of the particles escape?

Much like “finding the average” homeworks (of which this is a variation), you must keep track of the smallest wall thickness and the largest wall thickness while you are running your outer loop.

Then, AFTER you are done all the looping, figure out where between thinnest and thickest only 2% of the particles escape.

(This assumes, as does the homework, a normal fit to particle generation. The more particles you simulate, the more accurate the answer.)

Hope this helps.

[edit]
I like dutch’s solution too. Much nicer fit for the estimation.

[cash]

I tried redesigning dutch's program using syntax that I'm allowed to use,

#include <iostream>
#include <vector>
#include <ctime>
#include <algorithm>
#include <cmath>
using namespace std;
unsigned int seed = (unsigned)time(0);
double random(unsigned int &seed);
const int NumberOfNeutrons = 1000000;
const int CollisionsUntilAbsorbed = 7;
const double ToRad = acos(0) / 90.0;

int main() {
    vector<double> dist;
    
    double distrib_entry     = -70+(141*random(seed));
    double distrib_collision = -160+(321*random(seed));
    double distrib_distance  = .6+(1.7*random(seed));
    
    for (int n = 0; n < NumberOfNeutrons; ++n) {
        
        double total_dist = 0.0;
        double angle = distrib_entry;
        
        for (int i = 0; i < CollisionsUntilAbsorbed; ++i) {
            
            double ad = distrib_collision;
            if (ad < 0)
                ad -= 10.0 * ToRad;
            else
                ad += 10.0 * ToRad;
            
            angle += ad;
            
            total_dist += distrib_distance * cos(angle);
        }
        
        dist.push_back(total_dist);
    }
    
    sort(dist.begin(), dist.end());
    
    cout << dist[NumberOfNeutrons * 98 / 100] << '\n';
}

double random(unsigned int &seed)
{
    const int MODULUS = 15749;
    const int MULTIPLIER = 69069;
    const int INCREMENT = 1;
    seed = ((MULTIPLIER*seed)+INCREMENT)%MODULUS;
    return double (seed)/MODULUS;
}

However, I'm not getting good numbers for the thickness.

[Duthomhas]

You have made significant changes to it —primary of which is you have removed all randomness from the simulation. Do you see how?

I am also a little concerned about that random() function you are using. Is it being required by your professor?

[cash]

yes, I'm required to use the random() function, however, I don't see how I removed all randomness from the simulation. How can I use my random function and still make it random?

[Duthomhas]

Hmm, I am trying to think about this in a short way...

How can I use my random function and still make it random?

There is nothing wrong with your random function. It is working fine.

To explain, let us look at one of the original pseudo-random number generators: the John von Neumann’s Middle Squared method.^[1]

Start by taking a number. We will call this the seed.

    42

To get the next random number, square it, then take the middle two digits of the resulting four digit integer.

    42² = 1764
    Middle two digits are 76
    Answer: 76

There is your first pseudo-random number.

To get the next one, you just plug the number back in:

    76² = 5776
    Middle two digits are 77
    Answer: 77

Von Neumann used significantly larger numbers, but hopefully you should already see that this has some significant problems.^[2] (There are two. Can you figure out what they are?)


Even though the Middle-Squared Method has problems, it is the way that all PRNGs work. Except for the seed, which the user must supply, each successive generation is computed from the previous generation.

This is the case with that random() function you were given. The “seed” is initially a value obtained from the clock (via time(0)). After that, every time you want a new random number you use the last one you got out of it. That is why the seed argument is a reference.

I don't see how I removed all randomness from the simulation.

Knowing now how you compute a random number using successive permutations of a value, you must see that you cannot get random numbers without performing this permutation.^[3]

In other words, how often does your code call random(). The assignment says you must compute a random vector every time you start, and every time you have a collision.

If you pre-compute these vectors and reuse them, then all you are doing is sending the same particle through your simulation a thousand times.

But the whole point of the simulation is to send different particles through and tally what happens.

Look to see where you are computing randomness.
Hope this helps.



1. This was created back in the 1940s when developing computer simulations on nuclear weapon deployment. Interested people wanted to make sure detonating a bomb wouldn’t destroy all life on earth. (And then they just wanted to know how destructive their shiny new weapons could be.)

2. Von Neumann liked the problems because it became pretty obvious when the simulation failed.

3. Yes, this is a tautology.

[cash]

So basically in the for loop I should separately assign double angle, ad, and distrib_distance to random values.