Gravity Sim

COM-exclusive.c

#include <stdio.h>
//#pragma warning(disable: 4996)
#include <stdlib.h>
#include <math.h>
#define NUMPOINTS 3
#define STEPS 5000
//#define G 6.673E-11
#define G 5E-5 //bigger
#define CLOUDDIAMETER 1000
#define MINMASS 1000
#define MAXMASS 10000
#define PROGRESS 1000 // interval
#define FILENAME "COM exclusive.csv"

struct pointMass {
	double mass;
	double x; // position
	double y;
	double vx; // velocity
	double vy;
};

struct point {
	double x;
	double y;
};

int main(int argc, char *argv[]) {
	FILE *fp;

	struct pointMass point[NUMPOINTS];

	fp = fopen(FILENAME, "w");

	for (int np = 0; np < NUMPOINTS; np++) { // populate 100 point masses with data
		struct pointMass *p;
		p = &point[np];
		(*p).mass = (MINMASS + rand()) % (MAXMASS - MINMASS);
		(*p).x = rand() % CLOUDDIAMETER;
		(*p).y = rand() % CLOUDDIAMETER;
		(*p).vx = 0;
		(*p).vy = 0;
		fprintf(fp, "X%i,Y%i,", np, np);
	}
	fprintf(fp, "\n");
	// end setup

	for (int i = 0; i < STEPS; i++) { //foreaver ish
		struct point scratchPoint[NUMPOINTS];

		if (i % PROGRESS == 0) printf("%d ", i / PROGRESS); // progress bar

		for (int np = 0; np < NUMPOINTS; np++) { // for each point
			double forceX;
			double forceY;
			double massX = 0;
			double massY = 0;
			double dist;
			double force;
			double dir;
			struct pointMass *p;
			struct point *sp;
			struct pointMass com;
			com.mass = 0;
			p = &point[np];
			for (int nip = 0; nip < NUMPOINTS; nip++) { // check all influencing points
				struct pointMass *ip;
				ip = &point[nip]; // exerts force on p

				if (nip == np) continue; // we don't influence ourselves (technicly the calcs would agree, so this is optional)
				massX += (*ip).mass * (*ip).x; //http://hyperphysics.phy-astr.gsu.edu/hbase/cm.html
				massY += (*ip).mass * (*ip).y;
				com.mass += (*ip).mass;
			}

			com.x = massX / com.mass;
			com.y = massY / com.mass;

			dist = sqrt(
				pow((*p).x - com.x, 2) +
				pow((*p).y - com.y, 2)
				);
			force = G*(*p).mass*com.mass / pow(dist, 2);
			dir = atan2(com.x - (*p).x, com.y - (*p).y);
			forceX = force*cos(dir);
			forceY = force*sin(dir);

			//now we have forceX and forceY

			sp = &scratchPoint[np];

			// update velocity (we can to this live b.c velocity dosent affecy gravity)
			(*p).vx += forceX / (*p).mass; // f=ma a=delta v
			(*p).vy += forceY / (*p).mass;

			// update position (cant do this live b.c. it would affect later points)
			(*sp).x = (*p).x + (*p).vx;
			(*sp).y = (*p).y + (*p).vy;

			fprintf(fp, "%f,%f,", (*sp).x, (*sp).y);
		}

		for (int np = 0; np < NUMPOINTS; np++) { // move scratch points to points for next time
			struct pointMass *p;
			struct point *sp;
			p = &point[np];
			sp = &scratchPoint[np];
			(*p).x = (*sp).x;
			(*p).y = (*sp).y;
		}

		fprintf(fp, "\n");
	}
	puts("\n"); // it's just good form
}

COM-inclusive.c

#include <stdio.h>
//#pragma warning(disable: 4996)
#include <stdlib.h>
#include <math.h>
#define NUMPOINTS 3
#define STEPS 5000
//#define G 6.673E-11
#define G 5E-5 //bigger
#define CLOUDDIAMETER 1000
#define MINMASS 1000
#define MAXMASS 10000
#define PROGRESS 1000 // interval
#define FILENAME "COM inclusive.csv"
// future ideas http://www.sciencedirect.com/science/article/pii/S0019103513001942
struct pointMass {
	double mass;
	double x; // position
	double y;
	double vx; // velocity
	double vy;
};

struct point {
	double x;
	double y;
};

