subroutine nubeam_maxwell_plot(il_dep)
  use nbi_com
  use nbi_mpi_com
  implicit none
  logical, intent(in) :: il_dep
  real(fp),allocatable,dimension(:,:,:) :: wrk,zfbm
  real(fp),allocatable,dimension(:) :: v_wrk, v_wrk_pll ! v_perp, v_pll, grids
  integer, parameter ::  i_grid=51, i_grid_pll=101 ! v_perp, v_pll, grid size
  real(fp),parameter :: v_grid=2.5_fp ! max abs value of normalized v_perp and in normalized v_pll
                                    ! [-v_grid,v_grid]
  real(fp) :: vperp_norm,vpll_norm,hv,hv_pll,deltay
  integer :: ind,jp,j,i_mj,isb,i_vperp,i_vpll,inum,mpierr
  integer find_ptcl_ind

  isb=lsbeam
  if(allocated(wrk)) deallocate(wrk)
  if(allocated(v_wrk)) deallocate(v_wrk)
  if(allocated(v_wrk_pll)) deallocate(v_wrk_pll)
  i_mj=ledge-lcentr+1
  allocate(wrk(i_grid_pll,i_grid,i_mj))
  allocate(v_wrk(i_grid))
  allocate(v_wrk_pll(i_grid_pll))
  wrk=0.0_fp
  v_wrk=0.0_fp
  v_wrk_pll=0.0_fp
! creat a grid, same in normalized v_perp and in normalized v_pll,  [-v_grid,v_grid]
  v_wrk(1)=0.0_fp
  v_wrk(i_grid)=v_grid
  hv=v_grid/(i_grid-1)
  do j=2,i_grid-1
    v_wrk(j)=v_wrk(j-1)+hv
  enddo
  v_wrk_pll(1)=-v_grid
  v_wrk_pll(i_grid_pll)=v_grid
  hv_pll=2.0_fp*v_grid/(i_grid_pll-1)
  do j=2,i_grid_pll-1
    v_wrk_pll(j)=v_wrk_pll(j-1)+hv_pll
  enddo
!go over all mc ptcls
  do jp=1,nuse(lsbeam)
     ind=nbndex(jp,isb)
     ind=find_ptcl_ind(numprocd,myidd,ind) !find proper ptcl.id#
     deltay=abs(tbm2-tdepay(ind,isb))
     nflr=min(lep1,int(lcentr+xiay(ind,isb)*nzones))
     if(ebavps(nflr,isb).le.0.0_fp.or.dvolzc(nflr).le.0.0_fp.or.wghtay(ind,isb).eq.0.0_fp) then
        cycle
     endif
     if(abs(xksidy(ind,isb)).eq.1.0_fp) then
        vperp_norm=0.0_fp
     else
        vperp_norm=sqrt(vay(ind,isb)**2*(1.0_fp-xksidy(ind,isb)**2))/sqrt(ebavps(nflr,isb)/v2toev)
     endif
     i_vperp=int(1+vperp_norm/hv) !find index for v_perp
     vpll_norm=xksidy(ind,isb)*vay(ind,isb)/sqrt(ebavps(nflr,isb)/v2toev)
     i_vpll=int(1+(vpll_norm+v_grid)/hv_pll)   !find index for v_pll
     if(i_vperp.gt.i_grid.or.i_vpll.gt.i_grid_pll) then
        cycle
     endif
     if(i_vperp.lt.1.or.i_vpll.lt.1) then
        cycle
     endif
!fill distribution function
     nflr=nflr-lcentr+1 ! as my grid from 1 not lcenter
     if(vperp_norm.ne.0.0_fp) then
        wrk(i_vpll,i_vperp,nflr)=wrk(i_vpll,i_vperp,nflr)+wghtay(ind,isb)*deltay/delt/&
             dvolzc(nflr+1) /hv_pll/hv/vperp_norm/2.0_fp
     endif
  enddo
! mpi reduce
#ifdef __MPI
  if(allocated(zfbm)) deallocate(zfbm)
  allocate(zfbm(i_grid_pll,i_grid,i_mj))
  zfbm=0.0_fp
  inum = i_grid_pll*i_grid*i_mj
  
  mpierr=0
      call mpi_reduce(wrk(1,1,1),zfbm,inum,MPI_REAL8,MPI_SUM,0,nbi_comm,mpierr)
      if(myidd.eq.0) then
         wrk=zfbm 
      endif
      if(allocated(zfbm)) deallocate(zfbm)
#endif
! plot
      if(myidd.eq.0) then
         call r8_grf4f1(wrk,i_grid_pll,i_grid,v_wrk_pll,i_grid_pll,v_wrk,i_grid,xizc(lcentr:ledge),i_mj,&
             'max dis func','norm vpll','norm vperp','sqrt norm tf',' ',' ',' ',' ','maxwell test',0)
      endif
      if(allocated(wrk)) deallocate(wrk)
      if(allocated(v_wrk)) deallocate(v_wrk)
      if(allocated(v_wrk_pll)) deallocate(v_wrk_pll)
#ifdef __MPI
      call mpi_barrier(nbi_comm,mpierr) 
#endif
  return

end subroutine nubeam_maxwell_plot
