/*
 * Copyright (c) 1994 Pangfeng Liu.  All rights reserved.
 */


#include <stdio.h>
#include <math.h>
#include "particle.h"

PARTICLE        sample_particles[MAX_SAMPLE_POINTS];
int 		N; 		/* number of particles 	*/
int		seed;
char		*infile_name, *outfile_name;
FILE		*infile, *outfile, *unit_file;
int             debug;

double ran1();
double sqrt();

int shape;
int number_of_galaxy;
VECTOR center[MAX_NUM_GALAXY];
float radius[MAX_NUM_GALAXY];

char *pos_string(pos)
     VECTOR pos;
{
  char *string;
  char *malloc();
  
  string = malloc(128);
  
  sprintf(string, "%25.22lf %25.22lf %25.22lf", pos.cord[0], pos.cord[1], pos.cord[2]);
  
  return(string);
}


char *getparam(par)
     char *par;                      /* name of parameter */
{
  char buf[128], *malloc(), *strncpy();
  int len;
  /* extern int strlen(); */
  
  printf("Enter %s: ", par);
  gets(buf);
  len = strlen(buf) + 1;
  return (strncpy(malloc(len), buf, len));
}


void set_vs(v, s)
     VECTOR *v;
     REAL s;
{
  int i;
  
  for (i = 0; i < DIMENSION; i++)
    v->cord[i] = s;
}

VECTOR sub_vv(u, v)
     VECTOR u, v;
{
  VECTOR w;
  int i;
  
  for (i = 0; i < DIMENSION; i++)
    w.cord[i] = u.cord[i] - v.cord[i];
  
  return(w);
}

VECTOR add_vv(u, v)
     VECTOR u, v;
{
  VECTOR w;
  int i;
  
  for (i = 0; i < DIMENSION; i++)
    w.cord[i] = u.cord[i] + v.cord[i];
  
  return(w);
}

VECTOR div_vs(v, s)
     VECTOR v;
     REAL s;
{
  VECTOR w;
  int i;
  
  for (i = 0; i < DIMENSION; i++)
    w.cord[i] = v.cord[i] / s;
  
  return(w);
}

VECTOR mul_vs(v, s)
     VECTOR v;
     REAL s;
{
  VECTOR w;
  int i;
  
  for (i = 0; i < DIMENSION; i++)
    w.cord[i] = v.cord[i] * s;
  
  return(w);
}

/* compute the inner product of two vectors */

REAL dot(u, v)
     VECTOR u, v;
{
  REAL d = 0.0;
  int i;
  
  for (i = 0; i < DIMENSION; i++)
    d += u.cord[i] * v.cord[i];
  return(d);
}


VECTOR random_vector(distribution)
     int distribution;
{
  register int k;
  REAL rsq, xrand();
  double sqrt();
  VECTOR x, unit_vector;
  
  do {
    rsq = 0.0;
    for (k = 0; k < DIMENSION; k++) 
      x.cord[k] = xrand(-RADIUS, RADIUS);
    
    rsq = dot(x, x);
  } while (rsq > RADIUS * RADIUS);
  
  if (debug)
    printf("x = %s\n", pos_string(x));
  return(x);
}


REAL xrand(xl, xh)
     REAL xl, xh;				/* lower, upper bounds on number */
{
  double ran1();
  
  return(xl + (xh - xl) * (ran1(&seed)));
}

