! fdist3d class for QuickPIC Open Source 1.0
! update: 04/18/2016

      module fdist3d_class

      use perrors_class
      use parallel_pipe_class
      use spect3d_class
      use ufield3d_class
      use part3d_lib
      use input_class
         
      implicit none

      private

      public :: fdist3d, fdist3d_000

      type, abstract :: fdist3d

         private

         class(spect3d), pointer, public :: sp => null()
         class(perrors), pointer, public :: err => null()
         class(parallel_pipe), pointer, public :: p => null()
!
! ndprof = profile type
         integer :: npf, npmax
                         
         contains
         
         generic :: new => init_fdist3d         
         generic :: del => end_fdist3d
         generic :: dist => dist3d
         procedure(ab_init_fdist3d), deferred, private :: init_fdist3d
         procedure, private :: end_fdist3d
         procedure(ab_dist3d), deferred, private :: dist3d
         procedure :: getnpf, getnpmax
                  
      end type 

      abstract interface
!
      subroutine ab_dist3d(this,part3d,npp,fd)
         import fdist3d
         import ufield3d
         implicit none
         class(fdist3d), intent(inout) :: this
         real, dimension(:,:), pointer, intent(inout) :: part3d
         integer, intent(inout) :: npp
         class(ufield3d), intent(in), pointer :: fd         
      end subroutine ab_dist3d
!
      subroutine ab_init_fdist3d(this,input,i)
         import fdist3d
         import input_json
         implicit none
         class(fdist3d), intent(inout) :: this
         type(input_json), intent(inout), pointer :: input
         integer, intent(in) :: i
      end subroutine ab_init_fdist3d
!
      end interface

      type, extends(fdist3d) :: fdist3d_000

         private

         integer :: npx, npy, npz
         real :: qm, sigx, sigy, sigz
         real :: bcx, bcy, bcz, sigvx, sigvy, sigvz
         real :: cx1,cx2,cx3,cy1,cy2,cy3,gamma,np
         logical :: quiet

         contains
         procedure, private :: init_fdist3d => init_fdist3d_000
         procedure, private :: dist3d => dist3d_000

      end type fdist3d_000

      character(len=10), save :: class = 'fdist3d:'
      character(len=128), save :: erstr
      
      contains
!
      function getnpf(this)

         implicit none

         class(fdist3d), intent(in) :: this
         integer :: getnpf
         
         getnpf = this%npf

      end function getnpf
!      
      function getnpmax(this)

         implicit none

         class(fdist3d), intent(in) :: this
         integer :: getnpmax
         
         getnpmax = this%npmax

      end function getnpmax
