#include <stdio.h>
#include "psys2d.h"

/* each Cset should be inside the particle system itself I think */
Cset maxspeed = {"Max Speed",          110, 110,         1,           500, 10};
Cset gra_frc  = {"Graviational Force", 0.000000006, 0.000000006, 0,           1.0, 1e-10};
Cset ang_frc  = {"Angular Force",      30, 30,          1,           1};
Cset oth_frc  = {"Other Force",        0.35, 0.35,        0.0001,      100};
Cset spr_frc  = {"Spring Constant",    -0.000002, -0.000002,   -0.00000001, 1};
Cset spr_len  = {"Spring Rest Len",    90, 90,         90,          100};

void psys_reset_defaults(Particle_Sys *p) {
  maxspeed.cval = maxspeed.dflt;
  gra_frc.cval  = gra_frc.dflt;
  ang_frc.cval  = ang_frc.dflt;
  oth_frc.cval  = oth_frc.dflt;
  spr_frc.cval  = spr_frc.dflt;
  spr_len.cval  = spr_len.dflt;
}

void show_cvar(Particle_Sys *psys, int i) {
  if (psys->cvar[i]) {
    printf("%30s : ", psys->cvar[i]->control_name);
    printf("%e  ", psys->cvar[i]->cval);
    printf("(%e -> %e ->  %e by %e)\n",
	   psys->cvar[i]->min,
	   psys->cvar[i]->dflt,
	   psys->cvar[i]->max,
	   psys->cvar[i]->delta);
  }
}
    
void psys_show_cset(Particle_Sys *psys) {
  int i=0; for (i=0;i<4;i++) show_cvar(psys, i);
}

void psys_dec_var(Particle_Sys *psys, int vnum) {
  double temp;
  if (psys->cvar[vnum]) {
    temp=psys->cvar[vnum]->cval - psys->cvar[vnum]->delta;
    if (temp > psys->cvar[vnum]->min)
      psys->cvar[vnum]->cval=temp;
    else
      printf("%s is at min value %e\n",
	     psys->cvar[vnum]->control_name,
	     psys->cvar[vnum]->cval);
  }
#ifdef DEBUG
  psys_show_cset(psys);
#endif
}

void psys_inc_var(Particle_Sys *psys, int vnum) {
  double temp;
  if (psys->cvar[vnum]) {
    temp=psys->cvar[vnum]->cval + psys->cvar[vnum]->delta;
    if (temp < psys->cvar[vnum]->max)
      psys->cvar[vnum]->cval=temp;
    else
      printf("%s is at max value %e\n",
	     psys->cvar[vnum]->control_name,
	     psys->cvar[vnum]->cval);
  }
#ifdef DEBUG
  psys_show_cset(psys);
#endif
}

static double drandom(void);
static void normalise(double, double*,  double *);

static double drandom(void) { /* returns a number between 0 and 1 */
  return (random()/(double)RAND_MAX);
}

void normalise(double max, double *x, double *y) {  /* Normalises a 2d vector */
  /* need to normalise to a smaller value!! */
  double d=(25000/max) * sqrt((*x)*(*x) + (*y)*(*y));
  if (d>0) {
    (*x)=(*x)/d;
    (*y)=(*y)/d;
  }
}

void psys_stick(Particle *p1, Particle *p2) {
#ifdef DEBUG
  printf("STUCK\n");
#endif
  p1->xv   += p2->xv;
  p1->yv   += p2->yv;
}


/**********************************************************************/
/********************* PARTICLE FORCE FUNCTIONS ***********************/
/**********************************************************************/

void psys_set_stf(Particle_Sys *psys) {
  psys->frc = psys_frc_stf;
  psys->cvar[0] = &maxspeed;
  psys->cvar[1] = 0;
  psys->cvar[2] = 0;
  psys->cvar[3] = 0;
}

void psys_frc_stf(Particle_Sys *psys, Particle *p1, Particle *p2) {
/* FORCE FUNCTION : STARFIELD */
/* CVAR1 : maxspeed */
  /* This is a really awful, or worse, starfield simulation! */
  if ((p1->x > 0.9999) || ((p1->x) < 0.0001) ||
      (p1->y > 0.9999) || ((p1->y) < 0.0001)) {
    p1->x = drandom()/10000 + 0.4995;
    p1->y = drandom()/10000 + 0.4995;
  }

  p1->xv += (p1->x - 0.5)/54232;
  p1->yv += (p1->y - 0.5)/54232;
  normalise(psys->cvar[0]->cval, &(p1->xv), &(p1->yv) );
}

void psys_set_ang(Particle_Sys *psys) {
  psys->frc = psys_frc_ang;
  psys->cvar[0] = &maxspeed;
  psys->cvar[1] = &ang_frc;
  psys->cvar[2] = 0;
  psys->cvar[3] = 0;
}

