int main(){

  double bound[nsurf+1][4],bprimea[nsurf+1][4],bprimeb[nsurf+1][4];
  double vbound[nsurf+1][4],vbar[4],ubar[4];
  double delr[4],delrprime[4],nxprime[nsurf+1],nzprime[nsurf+1];
  double deladelt[nsurf+1];
  double ecmcern,sigmay=0.0,lorentz,n_lab[4],u[4],g[4][4];
  double t0,t,tau,delt=0.5;
  double sthet,cthet,deltat,rapidity,xbar,zbar,check,init_zsize;
  double p[4],r[4],pprime[4],rprime[4];
  double gamma,length,delx,delz,ftry,dndv,deladeltmax,wtot;
  double tauz,vx,vz,randy,pprime3,rprime3;
  double vzmax,uxmax,uzmax,init_time;
  double uz,ux,umag;
  double biggestx,biggestz,vol0;
  double etot=0.0,ettot=0.0,etmax=0.0;
  double stopping,pmag,pperp,ryyy,pyyy,cphi,sphi,phi;
  double detdy=0.0,ecoll=0.0,ekin=0.0;
  double dnpdy=0.0,dnbdy=0.0,dnady=0.0;
  double mbeforetest,maftertest;
  int nevents;
  int isurf,it,jnew,npart=0,nsuccess,success,itry,ntry,nt;
  int alpha,beta,ncheck_fail=0;
  int nwrite=0,ntrysurf,mix,nbaryon=0,ncentral=0;
  FILE *fptroutput;

  quantum=0;
#ifdef NOCASCADE
  quantum=1;
#endif
  evaporate=1;
#ifdef NOEVAPORATE
  evaporate=0;
#endif

  IDUM=-time(NULL);
  /* IDUM=iseed; */
  printf("How many events?\n",nevents);
  scanf("%d",&nevents);
  fptroutput=fopen("temp_data/explode.dat","w");

  vzmax=tanh(rapmax);
  uzmax=vzmax/sqrt(1.0-vzmax*vzmax);
  uxmax=vxmax/sqrt(1.0-vxmax*vxmax);
  t0=(xinitsize/uxmax);
  init_time=sqrt(t0*t0+xinitsize*xinitsize)-t0;
  init_zsize=vzmax*init_time;
  printf("init_time=%g  init_zsize=%g  vzmax=%g t0=%g\n",init_time,init_zsize,vzmax,t0);

  dndv=resonanceinit();
  printf("dN/dV=%f\n",dndv);
  j=(int)(jmax*ran2());
  biggestx=xinitsize;
  biggestz=init_zsize;

  for(alpha=0;alpha<4;alpha++) {
    n_lab[alpha]=0.0;
    for(beta=0;beta<4;beta++){
      g[alpha][beta]=0;
      if(alpha==beta&&beta==0) g[alpha][beta]=1;
      if(alpha==beta&&beta!=0) g[alpha][beta]=-1;
    }
  }
  n_lab[0]=1.0;
  for(isurf=0;isurf<=nsurf;isurf++){
    bound[isurf][0]=init_time;
    bound[isurf][2]=0.0;
    uz=uzmax*sin((pi/2.0)*(double)(isurf)/(double)nsurf);
    uz=fabs(uz);
    ux=uxmax*sqrt(fabs(1.0-uz*uz/(uzmax*uzmax)));
    umag=sqrt(ux*ux+uz*uz);
    gamma=sqrt(1.0+umag*umag);
    vz=uz/gamma;
    vx=ux/gamma;
    bound[isurf][3]=init_time*vz;
    tauz=sqrt(fabs(init_time*init_time-pow(bound[isurf][3],2)));
    vx=vxmax*sqrt(fabs(1.0-vz*vz/(vzmax*vzmax)));
    bound[isurf][1]=init_time*vx*(tauz+t0)/tauz;
    vbound[isurf][0]=1.0;
    vbound[isurf][1]=(bound[isurf][1]/init_time)*tauz/(tauz+t0);
    vbound[isurf][2]=0.0;
    vbound[isurf][3]=bound[isurf][3]/init_time;
  }
  printf("t0=%f  init_zsize=%f\n",t0,init_zsize);
  for(isurf=0;isurf<=nsurf;isurf++){
    for(alpha=0;alpha<4;alpha++){
      bound[isurf][alpha]=bound[isurf][alpha]-.5*delt*vbound[isurf][alpha];
    }
  }
  /* Now begin stepping over time */
  nt=2000;
  for(it=0;it<nt;it++){
    /* Now calculate the weightings for each area element */
    wtot=0.0;
    deladeltmax=0.0;
    for(isurf=0;isurf<=nsurf;isurf++){
      for(alpha=0;alpha<4;alpha++){
	bound[isurf][alpha]=bound[isurf][alpha]+delt*vbound[isurf][alpha];
      }
      t=bound[0][0];
      if(isurf>0){
	xbar=.5*(bound[isurf][1]+bound[isurf-1][1]);
	zbar=.5*(bound[isurf][3]+bound[isurf-1][3]);
	tauz=(fabs(t*t-zbar*zbar));
	if(tauz<0.0) printf("oops tauz^2=%f\n",tauz);
	tauz=sqrt(fabs(tauz));
	if(tauz+t0<fabs(xbar)) printf("oops (tauz+t0)^2<x^2\n");
	tau=-t0+sqrt(-xbar*xbar+pow(tauz+t0,2));
	if(tau>taumax) {
	  deladelt[isurf]=0.0;
	}
	else{
	  gamma=0.0;
	  for(alpha=0;alpha<4;alpha++) {
	    vbar[alpha]=.5*(vbound[isurf-1][alpha]+vbound[isurf][alpha]);
	    gamma=gamma+vbar[alpha]*vbar[alpha]*g[alpha][alpha];
	    delr[alpha]=bound[isurf][alpha]-bound[isurf-1][alpha];
	  }
	  gamma=1.0/sqrt(gamma);
	  ubar[0]=gamma;
	  for(alpha=1;alpha<4;alpha++) ubar[alpha]=-gamma*vbar[alpha];
	  for(alpha=0;alpha<4;alpha++){
	    delrprime[alpha]=0.0;
	    for(beta=0;beta<4;beta++){
	      lorentz=(-((ubar[alpha]+n_lab[alpha])*
		       (ubar[beta]+n_lab[beta])/(1.0+ubar[0]))
		+2.0*ubar[alpha]*n_lab[beta]+g[alpha][beta])*g[beta][beta];
	      delrprime[alpha]=delrprime[alpha]+lorentz*delr[beta];
	    }
	  }
	  length=sqrt(delrprime[1]*delrprime[1]+delrprime[3]*delrprime[3]);
	  /* These are unit vectors perp. to surface in frame of matter */
	  nxprime[isurf]=delrprime[3]/length;
	  nzprime[isurf]=-delrprime[1]/length;
	  /* if(nxprime[isurf]<0||nzprime[isurf]<0) { */
	  /*  printf("Curious backwards nhat=(%g,%g),isurf=%d,time=%g\n", */
	  /*  nxprime[isurf],nzprime[isurf],isurf,bound[isurf][0]); */
	  /* } */
	  length=sqrt(delr[1]*delr[1]+delr[3]*delr[3]);
	  delx=delr[1]/length;
	  delz=delr[3]/length;
	  deladelt[isurf]=2.0*pi*delt*xbar*length
	    /(gamma*sqrt(1.0-pow(vbar[3]*delx-vbar[1]*delz,2)));
	  if(xbar>biggestx) biggestx=xbar;
	  if(zbar>biggestz) biggestz=zbar;
	  if(deladelt[isurf]>deladeltmax) deladeltmax=deladelt[isurf];
	  wtot=wtot+deladelt[isurf];
	}
      }
    }
    if(wtot<1.0e-20&&it>0) {
      printf("Evaporation finished  time=%7.3f\n",t);
      goto ENDEVAP;
    }
    /* Now generate particles */
    ftry=dndv*wtot*nevents;
    ntry=(int)ftry;
    if(ran2()<(ftry-(double)ntry)) ntry=ntry+1;
    nsuccess=0;
    ntrysurf=0;
    for(itry=0;itry<ntry;itry++){
      for(mix=0;mix<2;mix++){
	jnew=(int)(jmax*ran2());
	if(jweight[jnew]>jweight[j]
	   ||ran2()<jweight[jnew]/jweight[j]) {
	  j=jnew;
	}
      }
      ires=iresfind[j];
      make_particle(pprime);
      pprime[3]=fabs(pprime[3]);
      if(pprime[3]/pprime[0]>ran2()){
	if(evaporate!=0) npart=npart+1;
	nsuccess=nsuccess+1;
      newsurfelement: 
	isurf=1+(int)(ran2()*nsurf);
	ntrysurf=ntrysurf+1;
	if(ran2()>deladelt[isurf]/deladeltmax) goto newsurfelement;
	pprime3=-pprime[1]*nxprime[isurf]+pprime[3]*nzprime[isurf];
	pprime[1]=pprime[1]*nzprime[isurf]+pprime[3]*nxprime[isurf];
	pprime[3]=pprime3;
	gamma=0.0;
	for(alpha=0;alpha<4;alpha++) {
	  vbar[alpha]=.5*(vbound[isurf-1][alpha]+vbound[isurf][alpha]);
	  gamma=gamma+vbar[alpha]*vbar[alpha]*g[alpha][alpha];
	  delr[alpha]=bound[isurf][alpha]-bound[isurf-1][alpha];
	}
	gamma=1.0/sqrt(gamma);
	ubar[0]=gamma;
	for(alpha=1;alpha<4;alpha++) ubar[alpha]=gamma*vbar[alpha];
	for(alpha=0;alpha<4;alpha++){
	  p[alpha]=0.0;
	  for(beta=0;beta<4;beta++){
	      lorentz=(-((ubar[alpha]+n_lab[alpha])*
		       (ubar[beta]+n_lab[beta])/(1.0+ubar[0]))
		+2.0*ubar[alpha]*n_lab[beta]+g[alpha][beta])*g[beta][beta];
	    p[alpha]=p[alpha]+lorentz*pprime[beta];
	  }
	}
	deltat=delt*(.5-ran2());
	r[0]=t+deltat;
	randy=ran2();
	r[1]=bound[isurf-1][1]*randy+bound[isurf][1]*(1.0-randy);
	r[2]=0.0;
	r[3]=bound[isurf-1][3]*randy+bound[isurf][3]*(1.0-randy);
	vbar[1]=vbound[isurf-1][1]*randy+vbound[isurf][1]*(1.0-randy);
	vbar[3]=vbound[isurf-1][3]*randy+vbound[isurf][3]*(1.0-randy);
	r[1]=r[1]+deltat*vbar[1];
	r[3]=r[3]+deltat*vbar[3];
	rapidity=.5*log((1.0+p[3]/p[0])/(1.0-p[3]/p[0]));
	/* Either rotate s.t. py=0 or randomly */
	phi=2.0*pi*ran2();
	cphi=cos(phi);
	sphi=sin(phi);
	pperp=sqrt(p[1]*p[1]+p[2]*p[2]);
	//cphi=p[1]/pperp;
	//sphi=-p[2]/pperp;
	ryyy=r[2]*cphi+r[1]*sphi;
	r[1]=r[1]*cphi-r[2]*sphi;
	r[2]=ryyy;
	pyyy=p[2]*cphi+p[1]*sphi;
	p[1]=p[1]*cphi-p[2]*sphi;
	p[2]=pyyy;
	pmag=sqrt(p[1]*p[1]+p[2]*p[2]+p[3]*p[3]);
	pperp=sqrt(p[1]*p[1]+p[2]*p[2]);
	etot=etot+p[0];
	ettot=ettot+p[0]*pperp/pmag;
	etmax=etmax+p[0]*pi/4.0;
	if(ran2()<.5){
	  r[3]=-r[3];
	  p[3]=-p[3];
	}
	if(fabs(baryon[ires]>.9)) nbaryon=nbaryon+(int)baryon[ires];
	if(evaporate!=0){
	  sigmay=sigmay+rapidity*rapidity;
	  if(nwrite!=0) fprintf(fptroutput,"\n");
	  fprintf(fptroutput,
		  "%9d  %9d  %9g %9g  %9g  %9g  %9g  %9g  %9g  %9g  %9g",
		  nwrite+1,itypefind[ires],
		  p[1]/1000.0,p[2]/1000.0,p[3]/1000.0,p[0]/1000.0,
		  mass[ires]/1000.0,r[1],r[2],r[3],r[0]);
	  nwrite=nwrite+1;
	  if(fabs(rapidity)<0.5){
	    if(baryon[ires]>.9) dnbdy=dnbdy+1.0;
	    if(itypefind[ires]==2212) dnpdy=dnpdy+1.0;
	    if(itypefind[ires]==2224) dnpdy=dnpdy+1.0;
	    if(itypefind[ires]==2214) dnpdy=dnpdy+0.6666667;
	    if(itypefind[ires]==2114) dnpdy=dnpdy+0.3333333;
	    if(baryon[ires]<-.9) dnady=dnady+1.0;
	    detdy=detdy+sqrt(p[0]*p[0]-p[3]*p[3]);
	    ecoll=ecoll+p[0]-pprime[0];
	    ekin=ekin+p[0]-mass[ires];
	    ncentral=ncentral+1;
	  }
	}
      }
    }
  }
ENDEVAP:
  printf("Total particles evaporated per event: %f\n",(double)npart/nevents);
  printf("biggest x,z of surface: %7.2f,%7.2f\n",biggestx,biggestz);
  biggestx=biggestx+delt*vxmax;
  biggestz=biggestz+delt*vzmax;
  /* test for vzmax */
  r[0]=sqrt(taumax*taumax+biggestz*biggestz);

  /* Now we begin inner disassociation */
  vol0=pi*biggestx*biggestx*2.0*biggestz;
  ftry=dndv*vol0*nevents;
  ntry=(int)ftry;
  if(ran2()<(ftry-(double)ntry)) ntry=ntry+1;
  nsuccess=0;
  for(itry=0;itry<ntry;itry++){
    r[1]=sqrt(biggestx*biggestx*ran2());
    r[2]=0.0;
    r[3]=biggestz*ran2();
    r[0]=r[3]*r[3]+pow( sqrt(pow(r[1],2)+pow(taumax+t0,2)) -t0,2);
    if(r[0]<0.0){
      printf("t^2<0, x=%f  z=%f t^2=%f\n",r[1],r[3],r[0]);
      exit(1);
    }
    r[0]=sqrt(r[0]);
    vz=r[3]/r[0];
    tauz=sqrt(r[0]*r[0]-r[3]*r[3]);
    vx=r[1]*tauz/((tauz+t0)*r[0]);
    gamma=1.0/sqrt(1.0-vx*vx-vz*vz);
    ux=gamma*vx;
    uz=gamma*vz;
    check=pow((ux/uxmax),2)+pow((uz/uzmax),2);
    if(check>1.0) ncheck_fail=ncheck_fail+1;
    if(check<=1.0&&ran2()<1.0/gamma){
      jnew=(int)(jmax*ran2());
      if(jweight[jnew]>jweight[j]
	 ||ran2()<jweight[jnew]/jweight[j]) {
	j=jnew;
      }
      ires=iresfind[j];
      make_particle(pprime);
      npart=npart+1;
      nsuccess=nsuccess+1;
      u[0]=gamma;
      u[1]=gamma*vx;
      u[2]=0.0;
      u[3]=gamma*vz;
      mbeforetest=0.0;maftertest=0.0;
      for(alpha=0;alpha<4;alpha++){
	mbeforetest=mbeforetest+pprime[alpha]*pprime[alpha]*g[alpha][alpha];
	p[alpha]=0.0;
	for(beta=0;beta<4;beta++){
	  lorentz=(-(u[alpha]+n_lab[alpha])
		   *(u[beta]+n_lab[beta])/(1.0+u[0])
		   +2.0*u[alpha]*n_lab[beta]+g[alpha][beta])*g[beta][beta];
	  p[alpha]=p[alpha]+lorentz*pprime[beta];
	  }
	maftertest=maftertest+pprime[alpha]*pprime[alpha]*g[alpha][alpha];
      }
      if(fabs(maftertest-mbeforetest)>1.0){
	printf("Lorentz Trans failed\n");
	exit(1);
      }
      rapidity=0.5*log((1.0+p[3]/p[0])/(1.0-p[3]/p[0]));
      /* Either rotate s.t. py=0, or randomly */
      phi=2.0*pi*ran2();
      cphi=cos(phi);
      sphi=sin(phi);
      pperp=sqrt(p[1]*p[1]+p[2]*p[2]);
      //cphi=p[1]/pperp;
      //sphi=-p[2]/pperp;
      ryyy=r[2]*cphi+r[1]*sphi;
      r[1]=r[1]*cphi-r[2]*sphi;
      r[2]=ryyy;
      pyyy=p[2]*cphi+p[1]*sphi;
      p[1]=p[1]*cphi-p[2]*sphi;
      p[2]=pyyy;
      pmag=sqrt(p[1]*p[1]+p[2]*p[2]+p[3]*p[3]);
      etot=etot+p[0];
      ettot=ettot+p[0]*pperp/pmag;
      etmax=etmax+p[0]*pi/4.0;
      if(ran2()<.5){
	  r[3]=-r[3];
	  p[3]=-p[3];
      }
      if(fabs(baryon[ires])>0.9) nbaryon=nbaryon+(int)baryon[ires];
      sigmay=sigmay+rapidity*rapidity;
      if(nwrite!=0) fprintf(fptroutput,"\n");
      fprintf(fptroutput,
	      "%9d  %9d  %9g %9g  %9g  %9g  %9g  %9g  %9g  %9g  %9g",
	      nwrite+1,itypefind[ires],
	      p[1]/1000.0,p[2]/1000.0,p[3]/1000.0,p[0]/1000.0,
	      mass[ires]/1000.0,r[1],r[2],r[3],r[0]);
      nwrite=nwrite+1;
      if(fabs(rapidity)<0.5){
	if(baryon[ires]>.9) dnbdy=dnbdy+1.0;
	if(itypefind[ires]==2212) dnpdy=dnpdy+1.0;
	if(itypefind[ires]==2224) dnpdy=dnpdy+1.0;
	if(itypefind[ires]==2214) dnpdy=dnpdy+0.66666667;
	if(itypefind[ires]==2114) dnpdy=dnpdy+0.33333333;
	if(baryon[ires]<-.9) dnady=dnady+1.0;
	detdy=detdy+sqrt(p[0]*p[0]-p[3]*p[3]);
	ecoll=ecoll+(p[0]-pprime[0]);
	ekin=ekin+p[0]-mass[ires];
	ncentral=ncentral+1;
      }
    }
  }
  fclose(fptroutput);
  printf("No. of tries= %d No. of Successes = %d\n",ntry,nsuccess); 
  printf("No. of failures due to being outside ellipse: %d\n",ncheck_fail);
  printf("Total Particles generated per event: %f\n",(double)npart/nevents);
  printf("Collective/Kinetic Energy per particle at mid.rap.= %g / %g\n",
	 ecoll/ncentral,ekin/ncentral);
  printf("Net baryons generated per event: %f\n",(double)nbaryon/nevents);
  printf("baryon dndy= %f\n",dnbdy/nevents);
  printf("antibaryon dndy= %f\n",dnady/nevents);
  printf(" ___________________________\n");
  printf("| proton dndy = %6.1f      |\n",dnpdy/nevents);
  printf("| dEt/dy(GeV) = %6.1f      |\n",
	 (double)(detdy-938.3*dnbdy+2.0*938.3*dnady)/(1000.0*nevents));
  ecmcern=sqrt(2.0*938.0*158.0*1000.0+4.0*938.0*938.0);
  printf("| %% of Beam E = %6.1f%%     |\n",
	 100.0*etot/(nevents*ecmcern*416.0));
  printf("|    sigma_y  = %6.2f      |\n",sqrt(sigmay/nwrite));
  printf("|___________________________|\n");
  return 0;
}