C ______________________________ Profiles for Galaxy Simulations
C		             
C		                     
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMcool.h'
      REAL         AMprevious(Np),AMlast(Np),Daverg(0:500)
      Real*8       Xc,Yc,Zc
      DIMENSION    Iaverg(0:500)
      REAL         INPUT
      Character    FileASCII*100
      Character  File1*100,File2*100,File3*100
      EQUIVALENCE (Box,extras(100))



C...................................................................
C			Read data and open files
       File3 ='cooling_part.dat' 
       write(*,*)' enter 3 digits of the snapshot'
       read(*,*) moment
         write(File1,'("FILES/PMcrda0.",i3.3,".DAT")')moment
         write(File2,'("FILES/PMcrs0a0.",i3.3,".DAT")')moment
         write(*,*) File1
         write(*,*) File2


c      File1 ='PMcrd.DAT'
c      File2 ='PMcrs0.DAT'


      CALL RDTAPE(File1,File2)
      iSnapshot =INT(AEXPN*10000.)
         write(FileASCII,'("RESULTS/Result.DM.",i5.5,".dat")')
     &    iSnapshot
      OPEN(37,FILE=FileASCII,STATUS='UNKNOWN') 
 
        write(FileASCII,'("RESULTS/Result.SPH.",i5.5,".dat")')
     &    iSnapshot
      OPEN(39,FILE=FileASCII,STATUS='UNKNOWN') 

         write(FileASCII,'("RESULTS/Result.CYL.",i5.5,".dat")')
     &    iSnapshot
      OPEN(41,FILE=FileASCII,STATUS='UNKNOWN') 

      OPEN(18,FILE='Stat.dat',STATUS='UNKNOWN',access='append')
      OPEN(20,FILE='particles.dat',STATUS='UNKNOWN')
      OPEN(30,FILE='particlesu.dat',STATUS='UNKNOWN',form='UNFORMATTED')
      OPEN(10,file=File3,form='UNFORMATTED',status='OLD')

       time      = (age(AEXPN)-age(0.6))/1.e9
      box          = extras(100)! simulation box
      xfrac        = extras(51) ! fraction of particles to be cooled
      cool_frac    = extras(52) ! total cooling
      cfrac        = extras(53) ! cooling per step
      aexpn_s      = extras(54) ! start cooling
      aexpn_e      = extras(55) ! and cooling
      i_cool_which = extras(57) ! cool every # step
      i_cool_all   = extras(58) ! total # of cooling steps
      xm_part      = extras(59) ! particle weight (h^-1 M_s)
      CALL ReadCool

      i_halo = 1   ! only one halo cooled
      n_cooled = iph_cool(i_halo,1)
      xm_part      =1.e12/lspecies(Nspecies)

      aMassDisk = xm_part/hubble*n_cooled
      write (*,*)  ' hubble =',hubble,' mass of disk=',aMassDisk
      write (*,*) ' RDTAPE is done Time=',time,' a=',AEXPN
c      Do i=37,41,2
c          WRITE (i,100) HEADER,AEXPN,AEXP0,AMPLT,
c     +                  ASTEP,ISTEP,PARTW,EKIN,
c     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,
c     +                  Ocurv,n_cooled,aMassDisk,time
c      EndDo 
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=',I8,/
     +            1x,' Omega_0=',F7.3,' OmLam_0=',F7.4,' Hubble=',f7.3,/
     +            1x,' Omega_curvature=',F7.3,/
     + 1x,' Ndisk =',i8,' Mass disk =',G11.3,' Age=',f7.3,'Gyr') 
c      Box  =INPUT(' Enter box size in comoving Mpc/h =')
c      Box =0.5 ! NG3109
c      Box = 2   ! 'bulge'

      BoxV   = Box*100.    ! Box size in km/s
      ScaleV = BoxV/AEXPN/NGRID      ! scale factor for Velocities (km/s)
      ScaleC = Box/NGRID*AEXPN/hubble  ! scale factor for Coordinates (real Mpc)
      ScaleH  = 100.*hubble*sqrt(Om0/AEXPN**3+Oml0) ! H(a)=km/s/Mpc
      ScaleM = xm_part/hubble

 
      CALL ReadPnts(N,ScaleC,ScaleV,ScaleM)

      Call ChangePnts(N,Box,ScaleH)
      Call DumpPnts(N,Box,ScaleM)
      Call Statist(N,Box,ScaleM)
c      Call AngMom(N,30000,3.,Box,ScaleM)
      Call Profile(N,Box,ScaleM)
c      Call Projection(N,Box,0.)
c      Call Projection(N,Box,5.)

      END
c-------------------------------------------------------------------- 
C                            Get basic statistics of the system
      SUBROUTINE Statist(N,Box,ScaleM)