void psys_frc_ang(Particle_Sys *psys, Particle *p1, Particle *p2) {
/* FORCE FUNCTION : ANGULAR */
/* CVAR1 : maxspeed */
/* CVAR2 : str_frc_ang */
  /* changes direction of p1 to point more towards p2 */
  /* velocity stays the same. */
  /* p2 stays still. (or at least, we don't alter it's velocity) */
  double a = p1->xv /psys->cvar[1]->cval;
  double b = p1->yv /psys->cvar[1]->cval;
  double dx =  p2->x - p1->x;
  double dy =  p2->y - p1->y;
  double dxy = b*dx - dy*a;

  if (dxy>0) {
    /* best route is anti-clockwise */
    p1->xv += b;
    p1->yv -= a;
    p2->xv -= b;
    p2->yv += a;
  } else {
    p1->xv -= b;
    p2->xv += b;
    p1->yv += a;
    p2->yv -= a;
  }
  normalise(psys->cvar[0]->cval, &(p1->xv), &(p1->yv));
  normalise(psys->cvar[0]->cval, &(p2->xv), &(p2->yv));
}


void psys_set_oth(Particle_Sys *psys) {
  psys->frc = psys_frc_oth;
  psys->cvar[0] = &maxspeed;
  psys->cvar[1] = &oth_frc;
  psys->cvar[2] = 0;
  psys->cvar[3] = 0;
}


void psys_frc_oth(Particle_Sys *psys, Particle *p1, Particle *p2) {
/* FORCE FUNCTION : OTHER*/
/* CVAR1 : str_frc_oth */

  double dx, dy, dxy;
  dx = p2->x - p1->x;
  dy = p2->y - p1->y;
  if (!dx) dx=1;
  if (!dy) dy=1;
  dxy=sqrt(p1->xv*p1->xv + p1->yv*p1->yv);
  p1->xv+=psys->cvar[1]->cval * dx * dxy;
  p1->yv+=psys->cvar[1]->cval * dy * dxy;
}


void psys_set_gra(Particle_Sys *psys) {
  psys->frc = psys_frc_gra;
  psys->cvar[0] = &maxspeed;
  psys->cvar[1] = &gra_frc;
  psys->cvar[2] = 0;
  psys->cvar[3] = 0;
}

void psys_frc_sgra(Particle_Sys *psys, Particle *p1, Particle *p2) {
/* FORCE FUNCTION : GRAVITATIONAL */
/* CVAR1 : str_frc_gra */
  double dx=p1->x - p2->x;
  double dy=p1->y - p2->y;
  double ddx, ddy;
  double dxy = psys->cvar[1]->cval /(dx*dx + dy*dy);

  ddx = dx * dxy;    ddy = dy * dxy;
  p1->xv -= ddx; p1->yv -= ddy;
  p2->xv += ddx; p2->yv += ddy;

  dxy=-(dxy*dxy*dxy/psys->cvar[1]->cval);
  ddx = dx * dxy;
  ddy = dy * dxy;
  p1->xv -= ddx; p1->yv -= ddy;
  p2->xv += ddx; p2->yv += ddy;
  

}

void psys_frc_gra(Particle_Sys *psys, Particle *p1, Particle *p2) {
/* FORCE FUNCTION : GRAVITATIONAL */
/* CVAR1 : str_frc_gra */
  double dx=p1->x - p2->x;
  double dy=p1->y - p2->y;
  double dxy = 4.0 /(dx*dx - dy*dy);

  p1->xv -= dy/dxy;
  p1->yv -= dx/dxy;
  p2->xv -= dy/dxy;
  p2->yv -= dx/dxy;
}

void psys_set_spr(Particle_Sys *psys) {
  psys->frc = psys_frc_spr;
  psys->cvar[0] = &maxspeed;
  psys->cvar[1] = &spr_frc;
  psys->cvar[2] = &spr_len;
  psys->cvar[3] = 0;
}

void psys_frc_spr(Particle_Sys *psys, Particle *p1, Particle *p2) {
/* FORCE FUNCTION : SPRING */
/* CVAR1 : maxspeed */
/* CVAR2 : kon_frc_spr */
/* CVAR3 : len_frc_spr */

 /* where konstant is the spring constant k in F=kx (where x is extension)
    and springlen is the rest length of the spring */
  /* attract particle 1 to Particle 2 using spring forces */
  double dx=p1->x - p2->x;
  double dy=p1->y - p2->y;
  double dxy = sqrt(dx*dx + dy*dy);
  /* dxy is now distance between the two particles */
  if (0.000001>dxy) {
    p1->xv=-p1->xv;
    p2->xv=-p2->xv;
  }
  //    psys_stick(p1,p2);  /* too close to risk FPU /0 etc. or over/under flows */
  else {
    dxy=psys->cvar[1]->cval * (dxy-psys->cvar[2]->cval);
    if (p1->life == -666)
      dxy=dxy/8;
    p1->xv -= dx*dxy;
    p1->yv -= dy*dxy;
    p2->xv += dx*dxy;
    p2->yv += dy*dxy; 
    normalise(psys->cvar[0]->cval, &(p1->xv),&(p1->yv));
  }
}