!      
      subroutine end_fdist3d(this)
          
         implicit none
         
         class(fdist3d), intent(inout) :: this
         character(len=18), save :: sname = 'end_fdist3d:'

         call this%err%werrfl2(class//sname//' started')
         call this%err%werrfl2(class//sname//' ended')
                  
      end subroutine end_fdist3d
!
      subroutine init_fdist3d_000(this,input,i)
      
         implicit none
         
         class(fdist3d_000), intent(inout) :: this
         type(input_json), intent(inout), pointer :: input
         integer, intent(in) :: i        
! local data
         integer :: npf,npx,npy,npz,npmax
         real :: qm,sigx,sigy,sigz,bcx,bcy,bcz,sigvx,sigvy,sigvz
         real :: cx1,cx2,cx3,cy1,cy2,cy3,gamma,np
         logical :: quiet
         real :: min, max, cwp, n0
         real :: alx, aly, alz, dx, dy, dz
         integer :: indx, indy, indz         
         character(len=20) :: sn,s1
         character(len=18), save :: sname = 'init_fdist3d_000:'
         
         this%sp => input%sp
         this%err => input%err
         this%p => input%pp

         call this%err%werrfl2(class//sname//' started')

         write (sn,'(I3.3)') i
         s1 = 'beam('//trim(sn)//')'

         call input%get('simulation.n0',n0)
         call input%get('simulation.indx',indx)
         call input%get('simulation.indy',indy)
         call input%get('simulation.indz',indz)

         cwp=5.32150254*1e9/sqrt(n0)
         call input%get('simulation.box.x(1)',min)
         call input%get('simulation.box.x(2)',max)
         call input%get(trim(s1)//'.center(1)',bcx)
         bcx = bcx - min
         alx = (max-min) 
         dx=alx/real(2**indx)
         call input%get('simulation.box.y(1)',min)
         call input%get('simulation.box.y(2)',max)
         call input%get(trim(s1)//'.center(2)',bcy)
         bcy = bcy -min
         aly = (max-min) 
         dy=aly/real(2**indy)
         call input%get('simulation.box.z(1)',min)
         call input%get('simulation.box.z(2)',max)
         call input%get(trim(s1)//'.center(3)',bcz)
         bcz = bcz -min
         alz = (max-min) 
         dz=alz/real(2**indz)

         call input%get(trim(s1)//'.profile',npf)
         call input%get(trim(s1)//'.np(1)',npx)
         call input%get(trim(s1)//'.np(2)',npy)
         call input%get(trim(s1)//'.np(3)',npz)
         call input%get(trim(s1)//'.q',qm)
         call input%get(trim(s1)//'.sigma(1)',sigx)
         call input%get(trim(s1)//'.sigma(2)',sigy)
         call input%get(trim(s1)//'.sigma(3)',sigz)
         call input%get(trim(s1)//'.sigma_v(1)',sigvx)
         call input%get(trim(s1)//'.sigma_v(2)',sigvy)
         call input%get(trim(s1)//'.sigma_v(3)',sigvz)
         call input%get(trim(s1)//'.centroid_x(1)',cx1)
         call input%get(trim(s1)//'.centroid_x(2)',cx2)
         call input%get(trim(s1)//'.centroid_x(3)',cx3)
         call input%get(trim(s1)//'.centroid_y(1)',cy1)
         call input%get(trim(s1)//'.centroid_y(2)',cy2)
         call input%get(trim(s1)//'.centroid_y(3)',cy3)
         call input%get(trim(s1)//'.quiet_start',quiet)
         call input%get(trim(s1)//'.gamma',gamma)
         call input%get(trim(s1)//'.peak_density',np)
         call input%get(trim(s1)//'.npmax',npmax)


         this%npf = npf
         this%npx = npx
         this%npy = npy
         this%npz = npz
         this%npmax = npmax
         qm = qm*np*(2*3.1415926535897932)**1.5*sigx*sigy*sigz
         qm = qm*(2**indz)
         qm = qm*(2**indx)
         qm = qm*(2**indy)/(npx*alx*aly*alz) 
         qm = qm/npy
         qm = qm/npz
         this%qm = qm
         this%bcx = bcx/dx
         this%bcy = bcy/dy
         this%bcz = bcz/dz
         this%sigx = sigx/dx
         this%sigy = sigy/dy
         this%sigz = sigz/dz
         this%sigvx = sigvx
         this%sigvy = sigvy
         this%sigvz = sigvz
         this%cx1 = cx1*dz*dz/dx
         this%cx2 = cx2*dz/dx
         this%cx3 = cx3/dx
         this%cy1 = cy1*dz*dz/dy
         this%cy2 = cy2*dz/dy
         this%cy3 = cy3/dy
         this%gamma = gamma
         this%np = np
         this%quiet = quiet


         call this%err%werrfl2(class//sname//' ended')

      end subroutine init_fdist3d_000
!
      subroutine dist3d_000(this,part3d,npp,fd)
      
         implicit none
         
         class(fdist3d_000), intent(inout) :: this
         real, dimension(:,:), pointer, intent(inout) :: part3d
         integer, intent(inout) :: npp
         class(ufield3d), intent(in), pointer :: fd
! local data1

! edges(1) = lower boundary in y of particle partition
! edges(2) = upper boundary in y of particle partition
! edges(3) = lower boundary in z of particle partition
! edges(4) = upper boundary in z of particle partition
         real, dimension(:,:), pointer :: pt => null()
         integer :: npx, npy, npz, nx, ny, nz, ipbc
         real :: vtx, vty, vtz, vdx, vdy, vdz
         real :: sigx, sigy, sigz, x0, y0, z0
         real, dimension(3) :: cx, cy
         real, dimension(4) :: edges
         integer, dimension(2) :: noff
         integer :: nps=1
         logical :: lquiet = .false.
         integer :: idimp, npmax, ierr = 0
         character(len=18), save :: sname = 'dist3d_000:'

         call this%err%werrfl2(class//sname//' started')
         
         npx = this%npx; npy = this%npy; npz = this%npz
         nx = fd%getnd1(); ny = fd%getnd2(); nz = fd%getnd3()
         ipbc = this%sp%getpsolver()
         pt => part3d
         vtx = this%sigvx; vty = this%sigvy; vtz = this%sigvz
         vdx = 0.0; vdy = 0.0; vdz = this%gamma
         sigx = this%sigx; sigy = this%sigy; sigz = this%sigz
         x0 = this%bcx; y0 = this%bcy; z0 = this%bcz
         cx = (/this%cx1,this%cx2,this%cx3/); cy = (/this%cy1,this%cy2,this%cy3/)
         lquiet = this%quiet
         idimp = size(part3d,1); npmax = size(part3d,2)
         noff = fd%getnoff()
         edges(1) = noff(1); edges(3) = noff(2)
         edges(2) = edges(1) + fd%getnd2p()
         edges(4) = edges(3) + fd%getnd3p()         
         
         call PRVDIST32_RANDOM(pt,this%qm,edges,npp,nps,vtx,vty,vtz,vdx,vdy,&
         &vdz,npx,npy,npz,nx,ny,nz,ipbc,idimp,npmax,1,1,4,sigx,sigy,sigz,&
         &x0,y0,z0,cx,cy,lquiet,ierr)

         if (ierr /= 0) then
            write (erstr,*) 'PRVDIST32_RANDOM error'
            call this%err%equit(class//sname//erstr)
         endif
         
         call this%err%werrfl2(class//sname//' ended')
         
      end subroutine dist3d_000
!
      end module fdist3d_class