X-Git-Url: https://git.ucc.asn.au/?a=blobdiff_plain;f=course%2Fsemester2%2Fpprog%2Fassignment1%2Fsingle-thread%2Fnbody.c;h=cb39b7feeb459b7b7efeace5bd995889b54e27ab;hb=HEAD;hp=e4b4bc701f1f1485231ee33c1b11867ae3d17431;hpb=f1878eefc82b0c465c0062b89ef39fdb4926e8cb;p=matches%2Fhonours.git

diff --git a/course/semester2/pprog/assignment1/single-thread/nbody.c b/course/semester2/pprog/assignment1/single-thread/nbody.c
index e4b4bc70..cb39b7fe 100644
--- a/course/semester2/pprog/assignment1/single-thread/nbody.c
+++ b/course/semester2/pprog/assignment1/single-thread/nbody.c
@@ -23,10 +23,12 @@ RUNSTATE runstate = RUN;
  * @param a - Body to print'
  * @param out - FILE* to print to
  */
-void Body_Print(Body * a, FILE * out)
+inline void Body_Print(Body * a, FILE * out)
 {
-	fprintf(out,"Body %p M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Vy=%f Vz=%f\n", 
-           (void*)a, a->mass, a->x[0], a->x[1], a->x[2], a->F[0], a->F[1], a->F[2], a->v[0], a->v[1], a->v[2]);
+	//fprintf(out,"Body %p M=%f X=%f Y=%f Z=%f Fx=%f Fy=%f Fz=%f Vx=%f Vy=%f Vz=%f\n", 
+           //(void*)a, a->mass, a->x[0], a->x[1], a->x[2], a->F[0], a->F[1], a->F[2], a->v[0], a->v[1], a->v[2]);
+	fprintf(out, "%f %f %f %f %f %f %f %f %f %f\n", 
+		a->mass, a->x[0], a->x[1], a->x[2], a->v[0], a->v[1], a->v[2], a->F[0], a->F[1], a->F[2]);
 }
 
 /**
@@ -35,11 +37,9 @@ void Body_Print(Body * a, FILE * out)
  * @param a - The Body
  * @param b - The System
  */
-void Body_Force(Body * a, System * s) 
+inline void Body_Forces(Body * a, System * s) 
 {
-	double distance;
-	double con;
-	double gd;
+
 
 	for (unsigned i = 0; i < DIMENSIONS; ++i) //Set Force to zero
 		a->F[i] = 0;
@@ -50,24 +50,38 @@ void Body_Force(Body * a, System * s)
 		if (b == a)
 			continue; //Otherwise we would have infinite force
 		
-		//Calculate distance between a and b
-		distance = 0.0;
-		for (unsigned i = 0; i < DIMENSIONS; ++i)
-			distance += square(b->x[i] - a->x[i]);
-		distance = sqrt(distance);
-		con = G * a->mass * b->mass / square(distance);
-		gd = con / distance;
-		for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
-			a->F[i] += gd * (b->x[i] - a->x[i]);
+		Body_Force(a, b);
 	}
 }
 
+/**
+ * @function Body_Force
+ * @purpose Compute force on body a due to body b
+ * @param a - Body a
+ * @param b - Body b
+ */
+inline void Body_Force(Body * a, Body * b)
+{
+	float distance;
+	float con;
+	float gd;
+	//Calculate distance between a and b
+	distance = 0.0;
+	for (unsigned i = 0; i < DIMENSIONS; ++i)
+		distance += square(b->x[i] - a->x[i]);
+	distance = sqrt(distance);
+	con = G * a->mass * b->mass / square(distance);
+	gd = con / distance;
+	for (unsigned i = 0; i < DIMENSIONS; ++i) //Add force from b to a's net force
+		a->F[i] += gd * (b->x[i] - a->x[i]);
+}
+
 /**
  * @function Body_Velocity
  * @purpose Compute velocity of a body
  * @param a - The Body
  */
