Computes the path of a neutral through the beam_grid
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=Float64), | intent(in), | dimension(3) | :: | rin | Initial position of particle |
|
| real(kind=Float64), | intent(in), | dimension(3) | :: | vin | Initial velocity of particle |
|
| type(ParticleTrack), | intent(inout), | dimension(:) | :: | tracks | Array of ParticleTrack type |
|
| integer(kind=Int32), | intent(out) | :: | ntrack | Number of cells that a particle crosses |
||
| logical, | intent(out), | optional | :: | los_intersect | Indicator whether particle intersects a LOS in spec_chords |
subroutine track(rin, vin, tracks, ntrack, los_intersect)
!+ Computes the path of a neutral through the [[libfida:beam_grid]]
real(Float64), dimension(3), intent(in) :: rin
!+ Initial position of particle
real(Float64), dimension(3), intent(in) :: vin
!+ Initial velocity of particle
type(ParticleTrack), dimension(:), intent(inout) :: tracks
!+ Array of [[ParticleTrack]] type
integer(Int32), intent(out) :: ntrack
!+ Number of cells that a particle crosses
logical, intent(out), optional :: los_intersect
!+ Indicator whether particle intersects a LOS in [[libfida:spec_chords]]
integer :: cc, i, j, ii, mind, ncross, id, k, adaptive, max_cell_splits
integer, dimension(3) :: ind
logical :: in_plasma1, in_plasma2, in_plasma_tmp, los_inter
real(Float64) :: dT, dt1, inv_50, dt2, n_cells, split_tol, inv_param, inv_tol, inv_N, inv_dl, param_sum, param1, param2
real(Float64), dimension(3) :: dt_arr, dr
real(Float64), dimension(3) :: vn, inv_vn, vp
real(Float64), dimension(3) :: ri, ri_tmp, ri_cell
type(LocalEMFields) :: fields
type(LocalProfiles) :: plasma1, plasma2
type(LOSInters) :: inter
real(Float64), dimension(n_stark) :: lambda
integer, dimension(3) :: sgn
integer, dimension(3) :: gdims
vn = vin ; ri = rin ; sgn = 0 ; ntrack = 0
adaptive = 0; max_cell_splits = 1; split_tol = 0.0
los_inter=.False.
if(.not.present(los_intersect)) then
los_inter = .True. !avoids computation if not needed
endif
if(dot_product(vin,vin).eq.0.0) then
return
endif
gdims(1) = beam_grid%nx
gdims(2) = beam_grid%ny
gdims(3) = beam_grid%nz
!! define actual cell
call get_indices(ri,ind)
ri_cell = [beam_grid%xc(ind(1)), &
beam_grid%yc(ind(2)), &
beam_grid%zc(ind(3))]
do i=1,3
if (vn(i).gt.0.0) sgn(i) = 1
if (vn(i).lt.0.0) sgn(i) =-1
if (vn(i).eq.0.0) vn(i) = 1.0d-3
enddo
!! Check adaptive switch, adaptive > 0 activates splitting
adaptive = inputs%adaptive
max_cell_splits = inputs%max_cell_splits
if(adaptive.gt.0) then
split_tol = inputs%split_tol
if(split_tol.eq.0.0) then
adaptive = 0
else
inv_tol = 1.0/split_tol
endif
endif
dr = beam_grid%dr*sgn
inv_vn = 1/vn
inv_50 = 1.0/50.0
cc=1
tracks%time = 0.d0
tracks%flux = 0.d0
ncross = 0
call in_plasma(ri,in_plasma1)
track_loop: do i=1,beam_grid%ntrack
if(cc.gt.beam_grid%ntrack) exit track_loop
dt_arr = abs(( (ri_cell + 0.5*dr) - ri)*inv_vn)
mind = minloc(dt_arr,1)
dT = dt_arr(mind)
ri_tmp = ri + dT*vn
!! Check if velocity intersects LOS and produces wavelength in the right region
inter = spec_chords%inter(ind(1),ind(2),ind(3))
if((.not.los_inter).and.(inter%nchan.ne.0))then
call get_fields(fields, pos=ri_tmp)
chan_loop: do j=1,inter%nchan
id = inter%los_elem(j)%id
vp = ri_tmp - spec_chords%los(id)%lens
call doppler_stark(vp, vn, fields, beam_lambda0, lambda)
los_inter = any((lambda.ge.inputs%lambdamin).and.(lambda.le.inputs%lambdamax))
if(los_inter) exit chan_loop
enddo chan_loop
endif
call in_plasma(ri_tmp,in_plasma2)
if(in_plasma1.neqv.in_plasma2) then
dt1 = 0.0
track_fine: do ii=1,50
dt1 = dt1 + dT*inv_50
ri_tmp = ri + vn*dt1
call in_plasma(ri_tmp,in_plasma_tmp)
if(in_plasma2.eqv.in_plasma_tmp) exit track_fine
enddo track_fine
tracks(cc)%pos = ri + 0.5*dt1*vn
tracks(cc+1)%pos = ri + 0.5*(dt1 + dT)*vn
tracks(cc)%time = dt1
tracks(cc+1)%time = dT - dt1
tracks(cc)%ind = ind
tracks(cc+1)%ind = ind
cc = cc + 2
ncross = ncross + 1
elseif(adaptive.eq.0) then
tracks(cc)%pos = ri + 0.5*dT*vn
tracks(cc)%time = dT
tracks(cc)%ind = ind
cc = cc + 1
elseif(adaptive.gt.0) then !! If splitting is activated, calculate n_cells based on gradient between entrance and exit points
dt2 = 0.0
call get_plasma(plasma1, ri)
call get_plasma(plasma2, ri_tmp)
!! Split cell according to gradient of parameter, parameter is selected according to value of adaptive
select case (adaptive)
case(1)
param1 = plasma1%dene
param2 = plasma2%dene
case(2)
param1 = sum(plasma1%denn(1,:))/size(plasma1%denn(1,:))
param2 = sum(plasma2%denn(1,:))/size(plasma2%denn(1,:))
case(3)
param1 = plasma1%denf
param2 = plasma2%denf
case(4)
param1 = sum(plasma1%deni(:))/size(plasma1%deni(:))
param2 = sum(plasma2%deni(:))/size(plasma2%deni(:))
case(5)
param1 = plasma1%denimp
param2 = plasma2%denimp
case(6)
param1 = plasma1%te
param2 = plasma2%te
case(7)
param1 = plasma1%ti
param2 = plasma2%ti
case default
param1 = plasma1%dene
param2 = plasma2%dene
end select
param_sum = param1 + param2
if(param_sum.le.0.0) then
inv_param = 1.0
else
inv_param = 2.0/param_sum
endif
inv_dl = 1/(dT*norm2(vn))
n_cells = ceiling(abs(param1 - param2)*inv_param*inv_dl*inv_tol)
if(n_cells.gt.max_cell_splits) then
n_cells = max_cell_splits
elseif(n_cells.lt.1.0) then
n_cells = 1.0
endif
inv_N = 1.0/n_cells
split_loop2: do k=1,int(n_cells) !! Split time step and position by n_cells
dt2 = dt2 + dT*inv_N
tracks(cc)%pos = ri + 0.5*dt2*vn
tracks(cc)%time = dT*inv_N
tracks(cc)%ind = ind
cc = cc + 1
enddo split_loop2
else
tracks(cc)%pos = ri + 0.5*dT*vn
tracks(cc)%time = dT
tracks(cc)%ind = ind
cc = cc + 1
endif
in_plasma1 = in_plasma2
ri = ri + dT*vn
ind(mind) = ind(mind) + sgn(mind)
ri_cell(mind) = ri_cell(mind) + dr(mind)
if (ind(mind).gt.gdims(mind)) exit track_loop
if (ind(mind).lt.1) exit track_loop
if (ncross.ge.2) then
cc = cc - 1 !dont include last segment
exit track_loop
endif
enddo track_loop
ntrack = cc-1
if(present(los_intersect)) then
los_intersect = los_inter
endif
end subroutine track