int main(int argc, char *argv[]) {
	FILE *fp;

	struct pointMass point[NUMPOINTS];

	fp = fopen(FILENAME, "w");

	for (int np = 0; np < NUMPOINTS; np++) { // populate 100 point masses with data
		struct pointMass *p;
		p = &point[np];
		(*p).mass = (MINMASS + rand()) % (MAXMASS - MINMASS);
		(*p).x = rand() % CLOUDDIAMETER;
		(*p).y = rand() % CLOUDDIAMETER;
		(*p).vx = 0;
		(*p).vy = 0;
		fprintf(fp, "X%i,Y%i,", np, np);
	}
	fprintf(fp, "\n");
	// end setup

	for (int i = 0; i < STEPS; i++) { //foreaver ish
		struct point scratchPoint[NUMPOINTS];

		if (i % PROGRESS == 0) printf("%d ", i / PROGRESS); // progress bar

		for (int np = 0; np < NUMPOINTS; np++) { // for each point
			double forceX;
			double forceY;
			double massX = 0;
			double massY = 0;
			double dist;
			double force;
			double dir;
			struct pointMass *p;
			struct point *sp;
			struct pointMass com;
			com.mass = 0;
			p = &point[np];
			for (int nip = 0; nip < NUMPOINTS; nip++) { // check all influencing points
				struct pointMass *ip;
				ip = &point[nip]; // exerts force on p


				//disabled as per http://ccl.northwestern.edu/netlogo/models/N-Bodies
				//if (nip == np) continue; // we don't influence ourselves (technicly the calcs would agree, so this is optional)
				massX += (*ip).mass * (*ip).x; //http://hyperphysics.phy-astr.gsu.edu/hbase/cm.html
				massY += (*ip).mass * (*ip).y;
				com.mass += (*ip).mass;
			}

			com.x = massX / com.mass;
			com.y = massY / com.mass;

			dist = sqrt(
				pow((*p).x - com.x, 2) +
				pow((*p).y - com.y, 2)
				);
			force = G*(*p).mass*com.mass / pow(dist, 2);
			dir = atan2(com.x - (*p).x, com.y - (*p).y);
			forceX = force*cos(dir);
			forceY = force*sin(dir);

			//now we have forceX and forceY

			sp = &scratchPoint[np];

			// update velocity (we can to this live b.c velocity dosent affecy gravity)
			(*p).vx += forceX / (*p).mass; // f=ma a=delta v
			(*p).vy += forceY / (*p).mass;

			// update position (cant do this live b.c. it would affect later points)
			(*sp).x = (*p).x + (*p).vx;
			(*sp).y = (*p).y + (*p).vy;

			fprintf(fp, "%f,%f,", (*sp).x, (*sp).y);
		}

		for (int np = 0; np < NUMPOINTS; np++) { // move scratch points to points for next time
			struct pointMass *p;
			struct point *sp;
			p = &point[np];
			sp = &scratchPoint[np];
			(*p).x = (*sp).x;
			(*p).y = (*sp).y;
		}

		fprintf(fp, "\n");
	}
	puts("\n"); // it's just good form
}

simple.c

#include <stdio.h>
//#pragma warning(disable: 4996)
#include <stdlib.h>
#include <math.h>
#define NUMPOINTS 3
#define STEPS 5000
//#define G 6.673E-11
#define G 5E-5 //bigger
#define CLOUDDIAMETER 1000
#define MINMASS 1000
#define MAXMASS 10000
#define PROGRESS 1000 // interval
#define FILENAME "simple.csv"

struct pointMass {
	double mass;
	double x; // position
	double y;
	double vx; // velocity
	double vy;
};

struct point {
	double x;
	double y;
};

int main(int argc, char *argv[]) {
	FILE *fp;

	struct pointMass point[NUMPOINTS];

	fp = fopen(FILENAME, "w");

	for (int np = 0; np < NUMPOINTS; np++) { // populate 100 point masses with data
		struct pointMass *p;
		p = &point[np];
		(*p).mass = (MINMASS + rand()) % (MAXMASS - MINMASS);
		(*p).x = rand() % CLOUDDIAMETER;
		(*p).y = rand() % CLOUDDIAMETER;
		(*p).vx = 0;
		(*p).vy = 0;
		fprintf(fp, "X%i,Y%i,", np, np);
	}
	fprintf(fp, "\n");
	// end setup

	for (int i = 0; i < STEPS; i++) { //foreaver ish
		struct point scratchPoint[NUMPOINTS];

		if (i % PROGRESS == 0) printf("%d ", i / PROGRESS); // progress bar

		for (int np = 0; np < NUMPOINTS; np++) { // for each point
			double forceX = 0;
			double forceY = 0;
			struct pointMass *p;
			struct point *sp;
			p = &point[np];
			for (int nip = 0; nip < NUMPOINTS; nip++) { // check all influencing points
				double dist;
				double force;
				double dir;
				struct pointMass *ip;
				ip = &point[nip]; // exerts force on p

				if (nip == np) continue; // we don't influence ourselves (technicly the calcs would agree, so this is optional)

				dist = sqrt(
					pow((*p).x - (*ip).x, 2) +
					pow((*p).y - (*ip).y, 2)
					);
				force = G*(*p).mass*(*ip).mass / pow(dist, 2);
				dir = atan2((*ip).x - (*p).x, (*ip).y - (*p).y);
				forceX += force*cos(dir);
				forceY += force*sin(dir);
			}
			//now we have forceX and forceY

			sp = &scratchPoint[np];

			// update velocity (we can to this live b.c velocity dosent affecy gravity)
			(*p).vx += forceX / (*p).mass; // f=ma a=delta v
			(*p).vy += forceY / (*p).mass;

			// update position (cant do this live b.c. it would affect later points)
			(*sp).x = (*p).x + (*p).vx;
			(*sp).y = (*p).y + (*p).vy;

			fprintf(fp, "%f,%f,", (*sp).x, (*sp).y);
		}

		for (int np = 0; np < NUMPOINTS; np++) { // move scratch points to points for next time
			struct pointMass *p;
			struct point *sp;
			p = &point[np];
			sp = &scratchPoint[np];
			(*p).x = (*sp).x;
			(*p).y = (*sp).y;
		}

		fprintf(fp, "\n");
	}
	puts("\n"); // it's just good form
}