/* Set particle system to use a particlar force system */
void psys_set_rnd(Particle_Sys *psys) {
  switch ((int)(random()%5)) {
  case 0 : psys_set_gra(psys);break;
  case 1 : psys_set_ang(psys);break;
  case 2 : psys_set_oth(psys);break;
  case 3 : psys_set_spr(psys);break;
  case 4 : psys_set_stf(psys);break;
  default : /* ONE WOULD HOPE THIS IS NEVER REACHED! */ break;
  }
}

/**********************************************************************/
/**********************************************************************/
/**********************************************************************/

void psys_set_region(Particle_Sys *psys, double xmin, double xmax, double ymin, double ymax) {
  psys->reg_xmin=xmin; psys->reg_xmax=xmax;
  psys->reg_ymin=ymin; psys->reg_ymax=ymax;
}

void psys_rand_edge_starts(Particle_Sys *psys) {
  switch (random()%4) {
    case (0) : psys_set_region(psys,   0,   1, 0.9,   1); break; /* top */
    case (1) : psys_set_region(psys,   0,   1,   0, 0.1); break; /* bottom */
    case (2) : psys_set_region(psys,   0, 0.1,   0,   1); break; /* left */
    case (3) : psys_set_region(psys, 0.9,   1,   0,   1); break; /* bottom */
  }
}

void psys_rand_burstpt(Particle_Sys *psys) {
  /* sets a random burst point for new particle creation */
  double x=(drandom() * 0.75)+0.1;
  double y=(drandom() * 0.75)+0.1;
  psys_set_region(psys,x,x+0.05,y,y+0.05);
}

void psys_toggle_gravity(Particle_Sys *psys) {
  psys->use_gravity=!psys->use_gravity;
}

void psys_set_defaults(Particle_Sys *psys) {
  /* SET PARTICLE SYSTEM DEFAULTS */
  /* See particlesys.h for explanatory notes */

  psys_set_region(psys,0,1,0,1);
    
  psys->loadimage=0;
  /*    psys->rep_frc_gra=1000; */
  
  /* PARTICLE SYSTEM DEFAULTS */
  psys->n_particles=10000;

  /* FRICTION */
  psys->use_friction=1;
  psys->friction=0.999999;

  /* lives will be set to between +/-50% of this value */
  psys->start_life=1000;

  /* AXIS-BASED GRAVITY */
  psys->use_gravity=0;
  psys->xgrav=0.0;
  psys->ygrav=0.0005;

  /* MOTION TOWARDS A CENTRE POINT */
  psys->tend_to_centpt=0;
  psys->xcent=0.00005;
  psys->ycent=0.00005;

  /* psys_set_ang */
  psys_set_spr(psys);
}

void psys_resetparticle(Particle_Sys *psys, Particle *d) {
  if (d->life!=-666) {
    d->y=psys->reg_ymin + drandom() * (psys->reg_ymax - psys->reg_ymin);
    d->x=psys->reg_xmin + drandom() * (psys->reg_xmax - psys->reg_xmin);
    d->life=random()%(psys->start_life) + psys->start_life;
    d->xv=(rand()%10)/100000.0 - 0.0000501;
    d->yv=(rand()%10)/100000.0 - 0.0000501;
    d->xv/=100000.0;
    d->yv/=100000.0;
  }
  else
    printf("Asked to reset a fixed particle ? Why!?\n");
}
  

void psys_setparticle_at(Particle_Sys *psys, Particle *d, double x, double y) {
  // fix a particle at a given location.
  d->y=0.1+y*0.8;
  d->x=0.1+x*0.8;
  d->xv=0;
  d->yv=0;
  d->life=-666;
}

