/*-------------------------------------------------------------

Copyright (C) 2000 Sebastian Will. 
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/


/************************************************************
 *
 *  program to determine the distribution of the 
 *  mean square distances between random walks
 *  either self avoiding or non self avoiding
 *
 ************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/* parameters */
#define LENGTH 30 /* chain length */
#define SAMPLES 100000 /* number of samples to generate */
#define INIT_RANDOM 0 /* initialize random number generator with this value */
#define SELF_AVOIDING 1 /* generate only self avoiding walks */
#define DISTANCE_MATRIX 1 /* compute mean square dist from distance matrix */

/* type definitions */
typedef int point_t[3];
typedef point_t walk_t[LENGTH];

typedef double dist_mat_t[LENGTH][LENGTH];
/* convention: distance matrix is symmetric, i.e.
   one could use triangular matrix.
   ==> We use only entries m[i][j], where i<j. */

/* ------------------------------------------------------------ */
/* general aux functions */

/* random integer */
inline
int rand_int(int n) {
  return random()%n;
}

inline
double square(double x) {
  return x*x;
}

/* ------------------------------------------------------------ */
/* functions for points */
void clear_point(point_t p) {
  int i;
  for (i=0; i<3; ++i) p[i] = 0;
}

void copy_point(point_t to, point_t from) {
  int i;
  for (i=0; i<3; ++i) to[i] = from[i];
}

int equal_points(point_t p, point_t q) {
  return (p[0]==q[0] &&
	  p[1]==q[1] &&
	  p[2]==q[2]);
}

/* Euclidian distance of two points */
double eucl_dist_points(point_t p, point_t q) {
  return sqrt(square(p[0]-q[0])+
	      square(p[1]-q[1])+
	      square(p[2]-q[2]));
}

/* ------------------------------------------------------------ */
/* functions for walks */
void print_walk(walk_t walk) {
  int i;
  for (i=0; i<LENGTH; ++i)
    printf("(%d,%d,%d) ",walk[i][0],walk[i][1],walk[i][2]);
  printf("\n");
}

/* generate point q as a random successor to point p,
   such that q and p are lattice neighbours */
void gen_random_succ(point_t p, point_t q) {
  int dim,add;
  dim = rand_int(3);     /* 0,1 or 2 */
  add = rand_int(2)*2-1; /* +1 or -1 */
  copy_point(q,p);
  q[dim] += add;
}

/* ------------------------------------------------------------ */
/* distance matrices */

/* compute the distance matrix d of walk w */
void gen_distance_matrix(dist_mat_t d, walk_t w) {
  int i,j;
  
  for (i=0; i<LENGTH; ++i)
    for (j=i+1; j<LENGTH; ++j)
      d[i][j] = eucl_dist_points(w[i],w[j]);
}

/* ------------------------------------------------------------ */
/* generate random walk */
void gen_walk(walk_t walk) {
  int i;
  
  clear_point(walk[0]);
  
  for (i=1; i<LENGTH; ++i) {
    gen_random_succ(walk[i-1],walk[i]);
  }
}

/* test if position p is already occupied by points of walk w 
   with index less than n */
int occupied(point_t p, walk_t w, int n) {
  int i;
  for (i=0;i<n;++i)
    if (equal_points(w[i],p)) return 1;

  return 0;
}

/* generate random self avoiding walk */
void gen_sa_walk(walk_t walk) {
  /* implementation of elementary simple sampling */
  int i;
  static point_t neighbours[6];
  int num;
  int dim,add;
  
  clear_point(walk[0]);

  do {
    for ( i=1 ; i<LENGTH ; i++ ) {
      /* compute list of unoccupied neighbour positions */
      
      num=0; /* num counts the number of unocc. neighbours, 
		i.e. the length of the neighbour list */
      for (dim=0; dim<3; ++dim)
	for (add=-1; add<=1; add=+2) {
	  copy_point(neighbours[num],walk[i-1]);
	  neighbours[num][dim] += add;
	  if (!occupied(neighbours[num],walk,i-1)) num++; 
	}
      
      if (num==0) /* case: no neighbours unoccupied */
	break; /* this breaks the for loop and thus restarts
		  the generation of a random walk */
      else
	/* select one neighbour from list */
	copy_point(walk[i],neighbours[rand_int(num)]);      
    }
    
  } while( i < LENGTH );
}

/* mean square distance of two walks */
double distance(walk_t p,walk_t q) {
  int i,j;
  double d = 0.0;
  
  /* distance of relative moves as in mail */
  for (i=1; i<LENGTH; ++i)
    for (j=0; j<3; ++j)
      d += square((p[i][j]-p[i-1][j])-(q[i][j]-q[i-1][j]));

  return d/LENGTH;
}


/* mean square distance of the distance matrices of two walks */
double distance_dist_mat(walk_t p,walk_t q) {
  int i,j;
  static dist_mat_t d_p;
  static dist_mat_t d_q;
  
  double d=0;

  gen_distance_matrix(d_p,p);
  gen_distance_matrix(d_q,q);
  
  for (i=0; i<LENGTH; ++i)
    for (j=i+1; j<LENGTH; ++j)
      d += square(d_p[i][j] - d_q[i][j]);

  return d/(LENGTH*(LENGTH-1)/2);
}

int
main() {
  int i,j;
  double dist;
  walk_t walk[2]; // the two walks p and q
    
  srandom(INIT_RANDOM);

  // foreach sample
  for (i=0;i<=SAMPLES;++i) {
  
    // generate walks
    for (j=0;j<2;++j) {
      if (SELF_AVOIDING)
	gen_sa_walk(walk[j]);
      else
	gen_walk(walk[j]);
      // print_walk(walk[j]);
    }
    
    // compute distance
    if (DISTANCE_MATRIX)
      dist = distance_dist_mat(walk[0],walk[1]);
    else
      dist = distance(walk[0],walk[1]);

    printf("%f\n",dist); 
  }
  
  exit(0);
}