/* Monte Carlo that creates data for Thompson atom scattering */
/* gcc thompsonfig3.c -L/usr/X11R6/lib -lplot -lXaw -lXmu -lXt -lSM -lICE -lXext -lX11 -lm   */
#include <stdlib.h>
#include<sys/types.h>
#include<unistd.h>
#include <math.h>
#include<stdio.h>
#include<strings.h>
#include<plot.h>
#define PI 3.1415926
int pl_handle;
FILE *fptr;
pid_t getpid(void);
main(int argc, char *argv[]){
    plPlotter *plotter;
       plPlotterParams *plotter_params;
    float m1,m2,test,max,b,theta,R_cl, E,E0,num,den, dth;
    int n, m, p, crap,det, type;
    double long  counter[181];
    char label[10];

    
    /* seed the random number generator, establish MAX=2^{31} */
    srand((unsigned int)getpid());
    max=pow(2.0,31.0);
    /* constants and the energy */
    if(argc!=3){
	printf("./scatter E (eV) Type (0=Thompson, 1=Rutherford, 2=crystalball)\n");
	exit(0);
    }
    E=atof(argv[1]);
    type=atoi(argv[2]);
    E0=1137.6;
    dth=0.05; /* width of a bin */
for(n=0;n<=180;n++){
	counter[n]=0;
    }


 pl_handle = pl_newpl ("ps", stdin, stdout, stderr);
  pl_selectpl (pl_handle);
  if (pl_openpl () < 0)         /* open Plotter */
    {
      fprintf (stderr, "Couldn't open Plotter\n");
      return 1;
    }
  pl_fspace (-1, -1, 10, 10);
  pl_fline(-1.0, 0.0, 9.5, 0.0);
  pl_fline(0.0, -1.0, 0.0, 9.5);
  pl_ffontsize(0.25);

    for(n=0;n<5000000;n++){
	/* here is the Monte Carlo step */
    generate:{
	m1=2.0*(float)rand()/max-1.0;
	m2=2.0*(float)rand()/max-1.0;
	b=sqrt(m1*m1+m2*m2);
	if(b>1.0) goto generate;
    }


	/* now create the event */
	if(type==0){  /* Thompson */
	num=2.0*b*sqrt(1.0-b*b);
	den=(E/E0)-(1.0-2.0*b*b);
	theta=atan2(num,den);
	/* print out detector logging event */
	det=floor(180.0*theta/PI);
	/*printf("%d\n", det); */
	counter[det]=counter[det]+1;
	}
	if(type==1){ /*Rutherford*/
	    theta=2.0*atan2(E0,2.0*E*b);
	    det=floor(180.0*theta/PI);
	/*printf("%d\n", det); */
	counter[det]=counter[det]+1;
	}

	if(type==2){ /*crystalball, V_0=2*E_0*/
	    theta=PI-2.0*asin(b);
	    if(b*b<1.0-2.0*E0/E)
		theta=theta-2.0*acos(sqrt(E/(E-2.0*E0))*b);
	    det=floor(180.0*theta/PI);
	/*printf("%d\n", det); */
	counter[det]=counter[det]+1;
	}

    }
	
	    pl_erase();
	    pl_fline(-0.9, 0.0, 9.5, 0.0);
	    pl_fline(0.0, -0.9, 0.0, 9.5);
	    pl_linemod("dotted");
	    for(p=1;p<8;p++){
		pl_fline(0.0, (double)p, 9.5, (double)p);
		pl_fmove(-0.5, (double)p);
		sprintf(label, "10\\sp%d\\ep", p);
		pl_label(label);
	    }
	    pl_textangle(90);
	    pl_fmove(-0.6, 4.0);
	    pl_label("Counts, N");
	    pl_textangle(0);
	    pl_linemod("solid");
	    for(p=1;p<=18;p++){
		pl_fline(10.0*(double)p*0.05, -0.2, 10.0*(double)p*0.05,0.0);
	    }
	    pl_fmove(4.45, -0.5);
	    pl_label("90");
	    pl_fmove(8.9, -0.5);
	    pl_label("180");
	    pl_fmove(5.0, -0.9);
	    pl_label("\\*H, degrees");
	    pl_filltype(1);
	    pl_fillcolorname("gray90");
	    sprintf(label, "E = %f eV", E);
	    pl_fmove(0.6, 9.0);
	    pl_label(label);
	     for(m=0;m<=180;m++){
		if(counter[m]>=1)
		pl_fbox((double)m*dth-0.5*dth, 0.0, (double)m*dth+0.5*dth, log(fabs((double)counter[m]))/2.302585); 
		    }
       


 pl_closepl ();
  pl_selectpl (0);
  pl_deletepl (pl_handle);

   
}