read_beam Subroutine

  1. fidasim.f90
  2. libfida
  3. read_beam

public subroutine read_beam()

Reads neutral beam geometry and stores the quantities in nbi

Arguments

None

Calls

proc~~read_beam~~CallsGraph proc~read_beam read_beam h5ltread_dataset_double_f h5ltread_dataset_double_f proc~read_beam->h5ltread_dataset_double_f h5ltread_dataset_int_f h5ltread_dataset_int_f proc~read_beam->h5ltread_dataset_int_f h5close_f h5close_f proc~read_beam->h5close_f h5gopen_f h5gopen_f proc~read_beam->h5gopen_f proc~h5ltread_dataset_double_scalar_f h5ltread_dataset_double_scalar_f proc~read_beam->proc~h5ltread_dataset_double_scalar_f h5ltpath_valid_f h5ltpath_valid_f proc~read_beam->h5ltpath_valid_f h5gclose_f h5gclose_f proc~read_beam->h5gclose_f proc~uvw_to_xyz uvw_to_xyz proc~read_beam->proc~uvw_to_xyz proc~h5ltread_dataset_int_scalar_f h5ltread_dataset_int_scalar_f proc~read_beam->proc~h5ltread_dataset_int_scalar_f h5open_f h5open_f proc~read_beam->h5open_f h5ltread_dataset_string_f h5ltread_dataset_string_f proc~read_beam->h5ltread_dataset_string_f proc~tb_zyx tb_zyx proc~read_beam->proc~tb_zyx h5fclose_f h5fclose_f proc~read_beam->h5fclose_f h5fopen_f h5fopen_f proc~read_beam->h5fopen_f proc~h5ltread_dataset_double_scalar_f->h5ltread_dataset_double_f proc~h5ltread_dataset_int_scalar_f->h5ltread_dataset_int_f

Contents

Source Code

read_beam

Source Code

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/(inputs%ab*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(*,*) ''
    endif

end subroutine read_beam