c  program for planar Couette flow
c     
      include "ParC.inc"
      call init
      do i=1,1000
        a=rf()
      end do
      jj=0
      time=0.
      do while (time<200.)
        time=time+dtm
        jj=jj+1
        if(jj/1000==jj/1000.)print*,time,Pxy
        call move 
        call inde
        call coll
        call sample
        if(time>100.and.nwt==0) call init1
      end do
      call out
      stop                                                              
      end 
c
      subroutine init               
      include "ParC.inc"
      T0=300.d0 ! equlibrium tem in K
      au=1.66053886d-7 ! atomic unit times 1.E20
      amu=22.669d-6 ! visc. of Ar at T0 
      amm=39.948 ! molecular mass of Ar
      dd=3.35741 ! distance in angstrom  when U=0, Ar
      vm=sqrt(2.d0*8314.4621*T0/amm) ! most prob. speed
      Tr=T0/143.123d0 ! reduced temperature for Ar
      pi=acos(-1.d0) ! pi number
      ANc=pi*(dd*bm)**2*amu*del*dtm*nc/(amm*au*vm*inm)
      nsmp=1                                                           
      do 101 m=1,nc  
        ccg(1,m)=4.d0 !the max relative speeds.
 101    ccg(2,m)=rf() ! remainder   
      do nm=1,inm  ! initial distribution
        pp(nm)=rf()-0.5   ! generation of coordinate
  103   vx=8.*rf()-4. ! generation of velocity
        if(exp(-vx*vx)<rf()) go to 103                             
        pv(1,nm)=vx                                                   
        a=sqrt(-log(rf()))                                            
        b=2.d0*pi*rf()                                           
        pv(2,nm)=a*sin(b)                                             
        pv(3,nm)=a*cos(b)                                             
      end do          
      amx=0.d0
      cs=0.d0   
      nwt=0
      open(10,file='chiAr.dat') ! reading of deflection angle
      do 102 j=1,ne
 102    read(10,*) (xchi(i,j),i=1,nb)
      return                                                            
      end                                                               
c                                                                       
      subroutine init1                                                  
      include "ParC.inc"                                             
      time=0.                                                           
      nsmp=1
      amx=0.d0
      cs=0.d0    
      nwt=1                                                             
      return                                                            
      end                                                               
c                                                                       
      subroutine move  !molecules are moved over time interval dtm
      include "ParC.inc"                                           
      do n=1,inm    
        xi=pp(n) ! old coordinate
        x=xi+pv(1,n)*dtm ! new coordinate
  100   if(x > 0.5) then ! molecule stricks plane x=1
          atr=(x-0.5)/pv(1,n) ! rest of dtm
c  calculation of velocity compnents after reflection
          pv(1,n)=-sqrt(-log(rf()))                          
          a=sqrt(-log(rf())) 
          b=2.d0*pi*rf()
          pv(2,n)=a*sin(b)-Uw/2.
          pv(3,n)=a*cos(b)
          x=0.5+pv(1,n)*atr
        end if !if 100
        if(x < -0.5) then !molecule stricks plane x=0
          atr=(x+0.5)/pv(1,n) ! rest of dtm
          pv(1,n)=sqrt(-log(rf()))                          
          a=sqrt(-log(rf())) 
          b=2.*pi*rf()                                          
          pv(2,n)=a*sin(b)+Uw/2.
          pv(3,n)=a*cos(b)
          x=-0.5+pv(1,n)*atr                                               
        end if
        if(x>0.5) go to 100        
c calculation of energy fluxes
        if(xi <0..and.x > 0.) amx=amx+pv(2,n)
        if(xi >0..and.x < 0.) amx=amx-pv(2,n)
        ipl(n)=(x+0.5d0)*Nc+1 ! new number of cell                              
        pp(n)=x                                                       
      end do          
      return                                                            
      end                                                               
c                                                                       
      subroutine inde                                                   
c  the nm molecule numbers are arranged in order of cells
      include "ParC.inc"                                      
      ic(2,:)=0                                                  
      do 101 n=1,inm ! number of molecules counted                                                       
        mc=ipl(n)
 101    ic(2,mc)=ic(2,mc)+1                                         
      m=0        
      do 102 n=1,nc            
        ic(1,n)=m                                                   
        m=m+ic(2,n)                                                 
 102    ic(2,n)=0 ! start address has been set for the cells                                                   
      do 103 n=1,inm
        mc=ipl(n)  
        ic(2,mc)=ic(2,mc)+1
        k=ic(1,mc)+ic(2,mc)
 103    ir(k)=n
      return                                                            
      end   
