Calculates the intersection arclength of a circle with the beam_grid
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=Float64), | intent(in), | dimension(3) | :: | r0 | Position of center enter of the circle in beam grid coordinates [cm] |
|
| real(kind=Float64), | intent(in), | dimension(3) | :: | e1 | Unit vector pointing towards (R, 0) (r,phi) position of the circle in beam grid coordinates |
|
| real(kind=Float64), | intent(in), | dimension(3) | :: | e2 | Unit vector pointing towards (R, pi/2) (r,phi) position of the circle in beam grid coordinates |
|
| real(kind=Float64), | intent(in) | :: | radius | Radius of circle [cm] |
||
| real(kind=Float64), | intent(out) | :: | beam_grid_phi_enter | Phi value where the circle entered the beam_grid [rad] |
||
| real(kind=Float64), | intent(out) | :: | beam_grid_phi_exit | Phi value where the circle exits the beam_grid [rad] |
subroutine circle_grid_intersect(r0, e1, e2, radius, beam_grid_phi_enter, beam_grid_phi_exit)
!+ Calculates the intersection arclength of a circle with the [[libfida:beam_grid]]
real(Float64), dimension(3), intent(in) :: r0
!+ Position of center enter of the circle in beam grid coordinates [cm]
real(Float64), dimension(3), intent(in) :: e1
!+ Unit vector pointing towards (R, 0) (r,phi) position of the circle in beam grid coordinates
real(Float64), dimension(3), intent(in) :: e2
!+ Unit vector pointing towards (R, pi/2) (r,phi) position of the circle in beam grid coordinates
real(Float64), intent(in) :: radius
!+ Radius of circle [cm]
real(Float64), intent(out) :: beam_grid_phi_enter
!+ Phi value where the circle entered the [[libfida:beam_grid]] [rad]
real(Float64), intent(out) :: beam_grid_phi_exit
!+ Phi value where the circle exits the [[libfida:beam_grid]] [rad]
real(Float64), dimension(3) :: i1_p,i1_n,i2_p,i2_n
real(Float64), dimension(4) :: d
real(Float64), dimension(6) :: p, gams
real(Float64), dimension(4,6) :: phi
logical, dimension(4,6) :: inter
integer, dimension(6) :: n
integer :: i
real(Float64) :: alpha,beta,delta,sinx1,cosx1,sinx2,cosx2,tmp
real(Float64) :: tol = 1.0d-10
logical :: r0_ing
p = [beam_grid%xmin, beam_grid%xmax, &
beam_grid%ymin, beam_grid%ymax, &
beam_grid%zmin, beam_grid%zmax ]
n = [1, 1, 2, 2, 3, 3]
inter = .False.
phi = 0.d0
r0_ing = in_grid(r0)
do i=1,6
alpha = e2(n(i))
beta = e1(n(i))
if((alpha.eq.0.0).and.(beta.eq.0.0)) cycle
gams(i) = (p(i) - r0(n(i)))/radius
delta = alpha**4 + (alpha**2)*(beta**2 - gams(i)**2)
if (delta.ge.0.0) then
cosx1 = (gams(i)*beta + sqrt(delta))/(alpha**2 + beta**2)
if ((cosx1**2).le.1.0) then
sinx1 = sqrt(1 - cosx1**2)
i1_p = r0 + radius*cosx1*e1 + radius*sinx1*e2
i1_n = r0 + radius*cosx1*e1 - radius*sinx1*e2
if (approx_eq(i1_p(n(i)),p(i),tol).and.in_grid(i1_p)) then
inter(1,i) = .True.
phi(1,i) = atan2(sinx1,cosx1)
endif
if (approx_eq(i1_n(n(i)),p(i),tol).and.in_grid(i1_n)) then
inter(2,i) = .True.
phi(2,i) = atan2(-sinx1,cosx1)
endif
endif
if(delta.gt.0.0) then
cosx2 = (gams(i)*beta - sqrt(delta))/(alpha**2 + beta**2)
if ((cosx2**2).le.1.0) then
sinx2 = sqrt(1 - cosx2**2)
i2_p = r0 + radius*cosx2*e1 + radius*sinx2*e2
i2_n = r0 + radius*cosx2*e1 - radius*sinx2*e2
if (approx_eq(i2_p(n(i)),p(i),tol).and.in_grid(i2_p)) then
inter(3,i) = .True.
phi(3,i) = atan2(sinx2,cosx2)
endif
if (approx_eq(i2_n(n(i)),p(i),tol).and.in_grid(i2_n)) then
inter(4,i) = .True.
phi(4,i) = atan2(-sinx2,cosx2)
endif
endif
endif
endif
enddo
beam_grid_phi_enter = 0.d0
beam_grid_phi_exit = 0.d0
if (count(inter).gt.2) then
write(*,'("CIRCLE_GRID_INTERSECT: Circle intersects grid more than 2 times: ",i2)') count(inter)
return
endif
if(any(inter)) then
beam_grid_phi_enter = minval(phi,inter)
beam_grid_phi_exit = maxval(phi,inter)
if(r0_ing.and.any(count(inter,1).ge.2)) then
if((beam_grid_phi_exit - beam_grid_phi_enter) .lt. pi) then
tmp = beam_grid_phi_enter
beam_grid_phi_enter = beam_grid_phi_exit
beam_grid_phi_exit = tmp + 2*pi
endif
else
if((beam_grid_phi_exit - beam_grid_phi_enter) .gt. pi) then
tmp = beam_grid_phi_enter
beam_grid_phi_enter = beam_grid_phi_exit
beam_grid_phi_exit = tmp + 2*pi
endif
endif
if(approx_eq(beam_grid_phi_exit-beam_grid_phi_enter,pi,tol).and.r0_ing) then
beam_grid_phi_enter = 0.0
beam_grid_phi_exit = 2*pi
endif
else
if(r0_ing) then
call grid_intersect(r0, e1, tmp, i1_n,i1_p)
call grid_intersect(r0, e2, tmp, i2_n,i2_p)
d(1) = norm2(r0 - i1_n)/radius
d(2) = norm2(r0 - i1_p)/radius
d(3) = norm2(r0 - i2_n)/radius
d(4) = norm2(r0 - i2_p)/radius
if(all(d.ge.1.0)) then
beam_grid_phi_enter = 0.d0
beam_grid_phi_exit = 2.d0*pi
endif
endif
endif
end subroutine circle_grid_intersect