c  Subroutine to read selected data from mock 2dF and SDSS catalogues. 
c  Compile such that the file is read assuming bigendian byte ordering
c  eg -convert big_endian   on a DEC-alpha
c     default on SUNS or SGI.
c
c  Options:
c       ic=0   -- read only np,zmax,sigma8,omega0 and lambda0
c       ic=3   -- also read redshift space x,y,z
c       ic=4   -- also read redshift space x,y,z and true distance zh
c       ic=5   -- also read redshift space x,y,z,d and apparent magnitude mag
c       ic=6   -- also read redshift space x,y,z,d,mag and zhmax
c       ic=7   -- also read redshift space x,y,z,d,mag,zhmax  and id
c  Note any array not read can be declared with unit dimension 
c  in the calling program. eg REAL d(1)
c
c
c  File Format:
c
c        0,ng
c        zmax
c        2
c        xc,yc,zc,basis(3,3),mockseed [,vxc,vyc,vzc]
c        sigma_8,omega0,lambda0
c        (x(ig),ig=1,ng)
c        (y(ig),ig=1,ng)
c        (z(ig),ig=1,ng)
c        (zh(ig),ig=1,ng)
c        (mag(ig),ig=1,ng)
c        (zhmax(ig),ig=1,ng)
c        (id(ig),ig=1,ng)
c
c  
c  
c Variables:
c
c Input
c     file  -- name of mock catalogue file eg 1.sdss
c     
c Output
c     ng     -- number of galaxies in the mock catalogue
c     zmax   -- highest galaxy redshift 
c     sigma8 -- rms mass fluctuation in spheres of radius=8Mpc/h
c     omega0 -- density parameter
c     lambda -- cosmological constant in units of 3H_0^2
c     x(),y(),z() -- cartesian redshift-space position of each galaxy
c     zh() --       Hubble flow redshift, the redshift each galaxy would have 
c                   if it had no peculiar velocity
c                    
c     mag() --   Apparent Bj magnitude
c     zhmax() -- The maximum Hubble flow redshift at which the galaxy
c                would be selected in the catalogue
c     id() --    Particle index in the orginal N-body simulation
c
c
c Coordinate Systems:
c     For the 2dF catalogues the cartesian x,y,z are related to redshift, zz
c     and RA and dec by
c     zz = sqrt(x^2 +y^2 +z^2)
c     sin(dec) = z/zz              [  dec= asin(z/zz) ]
c     tan(RA)  = y/x               [  ra = atan2(y/x) ]
c
c     For the SDSS catalogues the z-axis of coordinate system points to
c     the centre of the SDSS survey, which elliptical in shape. 
c     Redshift,zz is again zz = sqrt(x^2 +y^2 +z^2) but now
c     sin(lat)  = z/zz              [  lat = asin(z/zz) ]
c     tan(long) = y/x               [ long = atan2(y/x) ]c
c     where
c     lat = 90 deg is the centre of the SDSS survey and
c     long is measured relative to the long axis of the SDSS ellipse.
c
c     For the PSCZ catalogues the z-axis of coordinate system points to
c     the NGP.
c     Redshift,zz is again zz = sqrt(x^2 +y^2 +z^2) but now
c     sin(b)  = z/zz                [  b = asin(z/zz) ]
c     tan(l) = y/x                  [  l = atan2(y/x) ]
c
c 
c
c
c Conversion to real-space coordinates:
c     The effect of peculiar velocities can be removed by
c     rescaling each x,y,z by the factor zh/zz .
c
c
c
c

c-------------------------------------------------------------------------
      subroutine mockread(file,ng,zmax,sigma8,omega0,lambda0,
     &     x,y,z,zh,mag,zhmax,id,ic)
      implicit none
      integer ng,n,i,ic,idf,mockseed,id(*)
      real x(*),y(*),z(*),zh(*),mag(*),zhmax(*),
     & zmax,sigma8,omega0,lambda0
      character file*(*)
      real xc,yc,zc,basis(3,3),vxc,vyc,vzc
      common /mock_header/ xc,yc,zc,basis,mockseed,vxc,vyc,vzc 
