Calculates NPA weights
subroutine npa_weights
!+ Calculates NPA weights
type(LocalEMFields) :: fields
type(NPAProbability) :: phit
real(Float64) :: pitch
real(Float64) :: pcxa
integer(Int32) :: det
integer(Int32) :: ii, jj, kk, i, ic !!indices
integer,dimension(1) :: ipitch
real(Float64), dimension(3) :: vi,vi_norm
real(Float64) :: vabs, fbm_denf, dE, dP
real(Float64), dimension(nlevs) :: pcx !! Rate coefficiants for CX
real(Float64), dimension(nlevs) :: states, states_i ! Density of n-states
integer, dimension(5) :: neut_types=[1,2,3,4,5]
real(Float64), dimension(3) :: pos,dpos,r_gyro
integer(Int32) :: ichan
real(Float64), dimension(:), allocatable :: ebarr, ptcharr
!! define energy - array
allocate(ebarr(inputs%ne_wght))
do i=1,inputs%ne_wght
ebarr(i)=real(i-0.5)*inputs%emax_wght/real(inputs%ne_wght)
enddo
dE=abs(ebarr(2)-ebarr(1))
!! define pitch - array
allocate(ptcharr(inputs%np_wght))
do i=1,inputs%np_wght
ptcharr(i)=real(i-0.5)*2./real(inputs%np_wght)-1.
enddo
dP=abs(ptcharr(2)-ptcharr(1))
if(inputs%verbose.ge.1) then
write(*,'(T3,"Number of Channels: ",i3)') npa_chords%nchan
write(*,'(T3,"Nenergy: ",i3)') inputs%ne_wght
write(*,'(T3,"Npitch: ",i3)') inputs%np_wght
write(*,'(T3,"Maximum energy: ",f7.2)') inputs%emax_wght
write(*,*) ''
endif
!! define storage arrays
allocate(nweight%emissivity(beam_grid%nx, &
beam_grid%ny, &
beam_grid%nz, &
npa_chords%nchan))
allocate(nweight%attenuation(inputs%ne_wght, &
beam_grid%nx, &
beam_grid%ny, &
beam_grid%nz, &
npa_chords%nchan))
allocate(nweight%cx(inputs%ne_wght, &
beam_grid%nx, &
beam_grid%ny, &
beam_grid%nz, &
npa_chords%nchan))
allocate(nweight%weight(inputs%ne_wght, &
inputs%np_wght, &
npa_chords%nchan))
allocate(nweight%flux(inputs%ne_wght, npa_chords%nchan))
nweight%emissivity = 0.d0
nweight%attenuation = 0.d0
nweight%cx = 0.d0
nweight%weight = 0.d0
nweight%flux = 0.d0
loop_over_channels: do ichan=1,npa_chords%nchan
!$OMP PARALLEL DO schedule(guided) collapse(3) private(ii,jj,kk,fields,phit,&
!$OMP& ic,det,pos,dpos,r_gyro,pitch,ipitch,vabs,vi,pcx,pcxa,states,states_i,vi_norm,fbm_denf)
loop_along_z: do kk=1,beam_grid%nz
loop_along_y: do jj=1,beam_grid%ny
loop_along_x: do ii=1,beam_grid%nx
phit = npa_chords%phit(ii,jj,kk,ichan)
if(phit%p.gt.0.d0) then
pos = [beam_grid%xc(ii), beam_grid%yc(jj), beam_grid%zc(kk)]
call get_fields(fields,pos=pos)
if(.not.fields%in_plasma) cycle loop_along_x
!!Check if it hits a detector just to make sure
dpos = phit%eff_rd
vi_norm = phit%dir
call hit_npa_detector(pos,vi_norm,det)
if (det.ne.ichan) then
if(inputs%verbose.ge.0) then
write(*,'(a)') 'NPA_WEIGHTS: Missed detector'
endif
cycle loop_along_x
endif
!! Determine the angle between the B-field and the Line of Sight
pitch = phit%pitch
ipitch=minloc(abs(ptcharr - pitch))
loop_over_energy: do ic = istart, inputs%ne_wght,istep !! energy loop
vabs = sqrt(ebarr(ic)/(v2_to_E_per_amu*inputs%ab))
vi = vi_norm*vabs
!!Correct for gyro orbit
call gyro_step(vi,fields,r_gyro)
fbm_denf=0
if (inputs%dist_type.eq.1) then
!get dist at guiding center
call get_ep_denf(ebarr(ic),pitch,fbm_denf,pos=(pos+r_gyro))
endif
if (fbm_denf.ne.fbm_denf) cycle loop_over_energy
!! -------------- calculate CX probability -------!!
call get_beam_cx_rate([ii,jj,kk],pos,vi,beam_ion,neut_types,pcx)
if(sum(pcx).le.0) cycle loop_over_energy
!!Calculate attenuation
states = pcx*1.0d14 !!needs to be large aribitrary number so colrad works
states_i=states
call attenuate(pos,dpos,vi,states)
pcxa=sum(states)/sum(states_i)
!$OMP CRITICAL(npa_wght)
nweight%attenuation(ic,ii,jj,kk,ichan) = pcxa
nweight%cx(ic,ii,jj,kk,ichan) = sum(pcx)
nweight%weight(ic,ipitch(1),ichan) = nweight%weight(ic,ipitch(1),ichan) + &
2*sum(pcx)*pcxa*phit%p*beam_grid%dv/dP
nweight%flux(ic,ichan) = nweight%flux(ic,ichan) + &
2*beam_grid%dv*fbm_denf*sum(pcx)*pcxa*phit%p
!Factor of 2 above is to convert fbm to ions/(cm^3 dE (domega/4pi))
nweight%emissivity(ii,jj,kk,ichan)=nweight%emissivity(ii,jj,kk,ichan)+ &
2*fbm_denf*sum(pcx)*pcxa*phit%p*dE
!$OMP END CRITICAL(npa_wght)
enddo loop_over_energy
endif
enddo loop_along_x
enddo loop_along_y
enddo loop_along_z
!$OMP END PARALLEL DO
enddo loop_over_channels
#ifdef _MPI
call parallel_sum(nweight%weight)
call parallel_sum(nweight%flux)
call parallel_sum(nweight%cx)
call parallel_sum(nweight%attenuation)
call parallel_sum(nweight%emissivity)
#endif
do ichan=1,npa_chords%nchan
if(inputs%verbose.ge.1) then
write(*,'(T4,"Channel: ",i3)') ichan
write(*,'(T4,"Radius: ",f10.3)') npa_chords%radius(ichan)
write(*,'(T4,A,ES14.5)') 'Flux: ',sum(nweight%flux(:,ichan))*dE
write(*,'(T4,A,ES14.5)') 'Weight: ',sum(nweight%weight(:,:,ichan))*dE*dP
write(*,*) ''
endif
enddo
if(inputs%verbose.ge.1) then
write(*,*) 'write npa weights: ' , time(time_start)
endif
#ifdef _MPI
if(my_rank().eq.0) call write_npa_weights()
#else
call write_npa_weights()
#endif
end subroutine npa_weights