!####################################################################
!> @author Holger Grosshans
!> @brief solve E field
      subroutine solveElectrostatics
      use var
      real(kind=pr) err,err00,err2,syncMax
      integer it


      if (pnd.ne.0._pr.and.syncMax(maxval(q_el)).ne.0._pr) then

      do it=1,itmax
        call calcPhi_el(err)

        if (it.eq.1) err00=max(err,1e-19_pr)
        if (err/err00.lt.tol) exit
      enddo

      call calcE_el(err2)

      if (myid.eq.0) &
      write(*,'(a,es8.1e2,a,i3,a,es8.1e2)') &
      'E-pot.   |L2/L2_0 = ',err/err00,' |L2 #',it-1,' = ',err
       
      endif

      return
      end

!####################################################################
!> @author Holger Grosshans
!> @brief one sweep to relax phi_el field using Jacobi method
!> @param err        L2 error norm
      subroutine calcPhi_el(err)
      use var
      real(kind=pr),dimension(ii,jj,ll) :: phi_el1
      real(kind=pr) :: diag,dex,dey,dez,err,dphi_el,syncSum
      integer :: i,j,l

! phi_el:   electrostatic potential (V)
! rho_el:   charge density (C/m**3)
! eps_el:   permittivity vacuum/air (F/m)

      err=0._pr
      phi_el1=0._pr


      do 1 i=imin,imax,1
      do 1 j=jmin,jmax,1
      do 1 l=lmin,lmax,1

      if (celltype(i,j,l).ne.active) cycle

! solve the poisson equation:
        diag= (2._pr/((xc(i+1)-xc(i))*(xc(i)-xc(i-1))) &
              +2._pr/((yc(j+1)-yc(j))*(yc(j)-yc(j-1))) &
              +2._pr/((zc(l+1)-zc(l))*(zc(l)-zc(l-1))))

        dex= &
        (phi_el(i+1,j,l)*(xc(i)-xc(i-1))+phi_el(i-1,j,l)*(xc(i+1)-xc(i))) &
        /(0.5_pr*(xc(i+1)-xc(i-1))*(xc(i+1)-xc(i))*(xc(i)-xc(i-1)))
        dey= &
        (phi_el(i,j+1,l)*(yc(j)-yc(j-1))+phi_el(i,j-1,l)*(yc(j+1)-yc(j))) &
        /(0.5_pr*(yc(j+1)-yc(j-1))*(yc(j+1)-yc(j))*(yc(j)-yc(j-1)))
        dez= &
        (phi_el(i,j,l+1)*(zc(l)-zc(l-1))+phi_el(i,j,l-1)*(zc(l+1)-zc(l))) &
        /(0.5_pr*(zc(l+1)-zc(l-1))*(zc(l+1)-zc(l))*(zc(l)-zc(l-1)))

        phi_el1(i,j,l) = (rho_el(i,j,l)/eps_el + dex + dey + dez)/diag
        dphi_el=phi_el1(i,j,l)-phi_el(i,j,l)
        err=err+dphi_el*dphi_el

1     enddo


      phi_el=phi_el1
      call sync(phi_el)
      call bcPhi_el
      err=(syncSum(err)/dimgptot)**(0.5_pr)
      
      return
      end

!####################################################################
!> @author Holger Grosshans
!> @brief one sweep to relax E field using Jacobi method
!> @param err        L2 error norm
      subroutine calcE_el(err)
      use var
      real(kind=pr) :: dexel,deyel,dezel,temp,err,syncSum
      integer :: i,j,l

! phi_el:   electrostatic potential (V)
! rho_el:   charge density (C/m**3)
! E_el:     electrostatic field (V/m)
! eps_el:   permittivity vacuum/air (F/m)

      err=0._pr

      do 1 i=imin,imax,1
      do 1 j=jmin,jmax,1
      do 1 l=lmin,lmax,1

      if (celltype(i,j,l).ne.active) cycle

