Commits

Anonymous committed c261561 Draft

Simplify integrator.
Remove unused particle properties and functions.
Fix build.

Comments (0)

Files changed (4)

 all: $(TARG)
 
 $(TARG): $(OFILES) $(HFILES)
-	$(CC) $(LDFLAGS) -o $(TARG) $(OFILES)
+	$(CC) -o $(TARG) $(OFILES) $(LDFLAGS) 
 
 %.o: %.c
 	$(CC) $(CFLAGS) -c $<
 struct Part
 {
 	Vect	r;		/* current position */
-	Vect	oldr;	/* old position */
 	Vect	v;		/* current velocity */
-	Vect	f;		/* current forces */
 	Vect	a;		/* acceleration */
-	float	m;
-	float	k;
+	float	m;		/* mass */
+	float	k;		/* elastic constant */
 };
 
 struct Partsys
 		printparts(sys);
 	timestep(sys);
 	
-	usleep(sys->dt*1000000/2);
 	glutPostRedisplay();
 }
 
 #include <stdlib.h>
 #include <math.h>
 #include <sys/time.h>
+#include <time.h>
 #include "dat.h"
 #include "fns.h"
 
 {
 	Part *sp;
 	Partsys *sys;
-	int i;
 	
 	sys = emalloc(sizeof(Partsys));
 	sys->dt = dt;
 	sys->part = sp;
 	sys->nparts = n;
 	
-	/* 
-	 * verlet needs current r(t) and old r(t-dt) positions
-	 * to work. 'hardcode' them here.
-	 * comment to see it run wild.
-	 */
-	for(i=0; i < n; i++)
-		sys->part[i].oldr.x = sys->part[i].oldr.y = sys->part[i].oldr.z = 50;
-	
-	/* todo: parameterize in R^3 */
 	sys->box[0].x = sys->box[0].y = 0;
 	sys->box[1].x = sys->box[1].y = 780;
 	return sys;
 }
 
 void
-forces(Partsys *sys)
-{
-	int i;
-	
-	for(i=0; i < sys->nparts; i++)
-		vcopy(&sys->part[i].a, &sys->grav);
-}
-
-float
-verlet(float x, float oldx, float a, float dt)
-{
-	return 2*x - oldx + a*dt*dt;
-}
-
-void
 motion(Partsys *sys)
 {
 	int i;
-	Part *p;
-	Vect *a, *oldr, *r, *v, temp, *g;
-	float dt, vh, m;
+	Part *part;
+	Vect *a, *r, *v;
+	float dt;
  
 	dt = sys->dt;
-	g = &sys->grav;
 
 	for(i=0; i < sys->nparts; i++){
-		p = &sys->part[i];
-		r = &sys->part[i].r;
-		v = &sys->part[i].v;
-		oldr = &sys->part[i].oldr;
-		a = &sys->part[i].a;
-		m = sys->part[i].m;
-		vcopy(&temp, r);	
+		part = &sys->part[i];
+		r = &part->r;
+		v = &part->v;
+		a = &part->a;
 
 		if(debug)
-			printf("motion: r(%g,%g,%g) oldr(%g,%g,%g) a(%g,%g,%g)\n", 
-				r->x, r->y, r->z, oldr->x, oldr->y, oldr->z, a->x, a->y, a->z);
+			printf("motion: r(%g,%g,%g) a(%g,%g,%g)\n", r->x, r->y, r->z, a->x, a->y, a->z);
 	
 		/* velocity verlet */
 		r->x = r->x + v->x*dt + (1./2)*a->x*(dt*dt);
-		vh = v->x + (1./2)*a->x*dt;
-		a->x = (1/m)*m*g->x * r->x;
-		v->x = vh + (1./2)*a->x*dt;
+		v->x = v->x + a->x*dt;
 
 		r->y = r->y + v->y*dt + (1./2)*a->y*(dt*dt);
-		vh = v->y + (1./2)*a->y*dt;
-		a->y = (1/m)*m*g->y * r->y;
-		v->y = vh + (1./2)*a->y*dt;
-
-		vcopy(&sys->part[i].oldr, &temp);
+		v->y = v->y + a->y*dt;
 	}
 }
 
 void
 timestep(Partsys *sys)
 {
-	forces(sys);
+	/* we would apply forces here, but we only have gravity */
 	motion(sys);
 	constraints(sys);
 }
 	int i, n;
 	float xm, ym, k;
 
-	/* 
-	 * bug: when we change y to positive being up - 
-	 * and gravity to negative - every particle despite
-	 * of its initial position hits the border at the
-	 * same time.
-	 */
-	grav.x = 0;
-	grav.y = 0.5;
-	grav.z = 0;
-	
 	if(argc < 3){
 		fprintf(stderr, "usage: pys nparts k\n");
 		exit(1);
 		exit(1);
 	}
 
+	srand(time(NULL));
+
+	/* 
+	 * bug: when we change y to positive being up - 
+	 * and gravity to negative - every particle despite
+	 * of its initial position hits the border at the
+	 * same time.
+	 */
+	grav.x = 0;
+	grav.y = 0.5;
+	grav.z = 0;
+	
+
 	sys = newpartsys(n, &grav, 0.1);
-	if(!sys)
+	if(sys == NULL){
+		fprintf(stderr, "no memory for particle system\n");
 		exit(1);
+	}
 
 	xm = sys->box[1].x;
 	ym = sys->box[1].y;
 
 	for(i=0; i < sys->nparts; i++){
-		sys->part[i].r.x = rand()%(int)xm+Border;
-		sys->part[i].r.y = rand()%(int)ym+Border;
+		sys->part[i].r.x = (rand()%(int)xm)+Border;
+		sys->part[i].r.y = (rand()%(int)ym)+Border;
+		vcopy(&sys->part[i].a, &sys->grav);
 		sys->part[i].m = 1;
 		sys->part[i].k = k;
 		if(debug)