double triangle(double m0,double m1,double m2){
  double answer,m0sq,m1sq,m2sq;
#ifdef EXTRA_TESTS
  if(m0<m1+m2) {
    printf("Disaster with triangle\n");
    exit(1);
  }
#endif
  m0sq=m0*m0;m1sq=m1*m1;m2sq=m2*m2;
  answer=m0sq*m0sq+m1sq*m1sq+m2sq*m2sq;
  answer=answer-2.0*(m0sq*m1sq+m0sq*m2sq+m1sq*m2sq);
  answer=answer/(4.0*m0sq);
  return answer;
}

/* These functions perform lorentz transformations of p[4] to pprime[4], */
/* assuming the function is called, e.g. */
/* u[4] is the four velocity of the transformation */
/* Note: if one wanted to transform a vector measured in the frame noted */
/* by n[4] to a frame noted by uu[4]. */

void lorentz1(double *u,double *p1,double *p1prime){
  static int n[4]={1,0,0,0};
  int mu;
  double udotn=0.0,p1dotn=0.0,p1dotu=0.0;
  //double mtest=0.0,mtestprime=0.0;
  for(mu=0;mu<4;mu++){
    //mtest=mtest+p1[mu]*p1[mu]*g[mu];
    p1dotn=p1dotn+p1[mu]*n[mu]*g[mu];
    p1dotu=p1dotu+p1[mu]*u[mu]*g[mu];
    udotn=udotn+u[mu]*n[mu]*g[mu];
  }
  for(mu=0;mu<4;mu++){
    p1prime[mu]=-((p1dotu+p1dotn)/(1.0+udotn))*(n[mu]+u[mu])
      +2*p1dotn*u[mu]+p1[mu];
    //mtestprime=mtestprime+p1prime[mu]*p1prime[mu]*g[mu];
  }
  //if(fabs(mtestprime-mtest)>0.0000123) {
    //printf("FAILURE in lorentz1  mtest=%g mtestprime=%g\n",mtest,mtestprime);
    //printf("mass=%g\n",sqrt(mtest));
    //printf("u=%g %g %g %g\n",u[0],u[1],u[2],u[3]);
    //exit(1);
  //}
}

void lorentz2(double *u,double *p1,double *p1prime,double *p2,double *p2prime){
  static int n[4]={1,0,0,0};
  int mu;
  double udotn=0.0,p1dotn=0.0,p1dotu=0.0,p2dotn=0.0,p2dotu=0.0;
  //double mtest1=0.0,mtest1prime=0.0;
  //double mtest2=0.0,mtest2prime=0.0;
  for(mu=0;mu<4;mu++){
    //mtest1=mtest1+p1[mu]*p1[mu]*g[mu];
    //mtest2=mtest2+p2[mu]*p2[mu]*g[mu];
    p1dotn=p1dotn+p1[mu]*n[mu]*g[mu];
    p1dotu=p1dotu+p1[mu]*u[mu]*g[mu];
    p2dotn=p2dotn+p2[mu]*n[mu]*g[mu];
    p2dotu=p2dotu+p2[mu]*u[mu]*g[mu];
    udotn=udotn+u[mu]*n[mu]*g[mu];
  }
  for(mu=0;mu<4;mu++){
    p1prime[mu]=-((p1dotu+p1dotn)/(1.0+udotn))*(n[mu]+u[mu])
      +2*p1dotn*u[mu]+p1[mu];
    p2prime[mu]=-((p2dotu+p2dotn)/(1.0+udotn))*(n[mu]+u[mu])
      +2*p2dotn*u[mu]+p2[mu];
    //mtest1prime=mtest1prime+p1prime[mu]*p1prime[mu]*g[mu];
    //mtest2prime=mtest2prime+p2prime[mu]*p2prime[mu]*g[mu];
  }
  //if(fabs(mtest1prime-mtest1)>0.0000123
       //||fabs(mtest2prime-mtest2)>0.0000123) {
    //printf("FAILURE in lorentz2\n");
    //exit(1);
  //}
}