! calculate electrostatic field
        temp= Ex_el(i,j,l)
        Ex_el(i,j,l)= - (phi_el(i+1,j,l)-phi_el(i-1,j,l))/(xc(i+1)-xc(i-1))
        dexel=temp-Ex_el(i,j,l)

        temp= Ey_el(i,j,l)
        Ey_el(i,j,l)= - (phi_el(i,j+1,l)-phi_el(i,j-1,l))/(yc(j+1)-yc(j-1))
        deyel=temp-Ey_el(i,j,l)

        temp= Ez_el(i,j,l)
        Ez_el(i,j,l)= - (phi_el(i,j,l+1)-phi_el(i,j,l-1))/(zc(l+1)-zc(l-1))
        dezel=temp-Ez_el(i,j,l)
      
        err=err+dexel*dexel+deyel*deyel+dezel*dezel
                
1     enddo

      call sync(Ex_el); call sync(Ey_el); call sync(Ez_el)
      call bcE_el
      err=(syncSum(err)/dimgptot/3._pr)**(0.5_pr)
!      write(*,'(x,a,es9.2e2,a,3(es9.2e2))') &
!      'res. E =',err

      return
      end


!####################################################################
!> @author Holger Grosshans
!> @brief compute volumetric charge density
      subroutine chargeDensity
      use var
      real(kind=pr),dimension(3,3,3) :: weight
      real(kind=pr) :: volE
      integer       :: n,ibeg,iend,jbeg,jend,lbeg,lend,i,j,l, &
                       iw,jw,lw,ip,jp,lp

!      rho_el(imin:imax,jmin:jmax,lmin:lmax)=0._pr
      rho_el(:,jmin:jmax,lmin:lmax)=0._pr


      do 1 n=1,np
        if (celltype(ip(n),jp(n),lp(n)).ne.active) cycle

        call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,0)

        do 5 i=ibeg,iend
        do 5 j=jbeg,jend
        do 5 l=lbeg,lend
          iw=i+1-ibeg
          jw=j+1-jbeg
          lw=l+1-lbeg
          rho_el(i,j,l) = rho_el(i,j,l) + q_el(n)*partn(n)*weight(iw,jw,lw)/volE
5       enddo

1     enddo

      call bcRho_el

      return
      end


!####################################################################
!> @author Holger Grosshans
!> @brief compute electric forces on particles from E-field
      subroutine forcesGauss
      use var
      real(kind=pr),dimension(3,3,3) :: weight
      real(kind=pr) :: Extot,Eytot,Eztot,partmass,volE
      integer       :: n,ibeg,iend,jbeg,jend,lbeg,lend,i,j,l, &
                       iw,jw,lw,ip,jp,lp


      do 1 n=1,np
        if (celltype(ip(n),jp(n),lp(n)).ne.active) cycle

        Extot=0._pr   
        Eytot=0._pr   
        Eztot=0._pr   

        call weightLE8(weight,ibeg,iend,jbeg,jend,lbeg,lend,volE,n,0)

        do 5 i=ibeg,iend
        do 5 j=jbeg,jend
        do 5 l=lbeg,lend
          iw=i+1-ibeg
          jw=j+1-jbeg
          lw=l+1-lbeg
          Extot=Extot+Ex_el(i,j,l)*weight(iw,jw,lw)
          Eytot=Eytot+Ey_el(i,j,l)*weight(iw,jw,lw)
          Eztot=Eztot+Ez_el(i,j,l)*weight(iw,jw,lw)
5       enddo

        partmass=4._pr/3._pr*pi*rhop*partn(n)*radp(n)**3
        fx_el(n)= q_el(n)*Extot/partmass
        fy_el(n)= q_el(n)*Eytot/partmass
        fz_el(n)= q_el(n)*Eztot/partmass

1     enddo


      return
      end

!####################################################################
!> @author Holger Grosshans
!> @brief compute electric forces on particles from Coulomb's law
      subroutine forcesCoulomb
      use var
      real(kind=pr) :: dis,dis3,disx,disy,disz, &
          Extot,Eytot,Eztot,partmass
      integer :: n1,n2,ip,jp,lp,m1,m2,n1end


      select case (elforceScheme)
      case (2)
        fx_el=0._pr
        fy_el=0._pr
        fz_el=0._pr

        do 1 n1=2,np
        do 2 n2=1,(n1-1)
          if (n1.eq.n2) cycle
          call forcesCoulombN1N2(n1,n2)
2       enddo
1       enddo

      case (3)
        do 3 i=1,ii
        do 3 j=1,jj
        do 3 l=1,ll

        do 4 m1=2,npic(i,j,l)
        do 5 m2=1,(m1-1)
          n1=nic(i,j,l,m1)
          n2=nic(i,j,l,m2)
          call forcesCoulombN1N2(n1,n2)
