C ______________________________	 Selects and Scales PM particles 
C                       For setting refinement in initial conditions
C		             December 1998        A. Klypin (aklypin@nmsu.edu) 
C          Input: 1) catalog of halos, which define what regions should be
C                                 selected for mass refinement
C                     2) PMcrd.DAT and PMcrs0.DAT -- 
C                                 lowest resolution run with one specie
c                                 of particles covering all the computational Box
C
C          Output: File with lagrangian coordinates of selected particles.
C                     These coordinates are used by PMstartM to make mass refinement
C          Format of halo catalog:
C                        First 16 lines -- any text (ignored)
C                        Next line(s): one line per halo (=center):
C                        x,y,z,vx,vy,vz,Mass,Radius,v1,v2
C                            Only coordinates (x,y,z) and Radius are used.
C                            (x,y,z) = center of sphere of radius Radius
C                        Units: (x,y,z) are in Mpc/h; Radius is in kpc/h
C  Within the sphere(s) all particles are selected. In the lagrangian
C  space occupied by the particles mass resolution will be increased.
C  Spheres may overlap in any possible way    
       INCLUDE 'PMparameters.h'
       REAL     INPUT
       Character*70  FileASCII,Line,FileCatalog
C...................................................................
C			Read data and open files
       WRITE (*,'(A,$)') ' Enter Output File Name = '
      READ  (*,'(A)') FileASCII
      OPEN(17,FILE=FileASCII,STATUS='UNKNOWN')
       WRITE (*,'(A,$)') ' Enter File Name for Input Halo Catalog = '
      READ  (*,'(A)') FileCatalog
      OPEN(18,FILE=FileCatalog,STATUS='UNKNOWN')
      CALL RDTAPE
      CALL SetWeight
      write (*,*) ' RDTAPE is done'
      WRITE (17,100) HEADER,
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
     +                  Ocurv
100   FORMAT(1X,'Header=>',A45,/
     +            1X,' A=',F8.3,' A0=',F8.3,' Ampl=',F8.3,' Step=',F8.3,/
     +            1X,' I =',I4,' WEIGHT=',F8.3,' Ekin=',E12.3,/
     +            1X,' Nrow=',I4,' Ngrid=',I4,' Nrecl=',I6,/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,
     +            1x,' Omega_curvature=',F7.3,/
     +            1x,'   Npart   Euler Coordinates (Mpc/h)',
     +                 '      Velocities (km/s)',
     +                 '    Lagrangian Coords (limits: 1-N_1d)')
      If(extras(100).gt.0.)Then
         Box =extras(100)
         write(*,*) '  Box size =',Box,' Mpc/h'
      Else
         Box  =INPUT(' Enter box size in comoving Mpc/h =')
      EndIf

      IF(Nspecies.ne.1)Then
         write (*,*) '--- You can use this program only for runs with',
     &     ' one specie of particles at the lowest mass resolution.'
         write (*,*)' --- The number of species is ',Nspecies,' STOP'
         Stop
      Endif 
      BoxV    = Box*100.    ! Box size in km/s
      ScaleV = BoxV/AEXPN/NGRID  ! scale factor for Velocities
      ScaleC = Box/NGRID         ! scale factor for Coordinates
      Do i=1,17       ! skip header lines
         read(18,'(A)') Line
         write(*,'(A)') Line
      EndDo
 200  read(18,*,err=300,end=300) Xcentr,Ycentr,Zcentr,
     +                   Vvx0,Vvy0,Vvz0,Am,Radius,Vrms,Vcirc
c 200  read(18,*,err=300,end=300) Xcentr,Ycentr,Zcentr,
c     +                   Am,Vcir,Radius
         Radius =Radius/1000.  ! Scale radius to Mpc/h
         rad2   = Radius**2
         xxmax  =-1.e+9
         xxmin  = 1.e+9
         vmax   =-1.e+9
         vmin   = 1.e+9
         xrand  =0.
         Icount =0
         ifile  =1

            N_particles =lspecies(Nspecies)   ! Total number of particles
            Npages = (N_particles -1)/NPAGE +1
            N_in_last=N_particles -NPAGE*(Npages-1)
            N_line=INT((lspecies(1))**0.33333+0.5)
            N2          =N_line**2
      write (*,*) ' Pages=',Npages,' Species=',Nspecies
      write (*,*) ' Nparticles=',N_particles,' N_in_last=',N_in_last
      write (*,*) ' Particles_1d=',N_line


      DO  IROW=1, Npages         ! Loop over particle pages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
c            write (*,*)' Move page=',IROW,' file=',ifile,' N=',In_page
            iL = NPAGE*(IROW-1)
         CALL GETROW(IROW,ifile) ! read in a page of particles
	    DO  IN=1, In_page          ! Loop over particles
	        X  =ScaleC* (XPAR(IN)-1.)         ! scale coordinates
	        Y  =ScaleC* (YPAR(IN)-1.)
	        Z  =ScaleC* (ZPAR(IN)-1.)
                                                                     ! periodical conditions
              dX =min(ABS(X-Xcentr),Box-ABS(X-Xcentr))
              dY =min(ABS(Y-Ycentr),Box-ABS(Y-Ycentr))
              dZ =min(ABS(Z-Zcentr),Box-ABS(Z-Zcentr))
              dd =dx**2 +dY**2 +dZ**2
              IF(dd.lt.rad2)Then                !  get the particle
                 Ipart =IN+iL                     ! current particle number
                WPAR =iWeight(Ipart)   ! particles weight
                 Vxs=ScaleV* VX(IN)                    ! scale velocities
                 Vys=ScaleV* VY(IN)
                 Vzs=ScaleV* VZ(IN)
                 Icount =Icount +1                          ! count particles
                 xxmax =MAX(xxmax,X,Y,Z)
                 xxmin =MIN(xxmin,X,Y,Z)
                 vmax =MAX(vmax,Vxs,Vys,Vzs)
                 vmin =MIN(vmin,Vxs,Vys,Vzs)
                                                                    ! find lagrangian coords
c                   Jpart =N_particles -Ipart+1 ! reverse the order of particles
                   Jpart =Ipart ! direct  order of particles
                    k1 =INT((Jpart-1)/N2)+1
                    jj   =Jpart-(k1-1)*N2
                    j1  =INT((jj-1)/N_line)+1
                    i1  =jj-(j1-1)*N_line
                   write (17,250) Ipart,X,Y,Z,Vxs,Vys,Vzs,i1,j1,k1
     &                                                  ! write to a file
              EndIf
 250          Format(i9,3F9.4,3F9.2,3i5,3i9)
        ENDDO 
      ENDDO
      write (*,*) ' Scaled Coordinates were written to:',
     &            FileASCII
      write (*,*) ' Number of particles          =',Icount
      write (*,*) ' Center: ', Xcentr,Ycentr,Zcentr,' R=',Radius
      write (*,*) ' Min/Max of coordinates(Mpc/h)=',xxmin,xxmax
      write (*,*) ' Min/Max of velocities (km/s) =',vmin,vmax

      goto 200      ! get next center
 300  write (*,*) ' <======= the End ========> '
      END


