int colfind(int i,int j,int *scatcode,double *tevent,double *tau_event,double *delt){
  double bsquared,delr2,pidotpj,pidotdelr,pjdotdelr,ptot2,m2i,m2j;
  double a,a_i,a_j,denom,t_i,t_j;
  double zcoll,mcheck;
  double mass_i,mass_j,mass_r,sigpar,b2tot,delb2,width;
  double *p_i,*p_j,*r_i,*r_j;
  int sigptr;
  int mu,test,ires_r;
  test=0;
#ifdef NO_COLLISIONS
  goto GIVE_IT_UP;
#endif
  p_i=partinfo[i].p;
  p_j=partinfo[j].p;
  r_i=partinfo[i].r;
  r_j=partinfo[j].r;
  a=r_j[0]-r_i[0];
  pidotdelr=p_i[0]*a;
  pjdotdelr=p_j[0]*a;
  delr2=a*a;
  pidotpj=p_i[0]*p_j[0];
  m2i=p_i[0]*p_i[0];
  m2j=p_j[0]*p_j[0];
  for(mu=1;mu<4;mu++){
    m2i=m2i-p_i[mu]*p_i[mu];
    m2j=m2j-p_j[mu]*p_j[mu];
    a=r_j[mu]-r_i[mu];
    pidotdelr=pidotdelr-p_i[mu]*a;
    pjdotdelr=pjdotdelr-p_j[mu]*a;
    delr2=delr2-a*a;
    pidotpj=pidotpj-p_i[mu]*p_j[mu];
  }
  t_j=(m2i*pjdotdelr-pidotpj*pidotdelr);
  if(t_j<0.0) goto GIVE_IT_UP;
  t_i=(pidotpj*pjdotdelr-m2j*pidotdelr);
  if(t_i<0.0) goto GIVE_IT_UP;
  denom=pidotpj*pidotpj-m2i*m2j;
#ifdef EXTRA_TESTS
  if(denom<0.0) {
    printf("denom <0 I got it wrong, and denom can be negative\n");
    exit(1);
  }
#endif
  t_i=t_i/denom;
  t_j=t_j/denom;
#ifdef TURN_OFF_COLLS_AFTER
  *tevent=0.5*(r_i[0]+r_j[0]+t_i*p_i[0]+t_j*p_j[0]);
  zcoll=0.5*(r_i[3]+r_j[3]);
  *tau_event=sqrt(fabs(*tevent**tevent-zcoll*zcoll));
  if(*tau_event>max_coll_tau) goto GIVE_IT_UP;
#endif
  bsquared=-delr2-t_j*pjdotdelr+t_i*pidotdelr;
  if(bsquared>biggestb2[partinfo[i].ires][partinfo[j].ires]) 
    goto GIVE_IT_UP;

#ifdef EXTRA_TESTS
  if(bsquared<=0.0) {
    printf("DANGER: bsquared=%f\n",bsquared);
    printf("i=%d   j=%d    m2i=%g   m2j=%g\n",i,j,m2i,m2j);
    printf("i p: %-08.1e %-08.1e %-08.1e %-08.1e r: %-08.1e %-08.1e %-08.1e %-08.1e\n",
	   p_i[0],p_i[1],p_i[2],p_i[3],
	   r_i[0],r_i[1],r_i[2],r_i[3]);    
    printf("j p: %-08.1e %-08.1e %-08.1e %-08.1e r: %-08.1e %-08.1e %-08.1e %-08.1e\n",
	   p_j[0],p_j[1],p_j[2],p_j[3],
	   r_j[0],r_j[1],r_j[2],r_j[3]);
    exit(1);
  }
#endif
  b2tot=0.0;
  sigptr=firstsigptr[partinfo[i].ires][partinfo[j].ires];
  if(sigptr!=-1){
    ptot2=m2i+m2j+2.0*pidotpj;
#ifdef EXTRA_TESTS
    mass_i=sqrt(fabs(m2i));
    mass_j=sqrt(fabs(m2j));
    if(ptot2<pow((mass_i+mass_j),2)){
      printf("ptot2=%g  mass_r=%g  mass_i=%g  mass_j=%g\n",
	     ptot2,resinfo[siginfo[sigptr].ires0].mass,mass_i,mass_j);
      printf("p_i: %g %g %g %g    p_j: %g %g %g %g\n",
	     p_i[0],p_i[1],p_i[2],p_i[3],p_j[0],p_j[1],p_j[2],p_j[3]);
      printf("m2i=%g m2j=%g\n",m2i,m2j);
      printf("masses out of whack in colfind\n");
      exit(1);
    }
#endif
    while(sigptr!=-1){
      //mcheck=sqrt(m2i)+sqrt(m2j);
      mcheck=sqrt(ptot2);
      ires_r=siginfo[sigptr].ires0;
      sigpar=siginfo[sigptr].sigpar;
      mass_r=resinfo[ires_r].mass;
      width=resinfo[ires_r].width;
      delb2=sigpar/(.25*width*width+pow((mcheck-mass_r),2));
      if(resinfo[ires_r].minmass>mcheck) {
	delb2=0.0;
	printf("HOLY COW!!!  not enough oomph\n");
	printf("i=%d ires_i=%d     j=%d ires_j=%d\n",
	       i,partinfo[i].ires,j,partinfo[j].ires);
	printf("m2i=%g p_i=%g %g %g %g\n",m2i,p_i[0],p_i[1],p_i[2],p_i[3]);
	printf("m2j=%g p_j=%g %g %g %g\n",m2j,p_j[0],p_j[1],p_j[2],p_j[3]);
	printf("ires1=%d ires2=%d  ires_r=%d\n",
	       partinfo[i].ires,partinfo[j].ires,ires_r);
	printf("mi= %g mj=%g  mcheck=%g minmass=%g\n",
	       sqrt(m2i),sqrt(m2j),mcheck,resinfo[ires_r].minmass);
      }
      b2tot=b2tot+delb2;
      if(bsquared<b2tot){
	*scatcode=ires_r;
	//printf("Found resonant collision in colfind scatcode= %d\n",*scatcode);
	goto FOUND_RES_COLLISION;
      }
      sigptr=siginfo[sigptr].nextsigptr;
    }
  }
  b2tot=b2tot+b2elastic[partinfo[i].ires][partinfo[j].ires];
  *scatcode=0;
  if(bsquared>b2tot) goto GIVE_IT_UP;
FOUND_RES_COLLISION:
  zcoll=0.5*(r_i[3]+r_j[3]+p_i[3]*t_i+p_j[3]*t_j);
  t_i=t_i*p_i[0];
  t_j=t_j*p_j[0];
  *tevent=0.5*(r_i[0]+r_j[0]+t_i+t_j);
  *tau_event=sqrt(fabs(*tevent**tevent-zcoll*zcoll));
  *delt=r_j[0]-r_i[0]+t_j-t_i;
  test=1;
GIVE_IT_UP:
  return test;
}

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

int decayfind(int i,int *scatcode,double *tevent,
	      double *tau_event,double *delt){
  double gamma,tau,zcoll;
  int test,ires;
  test=0;
  ires=partinfo[i].ires;
  if(resinfo[ires].decay==1){
    gamma=partinfo[i].p[0]/partinfo[i].mass;
    tau=-(hbarc/resinfo[ires].width)*log(ran2());
    *tevent=partinfo[i].r[0]+tau*gamma;
    zcoll=partinfo[i].r[3]+((*tevent-partinfo[i].r[0])
      *partinfo[i].p[3]/partinfo[i].p[0]);
    *tau_event=sqrt(fabs(*tevent**tevent-zcoll*zcoll));
    *delt=0.0;
    *scatcode=-ires;
    test=1;
  }
  return test;
}

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