Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (breakdown) Browse

fennecs/src/neutcol_mod.f90
81171805e064breakdown

Last Change
Raw File

neutcol_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: neutcol
!
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!> Module responsible for handling the electron-neutral collisions and creating electrons
!> by ionisation Based on the paper by Birdsall 1991 and Sengupta et al.
!------------------------------------------------------------------------------
module neutcol
USE constants
IMPLICIT NONE
private
LOGICAL, SAVE :: nlcol=.false. !< Flag to activate or not electron neutral collisions
LOGICAL :: nlmaxwellio=.false. !< Flag to define how ionised electrons are created (physically or according to maxwellian)
INTEGER :: itcol = 1 !< number of dt between each evaluation of neutcol_step
Real(kind=db) :: neutdens=2.4e16 !< Neutral particle density in m-3
Real(kind=db) :: neuttemp=300 !< Neutral particle temperature in K
Real(kind=db) :: neutpressure !< Neutral particle pressure in mbar
Real(kind=db) :: scatter_fac = 24.2 !< Energy scattering factor for the considered gas (here for Ne) [eV] see Opal 1971 https://doi.org/10.1063/1.1676707
real(kind=db) :: Eion = 21.56 !< Ionisation energy (eV) (here for Ne)
Real(kind=db) :: E0 = 27.21 !< Atomic unit of energy used for calculation of deviation angles [eV]
Real(kind=db) :: dE = 1 !< resolution for the computation of the cross-sections
Real(kind=db) :: Emax = 1 !< resolution for the computation of the cross-sections
real(kind=db) :: collfactor !< Normalised collision factor (n_n \delta t)
INTEGER :: nb_io_cross=0
Real(kind=db), ALLOCATABLE :: io_cross_sec(:,:) !< Ionisation cross-section table
Real(kind=db), ALLOCATABLE :: io_cross_sec_lin(:,:) !< Ionisation cross-section table for linear interpolation
Real(kind=db), ALLOCATABLE :: io_growth_cross_sec(:) !< Ionisation exponential fitting factor
INTEGER :: nb_ela_cross=0
Real(kind=db), ALLOCATABLE :: ela_cross_sec(:,:) !< Elastic collision cross section table
Real(kind=db), ALLOCATABLE :: ela_cross_sec_lin(:,:) !< Elastic collision cross section table for linear interpolation
Real(kind=db), ALLOCATABLE :: ela_growth_cross_sec(:) !< Elastic collision exponential fitting factor
Real(kind=db) :: Escale !< Energy normalisation factor used to reduce computation costs
CHARACTER(len=128) :: io_cross_sec_file=''
CHARACTER(len=128) :: ela_cross_sec_file=''
Real(kind=db) :: etemp=22000 !< In case of nlmaxwelio, defines the temperature of created electrons [K]
Real(kind=db) :: vth !< In case of nlmaxwelio, defines the normalised thermal velocity of created electrons
LOGICAL :: nldragio=.true. !< Set if inpinging electrons are affected by ionising collisions
INTEGER :: species(2) !< species(1) contains the specie index in plist which stores the colliding particles, species(2) stores the specie index for the released ion.
LOGICAL :: isotropic = .false. !< is the scattering angle isotropic
type coll_process
Real(kind=db):: Emin !< minimum energy in EV for collision to happen
Real(kind=db):: Emax !< maximum energy in EV for linear interpolation
Real(kind=db):: deltaE !< energy step in ev for interpolation
Real(kind=db), allocatable:: cross_sec(:) !< cross-section as read from file
Real(kind=db), allocatable:: E_levels(:) !< energy values as read from file
Real(kind=db), allocatable:: growth_parameter(:) !< collision exponential fitting factor
Real(kind=db), allocatable:: lin_E(:) !< linear energy levels for linear interpolation
Real(kind=db), allocatable:: lin_cross(:) !< linear cross-sections for linear interpolation
end type
NAMELIST /neutcolparams/ neutdens, Eion, &
& scatter_fac, nlcol, io_cross_sec_file, ela_cross_sec_file, nlmaxwellio, etemp, &
& nldragio, itcol, species, isotropic
PUBLIC:: neutcol_init, neutcol_step, neutcol_diag, itcol, neutdens
PROCEDURE(rotate_vel), POINTER:: change_dir => NULL()!< Function evaluating the weight for Dirichelt boundary conditions
ABSTRACT INTERFACE
SUBROUTINE rotate_vel(Ur, Uthet, Uz, coschi, thet)
use constants
real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet
END SUBROUTINE
end interface
CONTAINS
subroutine neutcol_init(lu_in, p)
use mpi
Use basic, only: mpirank, dt, nlclassical,rnorm, vnorm
Use beam, only: particles
Use constants
implicit none
INTEGER, INTENT(IN) :: lu_in
TYPE(particles) :: p
INTEGER:: ierr, istat, i
character(len=1000) :: line
real(kind=db):: xsi
species(1)=1
species(2)=-1
Rewind(lu_in)
READ(lu_in, neutcolparams, iostat=istat)
if (istat.gt.0) then
backspace(lu_in)
read(lu_in,fmt='(A)') line
write(*,'(A)') &
'Invalid line in neutcolparams: '//trim(line)
call MPI_Abort(MPI_COMM_WORLD, -1, ierr)
stop
end if
IF(mpirank .eq. 0) THEN
WRITE(*, neutcolparams)
END IF
if(.not. nlcol) return
if(nlclassical)THEN
Escale=0.5*p%m/elchar*vlight**2
else
Escale=p%m*vlight**2/elchar
end if
if (nlmaxwellio) vth=sqrt(kb*etemp/p%m)/vnorm
if(io_cross_sec_file .ne.'') then
call read_cross_sec(io_cross_sec_file,io_cross_sec, nb_io_cross)
if(nb_io_cross .gt. 0) then
allocate(io_growth_cross_sec(nb_io_cross-1))
! Normalisations
io_cross_sec(:,2)=io_cross_sec(:,2)/rnorm**2
! Precomputing of exponential fitting factor for faster execution
io_growth_cross_sec=log(io_cross_sec(2:nb_io_cross,2)/io_cross_sec(1:nb_io_cross-1,2))/ &
& log(io_cross_sec(2:nb_io_cross,1)/io_cross_sec(1:nb_io_cross-1,1))
end if
end if
if(ela_cross_sec_file .ne.'') then
call read_cross_sec(ela_cross_sec_file,ela_cross_sec, nb_ela_cross)
if(nb_ela_cross .gt. 0) then
allocate(ela_growth_cross_sec(nb_ela_cross-1))
! Normalisations
ela_cross_sec(:,2)=ela_cross_sec(:,2)/rnorm**2
if(.not. isotropic) then
do i=1,nb_ela_cross
xsi=ela_cross_sec(i,1)/(0.25*E0+ela_cross_sec(i,1))
ela_cross_sec(i,2)=ela_cross_sec(i,2)*(2*xsi**2)/((1-xsi)*((1+xsi)*log((1+xsi)/(1-xsi))-2*xsi))
end do
end if
! Precomputing of exponential fitting factor for faster execution
ela_growth_cross_sec=log(ela_cross_sec(2:nb_ela_cross,2)/ela_cross_sec(1:nb_ela_cross-1,2))/ &
& log(ela_cross_sec(2:nb_ela_cross,1)/ela_cross_sec(1:nb_ela_cross-1,1))
end if
END IF
nlcol=nlcol .and. (allocated(io_cross_sec) .or. allocated(ela_cross_sec))
! Collision factor depending on neutral gas parameters
collfactor=neutdens*dt*rnorm**3*itcol
neutpressure=neutdens*kb*300/100
if (.not. isotropic)then
change_dir=> rotate
else
change_dir=> scatter
end if
end subroutine neutcol_init
Subroutine neutcol_diag(File_handle, str, vnorm)
use mpi
Use futils
Integer:: File_handle
Real(kind=db):: vnorm
Character(len=*):: str
CHARACTER(len=256):: grpname
Integer:: ierr, mpirank,i
Real(kind=db)::xsi
Real(kind=db),allocatable:: nonisotropic_rescale(:,:)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr)
IF(mpirank .eq. 0 .and. nlcol) THEN
Write(grpname,'(a,a)') trim(str),"/neutcol"
If(.not. isgroup(File_handle, trim(grpname))) THEN
CALL creatg(File_handle, trim(grpname))
END IF
Call attach(File_handle, trim(grpname), "neutdens", neutdens)
Call attach(File_handle, trim(grpname), "neuttemp", neuttemp)
Call attach(File_handle, trim(grpname), "neutpressure", neutpressure)
Call attach(File_handle, trim(grpname), "scatter_fac", scatter_fac)
Call attach(File_handle, trim(grpname), "Eion", Eion)
Call attach(File_handle, trim(grpname), "E0", E0)
Call attach(File_handle, trim(grpname), "Escale", Escale)
Call putarr(File_handle,trim(grpname)//"species", species)
if (allocated(io_cross_sec)) Call putarr(File_handle, trim(grpname)//"/io_cross_sec", io_cross_sec)
allocate(nonisotropic_rescale(nb_ela_cross,2))
nonisotropic_rescale=1
if(.not. isotropic) then
do i=1,nb_ela_cross
xsi=ela_cross_sec(i,1)/(0.25*E0+ela_cross_sec(i,1))
nonisotropic_rescale(i,2)=(2*xsi**2)/((1-xsi)*((1+xsi)*log((1+xsi)/(1-xsi))-2*xsi))
end do
end if
if (allocated(ela_cross_sec)) Call putarr(File_handle, trim(grpname)//"/ela_cross_sec", ela_cross_sec/nonisotropic_rescale)
END IF
End subroutine neutcol_diag
!-------------------------------------------------------------
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Simulates the elastic and ionising collisions for each particles in plist(species(1))
!
!> @param [inout] plist list of particle species considered in the code
!---------------------------------------------------------------------------
SUBROUTINE neutcol_step(plist)
!
USE random_mod
USE beam
USE omp_lib
USE basic, ONLY: nlclassical
USE distrib, ONLY: lodgaus
type(particles), TARGET::plist(:)
type(particles),pointer::p
INTEGER:: i, omp_thread, num_threads, j, nbcolls_ela, nbcolls_io, iend
real(kind=db):: Rand(5), rand_coll(0:24)
real(kind=db):: v2, v, ek, Everif, es, cosChi, thet, sig_io, sig_ela, vfact, xsi
type(linked_part_row):: ins_p
type(linked_part), POINTER:: created
real(kind=db):: collisionfact,nucol(3),vinit(3),vend(3)
p=>plist(species(1))
if(.not. nlcol .or. p%nploc .le. 0) return
num_threads=omp_get_max_threads()
nbcolls_ela=0
nbcolls_io=0
nucol=0
!!$OMP private(collisionfact,i,omp_thread,Rand,v2,ek,sig_io,sig_ela,es,coschi,thet,vfact, created, v, everif,xsi,vinit,vend)!, reduction(+:nbcolls_ela,nbcolls_io, nucol)
omp_thread=omp_get_thread_num()+1
!omp_thread=1
allocate(ins_p%start)
ins_p%n=0
created=>ins_p%start
!$OMP DO
DO j=1,p%Nploc,25
!for each particle
iend=min(p%Nploc,j+24)
CALL random_array(Rand_coll,iend-j+1,ran_index(omp_thread),ran_array(:,omp_thread))
! we calculate the kinetic energy and norm of the velocity
Do i=j,iend
v2=(p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2)
if(nlclassical) THEN
ek=v2*escale
v=sqrt(v2)
!vinit=p%U(:,i) ! (/p%UR(i),p%UTHET(i),p%UZ(i)/)
ELSE
ek=(p%gamma(i)-1)*escale
v=sqrt(v2)/p%gamma(i)
!vinit=p%U(:,i)/p%gamma(i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i)
end if
sig_io=0
sig_ela=0
! computes the ionisation and elastic collision cross-sections at this kinetic energy
! The ionisation event can only occur if the incoming electron energy is above the binding energy
if (ek .gt. Eion .and. nb_io_cross .gt. 1) then
sig_io=sig_fit(io_cross_sec,io_growth_cross_sec,ek, nb_io_cross)
end if
if (nb_ela_cross .gt. 1) then
sig_ela=sig_fit(ela_cross_sec,ela_growth_cross_sec,ek, nb_ela_cross)
end if
collisionfact=1-exp(-collfactor*(sig_io+sig_ela)*v)
! If we have a collision event
if (Rand_coll(i-j) .lt.collisionfact) THEN
CALL random_array(Rand,1,ran_index(omp_thread),ran_array(:,omp_thread))
! Check if elastic or ionising event is happening
IF(Rand(1).gt. sig_ela/(sig_io+sig_ela)) THEN ! An ionisation collision happened and we create the necessary electron
call io_coll(p,i,nbcolls_io,nucol,ins_p,created)
ELSE ! An elastic collision event happens
call elastic_coll(p,i,nbcolls_ela,nucol)
END IF
END IF
end do
END DO
!$OMP END DO NOWAIT
! clean up the memory after the loop
if(associated(created%prev)) then
created=>created%prev
ins_p%end=>created
deallocate(created%next)
else
deallocate(ins_p%start)
end if
! We collect all created particules into one linked list for easier insertion in plist
!Do i=1,num_threads
!$OMP CRITICAL (insertions)
if(species(2).gt.0.and.ins_p%n .gt.0) then
CALL add_created_part(plist(species(2)), ins_p, .false.,.true.)
end if
!$OMP END CRITICAL (insertions)
!$OMP CRITICAL(addcolls)
p%nbcolls=p%nbcolls+(/nbcolls_io, nbcolls_ela/)
p%nudcol=p%nudcol+nucol
!$OMP END CRITICAL(addcolls)
!$OMP BARRIER
!$OMP CRITICAL (insertelectrons)
if(ins_p%n .gt.0) then
CALL add_created_part(plist(species(1)),ins_p,.true.,.false.)
!exit
end if
!$OMP END CRITICAL (insertelectrons)
!end do
!Write(*,*)"mpirank: ", mpirank, " Nb colls ela, io: ",nbcolls_ela, nbcolls_io
!
END SUBROUTINE neutcol_step
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Simulates the elastic collision for particle p(i)
!
!> @param [inout] p structure of particle species considered in the code
!---------------------------------------------------------------------------
subroutine elastic_coll(p, i, nbcolls_ela, nucol)
use omp_lib
USE random_mod
USE beam
USE basic, ONLY: nlclassical
IMPLICIT NONE
type(particles), pointer::p
INTEGER:: i, omp_thread, nbcolls_ela
real(kind=db):: Rand(2)
real(kind=db):: ek, es, cosChi, thet, vfact, v, v2
real(kind=db):: nucol(3),vinit(3),vend(3)
! we calculate the kinetic energy and norm of the velocity
v2=(p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2)
if(nlclassical) THEN
ek=v2*escale
v=sqrt(v2)
vinit=p%U(:,i) ! (/p%UR(i),p%UTHET(i),p%UZ(i)/)
ELSE
ek=(p%gamma(i)-1)*escale
v=sqrt(v2)/p%gamma(i)
vinit=p%U(:,i)/p%gamma(i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i)
end if
omp_thread=omp_get_thread_num()+1
CALL random_array(Rand,2,ran_index(omp_thread),ran_array(:,omp_thread))
Es=Ek
vfact=1.0
nbcolls_ela=nbcolls_ela+1
! We calculate the scattered velocity angle for the scattered electron
if (isotropic) then
coschi=cos(Rand(1)*pi)
else
cosChi=1-2*Rand(1)/(1+8*Es/E0*(1-Rand(1)))
end if
thet=Rand(2)*2*pi
! Change the incident electron velocity direction and amplitude if necessary
p%U(:,i)=p%U(:,i)*vfact
call change_dir(p%U(1,i),p%U(2,i), p%U(3,i), coschi, thet)
if(nlclassical) THEN
vend=p%U(:,i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)
ELSE
p%gamma(i)=sqrt(1+p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2)
vend=p%U(:,i)/p%gamma(i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i)
END IF
nucol=nucol+1-vend/vinit
end subroutine
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief Simulates the ionisation collision for particle p(i)
!
!> @param [inout] p structure of particle species considered in the code
!---------------------------------------------------------------------------
subroutine io_coll(p, i, nbcolls_io, nucol, ins_p, created)
use omp_lib
USE random_mod
USE beam
USE basic, ONLY: nlclassical
USE distrib, ONLY: lodgaus
IMPLICIT NONE
type(particles), pointer::p
INTEGER:: i, omp_thread, nbcolls_io
real(kind=db):: Rand(3)
real(kind=db):: v, v2, ek, Everif, es, cosChi, thet, sig_io, sig_ela, vfact, xsi
real(kind=db):: collisionfact,nucol(3),vinit(3),vend(3)
type(linked_part_row):: ins_p
type(linked_part), POINTER:: created
! we calculate the kinetic energy and norm of the velocity
v2=(p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2)
if(nlclassical) THEN
ek=v2*escale
v=sqrt(v2)
vinit=p%U(:,i) ! (/p%UR(i),p%UTHET(i),p%UZ(i)/)
ELSE
ek=(p%gamma(i)-1)*escale
v=sqrt(v2)/p%gamma(i)
vinit=p%U(:,i)/p%gamma(i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i)
end if
omp_thread=omp_get_thread_num()+1
! prepare the memory for the released electron
ins_p%n=ins_p%n+1
allocate(created%next)
created%next%prev=>created
! Fill created particle new position
created%p%pos=p%pos(:,i)!(/p%R(i), p%THET(i), p%Z(i)/)
CALL random_array(Rand,3,ran_index(omp_thread),ran_array(:,omp_thread))
IF( nlmaxwellio ) THEN ! the new electron velocity is defined according to a Maxwellian
CALL lodgaus(0, Rand(1:3))
! get random velocity
created%p%U=vth*Rand(1:3)
ELSE
! Compute created electron energy
Es=scatter_fac*tan(Rand(1)*atan((Ek-Eion)/(2*scatter_fac)))
! Compute scattering angles for created electron
if (isotropic) then
coschi=cos(Rand(2)*pi)
else
cosChi=1-2*Rand(2)/(1+8*Es/E0*(1-Rand(2)))
end if
thet=Rand(3)*2*pi
if(nlclassical)THEN
! new velocity factor for created particle
vfact=sqrt(Es/Ek)
ELSE
! new velocity factor for created particle
vfact=sqrt(Es*(Es+2*Escale)/(Ek*(Ek+2*Escale)))
END IF
! Fill created particle velocity
created%p%U=vfact*p%U(:,i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)
! rotate the velocity vector due to the collision
call change_dir(created%p%U(1),created%p%U(2), created%p%U(3), coschi, thet)
END IF
vend=created%p%U
if(nlclassical)THEN
! Lorentz factor for created particle
created%p%gamma=1.0
ELSE
! Lorentz factor for created particle
created%p%gamma=sqrt(1+created%p%U(1)**2+created%p%U(2)**2+created%p%U(3)**2)
vend=vend/created%p%gamma
END IF
! We prepare the next created particle
ins_p%end=>created
created=>created%next
! We keep track of what changed
nbcolls_io=nbcolls_io+1
nucol=nucol-vend/vinit
! If we want the incoming electron to be scattered, we need to compute
! its new kinetic energy
if (nldragio) THEN
! We store the lossed energy in pot for keeping track of energy conservation
created%prev%p%pot=Eion+Es
CALL random_array(Rand,2,ran_index(omp_thread),ran_array(:,omp_thread))
Es=Ek-Eion-Es
if(nlclassical)THEN
! new velocity factor for scattered particle
vfact=sqrt(Es/Ek)
ELSE
! new velocity factor for scattered particle
vfact=sqrt(Es*(Es+2*Escale)/(Ek*(Ek+2*Escale)))
END IF
ELSE
return
END IF
! We calculate the scattered velocity angle for the scattered electron
if (isotropic) then
coschi=cos(Rand(1)*pi)
else
cosChi=1-2*Rand(1)/(1+8*Es/E0*(1-Rand(1)))
end if
thet=Rand(2)*2*pi
! Change the incident electron velocity direction and amplitude if necessary
p%U(:,i)=p%U(:,i)*vfact
!p%UTHET(i)=p%UTHET(i)*vfact
!p%UZ(i)=p%UZ(i)*vfact
call change_dir(p%U(1,i),p%U(2,i), p%U(3,i), coschi, thet)
if(nlclassical) THEN
vend=p%U(:,i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)
ELSE
p%gamma(i)=sqrt(1+p%U(1,i)**2+p%U(2,i)**2+p%U(3,i)**2)
vend=p%U(:,i)/p%gamma(i)!(/p%UR(i),p%UTHET(i),p%UZ(i)/)/p%gamma(i)
END IF
nucol=nucol+1-vend/vinit
end subroutine
FUNCTION sig_fit(sig_vec,growth_vec,ek,nb_cross)
use distrib, ONLY: closest
real(kind=db)::sig_fit, ek
real(kind=db):: sig_vec(:,:), growth_vec(:)
Integer:: k, nb_cross
sig_fit=0
k=closest(sig_vec(:,1),ek, nb_cross-1)
if(k.lt.1) return
!sig_fit=(sig_vec(k,1)-sig_vec(k-1,1))/(sig_vec(k,2)-sig_vec(k-1,2))*(sig_vec(k,2)-ek)+sig_vec(k-1,1)
! Exponential fitting relevant at high energies
!sig_fit=sig_vec(k,2)*(ek/sig_vec(k,1))**growth_vec(k)
sig_fit=sig_vec(k,2)*exp(log(ek/sig_vec(k,1))*growth_vec(k))
END FUNCTION sig_fit
FUNCTION sig_fit_lin(sig_vec,growth_vec,ek,nb_cross)
use distrib, ONLY: closest
real(kind=db)::sig_fit_lin, ek
real(kind=db):: sig_vec(:,:), growth_vec(:)
Integer:: k, nb_cross
sig_fit_lin=0
k=(k-sig_vec(k,1))/dE
if(k.lt.1) return
!sig_fit=(sig_vec(k,1)-sig_vec(k-1,1))/(sig_vec(k,2)-sig_vec(k-1,2))*(sig_vec(k,2)-ek)+sig_vec(k-1,1)
! Exponential fitting relevant at high energies
sig_fit_lin=sig_vec(k,2)+(ek-sig_vec(k,1))/(sig_vec(k+1,1)-sig_vec(k,1))*(sig_vec(k+1,2)-sig_vec(k,2))
END FUNCTION sig_fit_lin
SUBROUTINE rotate(Ur, Uthet, Uz, coschi, thet)
real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet
real(kind=db):: norm, perp(3), U(3), U0(3)
real(kind=db):: sinchi, sinthet, costhet
Integer :: iperp1,iperp2
U0=(/Ur,Uthet,Uz/)
norm=sqrt(sum(U0**2))
U=U0/norm
! Find a vector perpendicular to U for chi rotation
! find the direction with maximum amplitude
perp=(/1,1,1/)
iperp1=maxloc(abs(U),1)
! find second direction with next max amplitude
perp(iperp1)=0
iperp2=maxloc(abs(perp*U),1)
perp=0
perp(iperp2)=U(iperp1)
perp(iperp1)=-U(iperp2)
! Normalise the rotation vector
perp=perp/sqrt(sum(perp**2))
! Compute sinus and cosinus for rotation
sinchi=sqrt(1-coschi**2)
costhet=cos(thet)
sinthet=sin(thet)
! Rotation of angle chi around perp
Ur = (coschi+perp(1)**2*(1-coschi))*U0(1) + (perp(1)*perp(2)*(1-coschi)-perp(3)*sinchi)*U0(2) + (perp(1)*perp(3)*(1-coschi) + perp(2)*sinchi)*U0(3)
Uthet = (perp(1)*perp(2)*(1-coschi)+perp(3)*sinchi)*U0(1) + (coschi + perp(2)**2*(1-coschi))*U0(2) +(perp(2)*perp(3)*(1-coschi)-perp(1)*sinchi)*U0(3)
Uz = (perp(1)*perp(3)*(1-coschi)-perp(2)*sinchi)*U0(1) +(perp(3)*perp(2)*(1-coschi)+perp(1)*sinchi)*U0(2) +( coschi+perp(3)**2*(1-coschi))*U0(3)
U0 =(/Ur,Uthet,Uz/)
! second rotation according to uniform distribution
! Rotation of angle theta around U
Ur = (costhet+U(1)**2*(1-costhet))*U0(1) + (U(1)*U(2)*(1-costhet) - U(3)*sinthet)*U0(2) + (U(1)*U(3)*(1-costhet)+U(2)*sinthet)*U0(3)
Uthet = (U(2)*U(1)*(1-costhet)+U(3)*sinthet)*U0(1) + (costhet + U(2)**2*(1-costhet))*U0(2) + (U(2)*U(3)*(1-costhet)-U(1)*sinthet)*U0(3)
Uz = (U(3)*U(1)*(1-costhet) - U(2)*sinthet)*U0(1) + (U(3)*U(2)*(1-costhet)+U(1)*sinthet)*U0(2) + (costhet +U(3)**2*(1-costhet))*U0(3)
!normf=sqrt(Ur**2+Uthet**2+Uz**2)
!if(abs(norm-normf)/norm .gt. 1e-14) WRITE(*,*) "Error in rotate the norm of v changed"
END SUBROUTINE rotate
SUBROUTINE scatter(Ur, Uthet, Uz, coschi, thet)
real(kind=db), INTENT(INOUT):: Ur, uthet, uz, coschi, thet
real(kind=db):: norm
real(kind=db):: sinchi, sinthet, costhet
norm=sqrt(Ur**2+Uz**2+Uthet**2)
! Compute sinus and cosinus for rotation
sinchi=sqrt(1-coschi**2)
costhet=cos(thet)
sinthet=sin(thet)
Ur=norm*sinchi*costhet
Uthet=norm*sinchi*sinthet
Uz=norm*coschi
END SUBROUTINE scatter
SUBROUTINE read_cross_sec(filename,cross_sec, nb_cross)
CHARACTER(len=*) ::filename
Real(kind=db), ALLOCATABLE :: cross_sec(:,:)
INTEGER:: nb_cross
INTEGER :: lu_cross_sec=9999
INTEGER:: i, openerr, reason
CHARACTER(len=256) :: header
real(kind=db):: t1,t2
nb_cross=0
OPEN(UNIT=lu_cross_sec,FILE=trim(filename),ACTION='READ',IOSTAT=openerr)
header=' '
IF(openerr .ne. 0) THEN
CLOSE(unit=lu_cross_sec)
RETURN
END IF
! The cross section table is defined as a two column energy and cross_section
DO WHILE(.true.)
READ(lu_cross_sec,'(a)',IOSTAT=reason) header
header=adjustl(header)
if(reason .lt. 0 ) exit ! We reached end of file
if( header(1:1) .ne. '!') then
READ(header,*) t1, t2
if(t1 .ne. 0 .and. t2.ne. 0) nb_cross=nb_cross+1
end if
END DO
if (allocated(cross_sec)) deallocate(cross_sec)
allocate(cross_sec(nb_cross,2))
REWIND(lu_cross_sec)
! The cross section table is defined as a two column energy and cross_section
i=1
DO WHILE(i .le. nb_cross)
READ(lu_cross_sec,'(a)',IOSTAT=reason) header
header=adjustl(header)
if(reason .lt. 0 ) exit ! We reached end of file
if( header(1:1) .ne. '!') then
READ(header,*) cross_sec(i,1), cross_sec(i,2)
if(cross_sec(i,1) .ne. 0 .and. cross_sec(i,2).ne. 0) i=i+1
end if
END DO
CLOSE(unit=lu_cross_sec)
END subroutine read_cross_sec
subroutine init_cross_sec_type(cross_sec_data, nb, cross_sec_array, exp_growth_array, Emin, Emax, deltaE)
type(coll_process):: cross_sec_data
INTEGER:: nb
real(kind=db):: cross_sec_array(:,:), exp_growth_array(:)
real(kind=db):: Emin, Emax, deltaE
allocate(cross_sec_data%cross_sec(nb),cross_sec_data%E_levels(nb))
cross_sec_data%cross_sec=cross_sec_array(:,2)
cross_sec_data%E_levels=cross_sec_array(:,1)
allocate(cross_sec_data%growth_parameter(nb-1))
cross_sec_data%growth_parameter=exp_growth_array
cross_sec_data%Emin=Emin
cross_sec_data%Emax=Emax
end subroutine
end module neutcol
c4science · Help