-void Body_Velocity(Body * a) 
+inline void Body_Velocity(Body * a) 
 {
 	for (unsigned i = 0; i < DIMENSIONS; ++i)
 		a->v[i] += a->F[i] / a->mass * DELTA_T;
@@ -78,7 +92,7 @@ void Body_Velocity(Body * a)
  * @purpose Compute position of a body
  * @param a - The Body
  */
-void Body_Position(Body * a) 
+inline void Body_Position(Body * a) 
 {
 	for (unsigned i = 0; i < DIMENSIONS; ++i)
 		a->x[i] += a->v[i] * DELTA_T;
@@ -89,7 +103,7 @@ void Body_Position(Body * a)
  * @purpose Compute forces and then positions for bodies in System
  * NOTE: Only used in the single threaded version of the program
  */
-void System_Compute(System * s)
+inline void System_Compute(System * s)
 {
 	System_Forces(s, s);
 	System_Positions(s);
@@ -103,12 +117,12 @@ void System_Compute(System * s)
  *	 In multi threaded version, called with s2 == &universe, but s1 is constructed for each thread
  *		In nested multi-threaded version, s2 is constructed for the nested threads as well
  */
-void System_Forces(System * s1, System * s2) 
+inline void System_Forces(System * s1, System * s2) 
 {
 	//printf("Compute forces for %p - %d bodies...\n", (void*)s1, s1->N);
 	for (unsigned i = 0; i < s1->N; ++i)
 	{
-		Body_Force(s1->body+i, s2);
+		Body_Forces(s1->body+i, s2);
 		Body_Velocity(s1->body+i);
 	}
 	//printf("Compute positions for %p - Done!\n", (void*)s1);
@@ -119,13 +133,13 @@ void System_Forces(System * s1, System * s2)
  * @purpose Compute positions of all objects in System
  * @param s - The system
  */
-void System_Positions(System * s)
+inline void System_Positions(System * s)
 {
 	//printf("Compute positions for %p - %d bodies...\n", (void*)s, s->N);
 	for (unsigned i = 0; i < s->N; ++i)
 		Body_Position(s->body+i);
 	//printf("Compute positions for %p - Done!\n", (void*)s);
-	s->steps += 1;
+	
 }
 
 
@@ -137,13 +151,21 @@ void System_Positions(System * s)
  * @param s - The System
  * @param fileName - The input file
  */
-void System_Init(System * s, const char *fileName) 
+inline void System_Init(System * s, const char *fileName) 
 {
 	char line[LINE_SIZE];
 	char * token;
 	FILE * file;
 
 	file = fopen(fileName, "rt");
+
+	if (file == NULL)
+	{
+		fprintf(stderr, "%s\n", fileName);
+		perror("Couldn't open file");
+		exit(EXIT_FAILURE);
+	}
+
 	s->N = atoi(fgets(line, LINE_SIZE, file));
 	s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
 	s->steps = 0;
@@ -155,19 +177,23 @@ void System_Init(System * s, const char *fileName)
 		if (fgets(line, LINE_SIZE, file) != NULL)
 		{
 			Body * a = s->body+i;
-			token = strtok(line, ",; ");
+			token = strtok(line, ",; \t");
 			a->mass = atof(token);
 			
 			for (unsigned j = 0; j < DIMENSIONS; ++j)
 			{
-				token = strtok(NULL, ",; ");
+				token = strtok(NULL, ",; \t");
 				a->x[j] = atof(token);
 			}
 			for (unsigned j = 0; j < DIMENSIONS; ++j)
 			{
-				token = strtok(NULL, ",; ");
+				token = strtok(NULL, ",; \t");
 				a->v[j] = atof(token);
 			}
+
+			// Ignore any additional tokens
+			while (token != NULL)
+				token = strtok(NULL, ",; \t");
 			//Body_Print(a);
 		}
 	}
@@ -182,19 +208,30 @@ void System_Init(System * s, const char *fileName)
  * @purpose Cleans up the universe by freeing all memory
  *	 Also prints the bodies in the universe to a file specified in "options" if required.
  */
-void Universe_Cleanup()
+inline void Universe_Cleanup()
 {
+	//fprintf(stderr, "Called Universe_Cleanup()\n");
 	if (options.output != NULL) // An output file was specified...
 	{
-		FILE * save = fopen(options.output, "w");
+		FILE * save = NULL;
+		if (strcmp("stdout", options.output) == 0)
+			save = stdout;
+		else if (strcmp("stderr", options.output) == 0)
+			save = stderr;
+		else
+			save = fopen(options.output, "w");
 		if (save == NULL)
 			perror("Couldn't save universe to file.");
 		else
 		{
-			fprintf(save, "# Final field:\n");
+			// Print the output in the same format as the initial field file
+			fprintf(save, "%u\n", universe.N);
+				
 			for (unsigned i = 0; i < universe.N; ++i) // Print all bodies to the file
 			{
-				Body_Print(universe.body+i, save);
+				Body_Print(universe.body+i, save);			
+				
+								
 			}
 			fclose(save);
 		}
@@ -212,7 +249,7 @@ void Universe_Cleanup()
  *	- Clock cycles executed 
  *	- Equivelant time for a single thread to execute same number of clock cycles
  */
-void DisplayStatistics()
+inline void DisplayStatistics()
 {
 	clock_t end = clock();
 	struct timeval end_time;
@@ -232,11 +269,77 @@ void DisplayStatistics()
  * @purpose Helper to check whether the program is supposed to exit
  * 	Does not check for user triggered quit
  *	Checks for either a timeout, or computation of the required number of steps
+ * 
+ * The reason timeouts are integers is because this function is called a lot in lots of threads
+ *	So using floats (and working out microseconds every time) is expensive
+ *
+ */
+inline bool ExitCondition(void)
+{
+	bool result = (runstate != RUN || 
+			(options.num_steps >= 0 && universe.steps >= options.num_steps) ||
+			(options.timeout >= 0 && ((int)(time(NULL) - options.start_time.tv_sec) >= options.timeout)));
+	//fprintf(stderr,"runstate %d\n timeout %d\n steps %d\n", (int)(runstate != RUN), (int)(options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)), (int)(options.num_steps > 0 && universe.steps > options.num_steps));
+	return result;
+}
+
+/**
+ * @function Split_System
+ * @purpose Helper to divide one system into an array of systems
+ *	Each sub system will have N = (s->N / n) bodies in it
+ * @param s - The original system (typically &universe)
+ * @param n - The number of sub systems in the array
+ *
+ * WARNING: It is the caller's responsibility to free() the returned array
  */
-bool ExitCondition(void)
+System * Split_System(System * s, unsigned n)
 {
-	bool result = (runstate != RUN || (options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout))
-		|| (options.num_steps > 0 && universe.steps > options.num_steps));
-	//printf("runstate %d\n timeout %d\n steps %d\n", (int)(runstate != RUN), (int)(options.timeout > 0.00 && ((unsigned)(time(NULL) - options.start_time.tv_sec) >= options.timeout)), (int)(options.num_steps > 0 && universe.steps > options.num_steps));
+	System * result = (System*)(calloc(n, sizeof(System)));
+	if (result == NULL)
+	{
+		perror("Couldn't create array of sub systems");
+		exit(EXIT_FAILURE);
+	}
+
+	unsigned n_per_system = (s->N) / n;
+	unsigned remainder = (s->N) % n;
+
+	for (unsigned i = 0; i < n; ++i)	
+	{
+		result[i].N = n_per_system;
+		if (i == n-1)
+			result[i].N += remainder;
+		result[i].body = (s->body) + (n_per_system * i);
+		result[i].steps = 0;
+	}
 	return result;
 }
+
+/**
+ * @function System_Random
+ * @purpose Randomly generate initial body field
+ * @param s - The system
+ * @param n - Number of bodies
+ */
+void System_Random(System * s, int n)
+{
+	srand(0);
+	s->N = (unsigned)n;
+	s->body = (Body*) calloc((size_t)s->N, sizeof(Body));
+	s->steps = 0;
+
+	s->body[0].mass = 1e11;
+
+	for (unsigned a = 1; a < s->N; ++a)
+	{
+		s->body[a].mass = 1e7;
+		for (unsigned i = 0; i < DIMENSIONS; ++i)
+		{
+			s->body[a].x[i] = -10000 + rand() % 20000;
+			s->body[a].v[i] = -2 + rand() % 4; 
+		}
+		s->body[a].v[2] = 0;
+		s->body[a].x[2] = 0; //set z to zero; force in plane
+		
+	}
+}