c-------------------------------------------------------------------- 
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      Real*8 Lx,Ly,Lz,aL,V2,Lz1,Lz2,Lz3,Lz5,aM1,aM2,aM3,aM5
      Real*8 Lxh,Lyh,Lzh,aLh,dRh,dR2h,Wh,Lhalo
      Real*8 Xcc,Ycc,Zcc,Vcx,Vcy,Vcz,dR,dR2,X,Y,Z
      Real*8 Lzh1,Lzh2,Lzh3,Lzh5,aMh1,aMh2,aMh3,aMh5
      Real*8 R1,R2,R3,R5,Rh1,Rh2,Rh3,Rh5
      Real*8 Xx2,Yy2,Zz2

      DATA Vcx,Vcy,Vc,dR,dR2,V2/6*0.d+0/
      DATA Lx,Ly,Lz,Lxh,Lyh,Lzh,Wh,dRh,dR2h/9*0.d+0/
      DATA Lz1,Lz2,Lz3,Lz5,aM1,aM2,aM3,aM5/8*0.d+0/
      DATA Lzh1,Lzh2,Lzh3,Lzh5,aMh1,aMh2,aMh3,aMh5/8*0.d+0/
      DATA R1,R2,R3,R5,Rh1,Rh2,Rh3,Rh5/8*0.d+0/

      Rad1 =1.
      Rad2 =2.
      Rad3 =3.
      Rad5 =5. 
      Xcc = Box/2.
      Ycc = Box/2.
      Zcc = Box/2.
      Xx2=0.               ! find moments of inertia
      Yy2=0.
      Zz2=0.
      n_cooled = 0

      Do jp=1,N
         X = (Xp(jp)- Xcc) *1000.  ! scale to kpc inits
         Y = (Yp(jp)- Ycc) *1000.
         Z = (Zp(jp)- Zcc) *1000.
            aLx =Y*VZp(jp)  -Z*VYp(jp)
            aLy =Z*VXp(jp)   -X*VZp(jp)
            aLz =X*VYp(jp)   -Y*VXp(jp)
             dd = X**2 + Y**2 + Z**2
             ds = sqrt(dd)
                     If(iMask(jp).gt.0)Then
                        n_cooled =n_cooled +1
                         Lx =Lx +aLx
                         Ly =Ly +aLy
                         Lz =Lz +aLz
                        dR = dR + ds
                        dR2= dR2+ dd
                        If(ds.le.Rad5)Then
                           Lz5 =Lz5 +aLz
                           aM5 = aM5 + 1
                           R5  = R5 + ds
                           If(ds.le.Rad3)Then
                              Lz3 =Lz3 +aLz
                              aM3 = aM3 + 1
                              R3  = R3 + ds
                               If(ds.le.Rad2)Then
                                 Lz2 =Lz2 +aLz
                                 aM2 = aM2 + 1
                                 R2  = R2 + ds
                                  If(ds.le.Rad1)Then
                                    Lz1 =Lz1 +aLz
                                    aM1 = aM1 + 1
                                    R1  = R1 + ds
                                  EndIf
                               EndIf
                            EndIf
                         EndIf
                    Else
                        W = iWeight(jp)
                         Lxh =Lxh +aLx*W
                         Lyh =Lyh +aLy*W
                         Lzh =Lzh +aLz*W
                        dRh = dRh + ds*W
                        dR2h= dR2h+ dd*W
                        Vcx  = Vcx + VXp(jp) *W
                        Vcy  = Vcy + VYp(jp) *W
                        Vcz  = Vcz + VZp(jp) *W
                        Wh = Wh + iWeight(jp)
                        If(ds.le.Rad5)Then
                           Lzh5 =Lzh5 +aLz*W
                           aMh5 = aMh5 + W
                           Rh5  = Rh5 + ds*W
                           If(ds.le.Rad3)Then
                              Lzh3 =Lzh3 +aLz*W
                              aMh3 = aMh3 + W
                              Rh3  = Rh3 + ds*W
                               If(ds.le.Rad2)Then
                                 Lzh2 =Lzh2 +aLz*W
                                 aMh2 = aMh2 + W
                                 Rh2  = Rh2 + ds*W
                                  If(ds.le.Rad1)Then
                                    Lzh1 =Lzh1 +aLz*W
                                    aMh1 = aMh1 + W
                                    Rh1  = Rh1 + ds*W
                                  EndIf
                               EndIf
                            EndIf
                         EndIf
                     EndIf
         V2 = V2 +  VXp(jp)**2 +VYp(jp)**2+VZp(jp)**2
      EndDo
      
      Lx =Lx/n_cooled
      Ly =Ly/n_cooled
      Lz =Lz/n_cooled
      V2 =sqrt(V2/N)
      aL = sqrt(Lx**2 +Ly**2 +Lz**2)
      dR = dR/n_cooled
      Lhalo = sqrt(Lxh**2 +Lyh**2 +Lzh**2)
      
      write (*,110) 
     &           dR,sqrt(dR2/n_cooled),
     &           V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*n_cooled,
     &            aL,
     &            Lhalo,Lxh,Lyh,Lzh,
     &            Lhalo/Wh,Lxh/Wh,Lyh/Wh,Lzh/Wh
      write (37,110) 
     &           dR,sqrt(dR2/n_cooled),
     &           V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*n_cooled,
     &            aL,
     &            Lhalo,Lxh,Lyh,Lzh,
     &            Lhalo/Wh,Lxh/Wh,Lyh/Wh,Lzh/Wh
 110  format(
     &           5x,' Aver Distance Disk Particles(kpc)=',g11.4,
     &                  ' rms distance(kpc)=',g11.4,/
     &           5x,' rms Velocity=',g11.4,' km/s',/
     &           5x,' Disk Ang.Momentum: direction=',3g11.3,/
     &           5x,' Disk Ang.Momentum: L/<R> =',g11.4,' km/s',/
     &           5x,' Disk Ang.Momentum: L     =',g11.4,' kpc*km/s',/
     &           5x,' Disk Ang.Momentum: L/M   =',g11.4,' kpc*km/s',/
     &           5x,' Halo Ang.Momentum: L     =',4g11.4,' kpc*km/s',/
     &           5x,' Halo Ang.Momentum: L/M   =',4g11.4,' kpc*km/s')
      write (*,120) Vcx/Wh, Vcy/Wh,Vcz/Wh,Wh
      write (37,120) Vcx/Wh, Vcy/Wh,Vcz/Wh,Wh

 120  format(
     &           5x,' Aver Velocity of Dm Particles(km/s)=',3g11.3,
     &              ' Total Weight=',g11.3)
      write (37,140)      
      write (37,130) Rad1,
     &           Lz1,Lz1/max(aM1,1.d+0),Lz1/max(R1,1.d-10),
     &           R1/max(aM1,1.d+0),aM1,
     &           Lzh1,Lzh1/max(aMh1,1.d+0),Lzh1/max(Rh1,1.d-10),
     &           Rh1/max(aMh1,1.d+0),aMh1
      write (37,130) Rad2,
     &           Lz2,Lz2/max(aM2,1.d+0),Lz2/max(R2,1.d-10),
     &           R2/max(aM2,1.d+0),aM2,
     &           Lzh2,Lzh2/max(aMh2,1.d+0),Lzh2/max(Rh2,1.d-10),
     &           Rh2/max(aMh2,1.d+0),aMh2
      write (37,130) Rad3,
     &           Lz3,Lz3/max(aM3,1.d+0),Lz3/max(R3,1.d-10),
     &           R3/max(aM3,1.d+0),aM3,
     &           Lzh3,Lzh3/max(aMh3,1.d+0),Lzh3/max(Rh3,1.d-10),
     &           Rh3/max(aMh3,1.d+0),aMh3
      write (37,130) Rad5,
     &           Lz5,Lz5/max(aM5,1.d+0),Lz5/max(R5,1.d-10),
     &           R5/max(aM5,1.d+0),aM5,
     &           Lzh5,Lzh5/max(aMh5,1.d+0),Lzh5/max(Rh5,1.d-10),
     &           Rh5/max(aMh5,1.d+0),aMh5
           time = Age(AEXPN)-Age(0.60)
      write (18,150) AEXPN,time/1.e9,ISTEP,aL*n_cooled,
     &           bar_angle*180./3.1415926,
     &           bar_length,sqrt(1.-bar_ecc),
     &           Rad1,Lz1/max(aM1,1.d+0),aM1,
     &           Rad2,Lz2/max(aM2,1.d+0),aM2,
     &           Rad5,Lz5/max(aM5,1.d+0),aM5


 150  format(2f8.4,i4,g12.5,3f8.3,4(2x,f6.1,g11.3,g11.3))
 140  format(' Rad(kpc) Disk:Lz',4x,'Lz/M',2x,'<Vrotat>km/s',1x,
     &           '<R>kpc  Npart',3x,' Halo:  Lz',4x,'Lz/M',2x,
     &           '<Vrotat>km/s',1x,'<R>kpc  Npart')
 130  format(f7.2,2(2g11.3,2f7.2,g11.3,3x))
      Return
      End
C------------------------------------------
C                  write particles in a file       
      SUBROUTINE DumpPnts(N,Box,ScaleM)
C------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      Real*8 MassTot

      Xcc = Box/2.
      Ycc = Box/2.
      Zcc = Box/2.

      MassTot =0.
      Do jp=1,N
         MassTot = MassTot + iWeight(jp)
      EndDo
      MassTot = MassTot * ScaleM

      i_halo = 1   ! only one halo cooled
      n_cooled = iph_cool(i_halo,1)

      write(20,100)
      write(20,110) N,n_cooled,age(AEXPN)-age(0.6),ISTEP,
     &                  ScaleM, ScaleM*n_cooled, MassTot 
      write(20,115) 

c      CALL GetDensity(Box)
      Do jp=1,N
         If(iMask(jp).gt.0)Then   ! DM particle
         xx = (Xp(jp)-Xcc)*1000.
         yy = (Yp(jp)-Ycc)*1000.
         zz = (Zp(jp)-Zcc)*1000.


         write (20,130)  xx,yy,zz,
     &              VXp(jp) ,VYp(jp),VZp(jp)
c                   xx*VYp(jp)-yy*VXp(jp) 
c     &                          INT(iWeight(jp)+0.1)
         EndIf 
      EndDo

 100  format(' TotalPart  N: Disk ',
     &             T32,'Time',T42,'Step',
     &             T48,'MassParticle',T62,'MassDisk',
     &             T74,'MassTotal',/
     &             T32,'(yrs)',T48,'(Msun)',T62,'(Msun)',
     &             T74,'(Msun)')
 115  format(T3,'X(kpc)       Y          Z',
     &             T38,'Vx(kms)         Vy           Vz')
 110  format(i10,i7,T30,1P,g11.3,T41,i4,
     &              T48,g11.4,T62,g11.4,T74,g11.4)
 120  format(1P,3g13.5,2x,3g12.4,i3)
 130  format(3g13.5,2x,4g12.4,i3)
      Return
      End

C------------------------------------------
C                            Change positions: center, inclination
      SUBROUTINE ChangePnts(N,Box,ScaleH)
