! $Id: correlation.f90,v 1.8 2003/08/12 08:28:40 nilshau Exp $
!
!***********************************************************************
program correlation
!
! This program calculates the correlation function along the local magnetic
! field, in addition to along the normal and binormal to the local magnetic
! field. 
! The program reads in the data from ALL the different processors, 
! and saves it in a global array, so
! there is an upper limit on how large the runs can be.
!
!-----------------------------------------------------------------------
!   22-july-03/nils: code
!
use Cdata
use Cparam
use Sub
use Deriv
use General
!
implicit none
!
!
!  define parameters
!
character (len=130) ::  file,corr_result_file,file_ext,corr_in_file
character (len=30) :: direction,direction_
character (len=4) :: ch,csnap
integer :: cpu
real :: dummy
real :: x1p,y1p,z1p
real :: x1frac,x2frac,y1frac,y2frac,z1frac,z2frac
real, dimension (mx,my,mz,mvar) :: f
real, dimension (mx,my,mz,3) :: bb_loc,dir_loc
real, dimension (nxgrid,nygrid,nzgrid,3) :: dir_glob,bb_glob
real, dimension (nx,3) :: bb,BnB,b_parallel,normal,binormal
real, dimension (nx) :: b2,BnB2,l_d,Cx,Cy,Cz,b4
real, dimension (nx) :: xxx,xxy,xxz,yyx,yyy,yyz,zzx,zzy,zzz
real, dimension (nx) :: bxx,bxy,bxz,byx,byy,byz,bzx,bzy,bzz,a1,a2
double precision, dimension (nx/2) :: corr
real, dimension (3) :: B0,B_vec
integer :: snap,snap_start,snap_stop
integer :: l1_global,l2_global,m1_global,m2_global,n1_global,n2_global
integer :: xxp1,xxp2,yyp1,yyp2,zzp1,zzp2
integer :: ll,mmm,nnn,i,dir_nr,dir_count
!
! Setting some constants
!
iux=1 
iuy=2
iuz=3
ilnrho=4
iax=5
iay=6
iaz=7
!
ip=10
!
! Reading correlation.in
!
call safe_character_assign(corr_in_file,'correlation.in')
if (ip<8) print*,'Reading '//trim(corr_in_file)
open(1,FILE=corr_in_file,FORM='formatted')
read(1,*) direction_
read(1,*) snap_start,snap_stop
close(1)
!
! Doing all directions (parallel, normal and binormal) if direction_='all',
! otherwise doing only direction=direction_
!
if (direction_ .eq. 'all') then
   dir_count=3
else
   dir_count=1
endif
do dir_nr=1,dir_count
   if (direction_ .eq. 'all') then 
      if (dir_nr .eq. 1) direction='parallel'
      if (dir_nr .eq. 2) direction='normal'
      if (dir_nr .eq. 3) direction='binormal'
   else
      direction=direction_
   end if
