c--------------------------------------------------------------------
c Generates Random Galaxy positions. 
c
c
c  Input:
c     iseed   = random number seed (should be a negative integer)
c     zlimit  = redshift at which to truncate the catalogue eg z=0.5
c     selfile = name of file containing the tabulated selection function
c               eg om1:lam0:strc:a/sel1.2df_sgp
c                  (or 00002df_sgp for a volume limited catalogue)
c     vlim    = set .TRUE. if a volume limited cat is required
c
c     vran    = 1-dimensional random velocity dispersion in units
c               of km/s (set in parameter statement)
c
c  Output:
c     x,y,z  = Galaxy redshift space position
c     d      = Galaxy redshift if on pure Hubble Flow. NB Will equal
c               sqrt(x^2+y^2+z^2) if vran = 0.0 km/s .
c
c  Method:
c

      subroutine ran_cat(x,y,z,d,iseed,zlimit,selfile,vlim)
      implicit none
      integer iseed,ifirst,lnblnk
      character selfile*(*), cat*4
      real x,y,z,d,vx,vy,vz,zlimit,z_gal,z_hub,v_dot_rhat,vran,d2
      parameter (vran=300.0)
      logical in,vlim
      real ran3,gasdev
      save ifirst,cat
      data ifirst/1/
c
c On the first call read the specified selection function file
c and identify whether it is "sdss" or 2dF "_sgp" or "_ngp" so
c that the appropriate geometrical boudaries can be applied.
      if (ifirst.eq.1) then
         cat=selfile(lnblnk(selfile)-3:lnblnk(selfile))
         write(*,*) 'cat=',cat,' selfile=',selfile
         if ((cat.ne.'sdss').and.(cat.ne.'_sgp').and.
     &        (cat.ne.'_ngp').and.(cat.ne.'pscz') ) 
     &        stop 'Unrecognised selection function' 
         call read_sel(selfile,zlimit,vlim)
         ifirst=0
      end if
c
c
c We now loop until we find a galaxy in the catalogue
      in=.FALSE.
      do while (.not.in)
c        Generate a random real-space position in a cubic box
         x = (ran3(iseed)-0.5)*2.0
         y = (ran3(iseed)-0.5)*2.0
         z = (ran3(iseed)-0.5)*2.0
         d2= x**2 + y**2 + z**2
c        Apply constraint to retain only an isotropic distribution
         if (d2.lt.1.0) then
c           Apply Geometric constraint of specific survey
            d=sqrt(d2)
            if (cat.eq.'sdss') then
               call geometry_sdss(in,x,y,z,d)
            else if (cat.eq.'_sgp') then
               call geometry_2df_sgp(in,x,y,z,d)
            else if (cat.eq.'_ngp') then
               call geometry_2df_ngp(in,x,y,z,d)
            else if (cat.eq.'pscz') then
               call geometry_pscz(in,x,y,z,d)
            end if
            if (in) then       
c              Select a random velocity. 
               if (vran.gt.0.0) then
c                 Each component chosen from a Gaussian of 
c                 sigma=vran km/s 
                  vx=gasdev(iseed)*vran
                  vy=gasdev(iseed)*vran
                  vz=gasdev(iseed)*vran
               else
                  vx=0.0
                  vy=0.0
                  vz=0.0
               end if
               
c              Generate redshift randomly from the tabulated N(<z).
               call get_redshift(z_hub,iseed)
               
c              Compute velocity Component along the line of sight
               v_dot_rhat=(vx*x + vy*y + vz*z)/d
c              and compute redshift via  (1+z_obs)=(1+z_hub)(1+v/c)
               z_gal   = ((1.0+z_hub)*(1+v_dot_rhat/3.0E+05)) - 1.0
               x = x*z_gal/d     !scale to redshift-space positions
               y = y*z_gal/d
               z = z*z_gal/d
               d = z_hub           !save redshift if on Hubble Flow
            end if
         end if
      end do                     !try again
      return
      end
c --------------------------------------------------
c  Generate a redshift consistent with the tabulated distribution
c
c  Method:
c    1) Generate a random number, nran, in the interval [0,N(z<zlimit)].
c    2) Search the table to find the redshift at which nran=N(<z) and
c    return z
c
      subroutine get_redshift(z_hub,iseed)
      implicit none
      integer nrmax,ir,iseed,nr
      parameter (nrmax=2000)
      real nltz(nrmax),z(nrmax),nran,ran3,h,z_hub
      common/selfun/ nltz,z,nr