C------------------------------------------
      PARAMETER (pi =3.14159265) ! 
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      Real*8 xshift,yshift,zshift,Lx,Ly,Lz,aL,Xc,Yc,Zc,V2
      Real*8 Xcc,Ycc,Zcc,Vcx,Vcy,Vcz,dR,dR2,aCos,Ld
      Real*8 Xx,Yy,Zz,XY

      Ncenter = N/2
      If(Ncenter.eq.0)write (*,*) ' Error!!  No particles in ChangePnts'

      Xcc = Box/2.
      Ycc = Box/2.
      Zcc = Box/2.
      V2 =0.
      xmax =-1.d+9
      xmin = 1.d+9
      ymax =-1.d+9
      ymin = 1.d+9
      zmax =-1.d+9
      zmin = 1.d+9
      Rselect = 0.250
      Xc0       = Xcc
      Yc0       = Ycc
      Zc0       = Zcc
      R2        = Rselect**2
      Ncount    = 0
      dn        =0.
 300  nn          = 0
      dn0 =dn
      Xc = 0.
      Yc = 0.
      Zc = 0.
      Lx = 0.
      Ly = 0.
      Lz = 0.
      dR = 0.
      dR2 = 0.
      Vcx = 0.
      Vcy = 0.
      Vcz = 0.
c                                         Find center of mass and mean V
      Do jp=1,Ncenter
         dd = (Xp(jp) -Xc0)**2 +(Yp(jp) -Yc0)**2 +(Zp(jp) -Zc0)**2
         If(dd.lt.R2)Then
            nn = nn + 1
           Xc = Xc +Xp(jp)
           Yc = Yc +Yp(jp)
           Zc = Zc +Zp(jp)
           Vcx = Vcx + Vxp(jp)
           Vcy = Vcy + Vyp(jp)
           Vcz = Vcz + Vzp(jp)
         EndIf
      EndDo
      If(nn.eq.0)write (*,*) ' Error:  no particles found. R=',Rselect
      Xc = Xc/nn
      Yc = Yc/nn
      Zc = Zc/nn
      Vcx = Vcx/nn
      Vcy = Vcy/nn
      Vcz = Vcz/nn
      xshift = Xcc - Xc   ! displacements to bring
      yshift = Ycc - Yc   ! the center to {Xcc,Ycc,Zcc}
      zshift = Zcc - Zc
c      write (*,100) nn,Rselect*1000.,nn/Rselect**3,Xc,Yc,Zc,
c     &                   1000.*xshift,1000.*yshift,1000.*zshift,
c     &                    Vcx,Vcy,Vcz
 100  format('--- Center: n=',i8,f7.2,g11.2,
     &                 '  X(Mpc)=',3g12.6,/
     &           25x,'shift(kpc)             =',3g11.5,/
     &           25x,' average Velocity(km/s)=',3g11.4)
      Xc0       = Xc
      Yc0       = Yc
      Zc0       = Zc
      Ncount = Ncount +1
      Rselect = Rselect/1.1
      R2         = Rselect**2
      dn          = nn/Rselect**3
      If(Ncount.lt.180.and.Rselect.gt.0.2e-3.and.nn.gt.20)goto300
       write (37,100) nn,Rselect*1000.,nn/Rselect**3,Xc,Yc,Zc,
     &                   1000.*xshift,1000.*yshift,1000.*zshift,
     &                    Vcx,Vcy,Vcz
       n_cooled =0
      Do jp=1,N
         Xp(jp) =Xp(jp) +xshift
         Yp(jp) =Yp(jp) +yshift
         Zp(jp) = Zp(jp) +zshift
         VXp(jp) = VXp(jp) +ScaleH*(Xp(jp) -Xcc)-Vcx
         VYp(jp) = VYp(jp) +ScaleH*(Yp(jp) -Ycc)-Vcy
         VZp(jp) = VZp(jp) +ScaleH*(Zp(jp) -Zcc)-Vcz
                  IF(Xp(jp).lt.0.)Xp(jp)=Xp(jp) +Box 
                  IF(Yp(jp).lt.0.)Yp(jp)=Yp(jp) +Box 
                  IF(Zp(jp).lt.0.)Zp(jp)=Zp(jp) +Box
                  IF(Xp(jp).GE.Box)Xp(jp)=Xp(jp) -Box 
                  IF(Yp(jp).GE.Box)Yp(jp)=Yp(jp) -Box 
                  IF(Zp(jp).GE.Box)Zp(jp)=Zp(jp) -Box
                     If(iMask(jp).gt.0)Then   ! cooled particle
                        n_cooled = n_cooled +1
                         Lx =Lx +(Yp(jp)-Ycc)*VZp(jp)
     &                             -(Zp(jp)-Zcc)*VYp(jp)
                         Ly =Ly +(Zp(jp)-Zcc)*VXp(jp)
     &                             -(Xp(jp)-Xcc)*VZp(jp)
                         Lz =Lz +(Xp(jp)-Xcc)*VYp(jp)
     &                             -(Yp(jp)-Ycc)*VXp(jp)
                        dd = (Xp(jp)-Xcc)**2 +
     &                        (Yp(jp)-Ycc)**2 + (Zp(jp)-Zcc)**2
                        dR = dR + sqrt(dd)
                        dR2= dR2+ dd
                        V2 = V2 + VXp(jp)**2+VYp(jp)**2+VZp(jp)**2
                     EndIf
         xmin =min(xmin,Xp(jp))
         ymin =min(ymin,Yp(jp))
         zmin =min(zmin,Zp(jp))
         xmax =max(xmax,Xp(jp))
         ymax =max(ymax,Yp(jp))
         zmax =max(zmax,Zp(jp))
      EndDo


     
      Lx =Lx/n_cooled
      Ly =Ly/n_cooled
      Lz =Lz/n_cooled
      V2 =sqrt(V2/n_cooled)
      aL = sqrt(Lx**2 +Ly**2 +Lz**2)
      dR = dR/n_cooled
      write (*,110) xmin,xmax,ymin,ymax,zmin,zmax,
     &           1000.*dR,1000.*sqrt(dR2/n_cooled),
     &           V2,Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*1000.*n_cooled
 110  format(' --- All particles:  Xmin/max=',2g11.4,/
     &           30x,2g11.4,/30x,2g11.4/,
     &           25x,' Aver Distance Disk Particles(kpc)=',g11.4,
     &                  ' rms distance(kpc)=',g11.4,/
     &           25x,' rms Velocity=',g11.4,' km/s',/
     &           25x,' Disk Ang.Momentum: direction=',3g11.4,/
     &           25x,' Disk Ang.Momentum: L/<R> =',g11.4,' km/s',/
     &           25x,' Disk Ang.Momentum: L     =',g12.5,' kpc*km/s')
c                                      Rotate the system of coordinates
         return
c
      If(Lz.le.0.5*aL)Then
         write (*,*) '         z-component of angular momentum',
     &     '                is too small. Something is very wrong.'

c         Stop
      EndIf
      aCos = 0.99999d-0
      If(Lz.gt.aCos*aL)Then
         write (*,*) '     Tilt is too small (Cos>',aCos,'). No change' 
         Return
      EndIf
c                                 Find maxtrix for rotation
      Lx =Lx/aL
      Ly =Ly/aL
      Lz =Lz/aL
c      Ld =sqrt(Lx**2 +Ly**2)
c      If(Ly.gt.0.)Then
c         ex1 =  Ly/Ld
c         ey1 = -Lx/Ld
c         ex2 =  Lx*Lz/Ld
c         ey2 = Ly*Lz/Ld
c         ez2 = -Ld
c      Else
c         ex1 = -Ly/Ld
c         ey1 =  Lx/Ld
c         ex2 =  -Lx*Lz/Ld
c         ey2 = -Ly*Lz/Ld
c         ez2 =  Ld
c      EndIf
c      ez1 =0.
      Ld =sqrt(Lx**2 +Lz**2)
         ex1 =  Lz/Ld
         ey1 = 0.
         ez1 = -Lx/Ld
         ex2 = -Lx*abs(Ly)/Ld
         ey2 =  Ld
         ez2 = -Lz*abs(Ly)/Ld