c
      open(11,file=file,status='old',form='unformatted')
      read(11) n,ng !n=0  and ng=number of galaxies in the catalogue
      write(0,*) 'number of galaxies in the catalogue=',ng
      read(11) zmax     
      read(11) idf  ! idf=2
      if (idf.eq.2) then
      else if (idf.eq.5) then
         write(0,*) 'reading data from genuine 2dF catalogue'
      else
         stop 'idf/ismc value not recognized'
      end if
      read(11) xc,yc,zc,basis,mockseed,vxc,vyc,vzc !defines position, 
c                                   !orientation and velocity of the 
c                                   !observer in the N-body simulation
      read(11) sigma8,omega0,lambda0!mass fluctuation in 8Mpc/h, density
d      write(0,*) 'sigma8=',sigma8
d      write(0,*) 'omega0=',omega0,' lambda0=',,lambda0 !mass fluctuation in 8Mpc/h, density
c                                   !parameter and cosmological constant
c      write(0,*) 'n= ',n,' ng= ',ng
c      write(0,*) 'zmax= ',zmax
c      write(0,*) 'idf= ',idf
c      write(0,*) 'xc=',xc,' yc=',yc,' zc=',zc
c      write(0,*) 'vxc=',vxc,' vyc=',vyc,' vzc=',vzc
c      write(0,*) 'basis:',basis
c      write(0,*) 'mockseed=',mockseed      
c      write(0,*) 'sigma_8=',sigma8,' omega0=',omega0,' lambda0=',lambda0
      if (ic.eq.0) close(11) 
      if (ic.eq.0) return !stop after reading header information if ic=0
      read(11) (x(i),i=1,ng)
c      write(0,*) 'read x'
      read(11) (y(i),i=1,ng)
c      write(0,*) 'read y'
      read(11) (z(i),i=1,ng)
c      write(0,*) 'read z'
      call flush(0)
      write(0,*) 'Redshift-space positions read'
      if (idf.eq.2) then
         if (ic.ge.4) read(11) (zh(i),i=1,ng)  
         if (ic.ge.4) write(0,*) 'True distances read'
         if (ic.ge.5) read(11) (mag(i),i=1,ng)  
         if (ic.ge.5) write(0,*) 'Apparent magnitudes read'
         if (ic.ge.6) read(11) (zhmax(i),i=1,ng)  
         if (ic.ge.6) write(0,*) 'd_max values read'
         if (ic.eq.7) read(11) (id(i),i=1,ng)  
         if (ic.ge.7) write(0,*) 'Particle identifications read'
c        id() is the particle index in the original N-body simulation
      else if (idf.eq.5) then
         if (ic.ge.4) read(11) (zh(i),i=1,ng)  
         if (ic.ge.4) write(0,*) 'redshift distance read'
         if (ic.ge.5) read(11) (mag(i),i=1,ng)  
         if (ic.ge.5) write(0,*) 'Apparent magnitudes read'
         if (ic.ge.6) read(11) (zhmax(i),i=1,ng)  
         if (ic.ge.6) write(0,*) 'something else read!'
         if (ic.eq.7) read(11) (id(i),i=1,ng)  
         if (ic.ge.7) write(0,*) 'target flag read'
      end if
      close(11)
      return
      end
c ------------------------------------------------------------------
c------------------------------------------------------------------------------
c   Use the redshift space direction and observer (local group)
c   velocity to convert the redshift space position to that of 
c   the local group frame.
c
c   Note that detect_blueshift() should be called before using
c   this conversion subroutine (subroutine included below).
c
c   This subroutine also checks that the observer velocity is
c   is alligned with the direction of the CMB dipole. It can
c   cope automatically with the observer velocity being in the
c   simulation coordinate system or that of the mock catalogue.
c
c   Finally it reports the number of blue-shifted galaxies
c   that remain in the LG velocity frame.
c
c     
      subroutine  cmb_to_lg(ng,x,y,z,zh,vxc,vyc,vzc,basis)