5       enddo
4       enddo
3       enddo

      end select

      return
      end


!####################################################################
!> @author Holger Grosshans
!> @brief compute Coulomb force of n2 acting on n1
      subroutine forcesCoulombN1N2(n1,n2)
      use var
      real(kind=pr) :: dis,dis3,disx,disy,disz,partmass
      integer :: n1,n2,ip,jp,lp


      disx=xp(n1)-xp(n2)
      disy=yp(n1)-yp(n2)
      disz=zp(n1)-zp(n2)
      dis=(disx**2+disy**2+disz**2)**0.5
      dis3=dis**3
  
      partmass=4._pr/3._pr*pi*rhop*partn(n1)*radp(n1)**3

      fx_el(n1)= fx_el(n1) &
               + q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
      fy_el(n1)= fy_el(n1) &
               + q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
      fz_el(n1)= fz_el(n1) &
               + q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass

      partmass=4._pr/3._pr*pi*rhop*partn(n2)*radp(n2)**3

      fx_el(n2)= fx_el(n2) &
               - q_el(n1)*q_el(n2)*disx/4._pr/pi/eps_el/dis3/partmass
      fy_el(n2)= fy_el(n2) &
               - q_el(n1)*q_el(n2)*disy/4._pr/pi/eps_el/dis3/partmass
      fz_el(n2)= fz_el(n2) &
               - q_el(n1)*q_el(n2)*disz/4._pr/pi/eps_el/dis3/partmass


      return
      end


!####################################################################
!> @author Holger Grosshans
!> @brief compute particle-wall charging
      subroutine chargeParticleWall(n,direction)
      use var
      real(kind=pr) :: ut2,un2,alpha1,AoAtot,dqp
      integer :: n,i,direction
      character*70 :: filename


! pipe:
!        ut2=up(n)**2
!        un2=vp(n)**2+wp(n)**2

      if (direction.eq.2) then
        ut2=up(n)**2+wp(n)**2
        un2=vp(n)**2
      elseif (direction.eq.3) then
        ut2=up(n)**2+vp(n)**2
        un2=wp(n)**2
      endif

! John (1979)
      alpha1=(0.625_pr*pi*rhop*(1._pr+restRatio)*(un2+ut2)* &
            ((1._pr-nyw**2)/Ew+(1._pr-nyp**2)/Ep))**(0.4_pr)
      AoAtot=alpha1/4._pr
      dqp=Qaccfactor*AoAtot*(qpmax-q_el(n))
      q_el(n)=min(dqp+q_el(n),qpmax)
      

      write(filename,'(a,i3.3,a)') 'results/particlesWall_p',myid,'.dat'
      open(access='append',unit=10,file=filename)
      write(10,77) nt,t,sqrt(ut2),sqrt(un2),AoAtot,dqp,q_el(n)
77    format(i7,6(x,es11.2e2))
      close(10)
 
      return
      end


!####################################################################
!> @author Holger Grosshans
!> @brief compute particle-particle charging
      subroutine chargeParticleParticle(n1,n2)
      use var
      real(kind=pr) uax2,urad2,alpha1,dqp
      integer :: n1,n2
      character*70 :: filename

      uax2=(up(n1)-up(n2))**2
      urad2=(vp(n1)-vp(n2))**2+(wp(n1)-wp(n2))**2

      if(urad2.eq.0.) then
        dqp=0.5*(q_el(n2)-q_el(n1))
        goto 10
      endif

! Soo (1971)
      alpha1=radp(n1)*radp(n2)*(0.625*pi*rhop*(1+restratio)*(uax2+urad2) &
             *(radp(n1)+radp(n2))**(0.5)/(radp(n1)**3+radp(n2)**3) &
             *(1-nyp**(2.0))/Ep)**(0.4)

      dqp=alpha1/8._pr/radp(n1)*(q_el(n2)-q_el(n1))

10    q_el(n1)=q_el(n1)+dqp
      q_el(n2)=q_el(n2)-dqp


      write(filename,'(a,i3.3,a)') 'results/particlesPart_p',myid,'.dat'
      open(access='append',unit=10,file=filename)
      write(10,77) nt,t,q_el(n1),q_el(n2),dqp
77    format(i7,5(x,es11.2e2))
      close(10)

     
      return
      end
