subroutine ChargeAssign()  ! ChargeAssign PotFromPoisson
    use particles
    use scattering
    use simulation ! int- Algo_AssignCharge, num_Nodes // real- time_Poisson, NodeSize, ChargeForPoisson  real(dim)::loc_Node, pot_Node, num_Node, charge_Node  ! 
    integer :: i, j, MinNode
    real ( kind = 8 ) :: MinDist, tempDist, tempWeight
    num_Node(1:num_Nodes) = 0.0d0
    do i = 1,num_maxCarriers
      if (Par_use(i).eq.1) then ! used particle
        if (Algo_AssignCharge.eq.1) then ! 1:NGP, Nearest-Grid-Point 2: Cloud-in-Cell
          MinDist = box_len(1) ! num_Nodes
          do j = 1,num_Nodes ! 
            tempDist = abs(Par_pos(i,1) - loc_Node(j))
            if (tempDist.lt.MinDist) then
              MinDist = tempDist
              MinNode = j
            endif
          enddo
          num_Node(MinNode) = num_Node(MinNode) + 1.0
        elseif (Algo_AssignCharge.eq.2) then ! CIC 
          do j = 1,num_Nodes ! 
            tempDist = abs(Par_pos(i,1) - loc_Node(j))
            if (tempDist.lt.NodeSize) then
              tempWeight = 1 - (tempDist/NodeSize)  
              num_Node(j) = num_Node(j) + tempWeight  
            endif
          enddo
        endif
      endif
    enddo
! Charge Recalculation 
    num_Node(1) = num_Node(1)*2.0
    num_Node(num_Nodes) = num_Node(num_Nodes)*2.0
    do j = 1,num_Nodes !
      charge_Node(j) = -(num_Node(j)*ChargeForPoisson - n_d)*e_c
    enddo    
    return
end

subroutine PotFromPoisson() 
    use particles
    use scattering
    use simulation
    integer :: i, j, nAlBin
    real ( kind = 8 ) :: tempFactor,posX,LowAlBound,UpAlBound ! MinDist, tempDist, tempWeight
    real ( kind = 8 ), dimension(:), allocatable :: MatEle_a, MatEle_b, MatEle_c, MatEle_r, MatEle_bp, MatEle_rp, temp_field 
    allocate ( MatEle_a(num_Nodes) )
    allocate ( MatEle_b(num_Nodes) )
    allocate ( MatEle_c(num_Nodes) )
    allocate ( MatEle_r(num_Nodes) )
    allocate ( MatEle_bp(num_Nodes) )
    allocate ( MatEle_rp(num_Nodes) )
    allocate ( temp_field(num_Nodes) )
  ! Solving Poisson-Finite Difference Scheme..[b c 0 0... a b c 0 ..  0 0 a b] b = -2/ a_1~(n-1) = 1, a_n = 2 / MatEle_c_2~n = 1, MatEle_c_1 = 2
    
    do i = 1,num_Nodes
      posX = loc_Bin(i)
      call DecideAlBin(posX,nAlBin,LowAlBound,UpAlBound)
      MatEle_a(i) = 1.0
      MatEle_b(i) = -2.0
      MatEle_c(i) = 1.0
      MatEle_r(i) = charge_Node(i)*NodeSize*NodeSize*4*pi/nBin_eps_low(nAlBin) 
    enddo
   
    MatEle_c(1) = 0.0  
    MatEle_r(1) = 0.0  
    MatEle_b(1) = 1.0  
    MatEle_a(num_Nodes) = 0.0
    MatEle_r(num_Nodes) = 0.0
    MatEle_b(num_Nodes) = 1.0
 !   MatEle_a(num_Nodes) = 2.0
 !   MatEle_b(num_Nodes) = -2.0
    
    MatEle_bp(1) = MatEle_b(1)
    MatEle_rp(1) = MatEle_r(1) ! num_Nodes  charge_Node
    do i = 2,num_Nodes
      tempFactor = MatEle_a(i)/MatEle_bp(i-1)
      MatEle_bp(i) = MatEle_b(i) - tempFactor*MatEle_c(i-1)
      MatEle_rp(i) = MatEle_r(i) - tempFactor*MatEle_rp(i-1)
    enddo 
    pot_Node(num_Nodes) = MatEle_rp(num_Nodes)/MatEle_bp(num_Nodes)
        !The second pass (back-substition)
    do i = num_Nodes-1, 1, -1
      pot_Node(i) = (MatEle_rp(i) - MatEle_c(i)*pot_Node(i+1))/MatEle_bp(i)
    enddo
    do i = 1,(num_Nodes-1)
      temp_field(i) = (pot_Node(i+1)-pot_Node(i))/NodeSize !center between two nodes
    enddo
    do i = 2,(num_Nodes-1)
      field_Node(i) = (temp_field(i)+temp_field(i-1))*0.5
    enddo
    field_Node(1) = 0.0
    field_Node(num_Nodes) = 0.0
    return
end

function PoissonField(posX)
  use particles
  use scattering
  use simulation
  real ( kind = 8 ) :: posX
  real ( kind = 8 ) :: PoissonField
  real ( kind = 8 ) :: tempField, tempDist, tempWeight, MinDist
  integer :: j, MinNode
  if (PoissonSetting.eq.1) then 
    if (Algo_PoiFieldAssign.eq.1) then  ! NGP
      MinDist = box_len(1) ! num_Nodes
      do j = 1,num_Nodes ! 
        tempDist = abs(posX - loc_Node(j))
        if (tempDist.lt.MinDist) then
          MinDist = tempDist
          MinNode = j
        endif
      enddo
      tempField = field_Node(MinNode)
    elseif (Algo_PoiFieldAssign.eq.2) then ! CIC
      tempField = 0.0
      do j = 1,num_Nodes ! 
        tempDist = abs(posX - loc_Node(j))
        if (tempDist.lt.NodeSize) then
          tempWeight = 1 - (tempDist/NodeSize)  
          tempField = tempField + tempWeight*field_Node(j)  
        endif
      enddo
    endif
  else
    tempField = 0.0d0
  endif
  PoissonField = tempField
  return
end  
! logical, external :: moved_too_much  ListUpdateRequested = moved_too_much ( Skin )