c         If(Ly.lt.0.)ey2 = -Ld
      ex3 =Lx
      ey3 =Ly
      ez3 =Lz
c                     Check if the vectors are (1) unit (2) orthogonal
      prod1 = ex1*Lx +ey1*Ly +ez1*Lz
      prod2 = ex2*Lx +ey2*Ly +ez2*Lz
      prod3 = ex1*ex2 +ey1*ey2 +ez1*ez2
      a1       = Lx**2 +Ly**2 +Lz**2
      a2       = ex1**2 +ey1**2 +ez1**2
      a3       = ex2**2 +ey2**2 +ez2**2
      write (*,*) ' Vector 1=',ex1,ey1,ez1 
      write (*,*) ' Vector 2=',ex2,ey2,ez2 
      write (*,*) ' Vector 3=',ex3,ey3,ez3 
      write (*,*) ' Vector lengths=', a1,a2,a3,' must be 1'
      write (*,*) ' Vector products=', prod1,prod2,prod3,' must be 0'
      if(abs(a1-1.).gt.1.e-2)stop
      if(abs(a2-1.).gt.1.e-2)stop
      if(abs(a2-1.).gt.1.e-2)stop
      if(abs(prod1).gt.1.5e-1)stop
      if(abs(prod2).gt.1.5e-1)stop
      if(abs(prod3).gt.1.5e-1)stop
      xmax =-1.d+9
      xmin = 1.d+9
      ymax =-1.d+9
      ymin = 1.d+9
      zmax =-1.d+9
      zmin = 1.d+9
      dR = 0.
      dR2 = 0.
      V2 = 0.
      Lx = 0.
      Ly = 0.
      Lz = 0.
      Do jp=1,N
         Xnew = ex1*(Xp(jp)-Xcc) +ey1*(Yp(jp)-Ycc) +ez1*(Zp(jp)-Zcc)+Xcc
         Ynew = ex2*(Xp(jp)-Xcc) +ey2*(Yp(jp)-Ycc) +ez2*(Zp(jp)-Zcc)+Ycc
         Znew = ex3*(Xp(jp)-Xcc) +ey3*(Yp(jp)-Ycc) +ez3*(Zp(jp)-Zcc)+Zcc
         VXnew = ex1*Vxp(jp) +ey1*Vyp(jp) +ez1*Vzp(jp)
         VYnew = ex2*Vxp(jp) +ey2*Vyp(jp) +ez2*Vzp(jp)
         VZnew = ex3*Vxp(jp) +ey3*Vyp(jp) +ez3*Vzp(jp)
         Xp(jp) = Xnew
         Yp(jp) = Ynew
         Zp(jp) = Znew
                  IF(Xp(jp).lt.0.)Xp(jp)=Xp(jp) +Box 
                  IF(Yp(jp).lt.0.)Yp(jp)=Yp(jp) +Box 
                  IF(Zp(jp).lt.0.)Zp(jp)=Zp(jp) +Box
                  IF(Xp(jp).GE.Box)Xp(jp)=Xp(jp) -Box 
                  IF(Yp(jp).GE.Box)Yp(jp)=Yp(jp) -Box 
                  IF(Zp(jp).GE.Box)Zp(jp)=Zp(jp) -Box
         Vxp(jp) = VXnew 
         Vyp(jp) = VYnew 
         Vzp(jp) = VZnew 
                     xmax =MAX(xmax,Xp(jp))
                     xmin  =MIN(xmin,Xp(jp))
                     ymax =MAX(ymax,Yp(jp))
                     ymin  =MIN(ymin,Yp(jp))
                     zmax =MAX(zmax,Zp(jp))
                     zmin  =MIN(zmin,Zp(jp))
                     If(iMask(jp).gt.0)Then
                         Lx =Lx +(Yp(jp)-Ycc)*VZp(jp)
     &                             -(Zp(jp)-Zcc)*VYp(jp)
                         Ly =Ly +(Zp(jp)-Zcc)*VXp(jp)
     &                             -(Xp(jp)-Xcc)*VZp(jp)
                         Lz =Lz +(Xp(jp)-Xcc)*VYp(jp)
     &                             -(Yp(jp)-Ycc)*VXp(jp)
                        dd = (Xp(jp)-Xcc)**2 +
     &                        (Yp(jp)-Ycc)**2 + (Zp(jp)-Zcc)**2
                        dR = dR + sqrt(dd)
                        dR2= dR2+ dd
                     EndIf

         V2 = V2 +  VXp(jp)**2 +VYp(jp)**2+VZp(jp)**2         
      EndDo
      Lx =Lx/n_cooled
      Ly =Ly/n_cooled
      Lz =Lz/n_cooled
      V2 =sqrt(V2/N)
      aL = sqrt(Lx**2 +Ly**2 +Lz**2)
      dR = dR/n_cooled
      write (*,*) '--- After rotation of the coordinates'
      write (*,110) xmin,xmax,ymin,ymax,zmin,zmax,
     &           1000.*dR,1000.*sqrt(dR2/n_cooled),V2,
     &           Lx/aL,Ly/aL,Lz/aL,aL/dR,aL*1000.*n_cooled
      Return
      End
C------------------------------------------
C                    Distribution of the angular momentum
C                     of disk and halo particles
      SUBROUTINE AngMom(N,N1,R1,Box,ScaleM)
C------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      PARAMETER ( Nang = 50, Ntang = 2*Nang+1 )
      PARAMETER ( aMax = 4000. )
      Real*8 Lzd(-Nang:Nang),Lzd1(-Nang:Nang),Lzd2(-Nang:Nang)
      Integer Mzd(-Nang:Nang),Mzd1(-Nang:Nang),Mzd2(-Nang:Nang)
      Real*8 Lzh(-Nang:Nang),Lzh1(-Nang:Nang),Lzh2(-Nang:Nang)
      Integer Mzh(-Nang:Nang),Mzh1(-Nang:Nang),Mzh2(-Nang:Nang)

      DATA Lzd/Ntang*0.d+0/,Lzd1/Ntang*0.d+0/,Lzd2/Ntang*0.d+0/
      DATA Lzh/Ntang*0.d+0/,Lzh1/Ntang*0.d+0/,Lzh2/Ntang*0.d+0/
      DATA Mzd/Ntang*0/,Mzd1/Ntang*0/,Mzd2/Ntang*0/
      DATA Mzh/Ntang*0/,Mzh1/Ntang*0/,Mzh2/Ntang*0/

      Xcc = Box/2.
      Ycc = Box/2.
      Zcc = Box/2.
      hz = aMax/Nang
      Do jp=1,N
         X = (Xp(jp)- Xcc) *1000.  ! scale to kpc inits
         Y = (Yp(jp)- Ycc) *1000.
         Z = (Zp(jp)- Zcc) *1000.
            aLx =Y*VZp(jp)  -Z*VYp(jp)
            aLy =Z*VXp(jp)   -X*VZp(jp)
            aLz =X*VYp(jp)   -Y*VXp(jp)
            dd = X**2 + Y**2 + Z**2
            ds = sqrt(dd)
            ind = INT(aLz/hz+10000.) -10000
            ind = Min(Max(ind,-Nang),Nang)
            If(iMask(jp).gt.0)Then
                     Lzd(ind) =Lzd(ind) +aLz
                     Mzd(ind) =Mzd(ind) + 1
                     If(ds.le.R1)Then
                        Lzd2(ind) =Lzd2(ind) +aLz
                        Mzd2(ind) =Mzd2(ind) + 1
                     EndIf
                     If(jp.le.N1)Then
                        Lzd1(ind) =Lzd1(ind) +aLz
                        Mzd1(ind) =Mzd1(ind) + 1
                     EndIf
            Else
                     W = iWeight(jp)
                     Lzh(ind) =Lzh(ind) +aLz*W
                     Mzh(ind) =Mzh(ind) + W
                     If(ds.le.R1)Then
                        Lzh2(ind) =Lzh2(ind) +aLz*W
                        Mzh2(ind) =Mzh2(ind) + W
                     EndIf

            EndIf
      EndDo

      write (*,100)N1,R1,R1
      write (37,100)N1,R1,R1
 100  format(5x,'Disk',12x,'n<n1=',i5,5x,'r<r1=',f5.1,
     & 7x,'Halo',12x,'n<n1=',i5,5x,'r<r1=',f5.1,/
     &   ' Lz(kpckm/s)','<Lz> Npart',  
     &   '  <Lz>  Npart',
     &   '  <Lz>    Npart',
     &'    <Lz>    Npart',  
     &   '  <Lz>    Npart',
     &   '  <Lz>  Npart')
      Do i=-Nang,Nang
         a = i*hz
         aL = Lzd(i)/Max(Mzd(i),1)
         aL2= Lzd2(i)/Max(Mzd2(i),1)
         aL1= Lzd1(i)/Max(Mzd1(i),1)
         bL = Lzh(i)/Max(Mzh(i),1)
         bL2= Lzh2(i)/Max(Mzh2(i),1)
         write(*,110) a,aL,Mzd(i),aL1,Mzd1(i),aL2,Mzd2(i),
     &                        bL,Mzh(i),bL2,Mzh2(i)
         write(37,110) a,aL,Mzd(i),aL1,Mzd1(i),aL2,Mzd2(i),
     &                        bL,Mzh(i),bL2,Mzh2(i)
      EndDo
 110  format(f7.1,3(f7.1,i6),5x,f8.1,i8,2(f8.1,i6))
      Return
      End