*************************************variables*********************************
      implicit none
      integer ng,i,nblue
      real PI,EPS,RLIMIT
      parameter (PI=3.141592654,EPS=1.0,RLIMIT=2.0)
      real x(ng),y(ng),z(ng),zh(ng),vxc,vyc,vzc,zz,vdot,zp_zz,dir,
     & vxo,vyo,vzo,basis(3,3),vobs
      logical cmbdipole
*******************************************************************************
c     Rotate velocity from simulation coordinates to cartesian equatorial
      vxo=vxc*basis(1,1)+vyc*basis(2,1)+vzc*basis(3,1)
      vyo=vxc*basis(1,2)+vyc*basis(2,2)+vzc*basis(3,2)
      vzo=vxc*basis(1,3)+vyc*basis(2,3)+vzc*basis(3,3)
      vobs=sqrt(vxc**2 + vyc**2 +vzc**2)
c
c     Check direction compared with that of the CMB dipole and magnitude
      if (vobs.le.0) then
         write(0,*) 'Observer at rest. No LG transformation possible' 
         return
      else
         cmbdipole=.true.
         if (abs(asin(vzo/vobs)*180.0/PI+27.07).gt.EPS) cmbdipole=.false.
         if (abs(atan2(vyo,vxo)*180.0/PI-166.05).gt.EPS) cmbdipole=.false.
         if (.not.cmbdipole) then
            write(0,*) 'Trying unrotated velocities'
            vxo=vxc
            vyo=vyc
            vzo=vzc
            cmbdipole=.true.
            if (abs(asin(vzo/vobs)*180.0/PI+27.07).gt.EPS) cmbdipole=.false.
            if (abs(atan2(vyo,vxo)*180.0/PI-166.05).gt.EPS) cmbdipole=.false.
            if (.not.cmbdipole) stop 'observer velocity mis-aligned'
         end if
      end if
c
c     Map from CMB frame to LG frame
      do i=1,ng
         zz=sqrt(x(i)**2+y(i)**2+z(i)**2)
c        Component of observer's velocity in the direction
c        of the target galaxy in units of 3.0e+05 km/s
         vdot= ( vxo*x(i)+vyo*y(i)+vzo*z(i) ) / (3.0e+05*zz)
c        Detect encoded blue-shift
         dir=+1.0
         if (zh(i).lt.0.0) dir=-1.0   
c        (1+z^prime) = (1+z) * (1-vdot) 
c        but if blue-shift is encoded then z --> -z
c        ie (1+z^prime) = (1-z) * (1-vdot) 
c        Hence z^prime/abs(z) = zp_zz =
         zp_zz= ((1.0+dir*zz)*(1.0-vdot)-1.0)/zz
c        Rescale redshift space vector by the magnitude of this factor
         x(i)=abs(zp_zz)*x(i)
         y(i)=abs(zp_zz)*y(i)
         z(i)=abs(zp_zz)*z(i)
c        If and only if zp_zz is negative is this galaxy blue-shifted.
c        Set  the sign of zh(i) to encode this information
         if (zp_zz.lt.0.0) then
            zh(i)=-abs(zh(i))
         else
            zh(i)=abs(zh(i))
         end if
      end do
c
c     Count number of blue-shifted galaxies after the transformation
      nblue=0
      do i=1,ng
         if (zh(i).lt.0.0) then
            nblue=nblue+1
         end if
      end do
      write(0,*) 'In LG frame: nblue=',nblue,' nred=',ng-nblue

      return
      end
c------------------------------------------------------------------------------
c  Given a unit vector (px,py,pz) pointing roughly to the pole
c  of the hemisphere containing the survey, find blue-shifted
c  galaxies and change the sign of their redshift space coordinates.
c  At the same time make zh() negative to encode the fact that the
c  galaxy is blue-shifted.
c
c  If this has already been done then just count the blueshifted
c  galaxies by detecting how many have negative zh()
c     
      subroutine  detect_blueshift(ng,x,y,z,zh,px,py,pz)