void psys_load_fixed(Particle_Sys *psys) {
  unsigned int width_return;
  unsigned int height_return;
  unsigned char *data_return;
  int x_hot_return;
  int y_hot_return;
  printf("Am here\n");
  int rdbmpret = XReadBitmapFileData(psys->loadimage, &width_return, &height_return, &data_return, &x_hot_return, &y_hot_return);
  printf("read that in okay\n");
  printf("Reading from file %s returns %d\n", psys->loadimage, rdbmpret);
  printf("Opened file is %d x %d, requiring a max of %d particles\n",width_return, height_return, width_return*height_return);

  int x,y;
  int pcount =0;
  int data;
  int mult=10;
  // work out how many particles are needed
  for (y=0; y<height_return;y++) {
    for (x=0; x<width_return;x++) {
      data = data_return[x/8+y*(width_return+width_return%8)/8];
      data = (data >> (x%8)) & 1;
      pcount+=data*mult;
    }
  }
  psys->n_particles=pcount+10;
  printf("Set psys count to %d\n", pcount);
  psys->p=(Particle *)malloc(sizeof(Particle)*psys->n_particles);

  pcount=0;
  for (y=0; y<height_return;y++) {
    for (x=0; x<width_return;x++) {
      data = data_return[x/8+y*(width_return+width_return%8)/8];
      data = (data >> (x%8)) & 1;
      if (data) {
        psys_setparticle_at(psys, psys->p+pcount, ((double)x)/width_return, ((double)y)/height_return);
        pcount++;
        int i=0;
        for (i=0;i<mult-1;i++){
          psys_resetparticle(psys, psys->p+pcount);
          pcount++;
        }
      }
    }
  }
  printf("Finished off having setup %d pixels from %d\n", pcount, psys->n_particles);
  XFree(data_return);
  printf("Even got here!\n");
}
  

void psys_reset_all_particles(Particle_Sys *psys) {
  int i;
  for (i=0;i<psys->n_particles;i++)
    psys_resetparticle(psys,psys->p+i);
}
  
int psys_create_particles(Particle_Sys *psys) {
  if (psys->loadimage!=(char*)0) {
    printf("Loading image from %s instead\n", psys->loadimage);
    psys_load_fixed(psys);
    return 1;
  }
  psys->p=0;
  psys->p=(Particle *)malloc(sizeof(Particle)*psys->n_particles);
  if (!psys->p)
    return 0;
  psys_reset_all_particles(psys);
  return 1;
}

/* void psys_set_spr_forces(Particle_Sys *psys, double springlen, double konstant) { 
   psys->len_frc_spr=springlen; 0.0001 
   psys->kon_frc_spr=konstant;  0.00001 
   } */

void psys_frc_to_all(Particle_Sys *psys) {
  /* apply the force system to all particles from all particles */
  int count1,count2;
  count1 = psys->n_particles;
  while (count1--)  {
    count2=count1;
    while (count2--)
      psys->frc(psys, psys->p+count1, psys->p+count2);
  }
}

void psys_frc_to_lead(Particle_Sys *psys) {
  /* apply the force system to act as a force from the
     lead particle on all other particles */
  int count1=psys->n_particles;
  while (count1--)
    psys->frc(psys, psys->p+count1, psys->p);
}

void psys_frc_to_next(Particle_Sys *psys) {
  /* apply the force system to act as a chain where
     particle n is affected by particle n-1 */
  int count1=psys->n_particles;
  while (--count1)
    psys->frc(psys, psys->p+count1, psys->p+count1+1);
}

void psys_move(Particle_Sys *psys) {
    int count1 = psys->n_particles;
    Particle *cp;
    double newx, newy;

    while (count1--) {
	cp = psys->p + count1;

        if (cp->life!=-666) {
          /* particles with life set to -1 are fixed particles */
          
          if (!cp->life--) {
            psys_resetparticle(psys,cp);
          }
          
          if (psys->use_friction) {
            cp->xv *=psys->friction;
            cp->yv *=psys->friction;
          }
          
          if (psys->tend_to_centpt) {
            cp->xv += (cp->x < 0.5)? psys->xcent : -psys->xcent;
            cp->yv += (cp->y < 0.5)? psys->ycent : -psys->ycent;
          }
          
          normalise(psys->cvar[0]->cval, &cp->xv, &cp->yv);
          
          if (psys->use_gravity) {
            cp->yv+=psys->ygrav;
            cp->xv+=psys->xgrav;
          }
          
          newx = cp->x + cp->xv;
          newy = cp->y + cp->yv;
          
          /* WALL BOUNCING */
          if (newx > 1) { newx = 1; cp->xv = -fabs(cp->xv); }
          if (newx < 0) { newx = 0; cp->xv =  fabs(cp->xv); }
          if (newy > 1) { newy = 1; cp->yv = -fabs(cp->yv); }
          if (newy < 0) { newy = 0; cp->yv =  fabs(cp->yv); }

          cp->x=newx;
          cp->y=newy;
          
        }
    }
}

void psys_automovep(Particle_Sys *psys) {
    psys_frc_to_next(psys);
    psys_move(psys);
}