C------------------------------------------
C                                                      
      SUBROUTINE Projection(N,Box,Alpha)
c             Alpha is angle between plain of the disk and
c             plain of projection
c  Zold                                                                                Znew
c   `                      / Disk (Xold)                                            |
c        `                /                                                                |
c             `          /  Alpha                                                     |
c==>>> ------------------------------ Xnew    ==>    -------------------------------Xnew
c                      /                                                                   |
c                     /                                                                    |
c
c           Znew is the distance in the plain of projects
c           Vlos = V_x,new 
c            Average for two observers 90 degrees away
C------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      Parameter (pi180 = 180./3.1415926535)
      Parameter (pi = 3.1415926535)
      EQUIVALENCE (Ndisk,extras(90))

      drlog =0.05
      d_lim = 0.100
      cosa = COS(Alpha/pi180)
      sina = SIN(Alpha/pi180)
      write (*,*) ' Angle =',Alpha, ' cos/sin=',cosa,sina 
      Xcc = Box/2.
      Ycc = Box/2.
      Zcc = Box/2.
      Do i =-Nrad,Nrad
               Nradp(i) = 0
               Radp(i)  = 0.
               Vradp(i)  = 0.
               Vrad2p(i) = 0.
      EndDo
      
      Do jp=1,Ndisk
         X = (Xp(jp)- Xcc) *1000.  ! scale to kpc inits
         Y = (Yp(jp)- Ycc) *1000.
         Z = (Zp(jp)- Zcc) *1000.
c                        x-z coordinates
            z_new = X*sina +Z*cosa ! z-coord in projection
            If(abs(z_new).le.d_lim)Then     ! inside the strip
             ind = INT(LOG10(ABS(Y))/drlog+1000.) -1000
             ind = min(max(ind,-Nrad),Nrad)
             Vradial = -VXp(jp)*cosa +VZp(jp)*sina
             If(Y.le.0.) Vradial =-Vradial
               Nradp(ind) = Nradp(ind) + 1
               Radp(ind)  = Radp(ind)  + abs(Y)
               Vradp(ind)  = Vradp(ind)  + Vradial
               Vrad2p(ind) = Vrad2p(ind)  + Vradial**2
            EndIf
c                        y-z coordinayes
            z_new = Y*sina +Z*cosa ! z-coord in projection
            If(abs(z_new).le.d_lim)Then     ! inside the strip
             ind = INT(LOG10(ABS(X))/drlog+1000.) -1000
             ind = min(max(ind,-Nrad),Nrad)
             Vradial = VYp(jp)*cosa -VZp(jp)*sina
             If(X.le.0.) Vradial =-Vradial
               Nradp(ind) = Nradp(ind) + 1
               Radp(ind)  = Radp(ind)  + abs(X)
               Vradp(ind)  = Vradp(ind)  + Vradial
               Vrad2p(ind) = Vrad2p(ind)  + Vradial**2
            EndIf
      EndDo
      V1 =0
      write(37,90) Alpha
 90   format(1x,10('-'),'Projection along major axis',10('-'),/
     &    '  Rout(kpc)   Rdisk  ',
     &            '  Vrot   sig_v_rot   Npart Alpha=',f7.2)
      Do  i=-Nrad+1,Nrad-1
         rad_out = 10.**((i+1)*drlog)
         nn = Nradp(i)
         Din = Din + nn
         n0 = max(nn,1)
         Raddisk  = Radp(i)/n0
         If(nn.eq.0)Raddisk=rad_out 
         Vradial = Vradp(i)/n0
         Vradial2=sqrt(max(Vrad2p(i)/n0-Vradial**2,1.e-10))
         If(nn.gt.20)Then
            write(37,100)rad_out,Raddisk,
     &           Vradial/cosa,Vradial2/cosa,nn
         EndIf
      EndDo
 100  format(2g10.3,2g11.3,i7)

      Return
      End

C------------------------------------------
C                                                      
      SUBROUTINE Profile(N,Box,ScaleM)
C------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      Parameter (pi180 = 180./3.1415926535)
      Parameter (pi = 3.1415926535)

      drlog =0.075
      Xcc = Box/2.
      Ycc = Box/2.
      Zcc = Box/2.

c                             DM profile   --------------------- spherical shells
      Do jp=1,N
        If(iMask(jp).lt.0)Then  ! not cooled particle
         W = iWeight(jp)
         X = (Xp(jp)- Xcc) *1000.  ! scale to kpc inits
         Y = (Yp(jp)- Ycc) *1000.
         Z = (Zp(jp)- Zcc) *1000.
            dd        = X**2 + Y**2 + Z**2
            ds         = sqrt(dd)
            Vtot      = (VXp(jp)**2+VYp(jp)**2+VZp(jp)**2)
            Vradial = (VXp(jp)*X+VYp(jp)*Y+VZp(jp)*Z) /
     &                      max(ds,1.e-5)
             Vtan2    = Vtot -Vradial**2
             ind = INT(LOG10(ds)/drlog+1000.) -1000
             ind = min(max(ind,-Nrad),Nrad)
               Nradh(ind) = Nradh(ind) + W
               Radh(ind)  = Radh(ind)  + ds*W
               Vfih2(ind)  = Vfih2(ind)  + Vtan2*W
               Vradh(ind)  = Vradh(ind)  + Vradial*W
               Vradh2(ind) = Vradh2(ind)  + Vradial**2*W
        EndIf 
      EndDo 
      V1 = 0.
      Din = 0.
      write(37,90)
 90   format(1x,20('-'),' DM profile',/
     &    '  Rout     Rhalo     Vchalo ',
     &            ' Density  ',
     &            '   Vrad    sigRadial      V_tan     Beta     Nhalo')
      Do  i=-Nrad,Nrad
         rad_out  = 10.**((i+1)*drlog)
         V0          = V1
         V1          = 4.1888*rad_out**3
         wn          = Nradh(i)
         nn           =wn /iWeight(1)
         Din         = Din + wn
         w0          = max(wn,1.e-9)
         densD    = wn/(V1-V0)*ScaleM*1.e-9
         Radhalo = Radh(i)/w0
         If(wn.eq.0)Radhalo=rad_out 
         Vtan       = sqrt(Vfih2(i)/w0)
         Vradial   = Vradh(i)/w0
         Vradial2 =sqrt(max(Vradh2(i)/w0-Vradial**2,1.e-10))
         Vtot        = sqrt((Vradh2(i)+Vfih2(i))/w0)
         Vcirc = 2.08e-3*sqrt(max(Din*ScaleM/rad_out,1.e-10))
         beta  = 1. -Vtan**2/2./(Vradial2**2+Vradial**2)
         If(nn.gt.4)Then
         write(37,100)rad_out,Radhalo,Vcirc,densD,
     &           Vradial,Vradial2,Vtan,beta,nn,INT(Din/iWeight(1))
         EndIf
      EndDo
 100  format(2g10.3,f7.2,
     &        g11.3,1x,4g11.4,2x,i7,i9)
                            ! disk statistics: spherical and cylindrical
c      Do i=-Nrad,Nrad
c         Nradp(i) =0
c      EndDo
      i_halo = 1   ! only one halo cooled
      n_cooled = iph_cool(i_halo,1)
      Do i=1,n_cooled
         jp = ip_cool(i)
         X = (Xp(jp)- Xcc) *1000.  ! scale to kpc inits
         Y = (Yp(jp)- Ycc) *1000.
         Z = (Zp(jp)- Zcc) *1000.
            aLx =Y*VZp(jp)  -Z*VYp(jp)
            aLy =Z*VXp(jp)   -X*VZp(jp)
            aLz =X*VYp(jp)   -Y*VXp(jp)
            dd = X**2 + Y**2 + Z**2
            ds = sqrt(dd)
            rd = X**2 + Y**2 
            rs = sqrt(rd)
            Vrotat  = aLz/max(rs,1.e-5)
            Vradial = (VXp(jp)*X+VYp(jp)*Y+VZp(jp)*Z) /
     &                   max(ds,1.e-5)
             Vvert  = VZp(jp)
             ind = INT(LOG10(ds)/drlog+1000.) -1000
             ind = min(max(ind,-Nrad),Nrad)
                                ! sperical shells
               Nradd(ind) = Nradd(ind) + 1
               Radd(ind)  = Radd(ind)  + ds
               Vfi(ind)   = Vfi(ind)  + Vrotat
c               Vfi2(ind)  = Vfi2(ind)  + Vrotat**2
c               Vrad(ind)  = Vrad(ind)  + Vradial
c               Vrad2(ind) = Vrad2(ind)  + Vradial**2
c               Vzz(ind)   = Vzz(ind)  + Vvert
c               Vzz2(ind)  = Vzz2(ind)  + Vvert**2

                                ! projection on the plane
             ind = INT(LOG10(rs)/drlog+1000.) -1000
             ind = min(max(ind,-Nrad),Nrad)
               Nradp(ind) = Nradp(ind) + 1
               Radp(ind)  = Radp(ind)  + rs
               Vfip(ind)   = Vfip(ind)  + Vrotat
               Vfi2p(ind)  = Vfi2p(ind)  + Vrotat**2
               Vradp(ind)  = Vradp(ind)  + Vradial
               Vrad2p(ind) = Vrad2p(ind)  + Vradial**2
               Vzzp(ind)   = Vzzp(ind)  + Vvert
               Vzz2p(ind)  = Vzz2p(ind)  + Vvert**2
                               ! projection along wedges
               inf = -1
               If(abs(Y)/rs.lt.0.5) inf =1  ! along the bar
               If(abs(X)/rs.lt.0.5) inf =2  ! perpendic  bar
               IF(inf.gt.0)Then
                  Nradb(ind,inf) = Nradb(ind,inf) + 1
                  Radb(ind,inf)  = Radb(ind,inf)  + rs
                  Vfib(ind,inf)   = Vfib(ind,inf)  + Vrotat
                  Vfi2b(ind,inf)  = Vfi2b(ind,inf)  + Vrotat**2
                  Vradb(ind,inf)  = Vradb(ind,inf)  + Vradial
                  Vrad2b(ind,inf) = Vrad2b(ind,inf)  + Vradial**2
                  Vzzb(ind,inf)   = Vzzb(ind,inf)  + Vvert
                  Vzz2b(ind,inf)  = Vzz2b(ind,inf)  + Vvert**2
               EndIf
      EndDo
                       ! print statistics ------------------------------
      V1 = 0.
      Din = 0.
      write(39,92)
 92   format(1x,20('-'),'All components: Spherical shells',/
     &    '  Rout     Rdisk      Vchalo Vcdisk Vcbulge Vctot',
     &    '  Vrot_disk',6x,'Nhalo',11x,'Ndisk',14x,'Nbulge')
      Do  i=-Nrad,Nrad
         rad_out  = 10.**((i+1)*drlog)
         V0          = V1
         V1          = 4.1888*rad_out**3
         wn          = Nradh(i)
         nn           =wn /iWeight(1)
         Din         = Din + wn
         w0          = max(wn,1.e-9)
         Radhalo     = Radh(i)/w0
         If(wn.eq.0)Radhalo=rad_out 
         nd          = Nradd(i)
         Dind        = Dind + nd
         n0          = max(nd,1)
         Raddisk = Radd(i)/n0
         Vrotat =Vfi(i)/n0
         If(nd.eq.0)Raddisk=rad_out 
         nbl         = Nradg(i)
         Dinbl       = Dinbl + nbl
         n0          = max(nbl,1)
         Radbulge    = Radg(i)/n0
         If(nbl.eq.0)Radbulge=rad_out 

         Vcirc = 2.08e-3*sqrt(max((Din+Dind+Dinbl)*ScaleM/rad_out,
     &                                             1.e-10))
         VcD  =  2.08e-3*sqrt(max((    Dind)*ScaleM/rad_out,1.e-10))
         VcH  =  2.08e-3*sqrt(max((Din     )*ScaleM/rad_out,1.e-10))
         VcBl =  2.08e-3*sqrt(max((Dinbl   )*ScaleM/rad_out,1.e-10))

         If(nn.gt.10)Then
         write(39,105)rad_out,Raddisk,VcH,VcD,VcBl,Vcirc,
     &           Vrotat,INT(wn),INT(Din),
     &           nd,INT(Dind),nbl,INT(Dinbl)
         EndIf
      EndDo
 105  format(2g10.3,f8.2,4f7.2,i8,2x,i8,2x,i6,4(2x,i8))
                         ! cylindrical shells: print ---------------------
      V1 = 0.
      Din = 0.
      write(41,94)
 94   format(1x,20('-'),'Cylindrical shells profile: Disk and Bulge',/
     &    '  Rout     Rdisk     Vrot  sig_rotat sig_z sig_rad',
     &    ' SurfDisk SurfBulge   Bar: Vrot sig_rot  Surf',
     &    '  off bar: Vrot sig_rot  Surf')
      Do  i=-Nrad,Nrad
         rad_out  = 10.**((i+1)*drlog)
         V0          = V1
         V1          = 3.141592*rad_out**2
         nd          = Nradp(i)
         Dind       = Dind + nd
         n0          = max(nd,1)
         Raddisk = Radp(i)/n0
         If(nd.eq.0)Raddisk=rad_out 


         Vrotat =Vfip(i)/n0            ! all disk 
         V2p     =sqrt(max(Vfi2p(i)/n0-Vrotat**2,1.e-10))
         Vwz      =sqrt(Vzz2p(i)/n0)
         Vradial =sqrt(Vrad2p(i)/n0)
         Surf     = nd/(V1-V0)
         SurfBl   = Nradu(i)/(V1-V0)

         nb     = Nradb(i,1)           ! along bar +- 30 degrees
         n0b   = max(nb,1)
         VrotatB =Vfib(i,1)/n0b
         V2B     =sqrt(max(Vfi2b(i,1)/n0b-VrotatB**2,1.e-10))
         SurfB    = nb/(V1-V0)*3.     ! correction for 30degrees

         nb     = Nradb(i,2)             ! off bar +- 30 degrees
         n0b   = max(nb,1)
         VrotatN =Vfib(i,2)/n0b
         V2N     =sqrt(max(Vfi2b(i,2)/n0b-VrotatN**2,1.e-10))
         SurfN    = nb/(V1-V0)*3.   

         vvv =sqrt(Vfi2b(i,2)/n0b)
c         write (*,'(4i6,6g11.3)') nd,n0,Nradb(i,1),Nradb(i,2),
c     &    Radp(i),rad_out,VrotatB ,V2B,VrotatN ,V2N
         If(nd.gt.10)Then
         write(41,107)rad_out,Raddisk,
     &                      Vrotat ,V2p,Vwz,Vradial,Surf,SurfBl,
     &                      VrotatB,V2B,SurfB,
     &                      VrotatN,V2N,SurfN
         EndIf
      EndDo
 107  format(2g10.3,4f7.2,2g11.3,3x,2f7.2,g11.3,3x,2f7.2,g11.3)

      Return
      End

C------------------------------------------
C                      Read cooled particles, set Mask
      SUBROUTINE ReadCool
C------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      Real xmax(2),xmin(2),ymax(2),ymin(2),
     +          zmax(2),zmin(2),vmax(2),vmin(2)
      DATA xmax,ymax,zmax,vmax/8*0.e-0/
      DATA xmin,ymin,zmin,vmin/8*0.e+15/

      read(10) n_halo
      IF(n_halo .gt. nh_cool) THEN
        write(*,*) 'n_halo .gt. nh_cool', n_halo, nh_cool
        write(*,*) 'increase nh_cool in a_cooling.h'
        STOP
      ENDIF

      kkL = 0   ! last index of the previous halo
      DO k = 1, n_halo
        read(10) (iph_cool(k,i),i =1,3),(halo_prop(k,kk),kk = 1,nh_prop)
        write(*,*) 'halo # :',k, '  iph_cool(k,i) = ',
     +           (iph_cool(k,i),i =1,3)     
        IF(iph_cool(k,1) .gt. n_cool)  THEN
          write(*,*) 'iph_cool(k,1) .gt. nh_cool',iph_cool(k,1),n_cool
          write(*,*) 'increase n_cool in a_cooling.h'
          STOP
        ENDIF
        kkF = kkL + 1                 ! read First particle
        if (kkF .NE. iph_cool(k,2)) write(*,*) 'ERROR in kkF'
        kkL = kkF + iph_cool(k,1) - 1   ! read Last particle
        if (kkL .NE. iph_cool(k,3)) write(*,*) 'ERROR in kkL'
        read(10) (ip_cool(kk), kk = iph_cool(k,2),iph_cool(k,3))
      ENDDO
      close (10)

      i_halo = 1   ! only one halo cooled
      n_cooled = iph_cool(i_halo,1)
      write(*,*) 'n_cooled = ', n_cooled
      Do i=1,Np
         iMask(i) = -1
      EndDo 
      DO k = 1, n_cooled
        i = ip_cool(k)
        iMask(i) = k
       ENDDO

       Return
      End
C------------------------------------------
C                                                      Read particles
      SUBROUTINE ReadPnts(N,ScaleC,ScaleV,ScaleM)
C------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'

      Integer Xdist(0:25),Ydist(0:25),Zdist(0:25)
      Real*8 TotWeight

      Box     = ScaleC*NGRID
      Radius2 = Radius**2

      N       = 0                            ! Number of particles
      Nsm     = 0
      TotWeight =0.
      xmax =-1.e+9
      xmin = 1.e+9
      ymax =-1.e+9
      ymin = 1.e+9
      zmax =-1.e+9
      zmin = 1.e+9
      vmax =-1.e+9
      vmin = 1.e+9
c      Do i=0,25
c         Xdist(i) =0
c         Ydist(i) =0
c         Zdist(i) =0
c      Enddo 
      Icount   = 0
      Icount_p = Icount-1
      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
cc         N_particles =          N_particles/3 !!!!!

         Npages = (N_particles -1)/NPAGE +1
         N_in_last=N_particles -NPAGE*(Npages-1)
      EndIf
      write (*,*) ' Pages=',Npages,' Species=',Nspecies
c      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
c            write (*,*)' Read page=',IROW,' file=',ifile,' N=',N
            iL = NPAGE*(IROW-1)
	      CALL GETROW(IROW,ifile)
                  Icount_p=Icount
	      DO IN=1,In_page
	        X  =ScaleC* (XPAR(IN)-1.)
	        Y  =ScaleC* (YPAR(IN)-1.)
	        Z  =ScaleC* (ZPAR(IN)-1.)
C           write (*,'(i8,10g11.4)') IN+iL,
c     &              X,Y,Z,XPAR(IN),YPAR(IN),ZPAR(IN)
                Vxs=ScaleV* VX(IN)
                Vys=ScaleV* VY(IN)
                Vzs=ScaleV* VZ(IN)
                Icount =Icount +1
                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             N   = N +1
                  If(Ipart.le.lspecies(1))Nsm=Nsm +1
                  IF(N.gt.Np)Then    ! too many particles
                           write (*,*) ' Too many particles. Max=',Np
                           STOP
                  Endif 
                  Xp(N)  =X   
                  Yp(N)  =Y   
                  Zp(N)  =Z   
                  Vxp(N)  =Vxs   
                  Vyp(N)  =Vys   
                  Vzp(N)  =Vzs   
                  Ipart  =IN+iL                     ! current particle number
                  iWeight(N) = INT(W/wspecies(1)+0.1) !  particles weight   iWeight(Ipart) 
                  TotWeight = TotWeight + iWeight(N) 
                  If(iWeight(N).lt.1.001)Then
                     xmax =MAX(xmax,Xp(N))
                     xmin  =MIN(xmin,Xp(N))
                     ymax =MAX(ymax,Yp(N))
                     ymin  =MIN(ymin,Yp(N))
                     zmax =MAX(zmax,Zp(N))
                     zmin  =MIN(zmin,Zp(N))
c                     ii =Min(INT(X/Box*25.),25)
c                     jj =Min(INT(Y/Box*25.),25)
c                     kk =Min(INT(Z/Box*25.),25)
c                     Xdist(ii) =Xdist(ii) +1
c                     Ydist(jj) =Ydist(jj) +1
c                     Zdist(kk) =Zdist(kk) +1
                   EndIf
        ENDDO
        ENDDO 
        write (*,*) ' Number of particles selected:=',N,' Small=',Nsm
        write (*,*) ' Limits                              =',xmin,xmax       
        write (*,*) '                                     =',ymin,ymax
        write (*,*) '                                     =',zmin,zmax
        write (*,'(" Total weight         =",1P,g11.3)')TotWeight
        write (*,'(" Total Mass           =",1P,g11.3)')TotWeight*ScaleM

c        write (*,'(i3,2x,i9,2x,i9,2x,i9)')
c     &             (i,Xdist(i),Ydist(i),Zdist(i),i=0,25)
        Return
        End
c     -----------------------------
       Function Age(aexp)
c     -----------------------------
c                                  age for LCDM model 
       PARAMETER (Om0 =0.3)
       PARAMETER (tconst =9.31e+9)
       
       x = ((1.-Om0)/Om0)**0.333333*aexp
       Age = tconst /sqrt(1.-Om0)*(log(x**1.5+sqrt(1.+x**3)))

       Return
       End
C ______________________________ Assign initial values
C		             
C		                     
      BLOCK DATA
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMcool.h'

      DATA Nradd/Nrad2*0/,Vfi/Nrad2*0./,Vrad/Nrad2*0./,
     &          Vzz/Nrad2*0./,Vfi2/Nrad2*0./,Vrad2/Nrad2*0./,
     &           Vzz2/Nrad2*0./,Radd/Nrad2*0./
      DATA Nradb/Nbt*0/,Vfib/Nbt*0./,Vradb/Nbt*0./,
     &          Vzzb/Nbt*0./,Vfi2b/Nbt*0./,Vrad2b/Nbt*0./,
     &           Vzz2b/Nbt*0./,Radb/Nbt*0./
      DATA Nradp/Nrad2*0/,Vfip/Nrad2*0./,Vradp/Nrad2*0./,
     &          Vzzp/Nrad2*0./,Vfi2p/Nrad2*0./,Vrad2p/Nrad2*0./,
     &           Vzz2p/Nrad2*0./,Radp/Nrad2*0./
      DATA Nradu/Nrad2*0/
      DATA Nradh/Nrad2*0/,Vfih/Nrad2*0./,Vradh/Nrad2*0./,
     &         Vfih2/Nrad2*0./,Vradh2/Nrad2*0./,Radh/Nrad2*0./
      DATA Dens/Np*0./
      End
C---------------------------------------------------
      SUBROUTINE RDTAPE(File1,File2)
C---------------------------------------------------
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMcool.h'

      Character  File1*100,File2*100
  
C                                     Open file on a tape
      OPEN(UNIT=9,FILE=File1,
     +     FORM='UNFORMATTED', STATUS='UNKNOWN')
C                                 Read control information
C                                 and see whether it has proper
C                                 structure

      READ  (9,err=10,end=10) HEADER, 
     +                  AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +                  TINTG,EKIN,EKIN1,EKIN2,AU0,AEU0,
     +                  NROWC,NGRIDC,Nspecies,Nseed,Om0,Oml0,hubble,Wp5
     +                   ,Ocurv,extras
      Ocurv =0. 

      WRITE (*,100) HEADER, 
     +     AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW,
     +     EKIN,
     +     NROWC,NGRID,NRECL,Om0,Oml0,hubble, 
     +     Ocurv


      WRITE (16,100) HEADER,
     +                   AEXPN,AEXP0,AMPLT,ASTEP,ISTEP,PARTW, 
     +                  EKIN,
     +                  NROWC,NGRID,NRECL,Om0,Oml0,hubble,  
     +                  Ocurv
     
 100  FORMAT(1X,'Header=>',A45,/   
     +     1X,' A=',F8.5,' 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)


      IF(NROW.NE.NROWC) THEN
         WRITE (*,*)
     +            ' NROW in PARAMETER and in TAPE-FILE are different'
      ENDIF
      IF(NGRID.NE.NGRIDC) THEN
         WRITE (*,*)
     +           ' NGRID in PARAMETER and in TAPE-FILE are different:'
         write (*,*) ' Ngrid=',NGRID,' NgridC=',NGRIDC
      ENDIF
C                                         Open work file on disk
 10   NBYTE = NRECL*4
      NACCES= NBYTE!/4 !for athena               
      
      
c      write (*,*) ' Open direct access file'
         OPEN(UNIT=21,FILE=File2,ACCESS='DIRECT',
     +	               STATUS='UNKNOWN',RECL=NACCES)
 
      REWIND 9
      RETURN
      END
C--------------------------------------------
C                            Read current PAGE of particles  from disk
C                            NRECL - length of ROW block in words
C--------------------------------------------
      SUBROUTINE GETROW(IROW,Ifile)  
C--------------------------------------------
c      INCLUDE 'Statparam.h'
      INCLUDE 'PMparameters.h'
      INCLUDE 'PMcool.h'

         READ  (20+Ifile,REC=IROW) RECDAT
      RETURN
      END
C--------------------------------------------------------------
C                                           Make linker lists of particles in each cell
      SUBROUTINE GetDensity(Box)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      EQUIVALENCE (Ndisk,extras(90))

      Cell = Box/Nblink
      dStep = 15.e-6
      Do i=1,5
         RadDens(i) =dStep*i**2  
      EndDo
      write (*,*)  ' GetDens: Cell=',Cell,dStep
      CALL List(Box)
C$OMP PARALLEL DO DEFAULT(SHARED) 
C$OMP+PRIVATE (jp,Density)
      Do jp=1,Ndisk
         CALL Neib(Xp(jp),Yp(jp),Zp(jp),Density) 
         Dens(jp) = Density
      EndDo

      Return
      End
C--------------------------------------------------------------
C                                           Make linker lists of particles in each cell
      SUBROUTINE List(Box)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      EQUIVALENCE (Ndisk,extras(90)),
     &               (Nbulge,extras(97)),(aMassDisk,extras(93))
             Do i=1,Ndisk
                Lst(i)=-1
             EndDo
             Do k=Nmlink,Nblink
             Do j=Nmlink,Nblink
             Do i=Nmlink,Nblink
                Label(i,j,k)=0
             EndDo
             EndDo
             EndDo  
      Do jp=1,Ndisk
         i=INT(Xp(jp)/Cell)
         j=INT(Yp(jp)/Cell)
         k=INT(Zp(jp)/Cell)
         i=MIN(MAX(Nmlink,i),Nblink)
         j=MIN(MAX(Nmlink,j),Nblink)
         k=MIN(MAX(Nmlink,k),Nblink)
         Lst(jp)      =Label(i,j,k)
         Label(i,j,k) =jp
      EndDo

      Return
      End
C--------------------------------------------------------------
C                                Find all neibhours for a center
C                                Xc,Yc,Zc - center; a0 -its weight
C
      SUBROUTINE Neib(Xc,Yc,Zc,Density)
C--------------------------------------------------------------
      INCLUDE 'PMparameters.h'      
      INCLUDE 'PMcool.h'
      DIMENSION iDen(5)
C                                       Initiate counters
        Radius2= RadDens(5)**2
        Radius = RadDens(5)
        Do i=1,5
           iDen(i) =0
        EndDo 
c        write (*,*)  '==== Neib:',Xc,Yc,Zc,dStep
c                      limits for Label
         i1=INT((Xc-Radius)/Cell)
         j1=INT((Yc-Radius)/Cell)
         k1=INT((Zc-Radius)/Cell)
            i1=MIN(MAX(Nmlink,i1),Nblink)
            j1=MIN(MAX(Nmlink,j1),Nblink)
            k1=MIN(MAX(Nmlink,k1),Nblink)
         i2=INT((Xc+Radius)/Cell)
         j2=INT((Yc+Radius)/Cell)
         k2=INT((Zc+Radius)/Cell)
            i2=MIN(MAX(Nmlink,i2),Nblink)
            j2=MIN(MAX(Nmlink,j2),Nblink)
            k2=MIN(MAX(Nmlink,k2),Nblink)
c            write (*,'("    limits=",6i5)')i1,i2,j1,j2,k1,k2
C                                        Look for neibhours
         Do k=k1,k2
         Do j=j1,j2
         Do i=i1,i2
            jp =Label(i,j,k)
 10         If(jp.ne.0)Then
               dd =(Xc-Xp(jp))**2+(Yc-Yp(jp))**2+(Zc-Zp(jp))**2
               If(dd.lt.Radius2)Then
                  disd        = sqrt(sqrt(dd)/dStep)
                  ind       = min(INT(disd) +1,5)
                  iDen(ind) = iDen(ind) +1
c                  write (*,'(6x,i9,g11.4,3f9.5,2i5)')
c     &             jp,sqrt(dd)*1.e6,
c     &              (Xp(jp)-0.5)*1.e3,(Yp(jp)-0.5)*1.e3,
c     &              (Zp(jp)-0.5)*1.e3,
c     &              ind,iDen(ind) 
              EndIf  
                  jp =Lst(jp)
               GoTo10
            EndIf
         EndDo
         EndDo
         EndDo

         nn =0
         Do i=1,5
            nn =nn +iDen(i)
            If(nn.ge.5)Then
               Density = nn/RadDens(i)**3
               goto 20
            EndIf
         EndDo
         Density = nn/RadDens(5)**3
 20   Density =Density/4.188e9  ! number of particles per kpc**3
c      write(*,'(5x,"Dens=",g11.3,6i8)')Density,iDen,nn
      Return
      End
