Computes the path of a neutral through the pass_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_cylindrical(rin, vin, tracks, ntrack, los_intersect)
!+ Computes the path of a neutral through the [[libfida:pass_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]]
real(Float64), dimension(3,3) :: basis
real(Float64), dimension(3) :: dt_arr, dr !! Need to make this rzphi
real(Float64), dimension(3) :: vn, vn_cyl, vp
real(Float64), dimension(3) :: ri, ri_cyl, ri_tmp
real(Float64), dimension(3) :: p, nz
real(Float64), dimension(3) :: v_plane_cyl, v_plane
real(Float64), dimension(3) :: h_plane_cyl, h_plane
real(Float64), dimension(3) :: arc_cyl, arc
real(Float64), dimension(3) :: redge, tedge
real(Float64), dimension(3) :: redge_cyl, tedge_cyl
integer, dimension(3) :: sgn
integer, dimension(3) :: gdims
integer, dimension(1) :: minpos
integer, dimension(3) :: ind
real(Float64) :: dT, dt1, inv_50, t
real(Float64) :: s, c, phi
integer :: cc, i, j, ii, mind, ncross, id
type(LocalEMFields) :: fields
type(LOSInters) :: inter
real(Float64), dimension(n_stark) :: lambda
logical :: in_plasma1, in_plasma2, in_plasma_tmp, los_inter
integer :: ir, iz, iphi
vn = vin ; ri = rin ; sgn = 0 ; ntrack = 0
los_inter=.False.
if(.not.present(los_intersect)) then
los_inter = .True. !avoids computation if not needed
endif
if(dot_product(vn,vn).eq.0.0) then
return
endif
gdims(1) = pass_grid%nr
gdims(2) = pass_grid%nz
gdims(3) = pass_grid%nphi
phi = atan2(rin(2),rin(1))
s = sin(phi) ; c = cos(phi)
vn_cyl(1) = c*vn(1) + s*vn(2)
vn_cyl(3) = -s*vn(1) + c*vn(2)
vn_cyl(2) = vn(3)
do i=1,3
!! sgn is in R-Z-Phi coordinates
if (vn_cyl(i).gt.0.d0) then
sgn(i) = 1
else if (vn_cyl(i).lt.0.d0) then
sgn(i) =-1
end if
enddo
dr(1) = pass_grid%dr*sgn(1)
dr(2) = pass_grid%dz*sgn(2)
dr(3) = pass_grid%dphi*sgn(3)
!! Define actual cell
ri_cyl(1) = sqrt(ri(1)*ri(1) + ri(2)*ri(2))
ri_cyl(2) = ri(3)
ri_cyl(3) = atan2(ri(2), ri(1))
call get_passive_grid_indices(ri_cyl,ind)
arc_cyl(1) = pass_grid%r(ind(1))
arc_cyl(2) = pass_grid%z(ind(2))
arc_cyl(3) = pass_grid%phi(ind(3))
h_plane_cyl = arc_cyl
v_plane_cyl = arc_cyl
!! Define surfaces to intersect
if(sgn(1).gt.0.d0) arc_cyl(1) = pass_grid%r(ind(1)+1)
if(sgn(2).gt.0.d0) h_plane_cyl(2) = pass_grid%z(ind(2)+1)
if(sgn(3).gt.0.d0) v_plane_cyl(3) = pass_grid%phi(ind(3)+1)
! Special case of the particle being on the surace handled below
if((sgn(1).lt.0.d0).and.(arc_cyl(1).eq.ri_cyl(1))) then
arc_cyl(1) = pass_grid%r(ind(1)-1)
h_plane_cyl(1) = pass_grid%r(ind(1)-1)
v_plane_cyl(1) = pass_grid%r(ind(1)-1)
endif
if((sgn(2).lt.0.d0).and.(h_plane_cyl(2).eq.ri_cyl(2))) then
arc_cyl(2) = pass_grid%z(ind(2)-1)
h_plane_cyl(2) = pass_grid%z(ind(2)-1)
v_plane_cyl(2) = pass_grid%z(ind(2)-1)
endif
if((sgn(3).lt.0.d0).and.(v_plane_cyl(3).eq.ri_cyl(3))) then
arc_cyl(3) = pass_grid%phi(ind(3)-1)
h_plane_cyl(3) = pass_grid%phi(ind(3)-1)
v_plane_cyl(3) = pass_grid%phi(ind(3)-1)
endif
call cyl_to_uvw(arc_cyl,arc)
call cyl_to_uvw(h_plane_cyl,h_plane)
call cyl_to_uvw(v_plane_cyl,v_plane)
!! Normal vectors
nz(1) = 0.d0 ; nz(2) = 0.d0 ; nz(3) = 1.d0
redge_cyl(1) = v_plane_cyl(1) + pass_grid%dr
redge_cyl(2) = v_plane_cyl(2)
redge_cyl(3) = v_plane_cyl(3)
call cyl_to_uvw(redge_cyl,redge)
tedge_cyl(1) = v_plane_cyl(1)
tedge_cyl(2) = v_plane_cyl(2) + pass_grid%dz
tedge_cyl(3) = v_plane_cyl(3)
call cyl_to_uvw(tedge_cyl,tedge)
call plane_basis(v_plane, redge, tedge, basis)
!! Track the particle
inv_50 = 1.0/50.0
cc=1
tracks%time = 0.d0
tracks%flux = 0.d0
ncross = 0
call in_plasma(ri,in_plasma1,input_coords=1)
track_loop: do i=1,pass_grid%ntrack
if(cc.gt.pass_grid%ntrack) exit track_loop
call line_cylinder_intersect(ri, vn, arc, p, dt_arr(1))
call line_plane_intersect(ri, vn, h_plane, nz, p, dt_arr(2))
call line_plane_intersect(ri, vn, v_plane, -basis(:,3), p, dt_arr(3))
minpos = minloc(dt_arr, mask=dt_arr.gt.0.d0)
mind = minpos(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%cyl_inter(ind(1),ind(2),ind(3))
if((.not.los_inter).and.(inter%nchan.ne.0))then
call get_fields(fields, pos=ri_tmp, input_coords=1)
chan_loop: do j=1,inter%nchan
id = inter%los_elem(j)%id
vp = ri_tmp - spec_chords%los(id)%lens_uvw
call doppler_stark(vp, vn, fields, 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,input_coords=1)
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,input_coords=1)
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
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)
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
!! Particle advancement and basis update
arc_cyl(mind) = arc_cyl(mind) + dr(mind)
h_plane_cyl(mind) = h_plane_cyl(mind) + dr(mind)
v_plane_cyl(mind) = v_plane_cyl(mind) + dr(mind)
call cyl_to_uvw(arc_cyl,arc)
call cyl_to_uvw(h_plane_cyl,h_plane)
call cyl_to_uvw(v_plane_cyl,v_plane)
redge_cyl(1) = v_plane_cyl(1) + pass_grid%dr
redge_cyl(2) = v_plane_cyl(2)
redge_cyl(3) = v_plane_cyl(3)
call cyl_to_uvw(redge_cyl,redge)
tedge_cyl(1) = v_plane_cyl(1)
tedge_cyl(2) = v_plane_cyl(2) + pass_grid%dz
tedge_cyl(3) = v_plane_cyl(3)
call cyl_to_uvw(tedge_cyl,tedge)
call plane_basis(v_plane, redge, tedge, basis)
enddo track_loop
ntrack = cc-1
if(present(los_intersect)) then
los_intersect = los_inter
endif
end subroutine track_cylindrical