#include "cpgplot.h"

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

void main()
{
  int nmax=10000;
  float x[nmax], y[nmax];
  double pi, x1, x2, y1, y2, xc, dx, dy, xnew, ynew;
  int npar, nstep, seed;
  unsigned int i, istep;
  double a1, a, theta;
  double h;
  char STR[8];
  int L;

  printf("Input: Number of particles:\n");
  scanf("%u", &npar);
  printf("Input: Number of steps:\n");
  scanf("%u", &nstep);
  printf("Input: seed\n");
  scanf("%u", &seed);
  printf("Input: h\n");
  scanf("%lf", &h);

  pi= M_PI;

/* 初期条件の設定 */
  x1=0.0;
  x2=2.0;
  y1=0.0;
  y2=1.0;
  xc=1.0;

  if(cpgbeg(0, "?", 1, 1) != 1)
    exit(EXIT_FAILURE);
  cpgask(1);

/*  cpgqinf('HARDCOPY', &STR, L); */
  
  cpgpage();
  cpgvstd();
  cpgwnad(0.0, 2.0, 0.0, 1.0);

/*   seed=1; */
  srand(seed);

  for(i=0; i<npar; i++) {
    a1=rand();
    x[i]=a1/(RAND_MAX+1.0);
    a1=rand();
    y[i]=a1/(RAND_MAX+1.0);
  }

/* ランダムウォーク */
  for(istep=0; istep<nstep; istep++) {
    printf("istep=%u\n", istep);
    cpgbbuf();
    cpgeras();
    cpgsci(1);
    cpgbox("bcnst", 0.0, 0, "bcnst", 0.0, 0);

    for(i=0; i<npar; i++) {
      a1=rand();
      a=a1/(RAND_MAX+1.0);
      theta=2*pi*a;
      dx = h * cos(theta);
      dy = h * sin(theta);
      xnew = x[i]+dx;
      ynew = y[i]+dy;
      if(xnew < x1) xnew=2*x1-dx-x[i];
      if(ynew < y1) ynew=2*y1-dy-y[i];
      if(xnew > x2) xnew=2*x2-dx-x[i]; 
      if(ynew > y2) ynew=2*y2-dy-y[i];
      x[i]=xnew;
      y[i]=ynew;
    }

    cpgsci(3);
    cpgslw(5);
    cpgpt(npar, x, y, -1);
    cpgslw(1);
    sprintf(STR, "%u", istep);
    cpgmtxt("T", 2.0, 0.0, 0.0, STR);
    cpgebuf(); 
  }
  
  cpgend();

  return;
}