!------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: IIEE ! !> @author !> S. Guinchard - EPFL/SPC !> P. Giroud-Garampon - EPFL/SPC ! !> Last modif. !> 16.02.2024 ! ! DESCRIPTION: !> Module handling ion induced electron emissions (IIEE) !------------------------------------------------------------------------------ MODULE iiee USE particletypes USE constants USE basic USE materials !#include "mkl_vsl.f90" ! for random # generators using MKL intel library IMPLICIT NONE CONTAINS !--------------------------------------------------------------------------- ! SUBROUTINE ion_induced(pion, losthole, pelec, nblostparts) ! ! ! DESCRIPTION !> function to determine the number of electrons !> to add to a given species as a function of the number of lost ions !-------------------------------------------------------------------------- SUBROUTINE ion_induced(pion, losthole, pelec, nblostparts) USE geometry TYPE(particles), INTENT(INOUT):: pion, pelec ! ions and electrons lists REAL(KIND = db), DIMENSION(3) :: last_pos ! last position for lost ion REAL(KIND = db), DIMENSION(3) :: normal_dir, dir ! unit vector for directions INTEGER, DIMENSION(pion%Nploc):: losthole ! indices of lost ions INTEGER :: nblostparts, Nploc, Nploc_old ! number of lost ions INTEGER :: parts_size_increase, nbadded ! number of added electrons INTEGER :: i,j INTEGER :: neuttype_id, material_id INTEGER :: kmax ! max number of electrons generated REAL(KIND = db) :: Einc, thetainc ! angle of incidence of the impinging ion 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.4621 theta = 4.7778 Emax = 30 ! Max value for emitted electron energy neuttype_id = pion%neuttype_id material_id = pion%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 = compute_Ekin(pion%U(:,losthole(i)), pion%m) ! [MeV] last_pos = revert_push(pion, losthole(i)) normal_dir = find_normal(last_pos) thetainc = incidence(pion%U(:,losthole(i)), normal_dir) IF (pion%neuttype_id .eq. 1) THEN ! Yield for H2 lambda = 2.0*compute_yield(Einc/2.0, thetainc, neuttype_id, material_id) ! computed from yield by protons ELSE ! Yield for other neutral lambda = compute_yield(Einc, thetainc, neuttype_id, material_id) END IF nbadded = gen_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. .true.) THEN pelec%U(1,Nploc_old+j) = 0.0 pelec%U(2,Nploc_old+j) = 0.0 pelec%U(3,Nploc_old+j) = 0.0 ELSE Eem = gen_Eem(kappa, theta, Emax) ! generate an energy value following gamma distribution dir = gen_dir(normal_dir) ! generate a random direction following a cosine distribution Vem = compute_V(Eem, pelec%m) pelec%U(1,Nploc_old+j) = Vem * normal_dir(1) pelec%U(2,Nploc_old+j) = Vem * normal_dir(2) pelec%U(3,Nploc_old+j) = Vem * normal_dir(3) END IF END DO END DO END SUBROUTINE ion_induced !--------------------------------------------------------------------------- ! FUNCTION compute_Ekin(velocity, mass) ! ! ! DESCRIPTION !> Computes the kinetic energy in MeV of a particle given its 3 velocity components !-------------------------------------------------------------------------- FUNCTION compute_Ekin(velocity, mass) RESULT(Ekin) REAL(KIND = db), DIMENSION(3) :: velocity REAL(KIND = db) :: mass, Ekin Ekin = 5E-7 * mass * vlight**2 /elchar * (velocity(1)**2 + velocity(2)**2 + velocity(3)**2) END FUNCTION compute_Ekin !--------------------------------------------------------------------------- ! FUNCTION compute_V(Ekin, mass) ! ! ! DESCRIPTION !> Computes the normalized velocity of a particle with energy Ekin !-------------------------------------------------------------------------- FUNCTION compute_V(Ekin, mass) RESULT(V) REAL(KIND = db) :: Ekin, mass, V V = sqrt(2/mass * Ekin * elchar) / vlight END FUNCTION compute_V !--------------------------------------------------------------------------- ! FUNCTION find_normal(last_position) ! ! ! DESCRIPTION !> Computes the normal direction to the electrode surface !-------------------------------------------------------------------------- FUNCTION find_normal(last_position) RESULT(normal_dir) USE geometry REAL(KIND = db), DIMENSION(3) :: last_position ! last position of the ion REAL(KIND = db), DIMENSION(3) :: normal_dir ! normal direction vector REAL(KIND = db), DIMENSION(3) :: weight ! geom. weight at last position REAL(KIND = db) :: norm call geom_weight(last_position(1), last_position(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 ion !-------------------------------------------------------------------------- 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 revert_push(pion, partid) ! ! ! DESCRIPTION !> reverts Buneman algorithm over one time step !> to obtain ion position right before capture !-------------------------------------------------------------------------- FUNCTION revert_push(pion, partid) USE basic, ONLY: dt TYPE(particles), INTENT(INOUT):: pion REAL(KIND=db), DIMENSION(3):: revert_push INTEGER :: partid revert_push(1) = abs(pion%pos(1,partid) - pion%U(1,partid)*dt) revert_push(2) = pion%pos(2,partid) - 1/pion%pos(1,partid)* pion%U(2,partid)*dt revert_push(3) = pion%pos(3,partid) - pion%U(3,partid)*dt END FUNCTION revert_push !--------------------------------------------------------------------------- ! FUNCTION compute_yield(energy, angle, neuttype_id, material_id) ! ! ! DESCRIPTION !> Gives the value of the electron yield !> as a function of the energy of the incident ion, the angle of incidence, !> the type of neutral gas and the type of material !-------------------------------------------------------------------------- REAL(KIND = db) FUNCTION compute_yield(energy, angle, neuttype_id, material_id) REAL(KIND = db) :: energy, angle, gamma_zero INTEGER :: neuttype_id, material_id REAL(KIND = db) :: Lambda_exp Lambda_exp = 1E-3 ! Value for normal incidence SELECT CASE(material_id) CASE(1) !304 stainless steel SELECT CASE(neuttype_id) CASE(1) ! Hydrogen IF(energy.le. 1E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_H_SS_1, energy) ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_SS_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_SS_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_SS_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_SS_5, energy) END IF CASE DEFAULT ! Hydrogen IF(energy.le. 1E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_H_SS_1,energy) ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_SS_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_SS_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_SS_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_SS_5, energy) END IF END SELECT CASE(2) ! Copper SELECT CASE(neuttype_id) CASE(1) ! Hydrogen IF(energy.le. 1E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_H_Cu_1, energy) ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Cu_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Cu_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Cu_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Cu_5, energy) END IF CASE(2) ! Helium IF(energy.le. 3E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_He_Cu_1, energy) ELSE IF(energy.gt. 3E-3 .and. energy.le. 1E-2) THEN gamma_zero = eval_polynomial(iiee_yield_He_Cu_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = eval_polynomial(iiee_yield_He_Cu_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = eval_polynomial(iiee_yield_He_Cu_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = eval_polynomial(iiee_yield_He_Cu_5, energy) END IF CASE(3) ! Neon IF(energy.le. 2E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Cu_1, energy) ELSE IF(energy.gt. 2E-3 .and. energy.le. 1E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Cu_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Cu_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 5E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Cu_4, energy) END IF CASE(4) ! Argon IF(energy.le. 1.5E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Cu_1, energy) ELSE IF(energy.gt. 1.5E-3 .and. energy.le. 1E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Cu_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Cu_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Cu_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Cu_5, energy) END IF CASE DEFAULT ! Hydrogen IF(energy.le. 1E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_H_Cu_1, energy) ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Cu_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Cu_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Cu_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Cu_5, energy) END IF END SELECT CASE(3) ! Alumium SELECT CASE(neuttype_id) CASE(1) ! Hydrogen IF(energy.le. 1E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_H_Al_1, energy) ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Al_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Al_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Al_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Al_5, energy) END IF CASE(2) ! Helium IF(energy.le. 2E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_He_Al_1, energy) ELSE IF(energy.gt. 2E-3 .and. energy.le. 1E-2) THEN gamma_zero = eval_polynomial(iiee_yield_He_Al_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = eval_polynomial(iiee_yield_He_Al_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = eval_polynomial(iiee_yield_He_Al_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = eval_polynomial(iiee_yield_He_Al_5, energy) END IF CASE(3) ! Neon IF(energy.le. 2E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Al_1, energy) ELSE IF(energy.gt. 2E-3 .and. energy.le. 1E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Al_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Al_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Al_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ne_Al_5, energy) END IF CASE(4) ! Argon IF(energy.le. 1E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Al_1, energy) ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Al_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Al_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Al_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = eval_polynomial(iiee_yield_Ar_Al_5, energy) END IF CASE DEFAULT ! Hydrogen IF(energy.le. 1E-3 ) THEN gamma_zero = eval_polynomial(iiee_yield_H_Al_1, energy) ELSE IF(energy.gt. 1E-3 .and. energy.le. 1E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Al_2, energy) ELSE IF(energy.gt. 1E-2 .and. energy.le. 2E-2 ) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Al_3, energy) ELSE IF(energy.gt. 2E-2 .and. energy.le. 3E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Al_4, energy) ELSE IF(energy.gt. 3E-2 .and. energy.le. 5E-2) THEN gamma_zero = Lambda_exp * eval_polynomial(iiee_yield_H_Al_5, energy) END IF END SELECT END SELECT ! Angular dependence IF(energy .le. 1E-3) THEN ! potential electron emission compute_yield = gamma_zero ELSE ! kinetic electron emission IF(angle .lt. 1.3963) THEN ! angle < 80 deg compute_yield = gamma_zero/COS(angle) ! inverse cosine dependence ELSE compute_yield = gamma_zero*5.7588 END IF END IF END FUNCTION compute_yield !--------------------------------------------------------------------------- ! FUNCTION gen_elec(lambda, kmax) ! ! ! DESCRIPTION !> Gives random values distributed following a Poisson distribution !> of parameter lambda = electron yield !-------------------------------------------------------------------------- INTEGER FUNCTION gen_elec(lambda, kmax) USE random_mod REAL(KIND = db), DIMENSION(kmax) :: proba, CDF ! probabilities, CDF values REAL(KIND = db) :: lambda REAL(KIND = db) :: nb_alea(1:1) INTEGER :: j, k, kmax ! Compute probabilities and CDF values DO k = 1,kmax proba(k) = exp(-lambda)*lambda**(k-1)/factorial_fun(k-1); CDF(k) = sum(proba(1:k)); END DO ! Generate the number of electrons according to Poisson distribution call random_array(nb_alea,1,ran_index(1),ran_array(:,1)) ! uniform distribution in [0,1] DO j = 1,size(CDF)-1 IF (nb_alea(1) .lt. CDF(1)) THEN gen_elec = 0 ELSE IF ((CDF(j) .le. nb_alea(1)) .and. (nb_alea(1) .lt. CDF(j+1))) THEN gen_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_elec = r ! !------------------------------------------------------------ END FUNCTION gen_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) :: Emin, Emax, Eem REAL(KIND = db) :: kappa, theta REAL(KIND = db) :: nb_alea(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 distrib. call random_array(nb_alea,1,ran_index(1),ran_array(:,1)) diff = abs(nb_alea(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 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_fun(n) ! ! ! DESCRIPTION !> Gives the factorial of an integer !-------------------------------------------------------------------------- INTEGER FUNCTION factorial_fun(n) INTEGER :: ii INTEGER :: n factorial_fun = 1 IF (n .lt. 0) THEN RETURN ELSE IF (n .eq. 0) THEN factorial_fun = 1 ELSE DO ii=1,n factorial_fun = factorial_fun*ii END DO END IF END FUNCTION factorial_fun END MODULE iiee