Calculates energy grid and (untrapped) pitch angle weights for integrals which drop a point from the energy/lambda dimensions. These are only used for estimating the velocity space error in [[dist_fn::get_verr]] as a part of the adaptive collisionality algorithm.
subroutine init_weights
use file_utils, only: open_output_file, close_output_file
use constants, only: pi => dpi
use species, only: nspec, f0_values
implicit none
real, dimension (:), allocatable :: modzeroes, werrtmp ! (negrid-2)
real, dimension (:), allocatable :: lmodzeroes, wlerrtmp ! (ng2-1)
integer :: ipt, ndiv, divmax, is
logical :: eflag
if (init_weights_init) return
init_weights_init = .true.
allocate(lmodzeroes(ng2-1), wlerrtmp(ng2-1))
allocate(wlerr(ng2,ng2))
wlerr = 0.0; lmodzeroes = 0.0; wlerrtmp = 0.0
wdim = nesub
allocate(modzeroes(nesub-1), werrtmp(nesub-1))
allocate(werr(negrid-1,nesub,nspec))
werr = 0.0 ; modzeroes = 0.0 ; werrtmp = 0.0
! loop to obtain weights for energy grid points. negrid-1 sets
! of weights are needed because we want to compute integrals
! for negrid-1 sets of energy points (corresponding to negrid-1
! points that we can choose to drop from the guassian quadrature)
do ipt = 1, nesub
do is = 1, nspec
! drops the point corresponding to ipt from the energy grid
if (ipt /= 1) modzeroes(:ipt-1) = speed(:ipt-1,is)
if (ipt /= nesub) modzeroes(ipt:nesub-1) = speed(ipt+1:nesub,is)
! get weights for energy grid points
call get_weights (nmax,0.0,vcut,modzeroes,werrtmp,ndiv,divmax,eflag)
! a zero is left in the position corresponding to the dropped point
if (ipt /= 1) werr(:ipt-1,ipt,is) = werrtmp(:ipt-1)
if (ipt /= nesub) werr(ipt+1:nesub,ipt,is) = werrtmp(ipt:nesub-1)
werr(nesub+1:,ipt,is) = w(nesub+1:negrid-1,is)
! absorbing volume element into weights
werr(:nesub,ipt,is) = werr(:nesub,ipt,is)*energy(:nesub,is)*2.0*pi*f0_values(:nesub,is)
end do
end do
! same thing done here for lamdba as was
! done earlier for energy space
do ipt=1,ng2
if (ipt /= 1) lmodzeroes(:ipt-1) = xx(:ipt-1)
if (ipt /= ng2) lmodzeroes(ipt:ng2-1) = xx(ipt+1:)
call get_weights (nmax,1.0,0.0,lmodzeroes,wlerrtmp,ndiv,divmax,eflag)
if (ipt /= 1) wlerr(:ipt-1,ipt) = wlerrtmp(:ipt-1)
if (ipt /= ng2) wlerr(ipt+1:,ipt) = wlerrtmp(ipt:)
end do
deallocate(modzeroes,werrtmp,lmodzeroes,wlerrtmp)
eflag = .false.
end subroutine init_weights