fida_weights_los Subroutine

public subroutine fida_weights_los()

Calculates LOS averaged FIDA weights

Arguments

None

Calls

proc~~fida_weights_los~~CallsGraph proc~fida_weights_los fida_weights_los proc~get_fields get_fields proc~fida_weights_los->proc~get_fields proc~write_fida_weights write_fida_weights proc~fida_weights_los->proc~write_fida_weights proc~get_plasma get_plasma proc~fida_weights_los->proc~get_plasma proc~in_plasma in_plasma proc~fida_weights_los->proc~in_plasma proc~grid_intersect grid_intersect proc~fida_weights_los->proc~grid_intersect proc~track track proc~fida_weights_los->proc~track proc~colrad colrad proc~fida_weights_los->proc~colrad proc~bt_cx_rates bt_cx_rates proc~fida_weights_los->proc~bt_cx_rates interface~parallel_sum parallel_sum proc~fida_weights_los->interface~parallel_sum proc~pitch_to_vec pitch_to_vec proc~fida_weights_los->proc~pitch_to_vec proc~store_fw_photons_at_chan store_fw_photons_at_chan proc~fida_weights_los->proc~store_fw_photons_at_chan proc~bb_cx_rates bb_cx_rates proc~fida_weights_los->proc~bb_cx_rates proc~my_rank my_rank proc~fida_weights_los->proc~my_rank proc~get_fields->proc~in_plasma proc~calc_perp_vectors calc_perp_vectors proc~get_fields->proc~calc_perp_vectors proc~uvw_to_xyz uvw_to_xyz proc~get_fields->proc~uvw_to_xyz proc~xyz_to_uvw xyz_to_uvw proc~get_fields->proc~xyz_to_uvw h5fclose_f h5fclose_f proc~write_fida_weights->h5fclose_f h5ltmake_dataset_int_f h5ltmake_dataset_int_f proc~write_fida_weights->h5ltmake_dataset_int_f h5close_f h5close_f proc~write_fida_weights->h5close_f h5ltset_attribute_string_f h5ltset_attribute_string_f proc~write_fida_weights->h5ltset_attribute_string_f interface~h5ltmake_compressed_dataset_double_f h5ltmake_compressed_dataset_double_f proc~write_fida_weights->interface~h5ltmake_compressed_dataset_double_f h5fcreate_f h5fcreate_f proc~write_fida_weights->h5fcreate_f h5open_f h5open_f proc~write_fida_weights->h5open_f proc~get_plasma->proc~in_plasma proc~get_plasma->proc~uvw_to_xyz proc~get_position get_position proc~get_plasma->proc~get_position proc~get_plasma->proc~xyz_to_uvw proc~cyl_to_uvw cyl_to_uvw proc~in_plasma->proc~cyl_to_uvw proc~in_plasma->proc~xyz_to_uvw interface~interpol_coeff interpol_coeff proc~in_plasma->interface~interpol_coeff proc~in_passive_grid in_passive_grid proc~grid_intersect->proc~in_passive_grid proc~track->proc~get_fields proc~track->proc~in_plasma proc~get_indices get_indices proc~track->proc~get_indices proc~doppler_stark doppler_stark proc~track->proc~doppler_stark proc~linsolve linsolve proc~colrad->proc~linsolve proc~eigen eigen proc~colrad->proc~eigen proc~get_rate_matrix get_rate_matrix proc~colrad->proc~get_rate_matrix proc~bt_cx_rates->interface~interpol_coeff proc~parallel_sum_d4 parallel_sum_d4 interface~parallel_sum->proc~parallel_sum_d4 proc~parallel_sum_d5 parallel_sum_d5 interface~parallel_sum->proc~parallel_sum_d5 proc~parallel_sum_d2 parallel_sum_d2 interface~parallel_sum->proc~parallel_sum_d2 proc~parallel_sum_d3 parallel_sum_d3 interface~parallel_sum->proc~parallel_sum_d3 proc~parallel_sum_d0 parallel_sum_d0 interface~parallel_sum->proc~parallel_sum_d0 proc~parallel_sum_d1 parallel_sum_d1 interface~parallel_sum->proc~parallel_sum_d1 proc~parallel_sum_i1 parallel_sum_i1 interface~parallel_sum->proc~parallel_sum_i1 proc~parallel_sum_i0 parallel_sum_i0 interface~parallel_sum->proc~parallel_sum_i0 proc~parallel_sum_i2 parallel_sum_i2 interface~parallel_sum->proc~parallel_sum_i2 proc~spectrum spectrum proc~store_fw_photons_at_chan->proc~spectrum proc~bb_cx_rates->interface~interpol_coeff proc~cyl_to_xyz cyl_to_xyz proc~get_position->proc~cyl_to_xyz mpi_allreduce mpi_allreduce proc~parallel_sum_d4->mpi_allreduce proc~parallel_sum_d5->mpi_allreduce proc~parallel_sum_d2->mpi_allreduce proc~parallel_sum_d3->mpi_allreduce proc~parallel_sum_d0->mpi_allreduce proc~parallel_sum_d1->mpi_allreduce dgetrf dgetrf proc~linsolve->dgetrf proc~matinv matinv proc~linsolve->proc~matinv dgetrs dgetrs proc~linsolve->dgetrs proc~parallel_sum_i1->mpi_allreduce proc~parallel_sum_i0->mpi_allreduce proc~parallel_sum_i2->mpi_allreduce proc~elmtrans elmtrans proc~eigen->proc~elmtrans proc~balback balback proc~eigen->proc~balback proc~hqr2 hqr2 proc~eigen->proc~hqr2 proc~elmhes elmhes proc~eigen->proc~elmhes proc~balance balance proc~eigen->proc~balance proc~h5ltmake_compressed_dataset_double_f_7 h5ltmake_compressed_dataset_double_f_7 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_7 proc~h5ltmake_compressed_dataset_double_f_6 h5ltmake_compressed_dataset_double_f_6 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_6 proc~h5ltmake_compressed_dataset_double_f_5 h5ltmake_compressed_dataset_double_f_5 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_5 proc~h5ltmake_compressed_dataset_double_f_4 h5ltmake_compressed_dataset_double_f_4 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_4 proc~h5ltmake_compressed_dataset_double_f_3 h5ltmake_compressed_dataset_double_f_3 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_3 proc~h5ltmake_compressed_dataset_double_f_2 h5ltmake_compressed_dataset_double_f_2 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_2 proc~h5ltmake_compressed_dataset_double_f_1 h5ltmake_compressed_dataset_double_f_1 interface~h5ltmake_compressed_dataset_double_f->proc~h5ltmake_compressed_dataset_double_f_1 proc~get_rate_matrix->interface~interpol_coeff proc~interpol1d_coeff_arr interpol1D_coeff_arr interface~interpol_coeff->proc~interpol1d_coeff_arr proc~cyl_interpol3d_coeff_arr cyl_interpol3D_coeff_arr interface~interpol_coeff->proc~cyl_interpol3d_coeff_arr proc~interpol1d_coeff interpol1D_coeff interface~interpol_coeff->proc~interpol1d_coeff proc~interpol2d_coeff interpol2D_coeff interface~interpol_coeff->proc~interpol2d_coeff proc~cyl_interpol3d_coeff cyl_interpol3D_coeff interface~interpol_coeff->proc~cyl_interpol3d_coeff proc~interpol2d_coeff_arr interpol2D_coeff_arr interface~interpol_coeff->proc~interpol2d_coeff_arr proc~approx_ge approx_ge proc~in_passive_grid->proc~approx_ge proc~approx_le approx_le proc~in_passive_grid->proc~approx_le proc~uvw_to_cyl uvw_to_cyl proc~in_passive_grid->proc~uvw_to_cyl proc~cyl_to_xyz->proc~cyl_to_uvw proc~cyl_to_xyz->proc~uvw_to_xyz proc~approx_eq approx_eq proc~approx_ge->proc~approx_eq proc~interpol1d_coeff_arr->proc~interpol1d_coeff proc~cyl_interpol3d_coeff_arr->proc~interpol2d_coeff proc~cyl_interpol3d_coeff_arr->proc~cyl_interpol3d_coeff proc~approx_le->proc~approx_eq h5dwrite_f h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_7->h5dwrite_f h5dclose_f h5dclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5dclose_f h5sclose_f h5sclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5sclose_f h5ltmake_dataset_double_f h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_7->h5ltmake_dataset_double_f h5pset_chunk_f h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_chunk_f h5pcreate_f h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_7->h5pcreate_f proc~chunk_size chunk_size proc~h5ltmake_compressed_dataset_double_f_7->proc~chunk_size h5pset_deflate_f h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_deflate_f h5pset_shuffle_f h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_7->h5pset_shuffle_f h5screate_simple_f h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_7->h5screate_simple_f h5pclose_f h5pclose_f proc~h5ltmake_compressed_dataset_double_f_7->h5pclose_f h5dcreate_f h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_7->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_6->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_6->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_6->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_6->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_6->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_6->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_6->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_6->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_6->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_6->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_6->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_6->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_5->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_5->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_5->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_5->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_5->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_5->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_5->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_5->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_5->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_5->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_5->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_5->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_4->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_4->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_4->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_4->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_4->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_4->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_4->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_4->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_4->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_4->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_3->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_3->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_3->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_3->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_3->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_3->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_3->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_3->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_3->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_3->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_3->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_3->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_2->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_2->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_2->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_2->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_2->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_2->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_2->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_2->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_2->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_2->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_2->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_2->h5dcreate_f proc~h5ltmake_compressed_dataset_double_f_1->h5dwrite_f proc~h5ltmake_compressed_dataset_double_f_1->h5dclose_f proc~h5ltmake_compressed_dataset_double_f_1->h5sclose_f proc~h5ltmake_compressed_dataset_double_f_1->h5ltmake_dataset_double_f proc~h5ltmake_compressed_dataset_double_f_1->h5pset_chunk_f proc~h5ltmake_compressed_dataset_double_f_1->h5pcreate_f proc~h5ltmake_compressed_dataset_double_f_1->proc~chunk_size proc~h5ltmake_compressed_dataset_double_f_1->h5pset_deflate_f proc~h5ltmake_compressed_dataset_double_f_1->h5pset_shuffle_f proc~h5ltmake_compressed_dataset_double_f_1->h5screate_simple_f proc~h5ltmake_compressed_dataset_double_f_1->h5pclose_f proc~h5ltmake_compressed_dataset_double_f_1->h5dcreate_f proc~cyl_interpol3d_coeff->proc~interpol2d_coeff proc~hqrvec hqrvec proc~hqr2->proc~hqrvec proc~interpol2d_coeff_arr->proc~interpol2d_coeff proc~ludcmp ludcmp proc~matinv->proc~ludcmp proc~outerprod outerprod proc~ludcmp->proc~outerprod

