BiotSavart_VectorPotential Subroutine

SUBROUTINE BiotSavart_VectorPotential(n_positions, xyz, n_points, n_segments, coil_points, ehat, L, prefactor, A)
    ! MODULES
    USE MODgvec_Globals, ONLY: TWOPI,CROSS,wp
    ! INPUT VARIABLES -----------------------------------------------------------------------------------------------------------!
    INTEGER,  INTENT(IN) :: n_positions               !! number of positions where the magnetic field is evaluated
    REAL(wp), INTENT(IN) :: xyz(3,n_positions)        !! x,y,z positions where the magnetic field is evaluated
    INTEGER,  INTENT(IN) :: n_points                  !! number of points representing the segmented coil
    INTEGER,  INTENT(IN) :: n_segments                !! number of coil segments
    REAL(wp), INTENT(IN) :: coil_points(3,n_points)   !! x,y,z positions of the segment chain (=polygon), n_segments=n_points-1
    REAL(wp), INTENT(IN) :: prefactor                 !! factor applied to the final magnetic field
    REAL(wp), INTENT(IN) :: ehat(3,n_segments)         !! unit vector along segment
    REAL(wp), INTENT(IN) :: L(n_segments)             !! segment length
    ! OUTPUT VARIABLES ----------------------------------------------------------------------------------------------------------!
    REAL(wp), INTENT(OUT):: A(3,n_positions)          !! magnetic field evaluated at xyz positions.
    ! LOCAL VARIABLES -----------------------------------------------------------------------------------------------------------!
    INTEGER              :: iPosition,iSegment
    REAL(wp)             :: R_i(3)                    !! vector from current position to starting point of segment
    REAL(wp)             :: R_f(3)                    !! vector from current position to end point of segment (=starting point of next segment)
    REAL(wp)             :: norm_ehat
    REAL(wp)             :: mod_R_i,mod_R_f
    REAL(wp)             :: f_epsilon
    ! CODE ----------------------------------------------------------------------------------------------------------------------!
    A = 0.0_wp
    !$OMP PARALLEL DO     &
    !$OMP SCHEDULE(STATIC) DEFAULT(PRIVATE) SHARED(xyz, coil_points, n_positions, n_segments, L, ehat) &
    !$OMP REDUCTION(+:A)
    DO iPosition=1,n_positions
        R_f = xyz(:,iPosition)-coil_points(:,1)
        mod_R_f= SQRT(SUM(R_f*R_f))
        DO iSegment=1,n_segments
            R_i = R_f
            mod_R_i = mod_R_f
            R_f = xyz(:,iPosition)-coil_points(:,iSegment+1)
            mod_R_f = SQRT(SUM(R_f*R_f))
            f_epsilon = LOG( (mod_R_i + mod_R_f + L(iSegment)) / (mod_R_i + mod_R_f - L(iSegment))  )
            A(:,iPosition) = A(:,iPosition) + ehat(:,iSegment) * f_epsilon
        END DO !iSegment
    END DO !iPosition
    !$OMP END PARALLEL DO
    A = prefactor*A
END SUBROUTINE BiotSavart_VectorPotential