!
print*,'Calculating corr_',direction
!
!  Loop over snapshots
!
do snap=snap_start,snap_stop
  !
  !  Loop over number of CPU's to get global array (f_global)
  !  Due to the fact that we are reading in data from all processors there is 
  !  an upper limit on number of gridpoints that do fit in memory. This should 
  !  idealy have been done in parallel, but for the time beeing I do not think 
  !  it is worth the effort (it will be pretty messy and demand a lot of 
  !  communication between processors).
  !
  corr=0
  do cpu=0,ncpus-1
     !print*,'Reading processor',cpu+1,'of',ncpus,'processors'
     ipy = modulo(cpu, nprocy)
     ipz = cpu/(nprocy)
     !
     ! Read var.dat 
     !
     call chn(cpu,ch)
     call chn(snap,csnap)
     if (snap_stop .eq. 0) call safe_character_assign(file,trim(datadir) &
          //'/proc'//trim(ch)//'/var.dat')
     if (snap_stop .ne. 0) call safe_character_assign(file,trim(datadir) &
          //'/proc'//trim(ch)//'/VAR'//trim(csnap))
     if (ip<8) print*,'Reading '//trim(file)
     print*,'Reading ',trim(file)
     open(1,file=file,form='unformatted')
     read(1) f
     if (lshear) then
        read(1) t,x,y,z,dx,dy,dz,deltay
     else
        read(1) t,x,y,z,dx,dy,dz
     endif
     close(1)
     !
     ! Finding direction for calculation of correlation function
     !
     do m=m1,m2
        do n=n1,n2
           call curl(f,iax,bb)
           call dot2_mn(bb,b2)
           b4=b2*b2
           b_parallel(:,1)=bb(:,1)/sqrt(b2)
           b_parallel(:,2)=bb(:,2)/sqrt(b2)
           b_parallel(:,3)=bb(:,3)/sqrt(b2)
           bb_loc(l1:l2,m,n,:)=bb
           if (direction .eq. 'parallel') then
              dir_loc(l1:l2,m,n,:)=b_parallel
           else if ((direction .eq. 'normal') .or. &
                (direction .eq. 'binormal') ) then
              !
              ! The normal is in the direction of the curvature, i.e. pointing
              ! towards the center of the "circle".
              ! (normal = b_parallel . nabla b_parallel)
              ! The binormal is the cross product of the local magnetic field
              ! and the normal.
              ! (binormal = normal x b_parallel)
              !
              call derij(f,iaz,a1,2,1)
              call derij(f,iay,a2,3,1)
              bxx=a1-a2
              call der2(f,iaz,a1,2)
              call derij(f,iay,a2,3,2)
              bxy=a1-a2
              call derij(f,iaz,a1,2,3)
              call der2(f,iay,a2,3)
              bxz=a1-a2
              call der2(f,iaz,a1,1)
              call derij(f,iax,a2,3,1)
              byx=-a1+a2
              call derij(f,iaz,a1,1,2)
              call derij(f,iax,a2,3,2)
              byy=-a1+a2
              call derij(f,iaz,a1,1,3)
              call der2(f,iax,a2,3)
              byz=-a1+a2
              call der2(f,iay,a1,1)
              call derij(f,iax,a2,2,1)
              bzx=a1-a2
              call derij(f,iay,a1,1,2)
              call der2(f,iax,a2,2)
              bzy=a1-a2
              call derij(f,iay,a1,1,3)
              call derij(f,iax,a2,2,3)
              bzz=a1-a2
              !
              BnB(:,1)=(bb(:,1)*bxx+bb(:,2)*bxy+bb(:,3)*bxz)/b2
              BnB(:,2)=(bb(:,1)*byx+bb(:,2)*byy+bb(:,3)*byz)/b2
              BnB(:,3)=(bb(:,1)*bzx+bb(:,2)*bzy+bb(:,3)*bzz)/b2
              !
              ! Must add second part due to derivative of 1/|B|
              !
              Cx=bb(:,1)*bxx+bb(:,2)*byx+bb(:,3)*bzx
              Cy=bb(:,1)*bxy+bb(:,2)*byy+bb(:,3)*bzy
              Cz=bb(:,1)*bxz+bb(:,2)*byz+bb(:,3)*bzz
              !
              xxx=bb(:,1)*bb(:,1)*Cx/b4 
              yyx=bb(:,2)*bb(:,1)*Cy/b4 
              zzx=bb(:,3)*bb(:,1)*Cz/b4 
              !
              xxy=bb(:,1)*bb(:,2)*Cx/b4 
              yyy=bb(:,2)*bb(:,2)*Cy/b4 
              zzy=bb(:,3)*bb(:,2)*Cz/b4 
              !
              xxz=bb(:,1)*bb(:,3)*Cx/b4 
              yyz=bb(:,2)*bb(:,3)*Cy/b4 
              zzz=bb(:,3)*bb(:,3)*Cz/b4 
              !
              BnB(:,1)=BnB(:,1)-(xxx+yyx+zzx)
              BnB(:,2)=BnB(:,2)-(xxy+yyy+zzy)
              BnB(:,3)=BnB(:,3)-(xxz+yyz+zzz)
              !
              call dot2_mn(BnB,BnB2)
              normal(:,1)=BnB(:,1)/sqrt(BnB2)
              normal(:,2)=BnB(:,2)/sqrt(BnB2)
              normal(:,3)=BnB(:,3)/sqrt(BnB2)
              if (direction .eq. 'normal') then
                 dir_loc(l1:l2,m,n,:)=normal
              else
                 binormal(:,1)=+(bb(:,2)*normal(:,3)-bb(:,3)*normal(:,2))
                 binormal(:,2)=-(bb(:,1)*normal(:,3)-bb(:,3)*normal(:,1))
                 binormal(:,3)=+(bb(:,1)*normal(:,2)-bb(:,2)*normal(:,1))
                 call dot2_mn(binormal,l_d)
                 dir_loc(l1:l2,m,n,1)=binormal(:,1)/sqrt(l_d)
                 dir_loc(l1:l2,m,n,2)=binormal(:,2)/sqrt(l_d)
                 dir_loc(l1:l2,m,n,3)=binormal(:,3)/sqrt(l_d)
              endif
           else
              print*,'No such direction:',direction
           endif
        enddo
     enddo
     !
     ! Putting local variables in global array  
     !
     l1_global=1
     l2_global=nxgrid
     m1_global=1+ipy*ny
     m2_global=(1+ipy)*ny
     n1_global=1+ipz*nz
     n2_global=(1+ipz)*nz
     !
     dir_glob(l1_global:l2_global,m1_global:m2_global,n1_global:n2_global,:)= &
          dir_loc(l1:l2,m1:m2,n1:n2,:)
     bb_glob(l1_global:l2_global,m1_global:m2_global,n1_global:n2_global,:)= &
          bb_loc(l1:l2,m1:m2,n1:n2,:)
     !
     !  Finnishing loop over processors
     !
  enddo
  !
  !  Start calculation of correlation functions
  !
  print*,'Start calculation of correlation functions'
  do ll=1,nxgrid
     print*,'Doing ',ll,' of ',nxgrid
     do mmm=1,nygrid
        do nnn=1,nzgrid
           B0=bb_glob(ll,mmm,nnn,:)
           !
           !  Beginning loop over all separations
           !
           do i=0,nx/2-1
              !
              !  Finding fractional coordinates
              !
              x1p=ll+i*dir_glob(ll,mmm,nnn,1)
              y1p=mmm+i*dir_glob(ll,mmm,nnn,2)
              z1p=nnn+i*dir_glob(ll,mmm,nnn,3)
              !
              !  Finding integer coordinates to use for interpolation
              !
              xxp1=mod(int(x1p+3*nx-1),nx)+1
              xxp2=mod(xxp1+1+nx-1,nx)+1
              yyp1=mod(int(y1p+3*nx-1),nx)+1
              yyp2=mod(yyp1+1+nx-1,nx)+1
              zzp1=mod(int(z1p+3*nx-1),nx)+1
              zzp2=mod(zzp1+1+nx-1,nx)+1
              !
              !  Determining importance of each of the eight interpolation 
              !  points.
              !
              x2frac=mod(x1p+real(nx)-1,real(nx))+1 - real(xxp1)
              x1frac=1.-x2frac    
              y2frac=mod(y1p+real(nx)-1,real(nx))+1 - real(yyp1)
              y1frac=1.-y2frac
              z2frac=mod(z1p+real(nx)-1,real(nx))+1 - real(zzp1)
              z1frac=1.-z2frac
              !
              !  Using interpolation to find field at fractional
              !  grid point.
              !
              B_vec=  x1frac*y1frac*z1frac*bb_glob(xxp1,yyp1,zzp1,:)+ &
                   x1frac*y1frac*z2frac*bb_glob(xxp1,yyp1,zzp2,:)+ &
                   x1frac*y2frac*z1frac*bb_glob(xxp1,yyp2,zzp1,:)+ &
                   x1frac*y2frac*z2frac*bb_glob(xxp1,yyp2,zzp2,:)+ &
                   x2frac*y1frac*z1frac*bb_glob(xxp2,yyp1,zzp1,:)+ &
                   x2frac*y1frac*z2frac*bb_glob(xxp2,yyp1,zzp2,:)+ &
                   x2frac*y2frac*z1frac*bb_glob(xxp2,yyp2,zzp1,:)+ &
                   x2frac*y2frac*z2frac*bb_glob(xxp2,yyp2,zzp2,:)
              !
              !  Adding up the results
              !
              corr(i+1)=corr(i+1)+ &
                   B0(1)*B_vec(1)+B0(2)*B_vec(2)+B0(3)*B_vec(3)
           end do
        end do
     end do
  end do
  !
  ! Normalizing
  !
  corr=corr/corr(1)
  !
  !  Saving the results
  !
  if (direction .eq. 'parallel') then
     file_ext='parallel'
  else if (direction .eq. 'normal') then
     file_ext='normal'
  else if (direction .eq. 'binormal') then
     file_ext='binormal'
  else
     print*,'No such direction:',direction
  endif
  if (snap_stop .eq. 0) call safe_character_assign(corr_result_file, &
       trim(datadir)//'/corr_'//trim(file_ext)//'.dat')
  if (snap_stop .ne. 0) call safe_character_assign(corr_result_file, &
       trim(datadir)//'/corr_'//trim(file_ext)//'_VAR'//trim(csnap)//'.dat')
  close(1)
  if (ip<8) print*,'Writing '//trim(corr_result_file)
  open(1,file=corr_result_file)
  write(1,'(3e12.4)') t
  write(1,'(3e12.4)') corr
  close(1)
  !
  !  Finnishing loop over snapshots (if snap_stop .ne. 0)
  !
enddo
!
!  Finnishing loop over directions (if direction_ .eq. 'all')
!
enddo
!
end program correlation