Reads the NPA geometry and stores the quantities in npa_chords
subroutine read_npa
!+ Reads the NPA geometry and stores the quantities in [[libfida:npa_chords]]
integer(HID_T) :: fid, gid
integer(HSIZE_T), dimension(2) :: dims
logical :: path_valid
real(Float64), dimension(:,:), allocatable :: a_tedge,a_redge,a_cent
real(Float64), dimension(:,:), allocatable :: d_tedge,d_redge,d_cent
integer, dimension(:), allocatable :: a_shape, d_shape
character(len=20) :: system = ''
real(Float64), parameter :: inv_4pi = (4.d0*pi)**(-1)
real(Float64), dimension(3) :: xyz_a_tedge,xyz_a_redge,xyz_a_cent
real(Float64), dimension(3) :: xyz_d_tedge,xyz_d_redge,xyz_d_cent
real(Float64), dimension(3) :: eff_rd, rd, rd_d, r0, r0_d, v0
real(Float64), dimension(3,3) :: basis, inv_basis
real(Float64), dimension(50) :: xd, yd
type(LocalEMFields) :: fields
real(Float64) :: length,total_prob, hh, hw, dprob, dx, dy, r, pitch
integer :: ichan,i,j,k,ix,iy,d_index,nd,cnt
integer :: error
!!Initialize HDF5 interface
call h5open_f(error)
!!Open HDF5 file
call h5fopen_f(inputs%geometry_file, H5F_ACC_RDWR_F, fid, error)
!!Check if NPA group exists
call h5ltpath_valid_f(fid, "/npa", .True., path_valid, error)
if(.not.path_valid) then
if(inputs%verbose.ge.0) then
write(*,'(a)') 'NPA geometry is not in the geometry file'
write(*,'(a)') 'Continuing without NPA diagnostics'
endif
inputs%calc_npa = 0
inputs%calc_npa_wght = 0
call h5fclose_f(fid, error)
call h5close_f(error)
return
endif
!!Open NPA group
call h5gopen_f(fid, "/npa", gid, error)
call h5ltread_dataset_string_f(gid, "/npa/system", system, error)
call h5ltread_dataset_int_scalar_f(gid, "/npa/nchan", npa_chords%nchan, error)
if(inputs%verbose.ge.1) then
write(*,'(a)') "---- NPA settings ----"
write(*,'(T2,"NPA System: ", a)') trim(adjustl(system))
write(*,'(T2,"Number of channels: ",i3)') npa_chords%nchan
endif
allocate(a_tedge(3, npa_chords%nchan))
allocate(a_redge(3, npa_chords%nchan))
allocate(a_cent(3, npa_chords%nchan))
allocate(a_shape(npa_chords%nchan))
allocate(d_tedge(3, npa_chords%nchan))
allocate(d_redge(3, npa_chords%nchan))
allocate(d_cent(3, npa_chords%nchan))
allocate(d_shape(npa_chords%nchan))
allocate(npa_chords%radius(npa_chords%nchan))
allocate(npa_chords%det(npa_chords%nchan))
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.
dims = [3,spec_chords%nchan]
call h5ltread_dataset_double_f(gid,"/npa/radius", npa_chords%radius, dims(2:2), error)
call h5ltread_dataset_int_f(gid, "/npa/a_shape", a_shape, dims(2:2), error)
call h5ltread_dataset_double_f(gid, "/npa/a_tedge", a_tedge, dims, error)
call h5ltread_dataset_double_f(gid, "/npa/a_redge", a_redge, dims, error)
call h5ltread_dataset_double_f(gid, "/npa/a_cent", a_cent, dims, error)
call h5ltread_dataset_int_f(gid, "/npa/d_shape", d_shape, dims(2:2), error)
call h5ltread_dataset_double_f(gid, "/npa/d_tedge", d_tedge, dims, error)
call h5ltread_dataset_double_f(gid, "/npa/d_redge", d_redge, dims, error)
call h5ltread_dataset_double_f(gid, "/npa/d_cent", d_cent, dims, error)
!!Close NPA group
call h5gclose_f(gid, error)
!!Close file id
call h5fclose_f(fid, error)
!!Close HDF5 interface
call h5close_f(error)
! Define detector/aperture shape
npa_chords%det%detector%shape = d_shape
npa_chords%det%aperture%shape = a_shape
chan_loop: do ichan=1,npa_chords%nchan
! Convert to beam grid coordinates
call uvw_to_xyz(a_cent(:,ichan), xyz_a_cent)
call uvw_to_xyz(a_redge(:,ichan),xyz_a_redge)
call uvw_to_xyz(a_tedge(:,ichan),xyz_a_tedge)
call uvw_to_xyz(d_cent(:,ichan), xyz_d_cent)
call uvw_to_xyz(d_redge(:,ichan),xyz_d_redge)
call uvw_to_xyz(d_tedge(:,ichan),xyz_d_tedge)
! Define detector/aperture hh/hw
npa_chords%det(ichan)%detector%hw = norm2(xyz_d_redge - xyz_d_cent)
npa_chords%det(ichan)%aperture%hw = norm2(xyz_a_redge - xyz_a_cent)
npa_chords%det(ichan)%detector%hh = norm2(xyz_d_tedge - xyz_d_cent)
npa_chords%det(ichan)%aperture%hh = norm2(xyz_a_tedge - xyz_a_cent)
! Define detector/aperture origin
npa_chords%det(ichan)%detector%origin = xyz_d_cent
npa_chords%det(ichan)%aperture%origin = xyz_a_cent
! Define detector/aperture basis
call plane_basis(xyz_d_cent, xyz_d_redge, xyz_d_tedge, &
npa_chords%det(ichan)%detector%basis, &
npa_chords%det(ichan)%detector%inv_basis)
call plane_basis(xyz_a_cent, xyz_a_redge, xyz_a_tedge, &
npa_chords%det(ichan)%aperture%basis, &
npa_chords%det(ichan)%aperture%inv_basis)
v0 = xyz_a_cent - xyz_d_cent
v0 = v0/norm2(v0)
call grid_intersect(xyz_d_cent,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
cnt = 0
! For each grid point find the probability of hitting the detector given an isotropic source
!$OMP PARALLEL DO schedule(guided) collapse(3) private(i,j,k,ix,iy,total_prob,eff_rd,r0,r0_d, &
!$OMP& rd_d,rd,d_index,v0,dprob,r,fields)
do k=1,beam_grid%nz
do j=1,beam_grid%ny
do i=1,beam_grid%nx
cnt = cnt+1
total_prob = 0.d0
eff_rd = eff_rd*0.d0
r0 = [beam_grid%xc(i),beam_grid%yc(j),beam_grid%zc(k)]
r0_d = matmul(inv_basis,r0-xyz_d_cent)
do ix = 1, nd
do iy = 1, nd
rd_d = [xd(ix),yd(iy),0.d0]
rd = matmul(basis,rd_d) + xyz_d_cent
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_rd = eff_rd + dprob*rd
total_prob = total_prob + dprob
endif
enddo !yd loop
enddo !xd loop
if(total_prob.gt.0.0) then
eff_rd = eff_rd/total_prob
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 = total_prob
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
if(inputs%verbose.ge.2) then
WRITE(*,'(T4,"Channel: ",i5," ",f7.2,"% completed",a,$)') &
ichan, cnt/real(beam_grid%ngrid)*100,char(13)
endif
enddo !x loop
enddo !y loop
enddo !z loop
!$OMP END PARALLEL DO
total_prob = sum(npa_chords%phit(:,:,:,ichan)%p)
if(total_prob.le.0.d0) then
if(inputs%verbose.ge.0) then
WRITE(*,'("Channel ",i3," missed the beam grid")') ichan
endif
cycle chan_loop
endif
endif
enddo chan_loop
if(inputs%verbose.ge.1) write(*,'(50X,a)') ""
deallocate(a_shape,a_cent,a_redge,a_tedge)
deallocate(d_shape,d_cent,d_redge,d_tedge)
end subroutine read_npa