Reads neutral beam geometry and stores the quantities in nbi
subroutine read_beam
!+ Reads neutral beam geometry and stores the quantities in [[libfida:nbi]]
integer(HID_T) :: fid, gid
integer(HSIZE_T), dimension(1) :: dims
real(Float64), dimension(3) ::uvw_src, uvw_axis, pos
real(Float64) :: dis
logical :: path_valid
integer :: error
!!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, "/nbi", .True., path_valid, error)
if(.not.path_valid) then
if(inputs%verbose.ge.0) then
write(*,'(a)') 'READ_BEAM: NBI geometry is not in the geometry file'
endif
stop
endif
!!Open NBI group
call h5gopen_f(fid, "/nbi", gid, error)
!!Read in beam definitions
call h5ltread_dataset_string_f(gid, "/nbi/name",nbi%name, error)
dims(1) = 3
call h5ltread_dataset_double_f(gid, "/nbi/src", uvw_src, dims, error)
call h5ltread_dataset_double_f(gid, "/nbi/axis", uvw_axis, dims, error)
call h5ltread_dataset_double_f(gid, "/nbi/divy", nbi%divy, dims, error)
call h5ltread_dataset_double_f(gid, "/nbi/divz", nbi%divz, dims, error)
call h5ltread_dataset_int_scalar_f(gid, "/nbi/shape", nbi%shape, error)
call h5ltread_dataset_double_scalar_f(gid, "/nbi/focy", nbi%focy, error)
call h5ltread_dataset_double_scalar_f(gid, "/nbi/focz", nbi%focz, error)
call h5ltread_dataset_double_scalar_f(gid, "/nbi/widy", nbi%widy, error)
call h5ltread_dataset_double_scalar_f(gid, "/nbi/widz", nbi%widz, error)
!!Read in aperture definitions
!! Check for naperture for compatibility with old runs
call h5ltpath_valid_f(gid, "/nbi/naperture", .True., path_valid, error)
if(path_valid) then
call h5ltread_dataset_int_scalar_f(gid,"/nbi/naperture",nbi%naperture, error)
else
nbi%naperture = 0
endif
if(nbi%naperture.gt.0) then
allocate(nbi%ashape(nbi%naperture), nbi%adist(nbi%naperture), &
nbi%awidy(nbi%naperture), nbi%awidz(nbi%naperture), &
nbi%aoffy(nbi%naperture), nbi%aoffz(nbi%naperture) )
dims(1) = nbi%naperture
call h5ltread_dataset_int_f(gid, "/nbi/ashape", nbi%ashape, dims, error)
call h5ltread_dataset_double_f(gid, "/nbi/awidy", nbi%awidy, dims, error)
call h5ltread_dataset_double_f(gid, "/nbi/awidz", nbi%awidz, dims, error)
call h5ltread_dataset_double_f(gid, "/nbi/aoffy", nbi%aoffy, dims, error)
call h5ltread_dataset_double_f(gid, "/nbi/aoffz", nbi%aoffz, dims, error)
call h5ltread_dataset_double_f(gid, "/nbi/adist", nbi%adist, dims, error)
endif
!!Close NBI group
call h5gclose_f(gid, error)
!!Close file id
call h5fclose_f(fid, error)
!!Close HDF5 interface
call h5close_f(error)
!!Convert to beam grid coordinates
call uvw_to_xyz(uvw_src,nbi%src)
nbi%axis = matmul(beam_grid%inv_basis,uvw_axis)
nbi%vinj=sqrt(2.d0*nbi%einj*1.d3 *e0/(beam_mass*mass_u))*1.d2 !! [cm/s]
pos = nbi%src + 200.0*nbi%axis
dis = sqrt(sum((pos - nbi%src)**2))
nbi%beta = asin((nbi%src(3) - pos(3))/dis)
nbi%alpha = atan2(pos(2)-nbi%src(2),pos(1)-nbi%src(1))
call tb_zyx(nbi%alpha,nbi%beta,0.d0,nbi%basis,nbi%inv_basis)
if(inputs%verbose.ge.1) then
write(*,'(a)') '---- Neutral beam settings ----'
write(*,'(T2,"Beam: ",a)') nbi%name
write(*,'(T2,"Power: ",f5.2," [MW]")') nbi%pinj
write(*,'(T2,"Voltage: ",f6.2," [keV]")') nbi%einj
write(*,'(T2,"Current fractions: [",f5.2,",",f5.2,",",f5.2,"]")') nbi%current_fractions
write(*,*) ''
endif
end subroutine read_beam