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'
       INCLUDE 'PMrotation.h'
       REAL     INPUT
       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')
      WRITE (*,'(A,$)') 'Enter Name for input Centers     = '
      READ  (*,'(A)') FileASCII
      OPEN(30,FILE=FileASCII,STATUS='UNKNOWN')
      CALL RDTAPE


      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)

      Box  =25.
c      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
      ScaleC = Box/NGRID         ! scale factor for Coordinates

      Nlines =25
      CALL ReadMax(Ncent,Nlines)
      Radius =0.
      Do i=1,Ncent
         Radius =MAX(Rmc(i),Radius)
      EndDo

      CALL SetWeights
      CALL ReadPnt(N,ScaleC,ScaleV)                     ! particles
      CALL Points(Box,N,Ncp,Radius)            ! Make more points periodically
      CALL Rotation(Ncent,Ncp)

      Stop
      END
C--------------------------------------------------------------
C
      SUBROUTINE Rotation(Ncent,Ncp)
C             Lambda = J*sqrt(<v^2>/2)/[V_circ_vir*R_vir]
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMrotation.h'
      real*8 wm,V2,ax,ay,az,Vvx,Vvy,Vvz,
     &                ddx,ddy,ddz

      Do i=1,Ncent             ! loop over centers
         nn =0
         wm=0.
         V2 =0.
         ax =0.
         ay =0.
         az =0.
         x = Xm(i)
         y = Ym(i)
         z = Zm(i)
         wx =Vmx(i)
         wy =Vmy(i)
         wz =Vmz(i)
         R   =Rmc(i)
         R2   =R**2
         Vvx =0.
         Vvy =0.
         Vvz =0.
         ddx =0.
         ddy =0.
         ddz =0.
         Do j=1,Ncp             ! loop over particles
           dx = Xp(j)-x
           dy = Yp(j)-y
           dz = Zp(j)-z
           dd =dx**2+dy**2+dz**2
           If(dd.lt.R2)Then   ! inside halo
               nn = nn+1
               ww=iWeight(i)
               ux =VXp(j)
               uy =VYp(j)
               uz =VZp(j)
               Vvx =ux*ww +Vvx
               Vvy =uy*ww +Vvy
               Vvz =uz*ww +Vvz
               ddx =ddx + dx*ww
               ddy =ddy + dy*ww
               ddz =ddz + dz*ww
               wm=wm +ww
               V2 =V2 + (ux**2+uy**2+uz**2)*ww
               ax = ax + (dy*uz-dz*uy)
               ay = ay + (dz*ux-dx*uz)
               az = az + (dx*uy-dy*ux)
            Endif 
         EndDo
         If(nn.eq.0)Then
            write(*,*) ' Error: no particles'
            stop
         Else
            ax = (ax - (ddy*Vvz-ddz*Vvy)/wm)/wm*AEXPN
            ay = (ay - (ddz*Vvx-ddx*Vvz)/wm)/wm*AEXPN
            az = (az - (ddx*Vvy-ddy*Vvx)/wm)/wm*AEXPN
            aJm = sqrt(ax**2+ay**2+az**2)
            V2 = (V2 - (Vvx**2+Vvy**2+Vvz**2)/wm)/wm
            V      =sqrt(max(V2,1.e-10))
            Vcirc =6.581e-5*sqrt(Amc(i)/Rmc(i)/AEXPN)
            aLambda1(i) =aJm*V/sqrt(2.)/Vcirc**2/Rmc(i)/AEXPN
            aLambda2(i) =aJm/Vcirc/Rmc(i)/AEXPN
            aJ(i) =aJm
            Vvx =Vvx/wm
            Vvy =Vvy/wm
            Vvz =Vvz/wm
            write (*,100) i,nn,wm,Vcirc,V,m1(i),m2(i),m3(i),
     &            aJ(i),aLambda1(i) ,aLambda2(i), 
     &            Vvx,Vvy,Vvz,Vmx(i),Vmy(i),Vmz(i)
            write (17,120) AEXPN,x,y,z, Vvx,Vvy,Vvz,aMc(i),
     &                          Rmc(i)*1000.,Vrm(i),Vcmax(i),
     &                          m1(i),m2(i),m3(i),m4(i),
     &                          aLambda1(i),aJ(i),ax,ay,az
      
         EndIf
       EndDo
 100   format(i3,i8,g11.3,2f8.2,2i8,i5,3g11.3,/
     &        3x,3F8.2,3x,3F8.2,g11.3)
 120   format(F7.4,F8.3,2F8.3,3F7.1,g10.3,f8.2,2F7.2,2i8,i7,i5,
     &       f7.4,4g12.4 )
      Return
      End
