Reads the spectral geometry and stores the quantities in spec_chords
subroutine read_chords
!+ Reads the spectral geometry and stores the quantities in [[libfida:spec_chords]]
integer(HID_T) :: fid, gid
integer(HSIZE_T), dimension(2) :: dims
logical :: path_valid
real(Float64), dimension(:,:), allocatable :: lenses
real(Float64), dimension(:,:), allocatable :: axes
real(Float64), dimension(:,:,:), allocatable :: dlength
real(Float64), dimension(:), allocatable :: spot_size, sigma_pi
type(LOSElement), dimension(:), allocatable :: los_elem
real(Float64) :: r0(3), v0(3), r_enter(3), r_exit(3)
real(Float64) :: xyz_lens(3), xyz_axis(3), length
real(Float64), dimension(3,3) :: basis
real(Float64), dimension(2) :: randomu
real(Float64) :: theta, sqrt_rho
type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks
character(len=20) :: system = ''
integer :: i, j, ic, nc, ncell, ind(3), ii, jj, kk
integer :: error
if(inputs%verbose.ge.1) then
write(*,'(a)') '---- FIDA/BES settings ----'
endif
!!Initialize HDF5 interface
call h5open_f(error)
!!Open HDF5 file
call h5fopen_f(inputs%geometry_file, H5F_ACC_RDONLY_F, fid, error)
!!Check if SPEC group exists
call h5ltpath_valid_f(fid, "/spec", .True., path_valid, error)
if(.not.path_valid) then
if(inputs%verbose.ge.0) then
write(*,'(a)') 'FIDA/BES geometry is not in the geometry file'
write(*,'(a)') 'Continuing without spectral diagnostics'
endif
inputs%calc_spec = 0
inputs%calc_fida = 0
inputs%calc_bes = 0
inputs%calc_brems = 0
inputs%calc_fida_wght = 0
call h5fclose_f(fid, error)
call h5close_f(error)
return
endif
!!Open SPEC group
call h5gopen_f(fid, "/spec", gid, error)
call h5ltread_dataset_string_f(gid, "/spec/system", system, error)
call h5ltread_dataset_int_scalar_f(gid, "/spec/nchan", spec_chords%nchan, error)
allocate(lenses(3, spec_chords%nchan))
allocate(axes(3, spec_chords%nchan))
allocate(spot_size(spec_chords%nchan))
allocate(sigma_pi(spec_chords%nchan))
allocate(spec_chords%los(spec_chords%nchan))
allocate(spec_chords%radius(spec_chords%nchan))
allocate(dlength(beam_grid%nx, &
beam_grid%ny, &
beam_grid%nz) )
dims = [3,spec_chords%nchan]
call h5ltread_dataset_double_f(gid, "/spec/lens", lenses, dims, error)
call h5ltread_dataset_double_f(gid, "/spec/axis", axes, dims, error)
call h5ltread_dataset_double_f(gid, "/spec/spot_size", spot_size, dims(2:2), error)
call h5ltread_dataset_double_f(gid, "/spec/sigma_pi", sigma_pi, dims(2:2), error)
call h5ltread_dataset_double_f(gid, "/spec/radius", spec_chords%radius, dims(2:2), error)
!!Close SPEC group
call h5gclose_f(gid, error)
!!Close file id
call h5fclose_f(fid, error)
!!Close HDF5 interface
call h5close_f(error)
chan_loop: do i=1,spec_chords%nchan
call uvw_to_xyz(lenses(:,i),xyz_lens)
xyz_axis = matmul(beam_grid%inv_basis,axes(:,i))
spec_chords%los(i)%lens = xyz_lens
spec_chords%los(i)%axis = xyz_axis
spec_chords%los(i)%sigma_pi = sigma_pi(i)
spec_chords%los(i)%spot_size = spot_size(i)
r0 = xyz_lens
v0 = xyz_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.0) then
WRITE(*,'("Channel ",i5," missed the beam grid")') i
endif
cycle chan_loop
endif
if(spot_size(i).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, ncell, 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) = spot_size(i)*sqrt_rho*cos(theta)
r0(3) = spot_size(i)*sqrt_rho*sin(theta)
r0 = matmul(basis,r0) + xyz_lens
call grid_intersect(r0, v0, length, r_enter, r_exit)
call track(r_enter, v0, tracks, ncell)
track_loop: do j=1, ncell
ind = tracks(j)%ind
!inds can repeat so add rather than assign
!$OMP CRITICAL(read_chords_1)
dlength(ind(1),ind(2),ind(3)) = &
dlength(ind(1),ind(2),ind(3)) + tracks(j)%time/real(nc) !time == distance
!$OMP END CRITICAL(read_chords_1)
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
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, dlength(ii,jj,kk))
else
allocate(los_elem(nc))
los_elem(1:(nc-1)) = spec_chords%inter(ii,jj,kk)%los_elem
los_elem(nc) = LOSElement(i, dlength(ii,jj,kk))
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 chan_loop
if(inputs%verbose.ge.1) then
write(*,'(T2,"FIDA/BES System: ",a)') trim(adjustl(system))
write(*,'(T2,"Number of channels: ",i5)') spec_chords%nchan
write(*,*) ''
endif
deallocate(lenses,axes,spot_size,sigma_pi)
end subroutine read_chords