C ______________________________	 Convert PM format to ASCII 
C		             November 1997         A. Klypin (aklypin@nmsu.edu) 
C		                     
C          Scales internal PM coordinates and velocities 
C          to  coordinates in Mpc/h and km/s
C          In PM the coordinates are in the range 1 - (NGRID+1)
C                     velocities are P = a_expansion*V_pec/(x_0H_0)
C                     where     x_0 = comoving cell_size=Box/Ngrid
C                                    H_0 = Hubble at z=0
C                   
C	     NROW = number of particles in 1D
C	     NGRID= number of cells        in 1D
       INCLUDE 'PMparameters.h'
       REAL     INPUT
       logical inside
       Character  FileASCII*50
C...................................................................
C			Read data and open files
      WRITE (*,'(A,$)') 'Enter Name for output ASCII file = '
      READ  (*,'(A)') FileASCII
      OPEN(17,FILE=FileASCII,STATUS='UNKNOWN')
      CALL RDTAPE

      xsh = 0.
      ysh = 0.
      zsh = 0.
      write (*,*) ' RDTAPE is done'
      Fraction =INPUT(' Enter fraction of particles     =')
      If(Fraction.le.0.)Stop
      Ifraction = INT(1./Fraction+0.01)
      Radius =INPUT(' Enter radius of sphere or zero for all box    =')
      If(Radius.gt.0.)Then    
         write (*,'(A,$)') ' Enter Center: x,y,z (Mpc/h)    ='
         read  (*,*) xc,yc,zc
         xl =xc -Radius
         xr =xc +Radius
         yl =yc -Radius 
         yr =yc +Radius 
         zl =zc -Radius  
         zr =zc +Radius  
         Rad2 =Radius**2
      EndIf
      xrand     = 0.
      write (*,*) ' Every ',Ifraction,' particle will be selected'
      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)
      Box  =INPUT(' Enter box size in comoving Mpc/h =')
      BoxV =Box*100.    ! Box size in km/s
      ScaleV = BoxV/AEXPN/NGRID  ! scale factor for Velocities
c      ScaleC = Box/NGRID*AEXPN/hubble ! scale factor for Coordinates: real proper
      ScaleC = Box/NGRID ! scale factor for Coordinates: comoving h-1

c      xsh = -Box*AEXPN/hubble/2. 
c      ysh = -Box*AEXPN/hubble/2.
c      zsh = -Box*AEXPN/hubble/2.
      xsh = 0.
      ysh = 0.
      zsh = 0.

      xmax =-1.e+9
      xmin = 1.e+9
      vmax =-1.e+9
      vmin = 1.e+9
      Icount =0
      ifile      =1
      If(Nspecies.eq.0)Then ! old case: 1 complete set
         Npages =NROW          ! Number of pages
         N_in_last=NPAGE       ! Number of particles in last page
      Else                             ! multiple masses
         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      EndIf
      write (*,*) ' Pages=',Npages,' Species=',Nspecies
      write (*,*) ' N_in_last=',N_in_last
C				       Loop over particles
       DO IROW=1,Npages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            write (*,*)' Read page=',IROW,' file=',ifile,' N=',In_page
            iL = NPAGE*(IROW-1)
	      CALL GETROW(IROW,ifile)
	      DO IN=1,In_page
	        X  =ScaleC* (XPAR(IN)-1.)
	        Y  =ScaleC* (YPAR(IN)-1.)
	        Z  =ScaleC* (ZPAR(IN)-1.)
              X = X+ xsh
              Y = Y+ ysh
              Z = Z+ zsh
c              X = X*1000.
c              Y = Y*1000.
c              Z = Z*1000.
c              If(X.gt.Box)X=X-Box
c              If(X.le.0.)X=X+Box
c              If(Y.gt.Box)Y=Y-Box
c              If(Y.le.0.)Y=Y+Box
c              If(Z.gt.Box)Z=Z-Box
c              If(Z.le.0.)Z=Z+Box
             Vxs=ScaleV* VX(IN)
             Vys=ScaleV* VY(IN)
             Vzs=ScaleV* VZ(IN)
             Icount =Icount +1

             xmax =MAX(xmax,X,Y,Z)
             xmin =MIN(xmin,X,Y,Z)
             vmax =MAX(vmax,Vxs,Vys,Vzs)
             vmin =MIN(vmin,Vxs,Vys,Vzs)
             If(Nspecies.eq.0)Then   ! find weight of particle
                W =PARTW
             Else
                Ipart =IN+iL             ! current particle number
                Do i =1,Nspecies
                   if(Ipart.le.lspecies(i))Then
                      W =wspecies(i)
                      goto 50
                   endif 
                EndDo 
             EndIf
 50        If(Fraction.lt.0.999)xrand =RANDd(Nseed) ! random fraction
           inside =.true.
           If(Radius.gt.0.)Then
c              if(X.lt.xl.or.X.gt.xr)inside =.false.
c              if(Y.lt.yl.or.Y.gt.yr)inside =.false.
              if(Z.lt.zl.or.Z.gt.zr)inside =.false.
c              If(inside)dd =(X-xc)**2+(Y-yc)**2+(Z-zc)**2
c              If(dd.lt.Rad2.and.inside)Then
c                 inside=.true.
c              Else
c                 inside =.false.
c              EndIf
           EndIf
              If(xrand.le.Fraction.and.inside)
     &                  write (17,200) Ipart,X,Y,Z,Vxs,Vys,Vzs,sqrt(dd)
 200       Format(i8,3F10.4,3F9.2,f8.3)
c              If(X.lt.0..or.Y.lt.0..or.Z.lt.0.)
c     &              write(70,210)Ipart,X,Y,Z,Vxs,Vys,Vzs,W
c              If(X.gt.Box.or.Y.gt.Box.or.Z.gt.Box)
c     &              write(70,210)Ipart,X,Y,Z,Vxs,Vys,Vzs,W
              Vv =sqrt(Vxs**2+Vys**2+Vzs**2+0.001)
              If(VV.gt.10000.)
     &              write(70,210)Ipart,X,Y,Z,Vxs,Vys,Vzs,W
              
 210          format(i8,3G11.3,3G11.3,f8.3)
	      ENDDO
        ENDDO 

      write (*,*) ' Scaled Coordinates were written to:',
     &            FileASCII
      write (*,*) ' Number of particles          =',Icount
      write (*,*) ' Min/Max of coordinates(Mpc/h)=',xmin,xmax
      write (*,*) ' Min/Max of velocities (km/s) =',vmin,vmax

      END