C--------------------------------------------------------------
C                       Add points around the Box to take into account periodical conditions
C  It replicates points periodically and keeps those which are at  distance
C  less than Radius from the Box. 
C                        N     = number of points inside the Box
C                        Ncp = total number of points
      SUBROUTINE Points(Box,N,Ncp,Radius)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMrotation.h'
      Logical    Inside
         B1  =-Radius
         B2  =Box+Radius
         B3  =Box-Radius
         Ncp =N
      Do i=1,N               ! Add dark matter particles
         x0 =xp(i)
         y0 =yp(i)
         z0 =zp(i)
         Inside =            (x0.gt.Radius).and.(x0.lt.B3)
         Inside =Inside.and.((y0.gt.Radius).and.(y0.lt.B3))
         Inside =Inside.and.((z0.gt.Radius).and.(z0.lt.B3))
         If(.not.Inside)THEN         
         ux =VXp(i)
         uy =VYp(i)
         uz =VZp(i)
         ww=iWeight(i)
            Do k1 = -1,1
               zr =z0+k1*Box
               If(zr.GT.B1.AND.zr.LT.B2)Then
               kk =k1**2
               Do j1 = -1,1
                  yr =y0+j1*Box
                  If(yr.GT.B1.AND.yr.LT.B2)Then
                     jj = kk +j1**2
                     Do i1 = -1,1
                        xr =x0+i1*Box
                        ii =jj +i1**2
                        If((xr.GT.B1.AND.xr.LT.B2).AND.ii.ne.0)
     .                     Then
                             Ncp =Ncp +1
                             If(Ncp.le.Np)Then
                               xp(Ncp)  =xr
                               yp(Ncp)  =yr
                               zp(Ncp)  =zr
                               VXp(Ncp)  =ux
                               VYp(Ncp)  =uy
                               VZp(Ncp)  =uz
                               iWeight(Ncp) =ww
                             EndIf
                        EndIf
                     EndDo
                  EndIf
               EndDo
               EndIf
            EndDo
         EndIf
      EndDo
      If(Ncp.Gt.Np)Then
         write(*,*)' Too many points:',Ncp,' Maximum is =', Np
         STOP
      EndIf
      Return
      End
C--------------------------------------------------
C                             Set Weights of particles for
C                             fast access
      SUBROUTINE SetWeights
C--------------------------------------------------
      INCLUDE 'PMparameters.h'

      Wpr =NGRID**3/FLOAT(NROW**3)
      If(Nspecies.eq.0)Then  ! old  constant weights
         N_particles =lspecies(1)
         Ww  = wspecies(1)/Wpr
         Do i=1,N_particles
            iWeight(i) =Ww
         EndDo
      Else
         N_particles =lspecies(Nspecies)
         jstart =1
         write (*,*) ' N_Particles=',N_particles
         If(N_particles.le.0)Then
            write (*,*) ' Wrong number of particles:',N_particles
            write (*,*) ' Nspecies=',Nspecies,' N=',lspecies
            STop
         endif 
         Do j=1,Nspecies
            jend =lspecies(j)
            Do k=jstart ,jend
               iWeight(k) =wspecies(j)/Wpr
            EndDo
            jstart =jend
         EndDo
      EndIf
      write (*,*) ' Set Weights for ',N_particles,' particles:',
     &       ' w_first=' ,iWeight(1),' w_last=',  iWeight(N_particles)
      RETURN
      END

C---------------------------------------------------------
C                                        Read particles
      SUBROUTINE ReadPnt(N,ScaleC,ScaleV)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMrotation.h'

      N =0
         N_particles =lspecies(Nspecies)   ! Total number of particles
         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      write (*,*) ' Pages=',Npages,' Species=',Nspecies
c      write (*,*) ' N_in_last=',N_in_last

      DO  IROW=1, Npages         ! Loop over particle pages
            In_page =NPAGE
            If(IROW.eq.Npages)In_page =N_in_last
            iL = NPAGE*(IROW-1)
C				       Loop over particles
         READ  (21,REC=IROW) RECDAT
	   DO   IN=1,In_page
            N   = N +1
            Xp(N)  =ScaleC*(XPAR(IN)-1.)
            Yp(N)  =ScaleC*(YPAR(IN)-1.)
            Zp(N)  =ScaleC*(ZPAR(IN)-1.)
            VXp(N)=ScaleV*VX(IN)
            VYp(N)=ScaleV*VY(IN)
            VZp(N)=ScaleV*VZ(IN)
                Ipart =IN+iL                     ! current particle number
             iWeight(N) =iWeight(Ipart)   ! particles weight            
         ENDDO
      ENDDO 
      Return
      End

C---------------------------------------------------------
C                                        Read list of maxima
c                                     Xnew =X-1: coord. now are 0-NGRID       
      SUBROUTINE ReadMax(Ncent,Nlines)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMrotation.h'
      Character*100 Line

      Do i=1,Nlines
         read(30,'(A)') Line
         write(*,'(A)') Line
      EndDo
       Ncent =0
      xxm =1000.
      xxn =-1000.
 10   read (30,*,err=5,end=5) Aa,X0, Y0, Z0, VvX,VvY,VvZ,aM,
     &                          Rv,V1,V2,nn1,nn2,nn3,nn4
            Ncent=Ncent+1
            If(Ncent.le.Nc)Then
               Xm(Ncent)  =X0
               Ym(Ncent)  =Y0
               Zm(Ncent)  =Z0
               Vmx(Ncent)=VvX
               Vmy(Ncent)=VvY
               Vmz(Ncent)=VvZ
               Amc(Ncent)=aM
               Rmc(Ncent)=Rv / 1000.   ! scale to Mpc/h
               Vrm(Ncent)=V1
               Vcmax(Ncent)=V2
               m1(Ncent)=nn1
               m2(Ncent)=nn2
               m3(Ncent)=nn3
               m4(Ncent)=nn4
               xxm    =min(xxm,Xm(Ncent))
               xxn     =max(xxn,Xm(Ncent))
            EndIf
            GoTo 10
 5    write(*,*) ' Read Centers: X{Min,Max}=',xxm,xxn,' Nc=',Ncent
      If(Ncent.GT.Nc)Then
         write (*,*) ' Too many centers: Max=',Nc
         Stop
      EndIf
      Return
      End
