int main(){
  int ievent,n0,nevents,oldfirstevent;
  int alpha,beta,i,j,kk,kf[4],kkk,ipart,ipartcheck,prevpart,trun;
  int nextpart,partcount,ii;
  int firstnewevent,oldfirstemptyevent,newpart,inewevent,nnewevents,test,ires;
  int scatcode,newscatcode,code,nextevent,lastevent;
  int nelas[50]={0},nfusion[50]={0},nbaryonscat[50]={0};
  int netelas=0,netfusion=0,netbaryonscat=0,eventsprocessed=0;
  int itau;
  double t_event,tau_event,delt,pmag,m2;
  FILE *fptr;

#ifdef HOLLYWOOD
  int imovie,moviecode,nframes=56,movietest;
  double tmovie[100],deltmovie=0.5;
  double xmovie,ymovie,zmovie;
  FILE *fptrmovie;
  fptrmovie=fopen("results/xyz.dat","w+");
  for(imovie=0;imovie<nframes;imovie++) tmovie[imovie]=1.0+deltmovie*imovie;
  imovie=0;
  fprintf(fptrmovie,"158.0\n");
#endif

#ifdef LOCAL_TEMP_FIND
  int i_localtemp=0;
  int localtemp_tau[10]={6,8,10,12,14,16,18,20,22,24};
#endif

  /* INITIALIZE PROGRAM */
  printf("What is the event number?\n");
  scanf("%d",&event_number);
  trun=time(NULL);
  IDUM=-trun;
  resonanceinit();

  /*-------------------------------------------------------------------*/
  /* CREATE EMPTY LINKED LIST OF EVENTS */
  /*-------------------------------------------------------------------*/
  firstevent=-1;
  firstnewevent=-1;
  firstemptyevent=0;
  for(ievent=1;ievent<=neventmax;ievent++) {
    eventinfo[ievent-1].next=ievent;
  }
  /* LINKED EVENT LIST FINISHED */

  /*-------------------------------------------------------------------*/  
  /* CREATE EMPTY LINKED LIST OF PARTICLES */
  /*-------------------------------------------------------------------*/
  firstpart=-1;
  firstemptypart=0;
  for(ipart=1;ipart<=npartmax;ipart++) {
    partinfo[ipart-1].next=ipart;
  }
  /* LINKED PARTICLE LIST FINISHED */

  /*-------------------------------------------------------------------*/
  /* READ PHASE SPACE INFORMATION ABOUT PARTICLES, N0 IS THE NUMBER READ IN */
  /*-------------------------------------------------------------------*/
  read_particles(&n0);
  /* PARTICLES READ IN */

  /*-------------------------------------------------------------------*/
  /* TEST FOR ORIGINAL EVENTS AND ADD TO THE EVENT LIST */
  /*-------------------------------------------------------------------*/
  nnewevents=0;
  for(i=0;i<n0;i++){
    for(j=i;j<n0;j++){
      if(i!=j) {
	test=colfind(i,j,&scatcode,&t_event,&tau_event,&delt);
      }
      else{
	test=decayfind(i,&scatcode,&t_event,&tau_event,&delt);
      }
      if(test==1){
	nnewevents=nnewevents+1;
	addevent(&firstnewevent,t_event,tau_event,delt,scatcode,i,j);
	//printf("Event added %d\n",scatcode);
	n_newevent[nnewevents]=firstnewevent;
#ifdef tau_order
	time_newevent[nnewevents]=tau_event;
#else
	time_newevent[nnewevents]=t_event;	
#endif
      }
    }
  }

  /* SORT THE NEW EVENTS BY TIME */
  if(nnewevents>1) shell_sort(nnewevents,time_newevent,n_newevent);
  /* ADD THE NEW EVENTS TO THE EVENT LIST */
  ii=1;
  firstevent=n_newevent[ii];
  ievent=firstevent;
  for(ii=2;ii<=nnewevents;ii++){
    nextevent=n_newevent[ii];
    eventinfo[ievent].next=nextevent;
    ievent=nextevent;
  }
  eventinfo[nextevent].next=-1;

  printf("******  Initialization Completed, %d future events  *******\n",
	 nnewevents);  

  nevents=nnewevents;
  if(nevents>=neventmax){
    printf("Make neventmax parameter larger!\n");
    exit(1);
  }
  /*-------------------------------------------------------------------*/
  /* INITIALIZATION COMPLETED */
  /*-------------------------------------------------------------------*/

  while (nevents!=0) {
  /*-------------------------------------------------------------------*/
  /* PERFORM EVENTS ONE AT A TIME */
  /*-------------------------------------------------------------------*/

#ifdef CHAIN_TESTS
    event_chain_test(nevents);
#endif

    nevents=nevents-1;
    i=eventinfo[firstevent].i;
    j=eventinfo[firstevent].j;
    t_event=eventinfo[firstevent].time;
    tau_event=eventinfo[firstevent].tau;
    itau=(int)fabs(tau_event);

#ifdef LOCAL_TEMP_FIND
    if(itau>=localtemp_tau[i_localtemp]&&i_localtemp<10){
      localtemp(localtemp_tau[i_localtemp]);
      i_localtemp=i_localtemp+1;
    }
#endif

    scatcode=eventinfo[firstevent].scatcode;
    if(scatcode>0) {
      resonancemaker(firstevent,kf);
      if(itau<50){
	netfusion=netfusion+1;
	nfusion[itau]=nfusion[itau]+1;
	if(fabs(resinfo[partinfo[i].ires].code)>=1000){
	  netbaryonscat=netbaryonscat+1;
	  nbaryonscat[itau]=nbaryonscat[itau]+1;
	}
	if(fabs(resinfo[partinfo[j].ires].code)>=1000) {
	  netbaryonscat=netbaryonscat+1;
	  nbaryonscat[itau]=nbaryonscat[itau]+1;
	}
      }
    }
    if(scatcode<0) resonancedecay(firstevent,kf);
    if(scatcode==0) {
      elasticscatter(firstevent,kf);
      if(itau<50){
	netelas=netelas+1;
	nelas[itau]=nelas[itau]+1;
	if(fabs(resinfo[partinfo[i].ires].code)>=1000){ 
	  netbaryonscat=netbaryonscat+1;
	  nbaryonscat[itau]=nbaryonscat[itau]+1;
	}
	if(fabs(resinfo[partinfo[j].ires].code)>=1000){
	  netbaryonscat=netbaryonscat+1;
	  nbaryonscat[itau]=nbaryonscat[itau]+1;
	}
      }
    }
    /* DELETE THE EVENT JUST PERFORMED */
    oldfirstevent=firstevent;
    firstevent=eventinfo[firstevent].next;
    firstnewevent=-1;
    if(firstemptyevent==neventmax){
      printf("neventmax is too small\n");
      exit(1);
    }
    eventinfo[oldfirstevent].next=firstemptyevent;
    firstemptyevent=oldfirstevent;

#ifdef CHAIN_TESTS
    event_chain_test(nevents);
#endif

#ifdef CHAIN_TESTS
    part_chain_test();
#endif
    /* EVENT PERFORMED */
    
    /*-------------------------------------------------------------------*/
    /* DELETE EVENT WITH PARTICLES WHICH DISAPPEARED, FIND NEW EVENTS*/
    /*-------------------------------------------------------------------*/

    nnewevents=0;
    ipart=firstpart;
    prevpart=-1;
    while(ipart!=-1){
      nextpart=partinfo[ipart].next;
      /* FIRST SEE IF PARTICLE SHOULD BE DELETED */
      if(scatcode!=0&&(i==ipart||j==ipart)){
	if(prevpart!=-1) {
	  partinfo[prevpart].next=partinfo[ipart].next;
	}
	else{
	  firstpart=partinfo[ipart].next;
	}
	if(firstemptypart==npartmax){
	  printf("npartmax is too small\n");
	  exit(1);
	}
	partinfo[ipart].next=firstemptypart;
	firstemptypart=ipart;
	npart=npart-1;
	//printf("Deleted particle  ipart=%d\n",ipart);
      }
      /* OTHERWISE CHECK FOR NEW EVENTS (KF[0] IS THE NUMBER OF */
      /* CREATED PARTICLES */
      else{
	for(kk=1;kk<=kf[0];kk++){
	  test=0;
	  ipartcheck=1;
	  for(kkk=1;kkk<=kf[0];kkk++){
	    if(kf[kkk]==ipart) ipartcheck=0;
	  }
	  if(ipartcheck==1){
	    test=colfind(ipart,kf[kk],&newscatcode,&t_event,&tau_event,&delt);
	  }
	  if(ipart==kf[kk]){
	    test=decayfind(kf[kk],&newscatcode,&t_event,&tau_event,&delt);
	  }
	  if(test==1) {
	    nnewevents=nnewevents+1;
	    addevent(&firstnewevent,t_event,tau_event,delt,newscatcode,
		     ipart,kf[kk]);
	    //printf("Adding new event, scatcode=%d\n",newscatcode);
	    n_newevent[nnewevents]=firstnewevent;
#ifdef tau_order
	    time_newevent[nnewevents]=tau_event;
#else
	    time_newevent[nnewevents]=t_event;
#endif
	  }
	}
	prevpart=ipart;
      }
      ipart=nextpart;
    }

#ifdef CHAIN_TESTS
    part_chain_test();
#endif

    /*-------------------------------------------------------------------*/
    /* SORT THE NEW EVENTS BY TIME, */
    /* AND ADD THE NEW EVENTS TO THE NEW EVENT LIST */
    /*-------------------------------------------------------------------*/

    if(nnewevents>0){
      if(nnewevents>1) shell_sort(nnewevents,time_newevent,n_newevent);
      ii=1;
      firstnewevent=n_newevent[ii];
      inewevent=firstnewevent;
      nextevent=n_newevent[ii];
      for(ii=2;ii<=nnewevents;ii++){
	nextevent=n_newevent[ii];
	eventinfo[inewevent].next=nextevent;
	inewevent=nextevent;
      }
      eventinfo[nextevent].next=-1;

#ifdef CHAIN_TESTS
      event_chain_test(nevents);
#endif
    }
    /* NEW EVENTS FOUND AND PUT INTO LIST */

    /*-------------------------------------------------------------------*/
    /* THUMB THROUGH OLD EVENTS */
    /* ADD NEW EVENTS TO EVENT LIST, DELETE EVENTS WITH OLD PART.S */
    /*-------------------------------------------------------------------*/

    ievent=firstevent;
    lastevent=-1;
    nevents=nevents+nnewevents;
    if(nevents>=neventmax){
      printf("Make neventmax larger!\n");
      exit(1);
    }
    ii=1;
    inewevent=n_newevent[ii];

    do{

      /* FIRST CHECK TO SEE IF EVENTS NEED TO BE DELETED */
      while(ievent!=-1&&(eventinfo[ievent].i==i||eventinfo[ievent].j==j||
			 eventinfo[ievent].i==j||eventinfo[ievent].j==i)){
	if(firstemptyevent==neventmax){
	  printf("neventmax is too small\n");
	  exit(1);
	}
	oldfirstemptyevent=firstemptyevent;
	firstemptyevent=ievent;
	if(lastevent!=-1) {
	  eventinfo[lastevent].next=eventinfo[ievent].next;
	}
	else{
	  firstevent=eventinfo[ievent].next;
	}
	nevents=nevents-1;
	ievent=eventinfo[ievent].next;
	eventinfo[firstemptyevent].next=oldfirstemptyevent;
      }

      /* NEXT CHECK TO SEE WHICH EVENTS NEED TO BE INSERTED */
#ifdef tau_order
      while(nnewevents>0&&
	    (ievent==-1||time_newevent[ii]<eventinfo[ievent].tau)){
#else
	while(nnewevents>0&&
	    (ievent==-1||time_newevent[ii]<eventinfo[ievent].time)){
#endif
	if(lastevent>=0) {
	  eventinfo[lastevent].next=inewevent;
	}
	else{
	  firstevent=inewevent;
	}
	lastevent=inewevent;
	eventinfo[inewevent].next=ievent;
	firstnewevent=n_newevent[ii+1];
	ii=ii+1;
	nnewevents=nnewevents-1;
	inewevent=n_newevent[ii];
      }
      if(ievent!=-1) {
	lastevent=ievent;
	ievent=eventinfo[ievent].next;
      }
    } while(ievent!=-1||nnewevents>0);

    eventsprocessed=eventsprocessed+1;

#ifdef HOLLYWOOD
#include "hollywood.c"
#endif
  }

  /*-------------------------------------------------------------------*/
  /* CASCADE FINISHED  BEGIN WRITING RESULTS */
  /*-------------------------------------------------------------------*/

  printf("* %d events processed (%d elas. scatt.s, %d fusions), npart=%d *\n",
	 eventsprocessed,netelas,netfusion,npart);
  write_statements(netbaryonscat,nelas,nfusion,nbaryonscat);

  trun=time(NULL)-trun;
  printf("-------- Program Finishing Normally in %d seconds ---------\n",trun);

  return 0;
}