*************************************variables*********************************
      implicit none
      integer ng,i,nblue
      real x(ng),y(ng),z(ng),zh(ng),px,py,pz,dir
*******************************************************************************
      nblue=0
      do i=1,ng
         dir=(px*x(i)+py*y(i)+pz*z(i))
         if (dir.lt.0.0) then
            x(i)=-x(i)
            y(i)=-y(i)
            z(i)=-z(i)
            zh(i)=-abs(zh(i))
         end if
         if (zh(i).lt.0.0)   nblue=nblue+1
      end do
      write(0,*) 'nblue=',nblue,' nred=',ng-nblue

      return
      end
c----------------------------------------------------------------------------
c Read the samplefile which determines which galaxies from file catfile
c are to be included in the redshift sample.
c
c     
      subroutine sampleread(samplefile,ng,wsel,catfile)
*************************************variables*********************************
      implicit none
      integer ng,i
      real wsel(*)
      character samplefile*(*),catfile*(*)
*******************************************************************************
      open(10,file=samplefile,status='old',form='unformatted')
      read(10) ng
      read(10) (wsel(i),i=1,ng)
      read(10) catfile(1:80)
c     upgrade to, if necessary, alter catfile string if
c     directory name differs from that of samplefile
      close(10)
      return
      end
c----------------------------------------------------------------------------
c  Combined sampleread and mockread to actually return the 
c  coordinates etc of the chosen sample.
c
c  If ic<0 ignore the sampling and just behave like mockread but with 
c  ic=icp=abs(ic)
c
c
cc  Options:
c       ic=0   -- read only ng,zmax,sigma8,omega0 and lambda0
c       ic=3   -- also read wsel, redshift space x,y,z
c       ic=4   -- also read wsel, redshift space x,y,z and true distance zh
c       ic=5   -- also read wsel, redshift space x,y,z,d and magnitude mag
c       ic=6   -- also read wsel, redshift space x,y,z,d,mag and zhmax
c       ic=7   -- also read wsel, redshift space x,y,z,d,mag,zhmax  and id
c
c       ic=-5  -- read all the galaxies from the parent mock catalogue
c                 and not just those in the current redshift sample
c                 and as for ic=5 above 
c                 also read redshift space x,y,z,d and magnitude mag
c
c  Note any array not read can be declared with unit dimension 
c  in the calling program to avoid allocating storage that isn't used. 
c eg REAL d(1)
c
c
c
c     
      subroutine sampleread2(samplefile,ng,zlim,sigma8,omega0,lambda0,
     &     x,y,z,d,mag,zhmax,id,wsel,ic)
*************************************variables*********************************
      implicit none
      integer ng,i,ic,id(*),ii,icp
      real x(*),y(*),z(*),d(*),zhmax(*),mag(*),wsel(*)
      real zlim,sigma8,omega0,lambda0
      character samplefile*(*),catfile*80
*******************************************************************************
      icp=abs(ic)
      call sampleread(samplefile,ng,wsel,catfile)
      call mockread(catfile,ng,zlim,sigma8,omega0,lambda0,
     &     x,y,z,d,mag,zhmax,id,icp)
      if (icp.gt.0) then
         ii=0
         do i=1,ng
            if (wsel(i).gt.0.0 .or. ic.lt.0) then
               ii=ii+1
               x(ii)=x(i)
               y(ii)=y(i)
               z(ii)=z(i)
               wsel(ii)=wsel(i)
               if (ic.ge.4) d(ii)=d(i)
               if (ic.ge.5) mag(ii)=mag(i)
               if (ic.ge.6) zhmax(ii)=zhmax(i)
               if (ic.ge.7) id(ii)=id(i)
            end if
         end do   
         ng=ii
         write(0,*) 'number in redshift sample=',ng
      end if
      return
      end