Calculates halo emission from user supplied neutrals file
subroutine halo_spec
!+ Calculates halo emission from user supplied neutrals file
integer :: ic, i, j, k, it, is
real(Float64), dimension(3) :: ri, vhalo, random3, rc
integer,dimension(3) :: ind
!! Determination of the CX probability
type(LocalProfiles) :: plasma
real(Float64) :: halo_photons, wght
logical :: inp
integer :: n = 10000
!$OMP PARALLEL DO schedule(dynamic,1) private(i,j,k,ic,ind,is, &
!$OMP& halo_photons, wght, rc, ri, inp, vhalo, random3, plasma)
loop_over_cells: do ic = istart, spec_chords%ncell, istep
call ind2sub(beam_grid%dims,spec_chords%cell(ic),ind)
i = ind(1) ; j = ind(2) ; k = ind(3)
halo_photons = neut%halo%dens(3,i,j,k)*tables%einstein(2,3)
if(halo_photons.le.0.0) cycle loop_over_cells
rc = [beam_grid%xc(i), beam_grid%yc(j), beam_grid%zc(k)]
!Find a point in cell that is also in the plasma
ri = rc
call in_plasma(ri, inp)
do while (.not.inp)
call randu(random3)
ri = rc + beam_grid%dr*(random3 - 0.5)
call in_plasma(ri,inp)
enddo
call get_plasma(plasma, pos=ri)
do is=1, n_thermal
wght = plasma%deni(is)/sum(plasma%deni)
do it=1, n
!! Halo Spectra
call mc_halo(plasma, thermal_mass(is), vhalo)
call store_photons(ri, vhalo, thermal_lambda0(is), wght*halo_photons/n, &
&spec%halo(:,:,:,is),spec%halostokes(:,:,:,:,is))
enddo
enddo
enddo loop_over_cells
!$OMP END PARALLEL DO
#ifdef _MPI
!! Combine Spectra
call parallel_sum(spec%halo)
#endif
end subroutine halo_spec