/* Buffon needle program                       */
/* gcc -O4 -o buffon buffon.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM   -lICE -lXext -lX11 -lm */
/* The output reguires GNU plotutils-2.2 or above                   */
/* J. R. Schmidt 1999                                               */


#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<plot.h>
#include<unistd.h>
#include<sys/types.h>
#define PI 3.1415926
float ell, x0,y00,x1,y10,dx,dy,pi,theta,ind1, ind2, ind3, ind4;
int N,num,i,j;
int pl_handle;

main (int argc, char *argv[])
{
  /* Do some error checking to avoid segmentation faults */
 /* I hate segmentation faults.                          */

  if(argc < 2 || argc > 2){
    printf("usage: ./buffon N\n");
    exit(1);}
  N=atoi(argv[1]);
ell=50.0;
srand((unsigned int)getpid()); 

  pl_handle = pl_newpl ("X", stdin, stdout, stderr);
   pl_selectpl (pl_handle);
  if (pl_openpl () < 0)         /* open Plotter */
    {
      fprintf (stderr, "Couldn't open Plotter\n");
      return 1;
    }
  /* set scales for the size of the plotting window */
  pl_fspace (-0.1*ell, -0.1*ell, 5.1*ell, 5.1*ell);
  pl_pencolorname("blue");
  pl_fbox(0.0,0.0, 5.0*ell, 5.0*ell);
  pl_fline(0.0,ell, 5.0*ell, ell);
  pl_fline(0.0,2.0*ell, 5.0*ell, 2.0*ell);
  pl_fline(0.0,3.0*ell, 5.0*ell, 3.0*ell);
  pl_fline(0.0,4.0*ell, 5.0*ell, 4.0*ell);
  pl_pencolorname("black");
  num=0;
  for(i=0;i<N;i++){
      x0=5.0*ell*(float)rand()/(float)RAND_MAX;
      y00=5.0*ell*(float)rand()/(float)RAND_MAX;
      theta=2.0*PI*(float)rand()/(float)RAND_MAX;
      dx=0.99*ell*cos(theta);
      dy=0.99*ell*sin(theta);
      x1=x0+dx;
      y10=y00+dy;
      if(x1>5.0*ell || x1< 0.0)x1=x0-dx;
      if(y10>5.0*ell || y10< 0.0)y10=y00-dy;
      ind1=(y00-ell)*(y10-ell);
      ind2=(y00-2.0*ell)*(y10-2.0*ell);
      ind3=(y00-3.0*ell)*(y10-3.0*ell);
      ind4=(y00-4.0*ell)*(y10-4.0*ell);
      if(ind1*ind2*ind3*ind4<=0.0){
	  pl_pencolorname("red");
         pl_fline(x0,y00,x1,y10); 
	 num=num+1; 
	 pl_pencolorname("black");}
      else{        
	  pl_fline(x0,y00,x1,y10);}
  }
  pi=2.0*(float)N/(float)num;
  printf("The number of crossings is %d\n", num);
  printf("PI is then measured to be %f\n", pi);
       
  pl_closepl ();
  pl_selectpl (0);
  pl_deletepl (pl_handle);
}