c!*******************************************************************************************

c!*     calculate order parameter

c!*******************************************************************************************

      INCLUDE 'FGRAPH.FI'      

      program order()

      implicit none      

      integer i,j,k,m,ii,jj,kk,ix,iy,iz,Nx,Ny,Nz,ixx,iyy,row,remain

     +,izz,ccc,my,mz

      parameter(Nx=64,Ny=Nx/2,Nz=64)

      integer ab(1:2,1:Nx,1:Ny,1:Nz),vx1(1:4),vy1(1:4),vz1(1:4)

      integer vx2(1:4),vy2(1:4),vz2(1:4),vx3(1:4),vy3(1:4),vz3(1:4)

      integer Nk(1:2,1:Nx,1:Ny,1:Nz),Sk(1:2,1:Nx,1:Ny,1:Nz)

      integer x0(1:2,1:2,1:4),y0(1:2,1:2,1:4),z0(1:2,1:2,1:4)

      real*8  cc(1:2,1:Nx,1:Ny,1:Nz),nn(1:2,1:Nx,1:Ny,1:Nz)

     +,fa,fb,f,ra,lro

      character*4 file_a,file_b 

      

      common /dd1/ nn

      common /dd2/ ab

      ccc=0

      ra=0.0 



      open(2,file='ini.dat')

	do i=1,2

	 do j=1,2

	  do k=1,4

	   read(2,'(3I3)') x0(i,j,k),y0(i,j,k),z0(i,j,k)

	  enddo

	 enddo

	enddo

	do i=1,6

       read(2,'(3I3)') vx1(i),vy1(i),vz1(i)

	enddo

	do i=1,4

       read(2,'(3I3)') vx2(i),vy2(i),vz2(i)

	enddo

	do i=1,4

       read(2,'(3I3)') vx3(i),vy3(i),vz3(i)

	enddo

	close(2)

      

	fa=0.0

	fb=0.0 

	row=int(Nx/6)

	remain=mod(Nx,6) 

	

	print*,'Input the name for alpha lattice:' 

	

	read*,file_a 

      print*,file_a

	

	open(2,file=file_a,mode='read')



	do iz=1,Nz

	do iy=1,Ny



	

	ix=0

	do i=1,row

	ccc=ccc+6

	ix=ix+6

	read(2,'(6(i12))') ab(1,ix-5,iy,iz),ab(1,ix-4,iy,iz)

     +,ab(1,ix-3,iy,iz),ab(1,ix-2,iy,iz),ab(1,ix-1,iy,iz),ab(1,ix,iy,iz)  





	if(ab(1,ix,iy,iz).eq.1) ab(1,ix,iy,iz)=1

      if(ab(1,ix,iy,iz).eq.2) ab(1,ix,iy,iz)=-1  

	if(ab(1,ix-1,iy,iz).eq.1) ab(1,ix-1,iy,iz)=1

      if(ab(1,ix-1,iy,iz).eq.2) ab(1,ix-1,iy,iz)=-1  

	if(ab(1,ix-2,iy,iz).eq.1) ab(1,ix-2,iy,iz)=1

      if(ab(1,ix-2,iy,iz).eq.2) ab(1,ix-2,iy,iz)=-1  

	if(ab(1,ix-3,iy,iz).eq.1) ab(1,ix-3,iy,iz)=1

      if(ab(1,ix-3,iy,iz).eq.2) ab(1,ix-3,iy,iz)=-1  

	if(ab(1,ix-4,iy,iz).eq.1) ab(1,ix-4,iy,iz)=1

      if(ab(1,ix-4,iy,iz).eq.2) ab(1,ix-4,iy,iz)=-1  

	if(ab(1,ix-5,iy,iz).eq.1) ab(1,ix-5,iy,iz)=1

      if(ab(1,ix-5,iy,iz).eq.2) ab(1,ix-5,iy,iz)=-1  

	

	if(ab(1,ix,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix,iy,iz).eq.1) fb=fb+1.0 

	if(ab(1,ix-1,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix-1,iy,iz).eq.1) fb=fb+1.0 

	if(ab(1,ix-2,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix-2,iy,iz).eq.1) fb=fb+1.0 

	if(ab(1,ix-3,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix-3,iy,iz).eq.1) fb=fb+1.0 

	if(ab(1,ix-4,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix-4,iy,iz).eq.1) fb=fb+1.0 

	if(ab(1,ix-5,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix-5,iy,iz).eq.1) fb=fb+1.0 





	

	enddo

	



	read(2,'(4(i12))') ab(1,ix+1,iy,iz),ab(1,ix+2,iy,iz)

     +,ab(1,ix+3,iy,iz),ab(1,ix+4,iy,iz)



      

	if(ab(1,ix+1,iy,iz).eq.1) ab(1,ix+1,iy,iz)=1

	if(ab(1,ix+1,iy,iz).eq.2) ab(1,ix+1,iy,iz)=-1

	if(ab(1,ix+2,iy,iz).eq.1) ab(1,ix+2,iy,iz)=1

	if(ab(1,ix+2,iy,iz).eq.2) ab(1,ix+2,iy,iz)=-1

	if(ab(1,ix+3,iy,iz).eq.1) ab(1,ix+3,iy,iz)=1

	if(ab(1,ix+3,iy,iz).eq.2) ab(1,ix+3,iy,iz)=-1

	if(ab(1,ix+4,iy,iz).eq.1) ab(1,ix+4,iy,iz)=1

	if(ab(1,ix+4,iy,iz).eq.2) ab(1,ix+4,iy,iz)=-1





	

	if(ab(1,ix+1,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix+1,iy,iz).eq.1) fb=fb+1.0

	if(ab(1,ix+2,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix+2,iy,iz).eq.1) fb=fb+1.0	

	if(ab(1,ix+3,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix+3,iy,iz).eq.1) fb=fb+1.0	

	if(ab(1,ix+4,iy,iz).eq.-1) fa=fa+1.0

	if(ab(1,ix+4,iy,iz).eq.1) fb=fb+1.0	

	 

	ccc=ccc+4





	enddo

	enddo



      close(2)

      

	print*,'Input the name for beta lattice:' 

	

	read*,file_b

      print*,file_b

	

	open(2,file=file_b,mode='read')

	do iz=1,Nz

	do iy=1,Ny

	



	ix=0

	do i=1,row 

	ccc=ccc+6

	ix=ix+6

	read(2,'(6(i12))') ab(2,ix-5,iy,iz),ab(2,ix-4,iy,iz)

     +,ab(2,ix-3,iy,iz),ab(2,ix-2,iy,iz),ab(2,ix-1,iy,iz),ab(2,ix,iy,iz)  



      

	if(ab(2,ix,iy,iz).eq.1) ab(2,ix,iy,iz)=1

      if(ab(2,ix,iy,iz).eq.2) ab(2,ix,iy,iz)=-1  

	if(ab(2,ix-1,iy,iz).eq.1) ab(2,ix-1,iy,iz)=1

      if(ab(2,ix-1,iy,iz).eq.2) ab(2,ix-1,iy,iz)=-1  

	if(ab(2,ix-2,iy,iz).eq.1) ab(2,ix-2,iy,iz)=1

      if(ab(2,ix-2,iy,iz).eq.2) ab(2,ix-2,iy,iz)=-1  

	if(ab(2,ix-3,iy,iz).eq.1) ab(2,ix-3,iy,iz)=1

      if(ab(2,ix-3,iy,iz).eq.2) ab(2,ix-3,iy,iz)=-1  

	if(ab(2,ix-4,iy,iz).eq.1) ab(2,ix-4,iy,iz)=1

      if(ab(2,ix-4,iy,iz).eq.2) ab(2,ix-4,iy,iz)=-1  

	if(ab(2,ix-5,iy,iz).eq.1) ab(2,ix-5,iy,iz)=1

      if(ab(2,ix-5,iy,iz).eq.2) ab(2,ix-5,iy,iz)=-1  

	

	if(ab(2,ix,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix,iy,iz).eq.1) fb=fb+1.0 

	if(ab(2,ix-1,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix-1,iy,iz).eq.1) fb=fb+1.0 

	if(ab(2,ix-2,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix-2,iy,iz).eq.1) fb=fb+1.0 

	if(ab(2,ix-3,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix-3,iy,iz).eq.1) fb=fb+1.0 

	if(ab(2,ix-4,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix-4,iy,iz).eq.1) fb=fb+1.0 

	if(ab(2,ix-5,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix-5,iy,iz).eq.1) fb=fb+1.0 





	

	enddo

	



	read(2,'(4(i12))') ab(2,ix+1,iy,iz),ab(2,ix+2,iy,iz)

     +,ab(2,ix+3,iy,iz),ab(2,ix+4,iy,iz)



      

	if(ab(2,ix+1,iy,iz).eq.1) ab(2,ix+1,iy,iz)=1

	if(ab(2,ix+1,iy,iz).eq.2) ab(2,ix+1,iy,iz)=-1

	if(ab(2,ix+2,iy,iz).eq.1) ab(2,ix+2,iy,iz)=1

	if(ab(2,ix+2,iy,iz).eq.2) ab(2,ix+2,iy,iz)=-1

	if(ab(2,ix+3,iy,iz).eq.1) ab(2,ix+3,iy,iz)=1

	if(ab(2,ix+3,iy,iz).eq.2) ab(2,ix+3,iy,iz)=-1

	if(ab(2,ix+4,iy,iz).eq.1) ab(2,ix+4,iy,iz)=1

	if(ab(2,ix+4,iy,iz).eq.2) ab(2,ix+4,iy,iz)=-1





	

	if(ab(2,ix+1,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix+1,iy,iz).eq.1) fb=fb+1.0

	if(ab(2,ix+2,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix+2,iy,iz).eq.1) fb=fb+1.0	

	if(ab(2,ix+3,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix+3,iy,iz).eq.1) fb=fb+1.0	

	if(ab(2,ix+4,iy,iz).eq.-1) fa=fa+1.0

	if(ab(2,ix+4,iy,iz).eq.1) fb=fb+1.0	



	ccc=ccc+4





	enddo

	enddo

      close(2)

      

      print*,'fa+fb:', fa+fb

      

      fa=fa/(fa+fb)

      fb=1.0-fa



	print*,'A%=',fa,'B%=',fb 

      

      print*,ccc

	ccc=0

      do i=1,2

       do j=1,2 

        do k=1,4

	    ix=x0(i,j,k)-4

         do m=-2,Nx/2

          iy=y0(i,j,k)-4

	    ix=ix+2

	   do my=-2,Ny/2

	    iz=z0(i,j,k)-4

          iy=iy+2

	   do mz=-2,Nz/2

          iz=iz+2

          if((ix.le.Nx.and.ix.ge.1).and.(iy.le.Ny.and.iy.ge.1)

     +.and.(iz.le.Nz.and.iz.ge.1)) then

	     Nk(i,ix,iy,iz)=0

	     Sk(i,ix,iy,iz)=-1

	     if(ab(i,ix,iy,iz).eq.1) Nk(i,ix,iy,iz)=1

           if(i.eq.2.and.ab(i,ix,iy,iz).eq.1) Sk(i,ix,iy,iz)=1

	     if(i.eq.1.and.ab(i,ix,iy,iz).eq.-1) then

	      Sk(i,ix,iy,iz)=1

            ra=ra+1.0

	     endif

	ccc=ccc+1

          nn(i,ix,iy,iz)=0.0d0

	    cc(i,ix,iy,iz)=0.0d0

          endif

         enddo

	   enddo

	   enddo

        enddo 

       enddo  

      enddo 

	print*,ccc

	ccc=0



