[ create a new paste ] login | about

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

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

#include <stdlib.h>

#include <math.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 * this) {                        /* initalize neutron */

	this->r[0] = 0;

	this->r[1] = 0;

	double theta = 2 * pi * rand()/RAND_MAX;

    

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

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

    

	this->absorbed = false;

	this->escaped = false;

}



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

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

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

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

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

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

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

}



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

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

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

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

	double t_max;

	if (rv < 0 && fabs(rv) >= v2 * (r2 - a*a))

		t_max = fabs(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())/(RAND_MAX);

	if (t > t_max)

		t = t_max;

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

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

}



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

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

		this->absorbed = true;

	} else {                                                /* scatter from nucleus */

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

		double theta = (2 * pi * rand())/ RAND_MAX;

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

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

	}

}



bool step_succeeded(Neutron * this) {                        

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

		return false;

    

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

    

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

		this->escaped = true;

		return false;

	}

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

		move_to_a(this);

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

	interact(this);

	return true;

}







int main() {

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

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

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

	int n;

	scanf("%d", &n);

     
       (double) rand()/RAND_MAX;
   



	

	

 

    

    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=(Neutron *)malloc(sizeof(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: