calculate_rotation_angles Subroutine

  SUBROUTINE calculate_rotation_angles(X,bpa,apa)
    IMPLICIT NONE
    REAL(rp), DIMENSION(3), INTENT(IN) :: X
    LOGICAL, DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: bpa
    REAL(rp), DIMENSION(:,:,:), ALLOCATABLE, INTENT(INOUT) :: apa
    REAL(rp) :: R, D, psi
    REAL(rp) :: a, b, c, dis, xp, xn
    REAL(rp) :: xtmp, ytmp
    INTEGER :: ii,jj
    ! bpa(:,:,1) -- > xp
    ! bpa(:,:,2) -- > xn

    R = SQRT(SUM(X(1:2)**2))
    D = SQRT( (cam%position(1) - X(1))**2 + X(2)**2 )
    psi = -ATAN2(cam%position(2) - X(3),D)

    bpa = .TRUE.

    do ii=1_idef,cam%num_pixels(1)
       a = 1.0_rp + TAN(ang%beta(ii))**2
       b = -2.0_rp*TAN(ang%beta(ii))**2*cam%position(1)
       c = (TAN(ang%beta(ii))*cam%position(1))**2 - R**2
       dis = b**2 - 4.0_rp*a*c

       if (dis.GT.0.0_rp) then
          do jj=1_idef,cam%num_pixels(2)

             if ((psi.GE.ang%psi(jj)).AND.(psi.LT.ang%psi(jj+1_idef))) then
                xp = 0.5_rp*(-b + SQRT(dis))/a
                xn = 0.5_rp*(-b - SQRT(dis))/a

                xtmp = xp - cam%position(1)
                ytmp = SQRT(R**2 - xp**2)

                ! Check if particle is behind inner wall
                if ((ATAN2(ytmp,xtmp).GT.ang%threshold_angle).AND. &
                     (SQRT(xtmp**2+ytmp**2).GT.ang%threshold_radius)) then
                   bpa(ii,jj,1) = .FALSE.
                else
                   apa(ii,jj,1) = ATAN2(ytmp,xp)
                   if (apa(ii,jj,1).LT.0.0_rp) apa(ii,jj,1) = &
                        apa(ii,jj,1) + 2.0_rp*C_PI
                end if

                xtmp = xn - cam%position(1)
                ytmp = SQRT(R**2 - xn**2)

                ! Check if particle is behind inner wall
                if ((ATAN2(ytmp,xtmp).GT.ang%threshold_angle).AND. &
                     (SQRT(xtmp**2+ytmp**2).GT.ang%threshold_radius)) then
                   bpa(ii,jj,2) = .FALSE.
                else
                   apa(ii,jj,2) = ATAN2(ytmp,xn)
                   if (apa(ii,jj,2).LT.0.0_rp) apa(ii,jj,2) = &
                        apa(ii,jj,2) + 2.0_rp*C_PI
                end if
             else ! Not in pixel (ii,jj)
                bpa(ii,jj,:) = .FALSE.
             end if ! Check if in pixel (ii,jj)

          end do ! NY
       else ! no real solutions
          bpa(ii,:,:) = .FALSE.
       end if ! Checking discriminant
    end do !! NX
  END SUBROUTINE calculate_rotation_angles