c      open(2,file='nncc.dat')

      do i=1,2

       do j=1,2 

        do k=1,4

	   ix=x0(i,j,k)-2

         do m=-1,Nx/2

          iy=y0(i,j,k)-2

		ix=ix+2

	   do my=-1,Ny/2

	    iz=z0(i,j,k)-2

          iy=iy+2

	   do mz=-1,Nz/2

          iz=iz+2 

      if((ix.le.Nx.and.ix.ge.1).and.(iy.le.Ny.and.iy.ge.1)

     +.and.(iz.le.Nz.and.iz.ge.1)) then

           nn(i,ix,iy,iz)=(0.250d0*Sk(i,ix,iy,iz))+nn(i,ix,iy,iz)

	     cc(i,ix,iy,iz)=(0.250d0*Nk(i,ix,iy,iz))+cc(i,ix,iy,iz)

	     do ii=1, 6

	      ixx=mod(ix+vx1(ii)+Nx-1,Nx)+1

	      iyy=mod(iy+vy1(ii)+Ny-1,Ny)+1

	      izz=mod(iz+vz1(ii)+Nz-1,Nz)+1

       nn(i,ix,iy,iz)=nn(i,ix,iy,iz)+((1.0d0/24.0d0)*Sk(i,ixx,iyy,izz))

	 cc(i,ix,iy,iz)=cc(i,ix,iy,iz)+((1.0d0/24.0d0)*Nk(i,ixx,iyy,izz))

	     enddo

	     jj=mod(i,2)+1

	     kk=i-2

           do ii=1, 4

 	      ixx=mod(ix+vx2(ii)+kk+Nx-1,Nx)+1

	      iyy=mod(iy+vy2(ii)+kk+Ny-1,Ny)+1

	      izz=mod(iz+vz2(ii)+kk+Nz-1,Nz)+1

       nn(i,ix,iy,iz)=nn(i,ix,iy,iz)+((1.0d0/16.0)*Sk(jj,ixx,iyy,izz))

	 cc(i,ix,iy,iz)=cc(i,ix,iy,iz)+((1.0d0/16.0)*Nk(jj,ixx,iyy,izz))

	     enddo

	     do ii=1, 4

 	      ixx=mod(ix+vx3(ii)+kk+Nx-1,Nx)+1

	      iyy=mod(iy+vy3(ii)+kk+Ny-1,Ny)+1

	      izz=mod(iz+vz3(ii)+kk+Nz-1,Nz)+1

       nn(i,ix,iy,iz)=nn(i,ix,iy,iz)+((1.0d0/16.0)*Sk(jj,ixx,iyy,izz))

	 cc(i,ix,iy,iz)=cc(i,ix,iy,iz)+((1.0d0/16.0)*Nk(jj,ixx,iyy,izz))

	     enddo