void lorentz3(double *u,double *p1,double *p1prime,double *p2,double *p2prime,double *p3,double *p3prime){
  static int n[4]={1,0,0,0};
  int mu;
  double udotn=0.0,p1dotn=0.0,p1dotu=0.0,p2dotn=0.0,p2dotu=0.0,
    p3dotn=0.0,p3dotu=0.0;
  //double mtest1=0.0,mtest1prime=0.0;
  //double mtest2=0.0,mtest2prime=0.0;
  //double mtest3=0.0,mtest3prime=0.0;
  for(mu=0;mu<4;mu++){
    //mtest1=mtest1+p1[mu]*p1[mu]*g[mu];
    //mtest2=mtest2+p2[mu]*p2[mu]*g[mu];
    //mtest3=mtest3+p3[mu]*p3[mu]*g[mu];
    p1dotn=p1dotn+p1[mu]*n[mu]*g[mu];
    p1dotu=p1dotu+p1[mu]*u[mu]*g[mu];
    p2dotn=p2dotn+p2[mu]*n[mu]*g[mu];
    p2dotu=p2dotu+p2[mu]*u[mu]*g[mu];
    p3dotn=p3dotn+p3[mu]*n[mu]*g[mu];
    p3dotu=p3dotu+p3[mu]*u[mu]*g[mu];
    udotn=udotn+u[mu]*n[mu]*g[mu];
  }
  for(mu=0;mu<4;mu++){
    p1prime[mu]=-((p1dotu+p1dotn)/(1.0+udotn))*(n[mu]+u[mu])
      +2*p1dotn*u[mu]+p1[mu];
    p2prime[mu]=-((p2dotu+p2dotn)/(1.0+udotn))*(n[mu]+u[mu])
      +2*p2dotn*u[mu]+p2[mu];
    p3prime[mu]=-((p3dotu+p3dotn)/(1.0+udotn))*(n[mu]+u[mu])
      +2*p3dotn*u[mu]+p3[mu];
    //mtest1prime=mtest1prime+p1prime[mu]*p1prime[mu]*g[mu];
    //mtest2prime=mtest2prime+p2prime[mu]*p2prime[mu]*g[mu];
    //mtest3prime=mtest3prime+p3prime[mu]*p3prime[mu]*g[mu];
  }
  //if(fabs(mtest1prime-mtest1)>0.0000123){ 
    //printf("failure of test1 in lorentz3\n");
    //exit(1);
  //}
  //if(fabs(mtest2prime-mtest2)>0.0000123) {
    //printf("failure of test2 in lorentz3\n");
    //exit(1);
  //}
  //if(fabs(mtest3prime-mtest3)>0.0000123) {
    //printf("failure of test3 in lorentz3\n");
    //exit(1);
  //}
  //printf("exiting lorentz3 successfully\n");
}

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

/* This the Shell method which works very well for n <= 50 */

void shell_sort(int n, double *a, int *b){
  long i,j,inc;
  double va;
  int vb;
  inc=1;
  do {
    inc *=3;
    inc++;
  } while (inc <= n);
  do {
    inc /= 3;
    for (i=inc+1;i<=n;i++) {
      va=a[i];
      vb=b[i];
      j=i;
      while (a[j-inc]>va) {
	a[j]=a[j-inc];
	b[j]=b[j-inc];
	j-=inc;
	if(j<=inc) break;
      }
      a[j]=va;
      b[j]=vb;
    }
  } while (inc>1);
}

void addevent(int *firstn,double t,double tau,double delt,
	      int scatcode,int i,int j){
  int oldfirst,alpha;
  oldfirst=*firstn;
  *firstn=firstemptyevent;
  if(firstemptyevent==neventmax){
    printf("neventmax is too small\n");
    exit(1);
  }
  firstemptyevent=eventinfo[firstemptyevent].next;
  eventinfo[*firstn].time=t;
  eventinfo[*firstn].tau=tau;
  eventinfo[*firstn].delt=delt;
  eventinfo[*firstn].scatcode=scatcode;
  eventinfo[*firstn].next=oldfirst;
  eventinfo[*firstn].i=i;
  if(scatcode<0&&i!=j) {
    printf("DANGER event needs i=j\n");
    exit(1);
  }
  eventinfo[*firstn].j=j;
}



#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double ran2() {
  long J;
  long K;
  static long IDUM2=123456789;
  static long IY=0;
  static long IV[NTAB];
  double randy;
  if(IDUM <=0){
    if(-(IDUM)<1) IDUM=1;
    else IDUM = -(IDUM);
    IDUM2=(IDUM);
    for (J=NTAB+7;J>=0;J--) {
      K=(IDUM)/IQ1;
      IDUM=IA1*(IDUM-K*IQ1)-K*IR1;
      if(IDUM<0) IDUM +=IM1;
      if(J<NTAB) IV[J]=IDUM;
    }
    IY=IV[0];
    //printf("Initialization Completed %d\n",IDUM);
  }
  K=(IDUM)/IQ1;
  IDUM=IA1*(IDUM-K*IQ1)-K*IR1;
  if(IDUM<0) IDUM+=IM1;
  K=IDUM2/IQ2;
  IDUM2=IA2*(IDUM2-K*IQ2)-K*IR2;
  if(IDUM2<0) IDUM2+=IM2;
  J=IY/NDIV;
  IY=IV[J]-IDUM2;
  IV[J]=IDUM;
  if(IY<1) IY+=IMM1;
  randy=AM*IY;
  if(randy>RNMX) return RNMX;
  else return randy;
}

/* These are only compiled if CHAIN_TESTS is defined */
#ifdef CHAIN_TESTS
/* Test the chain */
void event_chain_test(int nevents){
  int ievent,ii;
  ievent=firstevent;
  for(ii=1;ii<=nevents;ii++){
    ievent=eventinfo[ievent].next;
  }
  if(ievent!=-1) exit(1);
}

void part_chain_test(){
  int partcount,prevpart,ipart,kk,nextpart;
  ipart=firstpart;
  prevpart=-1;
  partcount=0;
  while(ipart!=-1){
    nextpart=partinfo[ipart].next;
    partcount=partcount+1;
    if(partcount>npart) {
      ipart=firstpart;
      prevpart=-1;
      partcount=0;
      for(kk=0;kk<=npart;kk++){
	nextpart=partinfo[ipart].next;
	partcount=partcount+1;
	prevpart=ipart;
	ipart=nextpart;
      }
      printf("part_chain_test failed!");
      exit(1);
    }
    prevpart=ipart;
    ipart=nextpart;
  }
  if(npart!=partcount) {
    printf("FAILURE: npart=%d partcount=%d\n",
	   npart,partcount);
    exit(1);
  }
}

#endif