Homec4science
Diffusion FENNECS lite (3Dmpimode) Browse

fennecslite/src/secondary_electrons_mod.f90
b8eece80d5a53Dmpimode

Last Change
Raw File

secondary_electrons_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: SEE
!
!> @author
!> P. Giroud-Garampon - EPFL/SPC
!
!> Last modif.
!> 16.02.2024
!
! DESCRIPTION:
!> Module handling secondary electron emissions (SEE)
!------------------------------------------------------------------------------
MODULE see
USE particletypes
USE constants
USE basic
USE materials
IMPLICIT NONE
CONTAINS
!---------------------------------------------------------------------------
! SUBROUTINE secondary_electrons(pelec, losthole, pelec, nblostparts)
!
!
! DESCRIPTION
!> function to determine the number of electrons
!> to add as a function of the number of lost electrons
!--------------------------------------------------------------------------
SUBROUTINE secondary_electrons(pelec, losthole, nblostparts)
USE geometry
TYPE(particles), INTENT(INOUT):: pelec ! electrons list
REAL(KIND = db), DIMENSION(3) :: last_pos ! last position for lost electron
REAL(KIND = db), DIMENSION(3) :: normal_dir, dir ! unit vector for directions
INTEGER, DIMENSION(pelec%Nploc):: losthole ! indices of lost electrons
INTEGER :: nblostparts, Nploc, Nploc_old ! number of lost electrons
INTEGER :: parts_size_increase, nbadded ! number of added electrons
INTEGER :: i,j
INTEGER :: material_id
INTEGER :: kmax ! max number of electrons generated
REAL(KIND = db) :: Einc, thetainc ! angle of incidence of the impinging electron
REAL(KIND = db) :: eta ! backscattering coefficient
REAL(KIND = db) :: randnb(1:1) ! random number
REAL(KIND = db) :: lambda ! Poisson parameter
REAL(KIND = db) :: kappa, theta, Emax ! gamma distribution parameters
REAL(KIND = db) :: Eem, Vem ! energy of the emitted electron
kmax = 20
kappa = 1.4444
theta = 4.0606
Emax = 30 ! Max value for secondary electron energy
material_id = pelec%material_id
! increase the size of the parts array
IF(pelec%Nploc + 10*nblostparts .gt. size(pelec%pos,2)) THEN
parts_size_increase=Max(floor(0.1*size(pelec%pos,2)),10*nblostparts)
CALL change_parts_allocation(pelec, parts_size_increase)
END IF
DO i=1,nblostparts
Einc = comp_Ekin(pelec%U(:,losthole(i)), pelec%m) ! [eV]
last_pos = revert_push(pelec, losthole(i))
normal_dir = find_normal(last_pos)
thetainc = incidence(pelec%U(:,losthole(i)), normal_dir)
eta = comp_back_coef(Einc, thetainc, material_id)
call random_array(randnb,1,ran_index(1),ran_array(:,1)) ! uniform in [0,1[
IF(randnb(1) .lt. eta) THEN ! the electron is backscattered
pelec%nbadded = pelec%nbadded + 1
pelec%Nploc = pelec%Nploc + 1
Nploc = pelec%Nploc
! return to its former position
pelec%pos(1,Nploc) = last_pos(1)
pelec%pos(2,Nploc) = last_pos(2)
pelec%pos(3,Nploc) = last_pos(3)
pelec%newindex = pelec%newindex + 1
pelec%partindex(Nploc) = pelec%newindex
Eem = Einc ! keep the same energy [eV]
Vem = comp_V(Eem, pelec%m)
dir = gen_dir(normal_dir)
pelec%U(1,Nploc) = Vem*dir(1)
pelec%U(2,Nploc) = Vem*dir(2)
pelec%U(3,Nploc) = Vem*dir(3)
END IF
! Secondary electrons are generated
lambda = comp_yield(Einc, thetainc, material_id)
nbadded = gen_nb_elec(lambda, kmax)
Nploc_old = pelec%Nploc
pelec%Nploc = pelec%Nploc + nbadded
Nploc = pelec%Nploc
pelec%nbadded = pelec%nbadded + nbadded
DO j=1,nbadded
pelec%pos(1,Nploc_old+j) = last_pos(1)
pelec%pos(2,Nploc_old+j) = last_pos(2)
pelec%pos(3,Nploc_old+j) = last_pos(3)
pelec%newindex = pelec%newindex + 1
pelec%partindex(Nploc_old+j) = pelec%newindex
IF(pelec%zero_vel .eqv. .false.) THEN
Eem = gen_Eem(kappa, theta, Emax) ! generate an energy value following Gamma distribution
Vem = comp_V(Eem, pelec%m)
dir = gen_dir(normal_dir) ! generate a random direction following a cosine distribution
pelec%U(1,Nploc_old+j) = Vem*dir(1)
pelec%U(2,Nploc_old+j) = Vem*dir(2)
pelec%U(3,Nploc_old+j) = Vem*dir(3)
ELSE
pelec%U(1,Nploc_old+j) = 0.0
pelec%U(2,Nploc_old+j) = 0.0
pelec%U(3,Nploc_old+j) = 0.0
END IF
END DO
END DO
END SUBROUTINE secondary_electrons
!---------------------------------------------------------------------------
! FUNCTION comp_Ekin(velocity, m)
!
!
! DESCRIPTION
!> Computes the kinetic energy in eV of a particle given its 3 velocity components
!--------------------------------------------------------------------------
FUNCTION comp_Ekin(velocity, m) RESULT(Ekin)
REAL(KIND = db), DIMENSION(3) :: velocity
REAL(KIND = db) :: m, Ekin
Ekin = 0.5 * m * vlight**2 /elchar * (velocity(1)**2 + velocity(2)**2 + velocity(3)**2)
END FUNCTION comp_Ekin
!---------------------------------------------------------------------------
! FUNCTION find_normal(last_pos)
!
!
! DESCRIPTION
!> Computes the normal direction to the electrode surface
!--------------------------------------------------------------------------
FUNCTION find_normal(last_pos) RESULT(normal_dir)
USE geometry
REAL(KIND = db), DIMENSION(3) :: last_pos ! last position of the electron
REAL(KIND = db), DIMENSION(3) :: normal_dir ! unit normal direction vector
REAL(KIND = db), DIMENSION(3) :: weight ! geom. weight at last position
REAL(KIND = db) :: norm
call geom_weight(last_pos(1), last_pos(3), weight)
norm = sqrt(weight(2)**2 + weight(3)**2)
normal_dir(1) = 1/norm * weight(3) ! normal along r
normal_dir(2) = 0.0 ! normal along theta
normal_dir(3) = 1/norm * weight(2) ! normal along z
END FUNCTION find_normal
!---------------------------------------------------------------------------
! FUNCTION incidence(velocity, normal_dir)
!
!
! DESCRIPTION
!> Computes the angle of incidence of the electron
!--------------------------------------------------------------------------
FUNCTION incidence(velocity, normal_dir) RESULT(angle)
REAL(KIND = db), DIMENSION(3) :: velocity, normal_dir
REAL(KIND = db) :: angle
angle = ACOS(DOT_PRODUCT(-1*normal_dir,velocity)/SQRT(DOT_PRODUCT(velocity,velocity)))
END FUNCTION incidence
!---------------------------------------------------------------------------
! FUNCTION comp_V(Ekin, mass)
!
!
! DESCRIPTION
!> Computes the normalized velocity of a particle with energy Ekin
!--------------------------------------------------------------------------
FUNCTION comp_V(Ekin, mass) RESULT(V)
REAL(KIND = db) :: Ekin, mass, V
V = sqrt(2/mass * Ekin * elchar) / vlight
END FUNCTION comp_V
!---------------------------------------------------------------------------
! FUNCTION revert_push(pelec, partid)
!
!
! DESCRIPTION
!> reverts Buneman algorithm over one time step
!> to obtain electron position right before capture
!--------------------------------------------------------------------------
FUNCTION revert_push(pelec, partid)
USE basic, ONLY: dt
TYPE(particles), INTENT(INOUT):: pelec
REAL(KIND=db), DIMENSION(3):: revert_push
INTEGER :: partid
revert_push(1) = abs(pelec%pos(1,partid) - pelec%U(1,partid)*dt)
revert_push(2) = pelec%pos(2,partid) - 1/pelec%pos(1,partid)* pelec%U(2,partid)*dt
revert_push(3) = pelec%pos(3,partid) - pelec%U(3,partid)*dt
END FUNCTION revert_push
!---------------------------------------------------------------------------
! FUNCTION comp_back_coef(energy, angle, material_id)
!
!
! DESCRIPTION
!> Gives the theoretical value for the backscattering coefficient
!> as a function of the energy of the incident electron, the angle
!> of incidence and the material
!
!--------------------------------------------------------------------------
REAL(KIND = db) FUNCTION comp_back_coef(energy, angle, material_id)
REAL(KIND = db) :: energy, angle, eta_zero
INTEGER :: material_id
SELECT CASE(material_id)
CASE(2) ! Copper
! Value for normal incidence
! Angular dependence
CASE(3) ! Aluminum
! Value for normal incidence
IF(energy .le. 5.0E1) THEN
eta_zero = back_coef_Al_1
ELSE IF(energy .gt. 5.0E1 .and. energy .le. 2.0E2) THEN
eta_zero = eval_polynomial(back_coef_Al_2, energy)
ELSE IF(energy .gt. 2.0E2) THEN
eta_zero = back_coef_Al_3
END IF
! Angular dependence
IF(angle .le. 1.0472) THEN ! angle < 60 deg
comp_back_coef = eta_zero/COS(angle) ! inverse cosine dependence
ELSE
comp_back_coef = eta_zero*2
END IF
CASE DEFAULT ! Aluminum
IF(energy .le. 5.0E1) THEN
eta_zero = back_coef_Al_1
ELSE IF(energy .gt. 5.0E1 .and. energy .le. 2.0E2) THEN
eta_zero = eval_polynomial(back_coef_Al_2, energy)
ELSE IF(energy .gt. 2.0E2) THEN
eta_zero = back_coef_Al_3
END IF
! Angular dependence
IF(angle .le. 1.0472) THEN ! angle < 60 deg
comp_back_coef = eta_zero/COS(angle) ! inverse cosine dependence
ELSE
comp_back_coef = eta_zero*2
END IF
END SELECT
END FUNCTION comp_back_coef
!---------------------------------------------------------------------------
! FUNCTION comp_yield(energy, angle, material_id)
!
!
! DESCRIPTION
!> Gives the value for the secondary electron yield
!> as a function of the energy of the incident electron, the angle
!> of incidence and the material type
!
!--------------------------------------------------------------------------
REAL(KIND = db) FUNCTION comp_yield(energy, angle, material_id)
REAL(KIND = db) :: energy, angle, x, delta_m, delta_zero
INTEGER :: material_id
SELECT CASE(material_id)
CASE(2) ! Copper
! Value for normal incidence
x = energy/E_m_Cu
delta_m = delta_m_Cu
delta_zero = delta_m*1.28*x**(-0.67)*(1-EXP(-1.614*x**1.67))
! Angular dependence
comp_yield = delta_zero
CASE(3) ! Aluminum
! Value for normal incidence
x = energy/E_m_Al
delta_m = delta_m_Al
delta_zero = delta_m*1.28*x**(-0.67)*(1-EXP(-1.614*x**1.67))
! Angular dependence
IF(energy > 100) THEN
IF(angle .lt. 1.0472) THEN ! angle < 60 deg
comp_yield = delta_zero/COS(angle) ! inverse cosine dependence
ELSE
comp_yield = delta_zero*2
END IF
ELSE
comp_yield = delta_zero
END IF
CASE DEFAULT ! Aluminum
! Value for normal incidence
x = energy/E_m_Al
delta_m = delta_m_Al
delta_zero = delta_m*1.28*x**(-0.67)*(1-EXP(-1.614*x**1.67))
! Angular dependence
IF(energy > 100) THEN
IF(angle .lt. 1.0472) THEN ! angle < 60 deg
comp_yield = delta_zero/COS(angle) ! inverse cosine dependence
ELSE
comp_yield = delta_zero*2
END IF
ELSE
comp_yield = delta_zero
END IF
END SELECT
END FUNCTION comp_yield
!---------------------------------------------------------------------------
! FUNCTION gen_nb_elec(lambda, kmax)
!
!
! DESCRIPTION
!> Gives random values distributed following a Poisson distribution
!> of parameter lambda = secondary electron yield
!--------------------------------------------------------------------------
INTEGER FUNCTION gen_nb_elec(lambda, kmax)
USE random_mod
REAL(KIND = db), DIMENSION(kmax) :: proba, CDF ! probabilities, CDF values
REAL(KIND = db) :: lambda
REAL(KIND = db) :: randnb(1:1)
INTEGER :: j, k, kmax
! Compute probabilities and CDF values
DO k = 1,kmax
proba(k) = exp(-lambda)*lambda**(k-1)/factorial(k-1)
CDF(k) = sum(proba(1:k))
END DO
! Generate the number of electrons according to poisson distribution
call random_array(randnb,1,ran_index(1),ran_array(:,1)) ! uniform distribution in [0,1]
DO j = 1,size(CDF)-1
IF (randnb(1) .lt. CDF(1)) THEN
gen_nb_elec = 0
ELSE IF ((randnb(1) .ge. CDF(j)) .and. (randnb(1) .lt. CDF(j+1))) THEN
gen_nb_elec = j
END IF
END DO
! Below: see Intel oneAPI Math Kernel Library - Fortran
! to optimise running speed when generating random numbers
! TO DO : change if enough time
!-----------------------------------------------------------
!
! INTEGER, INTENT(IN) :: method
! TYPE (VSL_STREAM_STATE), INTENT(IN) :: stream
! INTEGER, INTENT(IN) :: n = 1
! INTEGER, INTENT(IN) :: r
! status = virngpoisson( method, stream, n, r, lambda )
! gen_nb_elec = r
!
!------------------------------------------------------------
END FUNCTION gen_nb_elec
!---------------------------------------------------------------------------
! FUNCTION gen_Eem(kappa, theta, Emax)
!
!
! DESCRIPTION
!> Gives random values distributed following a Gamma distrib.
!> of parameters (kappa, theta) in [0, Emax] eV
!--------------------------------------------------------------------------
FUNCTION gen_Eem(kappa, theta, Emax) RESULT(Eem)
USE random_mod
USE incomplete_gamma
REAL(KIND = db) :: Eem, Emin, Emax
REAL(KIND = db) :: kappa, theta
REAL(KIND = db) :: randnb(1:1)
REAL(KIND = db), DIMENSION(1000) :: energy, CDF, diff
INTEGER :: i, ifault, posid
Emin = 0.00
DO i= 1,size(energy)
energy(i) = (Emin + (i-1)*(Emax-Emin)/(size(energy)-1))
CDF(i) = gamain(energy(i)/theta, kappa, ifault)
END DO
! Generate energy according to Gamma distribution
call random_array(randnb,1,ran_index(1),ran_array(:,1))
diff = abs(randnb(1) - CDF)
posid = minloc(diff, dim = 1)
Eem = energy(posid)
END FUNCTION gen_Eem
!---------------------------------------------------------------------------
! FUNCTION gen_dir(normal_dir)
!
!
! DESCRIPTION
!> Compute a random direction for the emitted electron
!> following a cosine distribution
!--------------------------------------------------------------------------
FUNCTION gen_dir(normal_dir) RESULT(dir)
USE random_mod
REAL(KIND = db) :: randnb1(1:1), randnb2(1:1)
REAL(KIND = db) :: cos_theta, sin_theta, phi
REAL(KIND = db), DIMENSION(3) :: normal_dir, vec_i, vec_j, dir
! Generate an orthonormal basis (i, j, n)
vec_i = (/0, 1, 0/)
vec_j(1) = -1*normal_dir(3)
vec_j(2) = 0.0
vec_j(3) = normal_dir(1)
! Generate a random polar angle in [0,pi/2]
call random_array(randnb1,1,ran_index(1),ran_array(:,1))
cos_theta = randnb1(1)
sin_theta = SIN(ACOS(cos_theta))
! Generate a random azimuthal angle in [0,2pi]
call random_array(randnb2,1,ran_index(1),ran_array(:,1))
phi = 2*pi*randnb2(1)
dir = cos_theta*normal_dir + sin_theta*COS(phi)*vec_i + sin_theta*SIN(phi)*vec_j
END FUNCTION gen_dir
!---------------------------------------------------------------------------
! FUNCTION eval_polynomial(coef, val)
!
!
! DESCRIPTION
!> Evaluate a polynomial at a given point
!> with its coefficients provided in an array
!
!--------------------------------------------------------------------------
REAL(KIND = db) FUNCTION eval_polynomial(coef, val)
REAL(KIND = db), DIMENSION(:) :: coef ! coefficients of the polyn
REAL(KIND = db) :: val ! point where the polyn is evaluated
INTEGER :: ii
eval_polynomial = 0
DO ii=1, size(coef)
eval_polynomial = eval_polynomial + coef(ii)*val**(ii-1)
END DO
END FUNCTION eval_polynomial
!---------------------------------------------------------------------------
! FUNCTION factorial(n)
!
!
! DESCRIPTION
!> Gives the factorial of an integer
!--------------------------------------------------------------------------
INTEGER FUNCTION factorial(n)
INTEGER :: ii
INTEGER :: n
factorial = 1
IF (n .lt. 0) THEN
RETURN
ELSE IF (n .eq. 0) THEN
factorial = 1
ELSE
DO ii=1,n
factorial = factorial*ii
END DO
END IF
END FUNCTION factorial
END MODULE see
c4science · Help