c           write(2,'(4(i3,x),2(f40.20),x)') i,ix,iy,iz,nn(i,ix,iy,iz),cc(i,ix,iy,iz) 

          ccc=ccc+1

		endif

         enddo

	   enddo

	   enddo

        enddo 

       enddo  

      enddo

	print*,ccc

c	close(2)

      

      f=fa

      if(fa.gt.0.5) then

      ra=1.0*Nx*Ny*Nz-(2.0*Nx*Ny*Nz*fa-ra)

      f=fb

      endif

      

      ra=ra/(Nx*Ny*Nz)

	print*, 'less%site=',ra

	lro=abs((ra-f)/(1.0-f))

	print*, 'lro=',lro

      call c_sro(fa,fb)

	call sh_plane(1,-1,0,0,0,0)

	call c_sk(file_a)

	       

      end







c!***************************************************************************************

c!*     get the order parameter of points on plane (i,j,k)

c!***************************************************************************************

	subroutine sh_plane(i,j,k,x0,y0,z0)

	integer i,j,k,x0,y0,z0,ii,jj,kk,mm,Nx,Ny,Nz

	integer*2 count,L

      parameter(Nx=64,Ny=Nx/2,Nz=64)

      real*8  nn(1:2,1:Nx,1:Ny,1:Nz),ix,iy,iz,origin

	common /dd1/ nn  

	count=0  

	L=Nx



	open(2,file='plane.dat',form='formatted')

	do ii=1,2

	 origin=(ii-1.0)*0.5

	 ix=origin

       do jj=1,Nx

        iy=origin

	  ix=ix+1.0

        do kk=1, Ny

         iz=origin

	   iy=iy+1.0

         do mm=1, Nz

          iz=iz+1.0

          if((ix.le.(Nx+origin).and.ix.ge.0.0).and.(iy.le.(Ny+origin)

     +.and.iy.ge.0.0).and.(iz.le.(Nz+origin).and.iz.ge.0.0)) then

	     if((1.0*(ix-x0)*i+1.0*(iy-y0)*j+1.0*(iz-z0)*k)

     +  .eq.0.0) then

            write(2,'(3(f30.15))') sqrt((ix-x0)**2.0+

     +(iy-y0)**2.0),(iz-z0)*1.0, nn(ii,jj,kk,mm)

      count=count+1 

	     endif

	    endif

	   enddo

	  enddo

	 enddo

	enddo

	close(2)  

	read(*,*)

	 

	call slice(count,L)

	

	end







