/* Read in  particle's phase space points */

void read_particles(int *n0){
  int newpart,ii,ires,code,alpha,ipart;
  double pmag,mass,px,py,pz,e;
  FILE *fptr;
  fptr=fopen("temp_data/explode.dat","r");
  *n0=0;
  while(feof(fptr)==0){
    *n0=*n0+1;
    if(firstemptypart==npartmax){
      printf("npartmax is too small\n");
      exit(1);
    }
    newpart=firstemptypart;
    fscanf(fptr,"%d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf",
	   &ipart,&code,&px,&py,&pz,&e,&mass,
	   &partinfo[newpart].r[1],&partinfo[newpart].r[2],
	   &partinfo[newpart].r[3],&partinfo[newpart].r[0]);

    partinfo[newpart].p[1]=1000.0*px;
    partinfo[newpart].p[2]=1000.0*py;
    partinfo[newpart].p[3]=1000.0*pz;
    partinfo[newpart].p[0]=1000.0*e;
    mass=1000.0*mass;

    pmag=pow(partinfo[newpart].p[1],2)+pow(partinfo[newpart].p[2],2)
      +pow(partinfo[newpart].p[3],2);
    if(pmag<10.0) {
      printf("peculiar pmag=%g   newpart=%d\n",pmag,newpart);
      printf("%g %g %g %g     %g %g %g\n",partinfo[newpart].r[0],
	     partinfo[newpart].r[1],partinfo[newpart].r[2],
	     partinfo[newpart].r[3],partinfo[newpart].p[1],
	     partinfo[newpart].p[2],partinfo[newpart].p[3]);
    }
    pmag=sqrt(pmag);
    for(ires=1;ires<=nres;ires++){
      if(code==resinfo[ires].code){
	ii=ires;
	goto CODE_MATCH;
      }
    }
    printf("can't match code!\n");
    printf("code=%d  n0=%d\n",code,*n0);
    exit(1);

  CODE_MATCH:

    partinfo[newpart].p[0]=sqrt(pmag*pmag+resinfo[ii].mass
				*resinfo[ii].mass);
    partinfo[newpart].mass=sqrt(fabs(partinfo[newpart].p[0]
				     *partinfo[newpart].p[0]-
				     pmag*pmag));
    for(alpha=0;alpha<4;alpha++) partinfo[newpart].b[alpha]=0.0;
    partinfo[newpart].ires=ii;
    firstemptypart=partinfo[firstemptypart].next;
    partinfo[newpart].next=firstpart;
    firstpart=newpart;
    npart=npart+1;
  }
  printf("Initial Particles Created, npart=%d n0=%d\n",npart,*n0);
}

/* ------------------------------------------------------------- */

/* Write statements */ 

void write_statements(int netbaryonscat,
		      int *nelas,int *nfusion,int *nbaryonscat){
  double et,rapidity,dnpdy=0.0,sigmabar=0.0;
  int nbaryons=0,partcount,ipart,code;
  int ires,i,unity=1;
  FILE *fptr,*fptrscatrate;
  fptr=fopen("results/phasespace.dat","a");
  fprintf(fptr,"%6d  %9d  %9g  %9g\n",event_number,npart,0.0,0.0);

  ipart=firstpart;
  partcount=0;
  et=0.0;
  while(ipart!=-1){
    partcount=partcount+1;
    ires=partinfo[ipart].ires;
    code=resinfo[ires].code;
    rapidity=0.5*log((partinfo[ipart].p[0]+partinfo[ipart].p[3])/
		     (partinfo[ipart].p[0]-partinfo[ipart].p[3]));

    fprintf(fptr,"%9d  %9d  %9g  %9g  %9g  %9g  %9g  %9g  %9g  %9g  %9g\n",
	    partcount,code,
	    partinfo[ipart].p[1]/1000.0,partinfo[ipart].p[2]/1000.0,
	    partinfo[ipart].p[3]/1000.0,partinfo[ipart].p[0]/1000.0,
	    resinfo[ires].mass/1000.0,
	    partinfo[ipart].r[1],partinfo[ipart].r[2],partinfo[ipart].r[3],
	    partinfo[ipart].r[0]);

    if(code!=22) sigmabar=sigmabar+rapidity*rapidity;

    if(fabs(rapidity)<0.5){
      et=et+sqrt(pow(resinfo[ires].mass,2)+pow(partinfo[ipart].p[1],2)
		 +pow(partinfo[ipart].p[2],2));
      if(resinfo[ires].baryon!=0) et=et-939.0*resinfo[ires].baryon;
      if(code==2112||code==2212) dnpdy=dnpdy+0.5;
    }

    if(abs(code)==2112||abs(code)==2212||abs(code)==3334||abs(code)==3122||
       abs(code)==3222||abs(code)==3212||abs(code)==3112||abs(code)==3322||
       abs(code)==3312) nbaryons=nbaryons+abs(code)/code;

    ipart=partinfo[ipart].next;
  }

  fclose(fptr);
  sigmabar=sigmabar/(double)partcount;
  et=et/1000.0;

  if(partcount!=npart) printf("npart=%d, should be partcount=%d\n",
			      npart,partcount);
  printf("There are %d baryons, and there were %d baryonic scatterings\n",
	 nbaryons,netbaryonscat);
  fptrscatrate=fopen("temp_data/scatvstime.dat","w");
  fprintf(fptrscatrate,"%g %g %g %d\n",et,dnpdy,sigmabar,unity);
  printf("et, dnpdy and sigma_y are: %g %g %g\n",et,dnpdy,sigmabar);
  for(i=0;i<50;i++) fprintf(fptrscatrate,"%3d  %4d %4d %4d %4d\n",
			    i,nelas[i],nfusion[i],nelas[i]+nfusion[i],
			    nbaryonscat[i]);
  fclose(fptrscatrate);
}