! $Id: pdf.f90,v 1.7 2004/02/05 10:17:21 nilshau Exp $
!
!***********************************************************************
program pdf
!
!
!-----------------------------------------------------------------------
!   24-july-03/nils: code
!
use Cdata
use Cparam
use Sub
use Deriv
use General
!
implicit none
!
!
!  define parameters
!
character (len=130) ::  file,pdf_result_file,file_ext,pdf_in_file
character (len=30) :: variable
character (len=4) :: ch,csnap
integer :: cpu
real :: dummy
double precision :: rKB=0,mean_K2B2=0,mean_K2=0,mean_B2=0
real :: x1p,y1p,z1p
real :: x1frac,x2frac,y1frac,y2frac,z1frac,z2frac
real, dimension (mx,my,mz,mvar+maux) :: 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,bp,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) :: ppdf
real, dimension (nx) :: K,K2,K2B2
integer, dimension (nx) :: K_int
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,i,dir_nr,dir_count,jj,ind2
!
! Setting some constants
!
iux=1 
iuy=2
iuz=3
ilnrho=4
iax=5
iay=6
iaz=7
!
ip=10
!
! Reading pdf.in
!
call safe_character_assign(pdf_in_file,'pdf.in')
if (ip<8) print*,'Reading '//trim(pdf_in_file)
open(1,FILE=pdf_in_file,FORM='formatted')
read(1,*) variable
read(1,*) snap_start,snap_stop
close(1)
!
!  Loop over snapshots
!
do snap=snap_start,snap_stop
  ppdf=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)
     open(1,file=file,form='unformatted')
     read(1) f(1:mx,1:my,1:mz,1:mvar)
     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)
     !
     ! Calculating PDFs
     !
     do m=m1,m2
        do n=n1,n2
           call curl(f,iax,bb)
           call dot2_mn(bb,b2)
           b4=b2*b2
           bp(:,1)=bb(:,1)/sqrt(b2)
           bp(:,2)=bb(:,2)/sqrt(b2)
           bp(:,3)=bb(:,3)/sqrt(b2)
           if (variable .eq. 'curvature') then
              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)
              K=sqrt(BnB2)
              !print*,minval(sqrt(BnB2)/b2),maxval(sqrt(BnB2)/b2)
              K_int=int(K)
              do jj=1,nx
                 if (K_int(jj)+1 .gt. nx) K_int(jj)=nx-1
                 ind2=K_int(jj)+1
                 ppdf(ind2)=ppdf(ind2)+1
              enddo
              K2=K*K
              K2B2=K2*b2
              mean_K2B2=mean_K2B2+sum(K2B2)/nwgrid
              mean_b2=mean_b2+sum(b2)/nwgrid
              mean_K2=mean_K2+sum(K2)/nwgrid
           else
              print*,'No such variable:',variable
           endif
        enddo
     enddo
     !
     !  Finnishing loop over processors
     !
  enddo

print*,'PDF=',ppdf
print*,' '
print*,'r_{K,B}=',mean_K2B2/(mean_b2*mean_K2)-1.

print*,mean_K2
print*,mean_B2
print*,mean_K2B2

  !
  !  Saving the results
  !
  if (variable .eq. 'curvature') then
     file_ext='curvature'
  else if (variable .eq. 'normal') then
     file_ext='normal'
  else if (variable .eq. 'binormal') then
     file_ext='binormal'
  else
     print*,'No such variable:',variable
  endif
  if (snap_stop .eq. 0) call safe_character_assign(pdf_result_file, &
       trim(datadir)//'/pdf_'//trim(file_ext)//'.dat')
  if (snap_stop .ne. 0) call safe_character_assign(pdf_result_file, &
       trim(datadir)//'/pdf_'//trim(file_ext)//'_VAR'//trim(csnap)//'.dat')
  close(1)
  if (ip<8) print*,'Writing '//trim(pdf_result_file)
  open(1,file=pdf_result_file)
  write(1,'(3e12.4)') t
  write(1,'(3e12.4)') ppdf
  close(1)
  !
  !  Finnishing loop over snapshots (if snap_stop .ne. 0)
  !
enddo
!
end program pdf