void make_particle(double *p){
  double pmag,cthet,sthet,phi;
  double f,f0,fmax,sign,mu,e,rootmt,test,g1,g2;
  int isuccess,ident;
  sign=1.0;
  if(abs(baryon[ires])>0) sign=-1.0;
  mu=mus*strangeness[ires]+mub*baryon[ires];
  fmax=fugacity*exp(-(mass[ires]-mu)/temp);
  if(mass[ires]/temp<10.0){
  RETRY1:
    pmag=-temp*log(ran2()*ran2()*ran2());
    e=sqrt(mass[ires]*mass[ires]+pmag*pmag);
    f=fugacity*exp(-(e-mu)/temp);
    f0=fugacity*exp(-(pmag-mu)/temp);
    f=f/(1.0-sign*quantum*f);
    test=f/f0;
    if(sign>0) test=test*(1.0-quantum*fmax);
    if(test>1.0) printf("f too big in make_particle\n");
    if(ran2()>test) goto RETRY1;
    phi=2.0*pi*ran2();
    cthet=1.0-2.0*ran2();
    sthet=sqrt(1.0-cthet*cthet);
    p[0]=e;
    p[3]=pmag*cthet;
    p[1]=pmag*sthet*cos(phi);
    p[2]=pmag*sthet*sin(phi);
  }
  else{
    rootmt=sqrt(mass[ires]*temp);
  RETRY2:
    gauss(&g1,&g2);
    p[1]=rootmt*g1;p[2]=rootmt*g2;
    gauss(&g1,&g2);
    p[3]=rootmt*g1;
    e=mass[ires]+.5*(p[1]*p[1]+p[2]*p[2]+p[3]*p[3])/mass[ires];
    p[0]=e;
    f0=fugacity*exp(-(e-mu)/temp);
    test=1.0/(1.0-quantum*sign*f0);
    if(sign>0) test=test*(1.0-quantum*fmax);
    if(ran2()>test) goto RETRY2;
  }
}

double resonanceinit(){
  int itype,degen[nresonances],ip,nj;
  int zz,bb,ss;
  double spin,mm,resweight[nresonances],dndv;
  double f,pmag,e,mu,sign,delp,constant;
  char cdummy[80];
  FILE *initinput,*initoutput;
  initinput=fopen("resonance_info/resonance_list.dat","r");
  initoutput=fopen("temp_data/initoutput.dat","w");
  for(ires=0;ires<nresonances;ires++){
    fscanf(initinput,"%d %lf %d %d %d %lf ",&itype,&mm,
	   &zz,&bb,&ss,&spin);
    fgets(cdummy,80,initinput);
    degen[ires]=(int)(2.0*spin+1.0);
    fprintf(initoutput,"%d %lf %d\n",itype,mm,degen[ires]);
    charge[ires]=zz;
    strangeness[ires]=ss;
    baryon[ires]=bb;
    itypefind[ires]=itype;
    mass[ires]=mm;
    resweight[ires]=0.0;
    /* don't generate photons */
    if(itype==22) degen[ires]=0;
    /* ---------------------------------------- */
  }
  fclose(initinput);
  delp=temp/50.0;
  constant=delp/(2.0*pi*pi*pow(197.323,3));
  for(ip=0;ip<1000;ip++){
    pmag=(ip+.5)*delp;
    for(ires=0;ires<nresonances;ires++){
      sign=1.0;
      if(abs((baryon[ires]))>0) sign=-1.0;
      mu=mus*strangeness[ires]+mub*baryon[ires];
      e=sqrt(mass[ires]*mass[ires]+pmag*pmag);
      f=fugacity*exp(-(e-mu)/temp);
      f=f/(1.0-sign*quantum*f);
      resweight[ires]=resweight[ires]+pmag*pmag*degen[ires]*f*constant;
    }
  }
  dndv=0.0;
  for(ires=0;ires<nresonances;ires++){
    dndv=dndv+resweight[ires];
  }
  jmax=0;
  for(ires=0;ires<nresonances;ires++){
    nj=(int)floor(0.5+(jlimit-100)*resweight[ires]/dndv);
    //if(nj==0) printf("This particle never made:  m=%7.2lf\n",mass[ires]);
    for(j=jmax;j<jmax+nj;j++){
      iresfind[j]=ires;
      jweight[j]=jlimit*resweight[ires]/(nj*dndv);
      /*  /printf("%lf   %d   %lf\n",mass[ires],j,jweight[j]);*/
    }
    if(nj!=0) fprintf(initoutput,"%lf %lf %d\n",mass[ires],jweight[jmax],nj);
    jmax=jmax+nj;
  }
  if(jmax>jlimit) printf("DANGER  jmax too big   %d\n");
  fclose(initoutput);
  printf("Initialization complete   jmax=%d\n",jmax);
  return dndv;
}

#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 TEMP;
  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;
  TEMP=AM*IY;
  if(TEMP>RNMX) return RNMX;
  else return TEMP;
}
void gauss(double *g1,double *g2){
  double x,y,r2,r;
  do{
    x=1.0-2.0*ran2();
    y=1.0-2.0*ran2();
    r2=x*x+y*y;
  }while (r2>1.0);
  r=sqrt(r2);
  *g1=(x/r)*sqrt(-2.0*log(r2));
  *g2=(y/x)**g1;
}