c                                  
      subroutine coll !calculates collisions  
      include "ParC.inc"   
      dimension vccm(3),vrc(3),vrcp(3)
      do n=1,nc 
        asel=ANc*ccg(1,n)*ic(2,n)*(ic(2,n)-1)+ccg(2,n) 
        nsel=asel !number of pairs to be selected.          
        ccg(2,n)=asel-nsel                                     
        do isel=1,nsel
          k=int(rf()*(ic(2,n)-0.001))+ic(1,n)+1                
          l=ir(k) !first molecule has been chosen 
  100     k=int(rf()*(ic(2,n)-0.001))+ic(1,n)+1                
          m=ir(k) !second molecule has been chosen                                                 
         if(l==m) go to 100 ! if the first is again chosen            
         vrc(:)=pv(:,l)-pv(:,m) ! relative velocity
         vr=sqrt(vrc(1)**2+vrc(2)**2+vrc(3)**2) ! relative spped 
         if (vr>ccg(1,n)) ccg(1,n)=vr !the maximum  is upgraded 
         if (rf()<vr/ccg(1,n)) then
          vccm(:)=(pv(:,l)+pv(:,m))*0.5d0 ! center mass velocity
          vrr=sqrt(vrc(2)*vrc(2)+vrc(3)*vrc(3))  
          eps=2.d0*pi*rf()
          Er=0.5d0*Tr*vr*vr ! energy of relative motion
          ib=Nb*rf()+1
          je=Er*Ne/Emax+1
          if(je>Ne) je=Ne
          x=xchi(ib,je) ! deflection angle is chosen
          if (vrr>1.d-15) then
            vrcp(1)=vrc(1)*cos(x)+sin(x)*sin(eps)*vrr
            vrcp(2)=vrc(2)*cos(x)+sin(x)*
     &              (vr*vrc(3)*cos(eps)-vrc(1)*vrc(2)*sin(eps))/vrr
            vrcp(3)=vrc(3)*cos(x)-sin(x)*
     &              (vr*vrc(2)*cos(eps)+vrc(1)*vrc(3)*sin(eps))/vrr
          else
            vrcp(1)=vrc(1)*cos(x)
            vrcp(2)=vrc(1)*sin(x)*cos(eps)
            vrcp(3)=vrc(1)*sin(x)*sin(eps)
          end if   
          pv(:,l)=vccm(:)+0.5d0*vrcp(:)       
          pv(:,m)=vccm(:)-0.5d0*vrcp(:)
         end if
        end do
      end do 
      return                                                            
      end                                                               
c                                                                       
      subroutine sample 
      include "ParC.inc"                                              
      nsmp=nsmp+1 ! number of samples 
      do 100 n=1,nc                                                   
        l=ic(2,n)                                                     
        do 100 j=1,l                                                      
          k=ic(1,n)+j                               
          m=ir(k)                                                     
          cs(1,n)=cs(1,n)+1
          cs(2,n)=cs(2,n)+pv(2,m)
 100      cs(3,n)=cs(3,n)+pv(1,m)**2+pv(2,m)**2+pv(3,m)**2
      an(:)=cs(1,:)*Nc/float(nsmp*inm) ! density                    
      u(:)=cs(2,:)/cs(1,:)
      temp(:)=2./3.*(cs(3,:)/cs(1,:)-u(:)**2) ! temperature
      Pxy=2.*amx/(nsmp*dtm*inm*Uw) ! heat flux       
      return                                                            
      end                                                               
c                                                                       
      subroutine out                                                    
      include "ParC.inc"                                           
      open (10,file="Res_Couette.dat")
      write(10,901) Pxy
      do 100 m=1,nc/2 ! writing of results
100     write(10,900) (m-0.5)/nc-0.5, an(m),u(m),temp(m)
901   format (1x,'Pxy=',f10.5/5x,'x',8x,'den',8x,'u',10x,'Temp')
900   format (f9.4,3(2x,f9.5))
      return                                                            
      end   
c                                                                      
      double precision function rf()                                    
      implicit double precision (a-h, o-z)                              
      data ix1,ix2,ix3,ix4,ix5,ix6/1500419,1400159,1364,1528,1,3/       
      rr1=1.0/float(ix1)                                                
      rr2=1.0/float(ix2)                                                
      ix5=mod(ix5*ix3,ix1)                                              
      ix6=mod(ix6*ix4,ix2)                                              
      rf=rr1*ix5+rr2*ix6                                                
      if(rf.ge.1.0)rf=rf-1.0                                            
      return                                                            
      end