Repository navigation
solution of LOS of passive grid, calling get_plasma_extrema #301
Description
Calling get_plasma_extrema is called after if conditions for various combinations of requests for NPA only, FIDA only or both. In case that one of them is not requested and well defined, the plasma extrema are searched for parameters of not requested diagnostic. That could cause issues if unrequested diagnostic geometry is not realistic. Therefore it is better to call get_plasma_extrema inside the if cycles.
The LOS of passive grid for NPA is defined as a norm of detector plane. Once the detector is tilted, it will lead to the selection of simulated area of plasma which is out the region of NPA signal origin. The better approach is to define the NPA LOS of passive grid as a vector between detector and aperture centers.
subroutine make_passive_grid
!+ Makes [[libfida:pass_grid] from user defined inputs
real(Float64), dimension(3,spec_chords%nchan+npa_chords%nchan) :: r0_arr, v0_arr, **a_cent**
real(Float64), dimension(2,spec_chords%nchan+npa_chords%nchan) :: xy_enter, xy_exit
logical, dimension(8+2*(spec_chords%nchan+npa_chords%nchan)) :: yle, ygt
logical, dimension(spec_chords%nchan+npa_chords%nchan) :: skip
real(Float64), dimension(:), allocatable :: xarr, yarr
real(Float64), dimension(3,8) :: vertices_xyz, vertices_uvw
real(Float64), dimension(2,3) :: extrema
real(Float64), dimension(3,3) :: xyz_axis
real(Float64), dimension(3) :: r0, vi
real(Float64), dimension(8) :: xarr_beam_grid, yarr_beam_grid
real(Float64) :: xmin, ymin, xmax, ymax, zmin, zmax, max_length
real(Float64) :: dlength = 3.0 !cm
integer :: i, iin, iout, dim_le, dim_gt, dim, **j**
logical :: inp, phi_pos
!! Get beam grid boundaries
! Convert vertices_xyz to vertices_uvw
! Note: vertices_xyz has the following coordinate definitions:
! 111 = 1, 222 = 2, 112 = 3, 211 = 4, 121 = 5, 212 = 6, 122 = 7, 221 = 8
xmin = beam_grid%xmin ; xmax = beam_grid%xmax
ymin = beam_grid%ymin ; ymax = beam_grid%ymax
zmin = beam_grid%zmin ; zmax = beam_grid%zmax
!Initialize minimum and maximum vertices
vertices_xyz(1,1) = xmin ; vertices_xyz(2,1) = ymin ; vertices_xyz(3,1) = zmin
vertices_xyz(1,2) = xmax ; vertices_xyz(2,2) = ymax ; vertices_xyz(3,2) = zmax
!Initialize
vertices_xyz(:,3) = vertices_xyz(:,1)
vertices_xyz(:,4) = vertices_xyz(:,1)
vertices_xyz(:,5) = vertices_xyz(:,1)
!Update
vertices_xyz(3,3) = zmax
vertices_xyz(1,4) = xmax
vertices_xyz(2,5) = ymax
!Initialize
vertices_xyz(:,7) = vertices_xyz(:,2)
vertices_xyz(:,8) = vertices_xyz(:,2)
vertices_xyz(:,6) = vertices_xyz(:,2)
!Update
vertices_xyz(1,7) = xmin
vertices_xyz(3,8) = zmin
vertices_xyz(2,6) = ymin
do i=1, 8
call xyz_to_uvw(vertices_xyz(:,i),vertices_uvw(:,i))
xarr_beam_grid(i) = vertices_uvw(1,i)
yarr_beam_grid(i) = vertices_uvw(2,i)
enddo
!! Next consider passive diagnostic extrema relative to the plasma
if((inputs%calc_pfida.gt.0).and.(inputs%calc_pnpa.gt.0)) then
do i=1,(spec_chords%nchan)
r0_arr(:,i) = spec_chords%los(i)%lens_uvw
v0_arr(:,i) = spec_chords%los(i)%axis_uvw
enddo
do i=1, npa_chords%nchan
**call xyz_to_uvw(npa_chords%det(i)%detector%origin, r0_arr(:,i))
call xyz_to_uvw(npa_chords%det(i)%aperture%origin, a_cent(:,i))**
xyz_axis = npa_chords%det(i)%detector%basis
**v0_arr(:,i) = a_cent(:,i) - r0_arr(:,i)**
enddo
**call get_plasma_extrema(r0_arr,v0_arr,extrema,xarr_beam_grid,yarr_beam_grid)**
else if(inputs%calc_pfida.gt.0) then
do i=1, spec_chords%nchan
r0_arr(:,i) = spec_chords%los(i)%lens_uvw
v0_arr(:,i) = spec_chords%los(i)%axis_uvw
enddo
**call get_plasma_extrema(r0_arr(:,:spec_chords%nchan),v0_arr(:,:spec_chords%nchan),extrema,xarr_beam_grid,yarr_beam_grid)**
else !pnpa>=1 case
do i=1, npa_chords%nchan
call xyz_to_uvw(npa_chords%det(i)%detector%origin, r0_arr(:,i))
**call xyz_to_uvw(npa_chords%det(i)%aperture%origin, a_cent(:,i))**
xyz_axis = npa_chords%det(i)%detector%basis
**v0_arr(:,i) = a_cent(:,i) - r0_arr(:,i)**
enddo
**call get_plasma_extrema(r0_arr(:,:npa_chords%nchan),v0_arr(:,:npa_chords%nchan),extrema,xarr_beam_grid,yarr_beam_grid)**
endif
!! Store the passive neutral grid
pass_grid%dr = inter_grid%dr
pass_grid%dz = inter_grid%dz
pass_grid%nr = inter_grid%nr
pass_grid%nz = inter_grid%nz
allocate(pass_grid%r(pass_grid%nr), pass_grid%z(pass_grid%nz))
pass_grid%r = inter_grid%r
pass_grid%z = inter_grid%z
pass_grid%da = pass_grid%dr*pass_grid%dz
pass_grid%dphi = 2*pi/100 !TODO: make this user input
pass_grid%nphi = int(ceiling((extrema(2,3)-extrema(1,3))/pass_grid%dphi))
allocate(pass_grid%phi(pass_grid%nphi))
do i=1, pass_grid%nphi
pass_grid%phi(i) = extrema(1,3) + (i-1)*pass_grid%dphi
enddo
pass_grid%dv = pass_grid%dr*pass_grid%dphi*pass_grid%dz
pass_grid%dims = [pass_grid%nr, pass_grid%nz, pass_grid%nphi]
pass_grid%ntrack = (pass_grid%nr+pass_grid%nz+pass_grid%nphi)*inputs%max_cell_splits
pass_grid%ngrid = pass_grid%nr*pass_grid%nz*pass_grid%nphi
if(inputs%verbose.ge.1) then
write(*,'(a)') "---- Passive grid settings ----"
write(*,'(T2,"Nr: ", i3)') pass_grid%nr
write(*,'(T2,"Nz: ", i3)') pass_grid%nz
write(*,'(T2,"Nphi: ", i3)') pass_grid%nphi
write(*,'(T2,"R range = [",f6.2,",",f6.2,"]")') &
pass_grid%r(1),pass_grid%r(pass_grid%nr)
write(*,'(T2,"Z range = [",f7.2,",",f6.2,"]")') &
pass_grid%z(1),pass_grid%z(pass_grid%nz)
write(*,'(T2,"Phi range = [",f5.2,",",f5.2,"]")') &
pass_grid%phi(1),pass_grid%phi(pass_grid%nphi)
write(*,'(T2,"dA: ", f5.2," [cm^3]")') pass_grid%da
write(*,*) ''
endif
end subroutine make_passive_grid
Activity
lstagner commentedon Jul 30, 2025
This seems like a reasonable solution. You would need to have a check to make sure the detector norm vector is in roughly the same direction as the new vector. But that shouldn't happen if the npa detector is inputted correctly...so you can probably ignore that. Why don't you implement the change, make sure it works and then open a PR
JoseRuedaRueda commentedon Aug 14, 2025
If this resolved?
I am doing passive INPA sims, and when I plot the initial position of the markers, it makes no sense, it cut in z (my grid is not cut there)
JoseRuedaRueda commentedon Aug 14, 2025
Update, indeed, this solves the issue. It also reduce like a factor 10 the number of particles which miss the detector (Issue #304 )
lstagner commentedon Aug 20, 2025
If this fixes the issue could either @JoseRuedaRueda or @KlaraBogar make a pull request.
JoseRuedaRueda commentedon Aug 25, 2025
Added to #299
prechelg commentedon Oct 22, 2025
Since this PR was merged in already, can we close this @KlaraBogar @JoseRuedaRueda