FIXME : Add documentation
subroutine xigridset
use theta_grid, only: ntgrid, bmag
implicit none
integer :: ig, je, ixi, il
if (.not. allocated(xi)) allocate (xi(2*nlambda, -ntgrid:ntgrid))
if (.not. allocated(ixi_to_il)) allocate (ixi_to_il(2*nlambda, -ntgrid:ntgrid))
if (.not. allocated(ixi_to_isgn)) allocate (ixi_to_isgn(2*nlambda, -ntgrid:ntgrid))
if (.not. allocated(wxi)) allocate (wxi(2*nlambda, -ntgrid:ntgrid))
! define array 'sgn' that returns sign of vpa associated with isgn=1,2
sgn(1) = 1
sgn(2) = -1
! xi is vpa / v and goes from 1 --> -1
xi = 0.0
do ig = -ntgrid, ntgrid
xi(:jend(ig), ig) = sqrt(max(1.0 - al(:jend(ig))*spread(bmag(ig),1,jend(ig)),0.0))
xi(jend(ig)+1:2*jend(ig)-1, ig) = -xi(jend(ig)-1:1:-1, ig)
end do
! get mapping from ixi (which runs between 1 and 2*nlambda) and il (runs between 1 and nlambda)
do ig = -ntgrid, ntgrid
je = jend(ig)
! if no trapped particles
if (je == 0) then
do ixi = 1, 2*nlambda
if (ixi > nlambda) then
ixi_to_isgn(ixi, ig) = 2
ixi_to_il(ixi, ig) = 2*nlambda - ixi + 1
else
ixi_to_isgn(ixi, ig) = 1
ixi_to_il(ixi, ig) = ixi
end if
end do
else
!CMR, 1/11/2013:
! Sketch of how ixi=>il mapping is arranged
!===============================================================================================
! ixi= 1, ... , je-1, je, je+1, ... , 2je-1, || 2je, 2je+1, ... , nl+je, nl+je+1, ... , 2nl
! il= 1, ... , je-1, je, je-1, ... , 1, || je, je+1, ... , nl, nl, ... , je+1
! isgn= 1, ... , 1, 2, 2, ... , 2, || 1, 1, ... , 1, 2, ... , 2
! particles passing through || je, + forbidden trapped velocity space
! nb only need one isigma for je, as v||=0 at the bounce point
!===============================================================================================
do ixi = 1, 2*nlambda
if (ixi >= nlambda + je + 1) then
ixi_to_isgn(ixi, ig) = 2
ixi_to_il(ixi, ig) = 2*nlambda + je + 1 - ixi
else if (ixi >= 2*je) then
ixi_to_isgn(ixi, ig) = 1
ixi_to_il(ixi, ig) = ixi - je
else if (ixi >= je) then
ixi_to_isgn(ixi, ig) = 2
ixi_to_il(ixi, ig) = 2*je - ixi
else
ixi_to_isgn(ixi, ig) = 1
ixi_to_il(ixi, ig) = ixi
end if
end do
end if
end do
wxi = 0.0
do ig = -ntgrid, ntgrid
do ixi = 1, 2*nlambda
il = ixi_to_il(ixi, ig)
wxi(ixi, ig) = wl(ig, il)
end do
end do
end subroutine xigridset