1 #include <stdio.h>
2 #include "psys2d.h"
3 
4 /* each Cset should be inside the particle system itself I think */
5 Cset maxspeed = {"Max Speed",          110, 110,         1,           500, 10};
6 Cset gra_frc  = {"Graviational Force", 0.000000006, 0.000000006, 0,           1.0, 1e-10};
7 Cset ang_frc  = {"Angular Force",      30, 30,          1,           1};
8 Cset oth_frc  = {"Other Force",        0.35, 0.35,        0.0001,      100};
9 Cset spr_frc  = {"Spring Constant",    -0.000002, -0.000002,   -0.00000001, 1};
10 Cset spr_len  = {"Spring Rest Len",    90, 90,         90,          100};
11 
psys_reset_defaults(Particle_Sys * p)12 void psys_reset_defaults(Particle_Sys *p) {
13   maxspeed.cval = maxspeed.dflt;
14   gra_frc.cval  = gra_frc.dflt;
15   ang_frc.cval  = ang_frc.dflt;
16   oth_frc.cval  = oth_frc.dflt;
17   spr_frc.cval  = spr_frc.dflt;
18   spr_len.cval  = spr_len.dflt;
19 }
20 
show_cvar(Particle_Sys * psys,int i)21 void show_cvar(Particle_Sys *psys, int i) {
22   if (psys->cvar[i]) {
23     printf("%30s : ", psys->cvar[i]->control_name);
24     printf("%e  ", psys->cvar[i]->cval);
25     printf("(%e -> %e ->  %e by %e)\n",
26 	   psys->cvar[i]->min,
27 	   psys->cvar[i]->dflt,
28 	   psys->cvar[i]->max,
29 	   psys->cvar[i]->delta);
30   }
31 }
32 
psys_show_cset(Particle_Sys * psys)33 void psys_show_cset(Particle_Sys *psys) {
34   int i=0; for (i=0;i<4;i++) show_cvar(psys, i);
35 }
36 
psys_dec_var(Particle_Sys * psys,int vnum)37 void psys_dec_var(Particle_Sys *psys, int vnum) {
38   double temp;
39   if (psys->cvar[vnum]) {
40     temp=psys->cvar[vnum]->cval - psys->cvar[vnum]->delta;
41     if (temp > psys->cvar[vnum]->min)
42       psys->cvar[vnum]->cval=temp;
43     else
44       printf("%s is at min value %e\n",
45 	     psys->cvar[vnum]->control_name,
46 	     psys->cvar[vnum]->cval);
47   }
48 #ifdef DEBUG
49   psys_show_cset(psys);
50 #endif
51 }
52 
psys_inc_var(Particle_Sys * psys,int vnum)53 void psys_inc_var(Particle_Sys *psys, int vnum) {
54   double temp;
55   if (psys->cvar[vnum]) {
56     temp=psys->cvar[vnum]->cval + psys->cvar[vnum]->delta;
57     if (temp < psys->cvar[vnum]->max)
58       psys->cvar[vnum]->cval=temp;
59     else
60       printf("%s is at max value %e\n",
61 	     psys->cvar[vnum]->control_name,
62 	     psys->cvar[vnum]->cval);
63   }
64 #ifdef DEBUG
65   psys_show_cset(psys);
66 #endif
67 }
68 
69 static double drandom(void);
70 static void normalise(double, double*,  double *);
71 
drandom(void)72 static double drandom(void) { /* returns a number between 0 and 1 */
73   return (random()/(double)RAND_MAX);
74 }
75 
normalise(double max,double * x,double * y)76 void normalise(double max, double *x, double *y) {  /* Normalises a 2d vector */
77   /* need to normalise to a smaller value!! */
78   double d=(25000/max) * sqrt((*x)*(*x) + (*y)*(*y));
79   if (d>0) {
80     (*x)=(*x)/d;
81     (*y)=(*y)/d;
82   }
83 }
84 
psys_stick(Particle * p1,Particle * p2)85 void psys_stick(Particle *p1, Particle *p2) {
86 #ifdef DEBUG
87   printf("STUCK\n");
88 #endif
89   p1->xv   += p2->xv;
90   p1->yv   += p2->yv;
91 }
92 
93 
94 /**********************************************************************/
95 /********************* PARTICLE FORCE FUNCTIONS ***********************/
96 /**********************************************************************/
97 
psys_set_stf(Particle_Sys * psys)98 void psys_set_stf(Particle_Sys *psys) {
99   psys->frc = psys_frc_stf;
100   psys->cvar[0] = &maxspeed;
101   psys->cvar[1] = 0;
102   psys->cvar[2] = 0;
103   psys->cvar[3] = 0;
104 }
105 
psys_frc_stf(Particle_Sys * psys,Particle * p1,Particle * p2)106 void psys_frc_stf(Particle_Sys *psys, Particle *p1, Particle *p2) {
107 /* FORCE FUNCTION : STARFIELD */
108 /* CVAR1 : maxspeed */
109   /* This is a really awful, or worse, starfield simulation! */
110   if ((p1->x > 0.9999) || ((p1->x) < 0.0001) ||
111       (p1->y > 0.9999) || ((p1->y) < 0.0001)) {
112     p1->x = drandom()/10000 + 0.4995;
113     p1->y = drandom()/10000 + 0.4995;
114   }
115 
116   p1->xv += (p1->x - 0.5)/54232;
117   p1->yv += (p1->y - 0.5)/54232;
118   normalise(psys->cvar[0]->cval, &(p1->xv), &(p1->yv) );
119 }
120 
psys_set_ang(Particle_Sys * psys)121 void psys_set_ang(Particle_Sys *psys) {
122   psys->frc = psys_frc_ang;
123   psys->cvar[0] = &maxspeed;
124   psys->cvar[1] = &ang_frc;
125   psys->cvar[2] = 0;
126   psys->cvar[3] = 0;
127 }
128 
psys_frc_ang(Particle_Sys * psys,Particle * p1,Particle * p2)129 void psys_frc_ang(Particle_Sys *psys, Particle *p1, Particle *p2) {
130 /* FORCE FUNCTION : ANGULAR */
131 /* CVAR1 : maxspeed */
132 /* CVAR2 : str_frc_ang */
133   /* changes direction of p1 to point more towards p2 */
134   /* velocity stays the same. */
135   /* p2 stays still. (or at least, we don't alter it's velocity) */
136   double a = p1->xv /psys->cvar[1]->cval;
137   double b = p1->yv /psys->cvar[1]->cval;
138   double dx =  p2->x - p1->x;
139   double dy =  p2->y - p1->y;
140   double dxy = b*dx - dy*a;
141 
142   if (dxy>0) {
143     /* best route is anti-clockwise */
144     p1->xv += b;
145     p1->yv -= a;
146     p2->xv -= b;
147     p2->yv += a;
148   } else {
149     p1->xv -= b;
150     p2->xv += b;
151     p1->yv += a;
152     p2->yv -= a;
153   }
154   normalise(psys->cvar[0]->cval, &(p1->xv), &(p1->yv));
155   normalise(psys->cvar[0]->cval, &(p2->xv), &(p2->yv));
156 }
157 
158 
psys_set_oth(Particle_Sys * psys)159 void psys_set_oth(Particle_Sys *psys) {
160   psys->frc = psys_frc_oth;
161   psys->cvar[0] = &maxspeed;
162   psys->cvar[1] = &oth_frc;
163   psys->cvar[2] = 0;
164   psys->cvar[3] = 0;
165 }
166 
167 
psys_frc_oth(Particle_Sys * psys,Particle * p1,Particle * p2)168 void psys_frc_oth(Particle_Sys *psys, Particle *p1, Particle *p2) {
169 /* FORCE FUNCTION : OTHER*/
170 /* CVAR1 : str_frc_oth */
171 
172   double dx, dy, dxy;
173   dx = p2->x - p1->x;
174   dy = p2->y - p1->y;
175   if (!dx) dx=1;
176   if (!dy) dy=1;
177   dxy=sqrt(p1->xv*p1->xv + p1->yv*p1->yv);
178   p1->xv+=psys->cvar[1]->cval * dx * dxy;
179   p1->yv+=psys->cvar[1]->cval * dy * dxy;
180 }
181 
182 
psys_set_gra(Particle_Sys * psys)183 void psys_set_gra(Particle_Sys *psys) {
184   psys->frc = psys_frc_gra;
185   psys->cvar[0] = &maxspeed;
186   psys->cvar[1] = &gra_frc;
187   psys->cvar[2] = 0;
188   psys->cvar[3] = 0;
189 }
190 
psys_frc_sgra(Particle_Sys * psys,Particle * p1,Particle * p2)191 void psys_frc_sgra(Particle_Sys *psys, Particle *p1, Particle *p2) {
192 /* FORCE FUNCTION : GRAVITATIONAL */
193 /* CVAR1 : str_frc_gra */
194   double dx=p1->x - p2->x;
195   double dy=p1->y - p2->y;
196   double ddx, ddy;
197   double dxy = psys->cvar[1]->cval /(dx*dx + dy*dy);
198 
199   ddx = dx * dxy;    ddy = dy * dxy;
200   p1->xv -= ddx; p1->yv -= ddy;
201   p2->xv += ddx; p2->yv += ddy;
202 
203   dxy=-(dxy*dxy*dxy/psys->cvar[1]->cval);
204   ddx = dx * dxy;
205   ddy = dy * dxy;
206   p1->xv -= ddx; p1->yv -= ddy;
207   p2->xv += ddx; p2->yv += ddy;
208 
209 
210 }
211 
psys_frc_gra(Particle_Sys * psys,Particle * p1,Particle * p2)212 void psys_frc_gra(Particle_Sys *psys, Particle *p1, Particle *p2) {
213 /* FORCE FUNCTION : GRAVITATIONAL */
214 /* CVAR1 : str_frc_gra */
215   double dx=p1->x - p2->x;
216   double dy=p1->y - p2->y;
217   double dxy = 4.0 /(dx*dx - dy*dy);
218 
219   p1->xv -= dy/dxy;
220   p1->yv -= dx/dxy;
221   p2->xv -= dy/dxy;
222   p2->yv -= dx/dxy;
223 }
224 
psys_set_spr(Particle_Sys * psys)225 void psys_set_spr(Particle_Sys *psys) {
226   psys->frc = psys_frc_spr;
227   psys->cvar[0] = &maxspeed;
228   psys->cvar[1] = &spr_frc;
229   psys->cvar[2] = &spr_len;
230   psys->cvar[3] = 0;
231 }
232 
psys_frc_spr(Particle_Sys * psys,Particle * p1,Particle * p2)233 void psys_frc_spr(Particle_Sys *psys, Particle *p1, Particle *p2) {
234 /* FORCE FUNCTION : SPRING */
235 /* CVAR1 : maxspeed */
236 /* CVAR2 : kon_frc_spr */
237 /* CVAR3 : len_frc_spr */
238 
239  /* where konstant is the spring constant k in F=kx (where x is extension)
240     and springlen is the rest length of the spring */
241   /* attract particle 1 to Particle 2 using spring forces */
242   double dx=p1->x - p2->x;
243   double dy=p1->y - p2->y;
244   double dxy = sqrt(dx*dx + dy*dy);
245   /* dxy is now distance between the two particles */
246   if (0.000001>dxy) {
247     p1->xv=-p1->xv;
248     p2->xv=-p2->xv;
249   }
250   //    psys_stick(p1,p2);  /* too close to risk FPU /0 etc. or over/under flows */
251   else {
252     dxy=psys->cvar[1]->cval * (dxy-psys->cvar[2]->cval);
253     if (p1->life == -666)
254       dxy=dxy/8;
255     p1->xv -= dx*dxy;
256     p1->yv -= dy*dxy;
257     p2->xv += dx*dxy;
258     p2->yv += dy*dxy;
259     normalise(psys->cvar[0]->cval, &(p1->xv),&(p1->yv));
260   }
261 }
262 
263 
264 /* Set particle system to use a particlar force system */
psys_set_rnd(Particle_Sys * psys)265 void psys_set_rnd(Particle_Sys *psys) {
266   switch ((int)(random()%5)) {
267   case 0 : psys_set_gra(psys);break;
268   case 1 : psys_set_ang(psys);break;
269   case 2 : psys_set_oth(psys);break;
270   case 3 : psys_set_spr(psys);break;
271   case 4 : psys_set_stf(psys);break;
272   default : /* ONE WOULD HOPE THIS IS NEVER REACHED! */ break;
273   }
274 }
275 
276 /**********************************************************************/
277 /**********************************************************************/
278 /**********************************************************************/
279 
psys_set_region(Particle_Sys * psys,double xmin,double xmax,double ymin,double ymax)280 void psys_set_region(Particle_Sys *psys, double xmin, double xmax, double ymin, double ymax) {
281   psys->reg_xmin=xmin; psys->reg_xmax=xmax;
282   psys->reg_ymin=ymin; psys->reg_ymax=ymax;
283 }
284 
psys_rand_edge_starts(Particle_Sys * psys)285 void psys_rand_edge_starts(Particle_Sys *psys) {
286   switch (random()%4) {
287     case (0) : psys_set_region(psys,   0,   1, 0.9,   1); break; /* top */
288     case (1) : psys_set_region(psys,   0,   1,   0, 0.1); break; /* bottom */
289     case (2) : psys_set_region(psys,   0, 0.1,   0,   1); break; /* left */
290     case (3) : psys_set_region(psys, 0.9,   1,   0,   1); break; /* bottom */
291   }
292 }
293 
psys_rand_burstpt(Particle_Sys * psys)294 void psys_rand_burstpt(Particle_Sys *psys) {
295   /* sets a random burst point for new particle creation */
296   double x=(drandom() * 0.75)+0.1;
297   double y=(drandom() * 0.75)+0.1;
298   psys_set_region(psys,x,x+0.05,y,y+0.05);
299 }
300 
psys_toggle_gravity(Particle_Sys * psys)301 void psys_toggle_gravity(Particle_Sys *psys) {
302   psys->use_gravity=!psys->use_gravity;
303 }
304 
psys_set_defaults(Particle_Sys * psys)305 void psys_set_defaults(Particle_Sys *psys) {
306   /* SET PARTICLE SYSTEM DEFAULTS */
307   /* See particlesys.h for explanatory notes */
308 
309   psys_set_region(psys,0,1,0,1);
310 
311   psys->loadimage=0;
312   /*    psys->rep_frc_gra=1000; */
313 
314   /* PARTICLE SYSTEM DEFAULTS */
315   psys->n_particles=10000;
316 
317   /* FRICTION */
318   psys->use_friction=1;
319   psys->friction=0.999999;
320 
321   /* lives will be set to between +/-50% of this value */
322   psys->start_life=1000;
323 
324   /* AXIS-BASED GRAVITY */
325   psys->use_gravity=0;
326   psys->xgrav=0.0;
327   psys->ygrav=0.0005;
328 
329   /* MOTION TOWARDS A CENTRE POINT */
330   psys->tend_to_centpt=0;
331   psys->xcent=0.00005;
332   psys->ycent=0.00005;
333 
334   /* psys_set_ang */
335   psys_set_spr(psys);
336 }
337 
psys_resetparticle(Particle_Sys * psys,Particle * d)338 void psys_resetparticle(Particle_Sys *psys, Particle *d) {
339   if (d->life!=-666) {
340     d->y=psys->reg_ymin + drandom() * (psys->reg_ymax - psys->reg_ymin);
341     d->x=psys->reg_xmin + drandom() * (psys->reg_xmax - psys->reg_xmin);
342     d->life=random()%(psys->start_life) + psys->start_life;
343     d->xv=(rand()%10)/100000.0 - 0.0000501;
344     d->yv=(rand()%10)/100000.0 - 0.0000501;
345     d->xv/=100000.0;
346     d->yv/=100000.0;
347   }
348   else
349     printf("Asked to reset a fixed particle ? Why!?\n");
350 }
351 
352 
psys_setparticle_at(Particle_Sys * psys,Particle * d,double x,double y)353 void psys_setparticle_at(Particle_Sys *psys, Particle *d, double x, double y) {
354   // fix a particle at a given location.
355   d->y=0.1+y*0.8;
356   d->x=0.1+x*0.8;
357   d->xv=0;
358   d->yv=0;
359   d->life=-666;
360 }
361 
psys_load_fixed(Particle_Sys * psys)362 void psys_load_fixed(Particle_Sys *psys) {
363   unsigned int width_return;
364   unsigned int height_return;
365   unsigned char *data_return;
366   int x_hot_return;
367   int y_hot_return;
368   printf("Am here\n");
369   int rdbmpret = XReadBitmapFileData(psys->loadimage, &width_return, &height_return, &data_return, &x_hot_return, &y_hot_return);
370   printf("read that in okay\n");
371   printf("Reading from file %s returns %d\n", psys->loadimage, rdbmpret);
372   printf("Opened file is %d x %d, requiring a max of %d particles\n",width_return, height_return, width_return*height_return);
373 
374   int x,y;
375   int pcount =0;
376   int data;
377   int mult=10;
378   // work out how many particles are needed
379   for (y=0; y<height_return;y++) {
380     for (x=0; x<width_return;x++) {
381       data = data_return[x/8+y*(width_return+width_return%8)/8];
382       data = (data >> (x%8)) & 1;
383       pcount+=data*mult;
384     }
385   }
386   psys->n_particles=pcount+10;
387   printf("Set psys count to %d\n", pcount);
388   psys->p=(Particle *)malloc(sizeof(Particle)*psys->n_particles);
389 
390   pcount=0;
391   for (y=0; y<height_return;y++) {
392     for (x=0; x<width_return;x++) {
393       data = data_return[x/8+y*(width_return+width_return%8)/8];
394       data = (data >> (x%8)) & 1;
395       if (data) {
396         psys_setparticle_at(psys, psys->p+pcount, ((double)x)/width_return, ((double)y)/height_return);
397         pcount++;
398         int i=0;
399         for (i=0;i<mult-1;i++){
400           psys_resetparticle(psys, psys->p+pcount);
401           pcount++;
402         }
403       }
404     }
405   }
406   printf("Finished off having setup %d pixels from %d\n", pcount, psys->n_particles);
407   XFree(data_return);
408   printf("Even got here!\n");
409 }
410 
411 
psys_reset_all_particles(Particle_Sys * psys)412 void psys_reset_all_particles(Particle_Sys *psys) {
413   int i;
414   for (i=0;i<psys->n_particles;i++)
415     psys_resetparticle(psys,psys->p+i);
416 }
417 
psys_create_particles(Particle_Sys * psys)418 int psys_create_particles(Particle_Sys *psys) {
419   if (psys->loadimage!=(char*)0) {
420     printf("Loading image from %s instead\n", psys->loadimage);
421     psys_load_fixed(psys);
422     return 1;
423   }
424   psys->p=0;
425   psys->p=(Particle *)malloc(sizeof(Particle)*psys->n_particles);
426   if (!psys->p)
427     return 0;
428   psys_reset_all_particles(psys);
429   return 1;
430 }
431 
432 /* void psys_set_spr_forces(Particle_Sys *psys, double springlen, double konstant) {
433    psys->len_frc_spr=springlen; 0.0001
434    psys->kon_frc_spr=konstant;  0.00001
435    } */
436 
psys_frc_to_all(Particle_Sys * psys)437 void psys_frc_to_all(Particle_Sys *psys) {
438   /* apply the force system to all particles from all particles */
439   int count1,count2;
440   count1 = psys->n_particles;
441   while (count1--)  {
442     count2=count1;
443     while (count2--)
444       psys->frc(psys, psys->p+count1, psys->p+count2);
445   }
446 }
447 
psys_frc_to_lead(Particle_Sys * psys)448 void psys_frc_to_lead(Particle_Sys *psys) {
449   /* apply the force system to act as a force from the
450      lead particle on all other particles */
451   int count1=psys->n_particles;
452   while (count1--)
453     psys->frc(psys, psys->p+count1, psys->p);
454 }
455 
psys_frc_to_next(Particle_Sys * psys)456 void psys_frc_to_next(Particle_Sys *psys) {
457   /* apply the force system to act as a chain where
458      particle n is affected by particle n-1 */
459   int count1=psys->n_particles;
460   while (--count1)
461     psys->frc(psys, psys->p+count1, psys->p+count1+1);
462 }
463 
psys_move(Particle_Sys * psys)464 void psys_move(Particle_Sys *psys) {
465     int count1 = psys->n_particles;
466     Particle *cp;
467     double newx, newy;
468 
469     while (count1--) {
470 	cp = psys->p + count1;
471 
472         if (cp->life!=-666) {
473           /* particles with life set to -1 are fixed particles */
474 
475           if (!cp->life--) {
476             psys_resetparticle(psys,cp);
477           }
478 
479           if (psys->use_friction) {
480             cp->xv *=psys->friction;
481             cp->yv *=psys->friction;
482           }
483 
484           if (psys->tend_to_centpt) {
485             cp->xv += (cp->x < 0.5)? psys->xcent : -psys->xcent;
486             cp->yv += (cp->y < 0.5)? psys->ycent : -psys->ycent;
487           }
488 
489           normalise(psys->cvar[0]->cval, &cp->xv, &cp->yv);
490 
491           if (psys->use_gravity) {
492             cp->yv+=psys->ygrav;
493             cp->xv+=psys->xgrav;
494           }
495 
496           newx = cp->x + cp->xv;
497           newy = cp->y + cp->yv;
498 
499           /* WALL BOUNCING */
500           if (newx > 1) { newx = 1; cp->xv = -fabs(cp->xv); }
501           if (newx < 0) { newx = 0; cp->xv =  fabs(cp->xv); }
502           if (newy > 1) { newy = 1; cp->yv = -fabs(cp->yv); }
503           if (newy < 0) { newy = 0; cp->yv =  fabs(cp->yv); }
504 
505           cp->x=newx;
506           cp->y=newy;
507 
508         }
509     }
510 }
511 
psys_automovep(Particle_Sys * psys)512 void psys_automovep(Particle_Sys *psys) {
513     psys_frc_to_next(psys);
514     psys_move(psys);
515 }
516 
517 
518 
519