read_npa Subroutine

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 proc~get_fields get_fields proc~read_npa->proc~get_fields h5ltread_dataset_string_f h5ltread_dataset_string_f proc~read_npa->h5ltread_dataset_string_f h5close_f h5close_f proc~read_npa->h5close_f h5gopen_f h5gopen_f proc~read_npa->h5gopen_f proc~plane_basis plane_basis proc~read_npa->proc~plane_basis h5ltpath_valid_f h5ltpath_valid_f proc~read_npa->h5ltpath_valid_f h5gclose_f h5gclose_f proc~read_npa->h5gclose_f proc~grid_intersect grid_intersect proc~read_npa->proc~grid_intersect proc~hit_npa_detector hit_npa_detector proc~read_npa->proc~hit_npa_detector proc~h5ltread_dataset_int_scalar_f h5ltread_dataset_int_scalar_f proc~read_npa->proc~h5ltread_dataset_int_scalar_f h5open_f h5open_f proc~read_npa->h5open_f h5ltread_dataset_double_f h5ltread_dataset_double_f proc~read_npa->h5ltread_dataset_double_f h5fclose_f h5fclose_f proc~read_npa->h5fclose_f h5fopen_f h5fopen_f proc~read_npa->h5fopen_f proc~in_plasma in_plasma proc~get_fields->proc~in_plasma proc~calc_perp_vectors calc_perp_vectors proc~get_fields->proc~calc_perp_vectors proc~cross_product cross_product proc~plane_basis->proc~cross_product proc~line_plane_intersect line_plane_intersect proc~hit_npa_detector->proc~line_plane_intersect proc~in_boundary in_boundary proc~hit_npa_detector->proc~in_boundary h5ltread_dataset_int_f h5ltread_dataset_int_f proc~h5ltread_dataset_int_scalar_f->h5ltread_dataset_int_f interface~interpol_coeff interpol_coeff proc~in_plasma->interface~interpol_coeff proc~xyz_to_uvw xyz_to_uvw proc~in_plasma->proc~xyz_to_uvw proc~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_coeff proc~interpol1d_coeff interpol1D_coeff interface~interpol_coeff->proc~interpol1d_coeff proc~interpol2d_coeff_arr interpol2D_coeff_arr interface~interpol_coeff->proc~interpol2d_coeff_arr proc~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~interpol1d_coeff_arr->proc~interpol1d_coeff

Called by

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

Contents

Source Code


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
    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