Makes [[libfida:pass_grid] from user defined inputs, and stores the quantities in spec_chords and npa_chords
subroutine make_diagnostic_grids
!+ Makes [[libfida:pass_grid] from user defined inputs, and stores the quantities in
!+ [[libfida:spec_chords]] and [[libfida:npa_chords]]
real(Float64), dimension(:,:,:), allocatable :: dlength
type(LOSElement), dimension(:), allocatable :: los_elem
type(ParticleTrack), dimension(:), allocatable :: tracks
real(Float64) :: r0(3), v0(3), r_enter(3), r_exit(3), r0_cyl(3), ri(3)
real(Float64), dimension(3,3) :: basis
real(Float64), dimension(2) :: randomu
real(Float64) :: theta, sqrt_rho, dl, length
character(len=20) :: system = ''
real(Float64), dimension(:), allocatable :: probs
real(Float64), dimension(:,:), allocatable :: eff_rds
real(Float64), parameter :: inv_4pi = (4.d0*pi)**(-1)
real(Float64), dimension(3) :: eff_rd, rd, rd_d, r0_d
real(Float64), dimension(3,3) :: inv_basis
real(Float64), dimension(:),allocatable :: xd, yd
type(LocalEMFields) :: fields
real(Float64) :: total_prob, hh, hw, dprob, dx, dy, r, pitch
integer :: ichan,k,id,ix,iy,d_index,nd,ind_d(2)
integer :: i, j, ic, nc, ntrack, ind(3), ii, jj, kk
integer :: error
if((inputs%calc_pfida+inputs%calc_pnpa+inputs%calc_cold).gt.0) then
if(inter_grid%nphi.gt.1) then
pass_grid = inter_grid
else
call make_passive_grid()
endif
endif
!! Spectral line-of-sight passive neutral grid intersection calculations
allocate(spec_chords%cyl_inter(pass_grid%nr,pass_grid%nz,pass_grid%nphi))
if((inputs%calc_pfida+inputs%calc_cold).gt.0) then
allocate(tracks(pass_grid%ntrack))
allocate(dlength(pass_grid%nr, &
pass_grid%nz, &
pass_grid%nphi) )
pass_grid_chan_loop: do i=1,spec_chords%nchan
r0 = spec_chords%los(i)%lens_uvw
v0 = spec_chords%los(i)%axis_uvw
v0 = v0/norm2(v0)
call line_basis(r0,v0,basis)
call grid_intersect(r0,v0,length,r_enter,r_exit,passive=.True.)
if(length.le.0.d0) then
if(inputs%verbose.ge.1) then
WRITE(*,'("Channel ",i5," missed the passive neutral grid or starts inside the plasma")') i
endif
cycle pass_grid_chan_loop
endif
if(spec_chords%los(i)%spot_size.le.0.d0) then
nc = 1
else
nc = 100
endif
dlength = 0.d0
!$OMP PARALLEL DO schedule(guided) private(ic,randomu,sqrt_rho,theta,r0, &
!$OMP& length, r_enter, r_exit, j, tracks, ntrack, ind)
do ic=1,nc
! Uniformally sample within spot size
call randu(randomu)
sqrt_rho = sqrt(randomu(1))
theta = 2*pi*randomu(2)
r0(1) = 0.d0
r0(2) = spec_chords%los(i)%spot_size*sqrt_rho*cos(theta)
r0(3) = spec_chords%los(i)%spot_size*sqrt_rho*sin(theta)
r0 = matmul(basis,r0) + spec_chords%los(i)%lens_uvw
call grid_intersect(r0,v0,length,r_enter,r_exit,passive=.True.)
call track_cylindrical(r_enter, v0, tracks, ntrack)
pass_grid_track_loop: do j=1, ntrack
ind = tracks(j)%ind
!inds can repeat so add rather than assign
!$OMP ATOMIC UPDATE
dlength(ind(1),ind(2),ind(3)) = &
dlength(ind(1),ind(2),ind(3)) + tracks(j)%time/real(nc) !time == distance
!$OMP END ATOMIC
enddo pass_grid_track_loop
enddo
!$OMP END PARALLEL DO
do kk=1,pass_grid%nphi
do jj=1,pass_grid%nz
rloop: do ii=1, pass_grid%nr
if(dlength(ii,jj,kk).ne.0.d0) then
dl = dlength(ii,jj,kk)
nc = spec_chords%cyl_inter(ii,jj,kk)%nchan + 1
if(nc.eq.1) then
allocate(spec_chords%cyl_inter(ii,jj,kk)%los_elem(nc))
spec_chords%cyl_inter(ii,jj,kk)%los_elem(nc) = LOSElement(i, dl)
else
allocate(los_elem(nc))
los_elem(1:(nc-1)) = spec_chords%cyl_inter(ii,jj,kk)%los_elem
los_elem(nc) = LOSElement(i, dl)
deallocate(spec_chords%cyl_inter(ii,jj,kk)%los_elem)
call move_alloc(los_elem, spec_chords%cyl_inter(ii,jj,kk)%los_elem)
endif
spec_chords%cyl_inter(ii,jj,kk)%nchan = nc
endif
enddo rloop
enddo
enddo
enddo pass_grid_chan_loop
spec_chords%cyl_ncell = count(spec_chords%cyl_inter%nchan.gt.0)
allocate(spec_chords%cyl_cell(spec_chords%cyl_ncell))
nc = 0
do ic=1,pass_grid%ngrid
call ind2sub(pass_grid%dims,ic,ind)
ii = ind(1) ; jj = ind(2) ; kk = ind(3)
if(spec_chords%cyl_inter(ii,jj,kk)%nchan.gt.0) then
nc = nc + 1
spec_chords%cyl_cell(nc) = ic
endif
enddo
deallocate(dlength, tracks)
endif
!! Spectral line-of-sight beam grid intersection calculations
if((inputs%tot_spectra+inputs%calc_fida_wght-inputs%calc_pfida).gt.0) then
allocate(dlength(beam_grid%nx, &
beam_grid%ny, &
beam_grid%nz) )
allocate(tracks(beam_grid%ntrack))
spec_chan_loop: do i=1,spec_chords%nchan
r0 = spec_chords%los(i)%lens
v0 = spec_chords%los(i)%axis
v0 = v0/norm2(v0)
call line_basis(r0,v0,basis)
call grid_intersect(r0,v0,length,r_enter,r_exit)
if(length.le.0.d0) then
if(inputs%verbose.ge.1) then
WRITE(*,'("Channel ",i5," missed the beam grid")') i
endif
cycle spec_chan_loop
endif
if(spec_chords%los(i)%spot_size.le.0.d0) then
nc = 1
else
nc = 100
endif
dlength = 0.d0
!$OMP PARALLEL DO schedule(guided) private(ic,randomu,sqrt_rho,theta,r0, &
!$OMP& length, r_enter, r_exit, j, tracks, ntrack, ind)
do ic=1,nc
! Uniformally sample within spot size
call randu(randomu)
sqrt_rho = sqrt(randomu(1))
theta = 2*pi*randomu(2)
r0(1) = 0.d0
r0(2) = spec_chords%los(i)%spot_size*sqrt_rho*cos(theta)
r0(3) = spec_chords%los(i)%spot_size*sqrt_rho*sin(theta)
r0 = matmul(basis,r0) + spec_chords%los(i)%lens
call grid_intersect(r0, v0, length, r_enter, r_exit)
call track(r_enter, v0, tracks, ntrack)
track_loop: do j=1, ntrack
ind = tracks(j)%ind
!inds can repeat so add rather than assign
!$OMP ATOMIC UPDATE
dlength(ind(1),ind(2),ind(3)) = &
dlength(ind(1),ind(2),ind(3)) + tracks(j)%time/real(nc) !time == distance
!$OMP END ATOMIC
enddo track_loop
enddo
!$OMP END PARALLEL DO
do kk=1,beam_grid%nz
do jj=1,beam_grid%ny
xloop: do ii=1, beam_grid%nx
if(dlength(ii,jj,kk).ne.0.d0) then
dl = dlength(ii,jj,kk)
nc = spec_chords%inter(ii,jj,kk)%nchan + 1
if(nc.eq.1) then
allocate(spec_chords%inter(ii,jj,kk)%los_elem(nc))
spec_chords%inter(ii,jj,kk)%los_elem(nc) = LOSElement(i, dl)
else
allocate(los_elem(nc))
los_elem(1:(nc-1)) = spec_chords%inter(ii,jj,kk)%los_elem
los_elem(nc) = LOSElement(i, dl)
deallocate(spec_chords%inter(ii,jj,kk)%los_elem)
call move_alloc(los_elem, spec_chords%inter(ii,jj,kk)%los_elem)
endif
spec_chords%inter(ii,jj,kk)%nchan = nc
endif
enddo xloop
enddo
enddo
enddo spec_chan_loop
spec_chords%ncell = count(spec_chords%inter%nchan.gt.0)
allocate(spec_chords%cell(spec_chords%ncell))
nc = 0
do ic=1,beam_grid%ngrid
call ind2sub(beam_grid%dims,ic,ind)
ii = ind(1) ; jj = ind(2) ; kk = ind(3)
if(spec_chords%inter(ii,jj,kk)%nchan.gt.0) then
nc = nc + 1
spec_chords%cell(nc) = ic
endif
enddo
endif
!! NPA probability calculations
allocate(xd(50),yd(50))
allocate(probs(beam_grid%ngrid))
allocate(eff_rds(3,beam_grid%ngrid))
allocate(npa_chords%phit(beam_grid%nx, &
beam_grid%ny, &
beam_grid%nz, &
npa_chords%nchan) )
allocate(npa_chords%hit(beam_grid%nx, &
beam_grid%ny, &
beam_grid%nz) )
npa_chords%hit = .False.
npa_chan_loop: do ichan=1,npa_chords%nchan
v0 = npa_chords%det(ichan)%aperture%origin - npa_chords%det(ichan)%detector%origin
v0 = v0/norm2(v0)
call grid_intersect(npa_chords%det(ichan)%detector%origin,v0,length,r0,r0_d)
if(length.le.0.0) then
if(inputs%verbose.ge.0) then
WRITE(*,'("Channel ",i3," centerline missed the beam grid")') ichan
endif
endif
if(inputs%calc_npa_wght.ge.1) then
hw = npa_chords%det(ichan)%detector%hw
hh = npa_chords%det(ichan)%detector%hh
nd = size(xd)
do i=1,nd
xd(i) = -hw + 2*hw*(i-0.5)/real(nd)
yd(i) = -hh + 2*hh*(i-0.5)/real(nd)
enddo
dx = abs(xd(2) - xd(1))
dy = abs(yd(2) - yd(1))
basis = npa_chords%det(ichan)%detector%basis
inv_basis = npa_chords%det(ichan)%detector%inv_basis
eff_rds = 0.d0
probs = 0.d0
! For each grid point find the probability of hitting the detector given an isotropic source
!$OMP PARALLEL DO schedule(guided) private(ic,i,j,k,ix,iy,total_prob,eff_rd,r0,r0_d, &
!$OMP& rd_d,rd,d_index,v0,dprob,r,fields,id,ind_d,ind)
do ic=istart,beam_grid%ngrid,istep
call ind2sub(beam_grid%dims,ic,ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
r0 = [beam_grid%xc(i),beam_grid%yc(j),beam_grid%zc(k)]
r0_d = matmul(inv_basis,r0-npa_chords%det(ichan)%detector%origin)
do id = 1, nd*nd
call ind2sub([nd,nd],id,ind_d)
ix = ind_d(1) ; iy = ind_d(2)
rd_d = [xd(ix),yd(iy),0.d0]
rd = matmul(basis,rd_d) + npa_chords%det(ichan)%detector%origin
v0 = rd - r0
d_index = 0
call hit_npa_detector(r0,v0,d_index,det=ichan)
if(d_index.ne.0) then
r = norm2(rd_d - r0_d)**2
dprob = (dx*dy) * inv_4pi * r0_d(3)/(r*sqrt(r))
eff_rds(:,ic) = eff_rds(:,ic) + dprob*rd
probs(ic) = probs(ic) + dprob
endif
enddo
enddo
!$OMP END PARALLEL DO
#ifdef _MPI
call parallel_sum(eff_rds)
call parallel_sum(probs)
#endif
do ic = 1, beam_grid%ngrid
if(probs(ic).gt.0.0) then
call ind2sub(beam_grid%dims, ic, ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
r0 = [beam_grid%xc(i),beam_grid%yc(j),beam_grid%zc(k)]
eff_rd = eff_rds(:,ic)/probs(ic)
call get_fields(fields,pos=r0)
v0 = (eff_rd - r0)/norm2(eff_rd - r0)
npa_chords%phit(i,j,k,ichan)%pitch = dot_product(fields%b_norm,v0)
npa_chords%phit(i,j,k,ichan)%p = probs(ic)
npa_chords%phit(i,j,k,ichan)%dir = v0
npa_chords%phit(i,j,k,ichan)%eff_rd = eff_rd
npa_chords%hit(i,j,k) = .True.
endif
enddo
total_prob = sum(probs)
if(total_prob.le.0.d0) then
if(inputs%verbose.ge.0) then
WRITE(*,'("Channel ",i3," missed the beam grid")') ichan
endif
cycle npa_chan_loop
endif
endif
enddo npa_chan_loop
deallocate(probs,eff_rds,xd,yd)
end subroutine make_diagnostic_grids