void init_particles(distribution)
     int distribution;
{
  int i, j, k, n;
  VECTOR cmr, cmv, offset, x;
  PARTICLE *p;
  int par_number_in_one_galaxy;
  REAL rsq, r2;
  
  seed = -1234;
  ran1(&seed);
  
  set_vs(&cmr, 0.0);                          
  set_vs(&cmv, 0.0);
  set_vs(&offset, 5.0);
  
  printf("generating particles...");
  fflush(stdout);
  
  switch (distribution)
    {
    case RANDOM:
      
      for (p = sample_particles, n = 0; n < N; n++, p++)
	{ 
	  r2 = RADIUS * RADIUS;
	  p->mass = 1.0 / N;  
	  do {
	    rsq = 0.0;
	    for (k = 0; k < DIMENSION; k++) 
	      x.cord[k] = xrand(-RADIUS, RADIUS);
	    rsq = dot(x, x);
	  } while (rsq > r2);
	  p->position = x;
	  
	  cmr = add_vv(cmr, p->position);
	  pickshell(&(p->velocity), xrand(0.0, 1.0));	
	  cmv = add_vv(cmv, p->velocity);	   /* add to running sum */ 
	}
      break;
      
    case PLUMMER:
      {
	REAL rsc, vsc, r, v, x, y;
	VECTOR cmr, cmv;
	
	p = sample_particles;
	rsc = 3 * PI / 16;		/* set length scale factor (0.589) */
	vsc = sqrt(1.0 / rsc);		/* and recip. speed scale (1.303) */
	
	for (j = 0; j < N; j++, p++)
	  {
	    p->mass = 1.0 / N;		/* set masses equal */
	    r = 1 / sqrt(pow(xrand(0.0, MFRAC), /* pick r in struct units */
			     -2.0/3.0) - 1);
	    
	    pickshell(&(p->position), rsc * r); /* pick scaled position */
	    cmr = add_vv(cmr, p->position);	/* add to running sum */
	    do {				/* select from fn g(x) */
	      x = xrand(0.0, 1.0);		/* for x in range 0:1   */
	      y = xrand(0.0, 0.1);		/* max of g(x) is 0.092 */
	    } while (y > x*x * pow(1 - x*x, 3.5));/* using von Neumann tech */
	    v = sqrt(2.0) * x/pow(1 + r*r, 0.25);/* find v in struct units */
	    pickshell(&(p->velocity), vsc * v);	   /* pick scaled velocity */
	    cmv = add_vv(cmv, p->velocity);	   /* add to running sum */
	  }
      }
      break;
      
    case TWO_GALAXIES:
      {
	REAL rsc, vsc, r, v, x, y;
	VECTOR cmr, cmv;
	
	p = sample_particles;
	rsc = 3 * PI / 16;              /* set length scale factor (0.589) */
	vsc = sqrt(1.0 / rsc);          /* and recip. speed scale (1.303) */
	
	for (j = 0; j < N; j++, p++)
	  {
	    p->mass = 1.0 / N;          /* set masses equal */
	    do
	      {
		r = 1 / sqrt(pow(xrand(0.0, MFRAC), -2.0/3.0) - 1);
		/* pick r in struct units */
	      }
	    while (r > 5.0);
	    
	    pickshell(&(p->position), rsc * r); /* pick scaled position */
	    if (j < N/2)
	      p->position = add_vv(p->position, offset);
	    else
	      p->position = sub_vv(p->position, offset);
	    cmr = add_vv(cmr, p->position);     /* add to running sum */
	    
	    do {                                /* select from fn g(x) */
	      x = xrand(0.0, 1.0);              /* for x in range 0:1   */
	      y = xrand(0.0, 0.1);              /* max of g(x) is 0.092 */
	    } while (y > x*x * pow(1 - x*x, 3.5));/* using von Neumann tech */
	    v = sqrt(2.0) * x/pow(1 + r*r, 0.25);/* find v in struct units */
	    pickshell(&(p->velocity), vsc * v);    /* pick scaled velocity */
	    cmv = add_vv(cmv, p->velocity);        /* add to running sum */
	  }
      }
      break;
      
      
    default:
      panic("wrong distribution\n");
    }
  
  printf("done\n");
  
  cmr = div_vs(cmr, (float) N);           
  cmv = div_vs(cmv, (float) N);
  
  printf("normalizing data...");
  fflush(stdout);
  for (p = sample_particles, n = 0; n < N; n++, p++)
    p->position = sub_vv(p->position, cmr);          
  
  printf("done\n");
}

void main()     
{
  int i, j;
  int dis;
  PARTICLE *p;
  
  N = atoi(getparam("particle number"));
  if (N > MAX_SAMPLE_POINTS)
    panic("particle number is larger than MAX_SAMPLE_POINTS\n");
  
  outfile_name = getparam("outfile name");
  debug = atoi(getparam("debug_mode"));
  
  if ((outfile = fopen(outfile_name, "w")) == NULL)
    panic("par_gen: can't open output file\n");
  
  shape = atoi(getparam("0: uniform sphere\n      1: Plummer model\n      2: Two galaxies\n      (0/1)?"));
  
  init_particles(shape);  
  
  if (debug && (N <= 2000))
    {
      FILE *fp;
      fp = fopen("pos", "w");
      printf("position\n");
      for (i = 0, p = sample_particles; i < N; i++, p++)
	fprintf(fp, "%s\n", pos_string(p->position));
      fclose(fp);
      
      printf("velocity\n");
      fp = fopen("vel", "w");
      for (i = 0, p = sample_particles; i < N; i++, p++)
	fprintf(fp, "%s\n", pos_string(p->velocity));
      fclose(fp);
    }
  
  print_particles();
  fclose(outfile);
}

panic(string)
     char *string;
{
  printf("panic: %s\n", string);
  exit(-1);
}

print_particles()
{
  int i, d;
  PARTICLE *p;
  REAL REAL_buf[MAX_SAMPLE_POINTS * DIMENSION], *REAL_buf_ptr;
  
  if (debug)
    printf("par_gen: begins print_particles\n");
  
  fwrite(&N, sizeof(int), 1, outfile);
  
  REAL_buf_ptr = REAL_buf;
  for (i = 0; i < N; i++)
    {
      p = &sample_particles[i];
      for (d = 0; d < DIMENSION; d++)
	*REAL_buf_ptr++ = p->position.cord[d];
    }
  fwrite(REAL_buf, sizeof(REAL), N * DIMENSION, outfile);
  
  REAL_buf_ptr = REAL_buf;
  for (i = 0; i < N; i++)
    {
      p = &sample_particles[i];
      for (d = 0; d < DIMENSION; d++)
	*REAL_buf_ptr++ = p->velocity.cord[d];
    }
  fwrite(REAL_buf, sizeof(REAL), N * DIMENSION, outfile);
  
  fflush(outfile);
  
  if (debug)
    printf("par_gen: finished print_particles\n");
}

/*
 * PICKSHELL: pick a random point on a sphere of specified radius.
 */

pickshell(vec, rad)
     VECTOR *vec;
     REAL rad;
{
  int k;
  REAL rsq, rsc;
  
  do {					/* pick point in NDIM-space */
    for (k = 0; k < DIMENSION; k++)		/* loop over dimensions   */
      vec->cord[k] = xrand(-1.0, 1.0);   	/* pick from unit cube  */
    rsq = dot(*vec, *vec);			/* compute radius squared */
  } while (rsq > 1.0);                	/* reject if outside sphere */
  rsc = rad / sqrt(rsq);			/* compute scaling factor   */
  *vec = mul_vs(*vec, rsc);			/* rescale to radius given  */
}