c
c     Generate random value in the interval [0,N(z<zlimit)]
      nran=ran3(iseed)*nltz(nr)
c     Locate the closest redshift
      call locate(nltz,nr,nran,ir)  !numerical recipes routine
      ir=max(ir,1)
c     Interpolate
      h=(nran-nltz(ir))/(nltz(ir+1)-nltz(ir))
      z_hub=z(ir+1)*h+z(ir)*(1.0-h)
      return
      end
c --------------------------------------------------
c  Read the Tabulated Selection Function File 
c  into common /selfun/
c
      subroutine read_sel(selfile,zlimit,vlim)
      implicit none
      character*80 selfile
      integer nrmax,nr
      parameter (nrmax=2000)
      real sel,r,z(nrmax),nltz(nrmax),zlimit
      logical vlim
      common/selfun/ nltz,z,nr
      open(10,file=selfile,status='old',form='formatted')
      nr=0
      read(10,*)
      read(10,*)
c Magnitude limited: N(<r)
      if (.not.vlim) then
         do while ((nr.lt.nrmax).and.(z(max(nr,1)).le.zlimit))
            nr=nr+1
            read(10,*,err=12) z(nr),r,nltz(nr),sel 
         end do
c volume limited: Vol(<r)
      else
         do while ((nr.lt.nrmax).and.(z(max(nr,1)).le.zlimit))
            nr=nr+1
            read(10,*,err=12),z(nr),r
            nltz(nr)=r**3
         end do
      end if
      close(10)
      return

 12   continue
      stop 'z_limit is higher than tabulated z!'
      end
c------------------------------------------------------------------------
      SUBROUTINE locate(xx,n,x,j)
      INTEGER j,n
      REAL x,xx(n)
      INTEGER jl,jm,ju
      jl=0
      ju=n+1
10    if(ju-jl.gt.1)then
        jm=(ju+jl)/2
        if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then
          jl=jm
        else
          ju=jm
        endif
      goto 10
      endif
      j=jl
      return
      END
c
      FUNCTION ran3(idum)
      INTEGER idum
      INTEGER MBIG,MSEED,MZ
C     REAL MBIG,MSEED,MZ
      REAL ran3,FAC
      PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
C     PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
      INTEGER i,iff,ii,inext,inextp,k
      INTEGER mj,mk,ma(55)
C     REAL mj,mk,ma(55)
      SAVE iff,inext,inextp,ma
      DATA iff /0/
      if(idum.lt.0.or.iff.eq.0)then
        iff=1
        mj=MSEED-iabs(idum)
        mj=mod(mj,MBIG)
        ma(55)=mj
        mk=1
        do 11 i=1,54
          ii=mod(21*i,55)
          ma(ii)=mk
          mk=mj-mk
          if(mk.lt.MZ)mk=mk+MBIG
          mj=ma(ii)
11      continue
        do 13 k=1,4
          do 12 i=1,55
            ma(i)=ma(i)-ma(1+mod(i+30,55))
            if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
12        continue
13      continue
        inext=0
        inextp=31
        idum=1
      endif
      inext=inext+1
      if(inext.eq.56)inext=1
      inextp=inextp+1
      if(inextp.eq.56)inextp=1
      mj=ma(inext)-ma(inextp)
      if(mj.lt.MZ)mj=mj+MBIG
      ma(inext)=mj
      ran3=mj*FAC
      return
      END
c
      FUNCTION gasdev(idum)
      INTEGER idum
      REAL gasdev
CU    USES ran3
      INTEGER iset
      REAL fac,gset,rsq,v1,v2,ran3
      SAVE iset,gset
      DATA iset/0/
      if (iset.eq.0) then
1       v1=2.*ran3(idum)-1.
        v2=2.*ran3(idum)-1.
        rsq=v1**2+v2**2
        if(rsq.ge.1..or.rsq.eq.0.)goto 1
        fac=sqrt(-2.*log(rsq)/rsq)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      else
        gasdev=gset
        iset=0
      endif
      return
      END