Called by

proc~~fida_weights_los~~CalledByGraph proc~fida_weights_los fida_weights_los program~fidasim fidasim program~fidasim->proc~fida_weights_los

Contents

Source Code


Source Code

subroutine fida_weights_los
    !+ Calculates LOS averaged FIDA weights
    type(LocalProfiles) :: plasma, plasma_cell
    type(LocalEMFields) :: fields, fields_cell
    real(Float64) :: denf
    real(Float64) :: wght, wght_tot
    real(Float64) :: photons !! photon flux
    real(Float64) :: length
    type(ParticleTrack), dimension(beam_grid%ntrack) :: tracks
    integer :: nwav
    integer(Int32) :: i, j, k, ienergy, cid, cind
    integer(Int32) :: ipitch, igyro, icell, ichan
    real(Float64), dimension(:), allocatable :: ebarr,ptcharr,phiarr
    real(Float64), dimension(:,:), allocatable :: mean_f
    real(Float64), dimension(3) :: vi, vi_norm, vp
    real(Float64), dimension(3) :: vnbi_f, vnbi_h, vnbi_t, vhalo
    real(Float64), dimension(3) :: r_enter, r_exit
    real(Float64) :: vabs, dE, dP
    !! Determination of the CX probability
    real(Float64), dimension(nlevs) :: fdens,hdens,tdens,dcxdens,halodens
    real(Float64), dimension(nlevs) :: rates
    real(Float64), dimension(nlevs) :: states ! Density of n-states
    real(Float64), dimension(nlevs) :: denn  ! Density of n-states
    !! COLRAD
    real(Float64) :: dt, max_dens, dlength, sigma_pi
    type(LOSInters) :: inter
    real(Float64) :: eb, ptch, phi
    !! Solution of differential equation
    integer, dimension(3) :: ind  !!actual cell
    real(Float64), dimension(3) :: ri
    integer(Int32) :: ntrack
    logical :: inp

    real(Float64):: etov2, dEdP

    nwav = inputs%nlambda_wght

    !! Define energy array
    allocate(ebarr(inputs%ne_wght))
    do i=1,inputs%ne_wght
        ebarr(i)=real(i-0.5)*inputs%emax_wght/real(inputs%ne_wght)
    enddo
    dE = abs(ebarr(2)-ebarr(1))

    !! Define pitch array
    allocate(ptcharr(inputs%np_wght))
    do i=1,inputs%np_wght
        ptcharr(i)=real(i-0.5)*2./real(inputs%np_wght)-1.
    enddo
    dP = abs(ptcharr(2)-ptcharr(1))
    dEdP = dE*dP

    !! define gyro - array
    allocate(phiarr(inputs%nphi_wght))
    do i=1,inputs%nphi_wght
        phiarr(i)=real(i-0.5)*2.d0*pi/real(inputs%nphi_wght)
    enddo

    !! allocate storage arrays
    allocate(fweight%mean_f(inputs%ne_wght,inputs%np_wght,spec_chords%nchan))
    allocate(fweight%weight(nwav,inputs%ne_wght,inputs%np_wght,spec_chords%nchan))

    allocate(mean_f(inputs%ne_wght,inputs%np_wght))
    !! zero out arrays
    fweight%weight = 0.d0
    fweight%mean_f = 0.d0
    mean_f = 0.d0

    if(inputs%verbose.ge.1) then
        write(*,'(T3,"Number of Channels: ",i5)') spec_chords%nchan
        write(*,'(T3,"Nlambda: ",i4)') nwav
        write(*,'(T3,"Nenergy: ",i3)') inputs%ne_wght
        write(*,'(T3,"Npitch: ",i3)') inputs%np_wght
        write(*,'(T3,"Ngyro: ", i3)') inputs%nphi_wght
        write(*,'(T3,"Maximum Energy: ",f7.2)') inputs%emax_wght
        write(*,'(T3,"LOS averaged: ",a)') "True"
        write(*,*) ''
    endif

    etov2 = 1.0/(v2_to_E_per_amu*inputs%ab)

    chan_loop: do ichan=1,spec_chords%nchan
        fdens = 0.d0 ; hdens = 0.d0 ; tdens = 0.d0
        halodens = 0.d0 ; dcxdens = 0.d0
        plasma = plasma*0.d0
        fields = fields*0.d0
        wght_tot = 0.d0
        mean_f = 0.d0
        do k=1,beam_grid%nz
            do j=1,beam_grid%ny
                x_loop: do i=1,beam_grid%nx
                    inter = spec_chords%inter(i,j,k)
                    cid = 0
                    cind = 0
                    do while (cid.ne.ichan.and.cind.lt.inter%nchan)
                        cind = cind + 1
                        cid = inter%los_elem(cind)%id
                    enddo
                    if(cid.eq.ichan) then
                        ind = [i,j,k]
                        ri = [beam_grid%xc(i), beam_grid%yc(j), beam_grid%zc(k)]
                        call in_plasma(ri,inp)
                        if(.not.inp) cycle x_loop
                        dlength = inter%los_elem(cind)%length
                        fdens = fdens + neut%full(:,i,j,k)*dlength
                        hdens = hdens + neut%half(:,i,j,k)*dlength
                        tdens = tdens + neut%third(:,i,j,k)*dlength
                        dcxdens = dcxdens + neut%dcx(:,i,j,k)*dlength
                        halodens = halodens + neut%halo(:,i,j,k)*dlength
                        wght = (neut%full(3,i,j,k) + neut%half(3,i,j,k) + &
                                neut%third(3,i,j,k) + neut%dcx(3,i,j,k) + &
                                neut%halo(3,i,j,k))*dlength

                        call get_plasma(plasma_cell,pos=ri)
                        call get_fields(fields_cell,pos=ri)
                        plasma = plasma + wght*plasma_cell
                        fields = fields + wght*fields_cell
                        if (inputs%dist_type.eq.1) then
                            do ipitch=1,inputs%np_wght
                                do ienergy=1,inputs%ne_wght
                                    call get_ep_denf(ebarr(ienergy),ptcharr(ipitch),denf,coeffs=fields_cell%b)
                                    mean_f(ienergy,ipitch) = mean_f(ienergy,ipitch) + wght*denf
                                enddo
                            enddo
                        endif
                        wght_tot = wght_tot + wght
                    endif
                enddo x_loop
            enddo
        enddo

        if(wght_tot.le.0) then
            if(inputs%verbose.ge.1) then
                write(*,'(T4,"Skipping channel ",i5,": Neutral density is zero")') ichan
            endif
            cycle chan_loop
        else
            plasma = plasma/wght_tot
            plasma%in_plasma = .True.
            fields = fields/wght_tot
            fields%in_plasma= .True.
            mean_f = mean_f/wght_tot
            if(inputs%verbose.ge.1) then
                write(*,'(T4,"Channel: ",i5)') ichan
                write(*,'(T4,"Radius: ",f7.2)') spec_chords%radius(ichan)
                write(*,'(T4,"Mean Fast-ion Density: ",ES14.5)') sum(mean_f)*dEdP
                write(*,*)''
            endif
        endif

        ri = plasma%pos
        vp = ri - spec_chords%los(ichan)%lens
        vnbi_f = ri - nbi%src
        vnbi_f = vnbi_f/norm2(vnbi_f)*nbi%vinj
        vnbi_h = vnbi_f/sqrt(2.d0)
        vnbi_t = vnbi_f/sqrt(3.d0)
        sigma_pi = spec_chords%los(ichan)%sigma_pi
        dlength = 1.d0

        !$OMP PARALLEL DO schedule(guided) collapse(3) private(eb,vabs,ptch,phi,vi,vi_norm, &
        !$OMP& r_enter,r_exit,length,max_dens,ind,tracks,ntrack,dt,icell,states,rates, &
        !$OMP& vhalo,denn,denf,photons,ienergy,ipitch,igyro)
        do ienergy=istart,inputs%ne_wght,istep
            do ipitch=1,inputs%np_wght
                do igyro=1,inputs%nphi_wght
                    eb = ebarr(ienergy)
                    vabs = sqrt(eb*etov2)
                    ptch = ptcharr(ipitch)
                    phi = phiarr(igyro)
                    call pitch_to_vec(ptch,phi,fields,vi_norm)
                    vi = vabs*vi_norm

                    call grid_intersect(ri,vi,length,r_enter,r_exit)
                    call track(r_enter, vi, tracks, ntrack)
                    max_dens = 0.d0
                    do icell=1,ntrack
                        ind = tracks(icell)%ind
                        tracks(icell)%flux = neut%full(3,ind(1),ind(2),ind(3))  + &
                                             neut%half(3,ind(1),ind(2),ind(3))  + &
                                             neut%third(3,ind(1),ind(2),ind(3)) + &
                                             neut%dcx(3,ind(1),ind(2),ind(3))   + &
                                             neut%halo(3,ind(1),ind(2),ind(3))
                        if(tracks(icell)%flux.gt.max_dens) max_dens=tracks(icell)%flux
                    enddo
                    dt = 0.d0
                    do icell=1,ntrack
                        if(tracks(icell)%flux.gt.(0.5*max_dens)) then
                            dt = dt + tracks(icell)%time
                        endif
                    enddo

                    states=0.d0
                    call bb_cx_rates(fdens,vi,vnbi_f,rates)
                    states = states + rates
                    call bb_cx_rates(hdens,vi,vnbi_h,rates)
                    states = states + rates
                    call bb_cx_rates(tdens,vi,vnbi_t,rates)
                    states = states + rates
                    call bt_cx_rates(plasma, dcxdens + halodens, vi, beam_ion, rates)
                    states = states + rates

                    call colrad(plasma,beam_ion,vi,dt,states,denn,photons)
                    denf = mean_f(ienergy,ipitch)*dEdP
                    photons = photons/real(inputs%nphi_wght)
                    call store_fw_photons_at_chan(ichan, ienergy, ipitch, &
                         vp, vi, fields, dlength, sigma_pi, denf, photons)

                enddo
            enddo
        enddo
        !$OMP END PARALLEL DO

    enddo chan_loop

    fweight%mean_f = fweight%mean_f/(dEdP)

    if(inputs%verbose.ge.1) then
        write(*,*) 'write fida weights:    ' , time(time_start)
    endif

#ifdef _MPI
    call parallel_sum(fweight%weight)
    call parallel_sum(fweight%mean_f)
    if(my_rank().eq.0) call write_fida_weights()
#else
    call write_fida_weights()
#endif

end subroutine fida_weights_los