c!**********************************************************************************************

c!*     show the order parameter on one plane

c!**********************************************************************************************      

      subroutine slice(count,L)      

      INCLUDE 'FGRAPH.FD'

      integer*2 status2,i, count,n_c,color,L

      real*8  c ,x,y ,norm1,norm2                  

      parameter(n_c=8) 

      

      norm1=1.0/(sqrt(3.0)*L)*800.0 

      norm2=1.0/(L)*600.0



      status2 = SETVIDEOMODE( $VRES16COLOR ) 

      open(2,file='plane.dat',mode='read')

      

      

      do i=1, count



      read(2,'(3(f30.15))') x,y,c 

      

      x=x*norm1+5.0

      y=y*norm2+5.0    



      color=int((c+1.0)/2.0*16+1)

      status2 = SETCOLOR( color )

      status2 = ELLIPSE_W( $GFILLINTERIOR, y, x, y+4.0, x+4.0 )   

      

      enddo



      close(2)

      

      read(*,*)

      

      CALL CLEARSCREEN( $GCLEARSCREEN )

      status2 = SETVIDEOMODE( $DEFAULTMODE )

           

      END





c!**********************************************************************************************

c!*     calculate short range order parameter

c!**********************************************************************************************

	subroutine c_sro(fa,fb)

      integer i,Nx,Ny,Nz,ix,iy,iz,ii,iix,iiy,iiz,x,y,z

      real*8 fa,fb,f,ccc,grab,pt,sro

      parameter(Nx=64,Ny=Nx/2,Nz=64)

	integer ab(1:2,1:Nx,1:Ny,1:Nz)

	common /dd2/ ab

	

	f=fa

	if(fb.le.fa) f=fb

      pt=Nx*Ny*Nz*8.0



	ccc=0.0

	grab=0.0



	do i=1,2

       do ix=1,Nx

	  do iy=1,Ny

	   do iz=1,Nz

	  

      

	 if(i.eq.1) then

	

	ii=2

      do x=ix-1,ix

	iix=mod(x+Nx,Nx)+1

	do y=iy-1,iy

	iiy=mod(y+Ny,Ny)+1

	do z=iz-1,iz

	iiz=mod(z+Nz,Nz)+1



	if(ab(ii,iix,iiy,iiz).ne.ab(i,ix,iy,iz)) grab=grab+0.5

	ccc=ccc+1.0



	enddo  

      enddo

      enddo



       else

	

	ii=1

      do x=ix,ix+1

	iix=mod(x+Nx,Nx)+1

	do y=iy,iy+1

	iiy=mod(y+Ny,Ny)+1

	do z=iz,iz+1

	iiz=mod(z+Nz,Nz)+1



	if(ab(ii,iix,iiy,iiz).ne.ab(i,ix,iy,iz)) grab=grab+0.5

	ccc=ccc+1.0



	enddo  

      enddo

      enddo



	 endif





	   enddo

	  enddo

	 enddo

	enddo

	

	sro=abs((grab/(pt*f)-(1-f))/(1.0+f))

	print*, 'sro=',sro,grab,pt*f*(1-f)



	

	end





