read_chords Subroutine

public subroutine read_chords()

Reads the spectral geometry and stores the quantities in spec_chords

Arguments

None

Calls

proc~~read_chords~~CallsGraph proc~read_chords read_chords h5ltread_dataset_string_f h5ltread_dataset_string_f proc~read_chords->h5ltread_dataset_string_f proc~randu randu proc~read_chords->proc~randu h5close_f h5close_f proc~read_chords->h5close_f h5gopen_f h5gopen_f proc~read_chords->h5gopen_f h5ltpath_valid_f h5ltpath_valid_f proc~read_chords->h5ltpath_valid_f h5gclose_f h5gclose_f proc~read_chords->h5gclose_f proc~grid_intersect grid_intersect proc~read_chords->proc~grid_intersect proc~track track proc~read_chords->proc~track proc~h5ltread_dataset_int_scalar_f h5ltread_dataset_int_scalar_f proc~read_chords->proc~h5ltread_dataset_int_scalar_f h5open_f h5open_f proc~read_chords->h5open_f h5ltread_dataset_double_f h5ltread_dataset_double_f proc~read_chords->h5ltread_dataset_double_f proc~line_basis line_basis proc~read_chords->proc~line_basis h5fclose_f h5fclose_f proc~read_chords->h5fclose_f h5fopen_f h5fopen_f proc~read_chords->h5fopen_f omp_get_thread_num omp_get_thread_num proc~randu->omp_get_thread_num proc~rng_uniform rng_uniform proc~randu->proc~rng_uniform proc~get_indices get_indices proc~track->proc~get_indices proc~in_plasma in_plasma proc~track->proc~in_plasma h5ltread_dataset_int_f h5ltread_dataset_int_f proc~h5ltread_dataset_int_scalar_f->h5ltread_dataset_int_f proc~tb_zyx tb_zyx proc~line_basis->proc~tb_zyx 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_chords~~CalledByGraph proc~read_chords read_chords program~fidasim fidasim program~fidasim->proc~read_chords

Contents

Source Code


Source Code

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