c     Version 2.0, November 1, 1994.
      program ctester
      implicit none
      double precision mass1,mass2,gammav,gamma,temp
      double precision pt1,pt2,rx,ry,rz,tau,lambda,b,bmin,bmax
      double precision bw0,bcut,x,y,z,t,e,px,py,pz,znew,z0,gauss
      integer*4 IDUM,IDUM2,IY,IV(32)
      common/RAN2JUNK/IDUM,IDUM2,IY,IV
      double precision randy,ran2
      double precision rperp
      double precision dely,rapid,rapid0,rapid1,rapid2
      integer ib,nb,i,ident,npts1,npts2,iseed
c     This program works very similarly to 'gaussmaker, only points
c     are generated according to a Bjorken geometry with no sidewards
c     expansion.  The user needs to enter the range of rapidities to
c     generate points for.  The program correctly finds the x,y,z,t assuming
c     the Bjorken source is infinite along the z-axis, only it only writes
c     momenta which appear in the specified rapidity range.
c     
      mass2=938.28d0
c     mass2=939.6d0
c     mass2=139.58d0
c     mass2=493.646d0
      mass1=mass2
c     mass1=1875.585d0
c     mass2=3727.323d0
c     mass1=938.28d0
c     Initialize random number generator
      IDUM=-1234
      nb=1
c     Make ident 0 if particles are not identical, 1 if they are identical
      ident=1
      open(20,form='formatted',file='phasedescription.dat',
     +     status='unknown')
      open(21,form='formatted',file='phasespace1.dat',
     +     status='unknown')
      write(20,*) nb
      write(6,*) 'what are lower and upper limits of the rapidities?'
      read(5,*) rapid1,rapid2
      write(6,*) 'what is the temperature?'
      read(5,*) temp
      pt1=dsqrt(mass1*temp)
      pt2=dsqrt(mass2*temp)
      write(6,*) 'what is rperp?'
      read(5,*) rperp
      write(6,*) 'what is bjorken tau?'
      read(5,*) tau
      write(6,*) 'what is lambda?'
      read(5,*) lambda
      lambda=dsqrt(lambda)
      do 111 ib=1,1
         b=0.d0
         bmin=0.d0
         bmax=.5d0
         bw0=1.d0
         bcut=1.d0
         npts1=10000
         write(20,*) b,bmin,bmax,npts1,npts2,bw0,bcut
         do 727 i=1,npts1
            call relpfind(px,py,pz,mass1,temp)
            e=dsqrt(mass1**2+px**2+py**2+pz**2)
            randy=ran2()
            if(randy.lt.lambda) then
               x=rperp*gauss()
               y=rperp*gauss()
               z0=0.d0
            else
               x=1000000.d0*gauss()
               y=1000000.d0*gauss()
               z0=1000000.0*gauss()
            endif
            t=dsqrt(tau**2+z0**2)
            rapid=rapid1+(rapid2-rapid1)*ran2()
            rapid0=.5d0*dlog((1.d0+pz/e)/(1.d0-pz/e))
            dely=rapid-rapid0
            gamma=cosh(dely)
            gammav=sinh(dely)
            pz=gamma*pz+gammav*e
            e=dsqrt(mass1**2+pz**2+px**2+py**2)
            z=gamma*z0+gammav*t
            t=gamma*t+gammav*z0
            write(21,*) px,py,pz,x,y,z,t
 727     continue
 111  continue
      close(20)
      close(21)
      stop
      end
c     __________________________________________
      double precision function gauss()
      implicit double precision(a-h,o-z)
 100  x=(1.-2.*ran2())
      y=(1.-2.*ran2())
      r2=x**2+y**2
      if(r2.gt.1.0) go to 100
      r=dsqrt(r2)
      c=x/r
      gauss=c*dsqrt(-2.0*dlog(r2))
      return
      end

c     
c     ---------Random Number generating Routine------------
c     This needs to go in the main program
c     integer*4 IDUM,IDUM2,IY,IV(32)
c     common/RAN2JUNK/IDUM,IDUM2,IY,IV
c     
      double precision function ran2()
      implicit none
      integer IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      real EPS,RNMX,AM
      parameter (IM1=2147483563,IM2=2147483399,IMM1=(IM1-1))
      parameter (IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211)
      parameter (IR2=3791,NTAB=32,NDIV=(1+IMM1/NTAB))
      parameter (EPS=1.2e-7,RNMX=(1.0-EPS),AM=(1.0/IM1))
c
      integer*4 IDUM
      integer*4 J
      integer*4 K
      integer*4 IDUM2
      integer*4 IY
      integer*4 IV(0:NTAB-1)
      double precision TEMP
      common/RAN2JUNK/IDUM,IDUM2,IY,IV
      if(IDUM.le.0) then
         IDUM2=123456789
         IY=0
         IDUM =-IDUM
         IDUM2=IDUM
         do J=NTAB+7,0,-1
            K=(IDUM)/IQ1
            IDUM=IA1*(IDUM-K*IQ1)-K*IR1
            if(IDUM.lt.0) IDUM =IDUM+IM1
            if(J.lt.NTAB) IV(J)=IDUM
         enddo
         IY=IV(0)
c         write(6,*) 'Initialization completed, IDUM=', IDUM
      endif
      K=(IDUM)/IQ1
      IDUM=IA1*(IDUM-K*IQ1)-K*IR1
      if(IDUM.lt.0) IDUM=IDUM+IM1
      K=IDUM2/IQ2
      IDUM2=IA2*(IDUM2-K*IQ2)-K*IR2
      if(IDUM2.lt.0) IDUM2=IDUM2+IM2
      J=IY/NDIV
      IY=IV(J)-IDUM2
      IV(J)=IDUM
      if(IY.lt.1) IY=IY+IMM1
      TEMP=AM*IY
      if(TEMP.gt.RNMX) then
         ran2=RNMX
      else
         ran2=TEMP
      endif
      return
      end
c
      subroutine relpfind(px,py,pz,mass,temp)
      implicit none
      double precision px,py,pz,mass,temp,e,randy,ran2
      double precision p,sth,cth,phi,w,arg
      integer i,idummy
 111  arg=ran2()*ran2()*ran2()
      if(arg.lt.0.0000001d0) go to 111
      p=-temp*dlog(arg)
      e=dsqrt(mass**2+p**2)
      w=dexp(-(e-p)/temp)
      randy=ran2()
      if(randy.gt.w) go to 111
      cth=1.d0-2.d0*ran2()
      sth=dsqrt(1.d0-cth**2)
      phi=2.d0*3.14159265358979323844d0*ran2()
      pz=p*cth
      px=p*sth*dcos(phi)
      py=p*sth*dsin(phi)
      return
      end