c!***************************************************************************************

c!*     get s of k

c!***************************************************************************************

      subroutine c_sk(file_a)

	integer i,j,k,ii,jj,kk,mm,Nx,Ny,Nz,r,b,Nk,NL

	real bin

      parameter(Nx=64,Ny=Nx/2,Nz=64,bin=0.02,Nk=8,b=250)

	integer ab(1:2,1:Nx,1:Ny,1:Nz),nsk(1:b),nsk1(1:b)

	real*8 di,rabi(0:Nk,0:Nk,0:Nk),sk(0:b),rabr(0:Nk,0:Nk,0:Nk)

	real*8 ccc,k1,k2,k3,kx,x1,x2,x3,sk1(0:b),m0,m1,n0,n2,NNN

     +,Inst1,Inst2

	character*4 file_a,file_b    

	character*7 fn

	character*1 filen(1:7) ,filea(1:4)

	      

	equivalence(filen,fn)

	equivalence(file_b,filea)

	

	

	common /dd2/ ab  



	do i=0,b

	     

	sk(i)=0.0

	sk1(i)=0.0

	   

	enddo



	do i=0,NK

      do j=0,NK

      do k=0,NK





      rabi(i,j,k)=0.0

      rabr(i,j,k)=0.0



	enddo

	enddo

	enddo



	

	file_b=file_a 

	filen(1)=filea(1) 

	filen(2)=filea(2) 

	filen(3)=filea(3) 

	filen(4)=filea(4)  

	filen(5)='.'

	filen(6)='s'

	filen(7)='k'

      kx=3.1415926/Nx

      NNN=1.0*Nx*Ny*Nz

	NL=int(NK/bin/2)

	ccc=0.0

	Inst1=0.0   

	Inst2=0.0



	do i=0,NK

	k1=kx*i



      do j=0,NK

	k2=kx*j



      do k=0,NK

	k3=kx*k



      

	do ii=1,2

	di=(ii-1.0)*0.5

      do jj=1, Nx

      x1=(jj-1.0)+di

      do kk=1, Ny

	x2=(kk-1.0)+di

      do mm=1, Nz

	x3=(mm-1.0)+di





      rabi(i,j,k)=rabi(i,j,k)+ab(ii,jj,kk,mm)

     +*sin(k1*x1+k2*x2+k3*x3)/(NNN)

      rabr(i,j,k)=rabr(i,j,k)+ab(ii,jj,kk,mm)

     +*cos(k1*x1+k2*x2+k3*x3)/(NNN)

	ccc=ccc+1.0



	enddo

	enddo

	enddo

	enddo

	

	r=int(sqrt(k1**2.0+k2**2.0+k3**2.0)/bin+1)

      

	if(mod((i+j+k),2).eq.0) then



	sk(r)=sk(r)+rabi(i,j,k)**2.0+rabr(i,j,k)**2.0

      nsk(r)=nsk(r)+1



	else

      

      sk1(r)=sk1(r)+rabi(i,j,k)**2.0+rabr(i,j,k)**2.0

	nsk1(r)=nsk1(r)+1



	endif



	enddo

	enddo

	enddo

	open(2,file=fn)

	do i=0, b

      write(2,'(f10.6,2(f40.20))') i*bin,sk(i),sk1(i)   

      Inst1=Inst1+sk(i)

      Inst2=Inst2+sk1(i)

	enddo

	close(2)



      m0=0.0

	m1=0.0

	n0=0.0

	n2=0.0

      do i=1,NL



	if(nsk(i).ne.0.and.nsk1(i).ne.0) then

      

      m0=m0+abs(1.0*sk(i)/nsk(i))

	m1=m1+abs(1.0*i*sk(i)/nsk(i))



	n0=n0+abs(1.0*sk1(i)/nsk1(i))

	n2=n2+abs(1.0*((i*bin)**2.0)*sk1(i)/nsk1(i))

      

	endif



	enddo

	

	print*,' m0:',m0,' m1:',m1,' n0:',n0,' n2:',n2



	print*,'R1=',m0/m1,'  R2=',sqrt(n0/n2)

	

	print*,'Inst1:',Inst1, 'Inst2:',Inst2



	print*,'calculation in s of k:',ccc





	end



      

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c      random number generator

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      

      

      double precision function ranno(ix)

      integer*4 a,p,ix,b15,b16,xhi,xalo,leftlo,fhi,k

      real*8 inverse

      parameter (inverse=4.65661287d-10)

      parameter (a=16807)

      parameter (b15=32768)

      parameter (b16=65536)

      parameter (p=2147483647)

      xhi=ix/b16

      xalo=(ix-xhi*b16)*a

      leftlo=xalo/b16

      fhi=xhi*a+leftlo

      k=fhi/b15

      ix=(((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k

        if (ix.lt.0) ix=ix+p

      ranno=dfloat(ix)*inverse

      continue

      return

      end  

      

c!**********************************************************************************************

c!*     calculate g of r

c!**********************************************************************************************

	subroutine c_gr(fa,fb)

      integer i,Nx,Ny,Nz,ix,iy,iz,ii,iix,iiy,iiz,r,b,NL,x,y,z,LN

      real*4 bin,x1,x2,y1,y2,z1,z2,dx,dy,dz,di,fa,fb

      parameter(Nx=64,Ny=Nx/2,Nz=64,bin=0.02,b=Nx*50/4+1)

	integer ab(1:2,1:Nx,1:Ny,1:Nz),gra(0:b),grb(0:b),grab(0:b),gr(0:b)

      real*8  ccc

	common /dd2/ ab

      

	NL=3

	NN=int(0.5*(3**(1.0/3.0))/bin)

	LN=int(NL/bin)

	print*,fa,fb



	ccc=0.0

	do i=1,2

       do ix=1,Nx

        x1=1.0*ix

	  do iy=1,Ny

	   y1=1.0*iy

	   do iz=1,Nz

          z1=1.0*iz

	  

	  do iix=ix-NL,ix+NL

	   x=mod(iix+Nx,Nx)+1

         x2=1.0*iix

	   do iiy=iy-NL,iy+NL

	    y=mod(iiy+Ny,Ny)+1

  	    y2=1.0*iiy

		do iiz=iz-NL,iz+NL

		 z=mod(iiz+Nz,Nz)+1	    

	     z2=1.0*iiz

	do ii=1,2

	if(i.eq.ii) then

	 di=0.0

	else 

	if(i.eq.1) then

	 di=-0.5

	else

	 di=0.5

	endif

	endif

	  dx=x2-x1+di

	  dy=y2-y1+di

	  dz=z2-z1+di

	  r=int(sqrt(dx**2.0+dy**2.0+dz**2.0)/bin)

	if(r.ge.NN.and.r.le.LN) then

	ccc=ccc+1.0

	  if(ab(i,ix,iy,iz).eq.ab(ii,x,y,z)) then

	    if(ab(i,ix,iy,iz).eq.1) then

           grb(r)=grb(r)+1

		   gr(r)=gr(r)+1

	    else

	     gra(r)=gra(r)+1

		   gr(r)=gr(r)+1

          endif

	  else

	     grab(r)=grab(r)+1

		   gr(r)=gr(r)+1

	  endif

	endif

	enddo



       

	 

	    enddo

	   enddo

	  enddo





	   enddo

	  enddo

	 enddo

	enddo	





	

	open(2,file='gr.dat',form='formatted')

	do i=1,b

       if(gr(i).eq.0.0) then

        write(2,'(f10.5,4(f35.15))') i*bin, 0.0, 0.0, 0.0,0.0

	 else

        write(2,'(f10.5,4(f35.15))') i*bin,1.0*gra(i)/(gr(i)*fa)

     +,1.0*grb(i)/(gr(i)*fb),1.0*grab(i)/(gr(i)),1.0*gr(i)

	 endif

	enddo

      close(2)

	print*,ccc

	end





      

      



        