[ create a new paste ] login | about

Link: http://codepad.org/GfZBpNyF    [ raw code | output | fork ]

C, pasted on Mar 11:
#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <time.h>



const double 	pi = 3.1415926; /* pi */

double 			a = 1; /*inner radius of wall */

double 			b = 1.5; /* outer radius of wall */

double 			tau = 0.2; /* mean free time for v=1 */

double 			p_a = 0.05; /* probability of absorption */

double 			f = 0.95; /* inelasticity parameter */

double 			v_0 = 1; /*initial speed */



double dotProduct(double * v1, double * v2) {

	return v1[0]*v2[0] + v1[1]*v2[1];

}



typedef int bool;

#define true 1

#define false 0



/* define Neutron object */

typedef struct neutron_t { 

	double r[2];

	double v[2];                      /* position and velocity */

	bool absorbed;

	bool escaped;                   /* final fate */

} Neutron;



void initNeutron(Neutron * self) {                        /* initalize neutron */

	self->r[0] = 2;

	self->r[1] = 0;

	double theta = 2 * pi * rand();

    

	self->v[0] = (v_0 * cos(theta));

	self->v[1] = (v_0 * sin(theta));

    

	self->absorbed = false;

	self->escaped = false;

}



void move_to_a(Neutron * self) {                                /* move to inner wall of shield*/

	double r2 = dotProduct(self->r, self->r);

	double v2 = dotProduct(self->v, self->v);

	double rv = dotProduct(self->r, self->v);

	double t = sqrt((a*a - r2)/v2 + rv*rv/v2/v2) - rv/v2;

	self->r[0] += t * self->v[0];

	self->r[1] += t * self->v[1];

}



void move_free(Neutron * self) {                              /* in shield*/

	double r2 = dotProduct(self->r, self->r);

	double v2 = dotProduct(self->v, self->v);

	double rv = dotProduct(self->r, self->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(rand());

	if (t > t_max)

		t = t_max;

	self->r[0] += t * self->v[0];

	self->r[1] += t * self->v[1];

}



void interact(Neutron * self) {                               /* interaction with nucleus in shield */

	if (rand() < p_a) {                                  /* absorbed by nucleus */

		self->absorbed = true;

	} else {                                                /* scatter from nucleus */

		double v_mag = f * sqrt(dotProduct(self->v, self->v));

		double theta = 2 * pi * rand();

		self->v[0] = v_mag * cos(theta);

		self->v[1] = v_mag * sin(theta);

	}

}



bool step_succeeded(Neutron * self) {                        

	if (self->absorbed || self->escaped)

		return false;

    

	double r_mag = sqrt(dotProduct(self->r, self->r));

    

	if (r_mag > b) {                              /* outside the shield */

		self->escaped = true;

		return false;

	}

	if (r_mag < a)                                     /* inside core */ 

		move_to_a(self);

	move_free(self);                                     /* inside the shield */

	interact(self);

	return true;

}







int main() {

	printf( "Neutron Reactor Shielding simultion\n");

	printf( " -------------------------------------\n");

	printf( " Number of particles to simulate:\n ");

	int n;

	scanf("%d", &n);

    



	srand(time(NULL));

	

    

    FILE *file; 

	file = fopen("neutron.data","w");

    

	int step_sum = 0; 

	int absorbed = 0; 

	int escaped = 0;

    

    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;

	}

	fclose(file);

	printf( " Number escaped = %i\n", escaped);

	printf( " Number absorbed = %i\n", absorbed);

	printf( " %i trajectories in file neutron.data\n", n);

}


Output:
1
2
3
4
5
6
7
8
In function `interact':
undefined reference to `cos'
undefined reference to `sin'
In function `initNeutron':
undefined reference to `cos'
undefined reference to `sin'
In function `move_free':
undefined reference to `log'


Create a new paste based on this one


Comments: