read_npa Subroutine

  1. fidasim.f90
  2. libfida
  3. read_npa

public subroutine read_npa()

Reads the NPA geometry and stores the quantities in npa_chords

Arguments

None

Calls

proc~~read_npa~~CallsGraph proc~read_npa read_npa h5ltpath_valid_f h5ltpath_valid_f proc~read_npa->h5ltpath_valid_f h5gopen_f h5gopen_f proc~read_npa->h5gopen_f h5fopen_f h5fopen_f proc~read_npa->h5fopen_f h5open_f h5open_f proc~read_npa->h5open_f h5fclose_f h5fclose_f proc~read_npa->h5fclose_f proc~h5ltread_dataset_int_scalar_f h5ltread_dataset_int_scalar_f proc~read_npa->proc~h5ltread_dataset_int_scalar_f h5ltread_dataset_int_f h5ltread_dataset_int_f proc~read_npa->h5ltread_dataset_int_f h5close_f h5close_f proc~read_npa->h5close_f h5ltread_dataset_string_f h5ltread_dataset_string_f proc~read_npa->h5ltread_dataset_string_f h5ltread_dataset_double_f h5ltread_dataset_double_f proc~read_npa->h5ltread_dataset_double_f h5gclose_f h5gclose_f proc~read_npa->h5gclose_f proc~uvw_to_xyz uvw_to_xyz proc~read_npa->proc~uvw_to_xyz proc~plane_basis plane_basis proc~read_npa->proc~plane_basis proc~h5ltread_dataset_int_scalar_f->h5ltread_dataset_int_f proc~cross_product cross_product proc~plane_basis->proc~cross_product

Called by

proc~~read_npa~~CalledByGraph proc~read_npa read_npa program~fidasim fidasim program~fidasim->proc~read_npa

Contents

Source Code

read_npa

Source Code

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
    character(len=20) :: system = ''

    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,3) :: basis, inv_basis
    real(Float64) :: hh, hw
    integer :: ichan
    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_pnpa = 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(d_tedge(3, npa_chords%nchan))
    allocate(d_redge(3, npa_chords%nchan))
    allocate(d_cent(3,  npa_chords%nchan))
    allocate(npa_chords%radius(npa_chords%nchan))
    allocate(npa_chords%det(npa_chords%nchan))

    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", npa_chords%det%aperture%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", npa_chords%det%detector%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)

    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)
    enddo chan_loop

    if(inputs%verbose.ge.1) write(*,'(50X,a)') ""

    deallocate(a_cent,a_redge,a_tedge)
    deallocate(d_cent,d_redge,d_tedge)

end subroutine read_npa