I've implemented code to find the polar coordinates of a point in 2D space. if the point lies in the 1st or 2nd Quadrant, 0<=theta<=pi and if it lies in the 3rd or 4th Quadrant, -pi <= theta <= 0.
module thetalib
contains
real function comp_theta( x1, x2)
implicit none
real , intent(in) :: x1, x2
real :: x1p, x2p
real :: x1_c=0.0, x2_c=0.0
real :: pi=4*atan(1.0)
x1p = x1 - x1_c
x2p = x2 - x2_c
! - Patch
!if ( x1p == 0 .and. x2p /= 0 ) then
! comp_theta = sign(pi/2.0, x2p)
!else
! comp_theta = atan ( x2p / x1p )
!endif
comp_theta = atan( x2p / x1p)
if ( x1p >= 0.0 .and. x2p >= 0.0 ) then
comp_theta = comp_theta
elseif ( x1p < 0 .and. x2p >= 0.0 ) then
comp_theta = pi + comp_theta
elseif( x1p < 0.0 .and. x2p < 0.0 ) then
comp_theta = -1* (pi - comp_theta)
elseif ( x1p >= 0.0 .and. x2p < 0.0 ) then
comp_theta = comp_theta
endif
return
end function comp_theta
end module thetalib
program main
use thetalib
implicit none
! Quadrant 1
print *, "(0.00, 1.00): ", comp_theta(0.00, 1.00)
print *, "(1.00, 0.00): ", comp_theta(1.00, 0.00)
print *, "(1.00, 1.00): ", comp_theta(1.00, 1.00)
! Quadrant 2
print *, "(-1.00, 1.00): ", comp_theta(-1.00, 1.00)
print *, "(-1.00, 0.00): ", comp_theta(-1.00, 0.00)
! Quadrant 3
print *, "(-1.00, -1.00): ", comp_theta(-1.00, -1.00)
! Quadrant 4
print *, "(0.00, -1.00): ", comp_theta(0.00, -1.00)
print *, "(1.00, -1.00): ", comp_theta(1.00, -1.00)
end program main
In the function thetalib::comp_theta, when there is a division by zero and the numerator is +ve, fortran evaluates it to be -infinity and when the numerator is -ve, it evaluates it to be +infinity ( see output )
(0.00, 1.00): -1.570796
(1.00, 0.00): 0.0000000E+00
(1.00, 1.00): 0.7853982
(-1.00, 1.00): 2.356194
(-1.00, 0.00): 3.141593
(-1.00, -1.00): -2.356194
(0.00, -1.00): 1.570796
(1.00, -1.00): -0.7853982
This baffled me. I've also implemented the patch you see to work around it. And to investigate it further, I setup a small test:
program main
implicit none
real :: x1, x2
x1 = 0.0 - 0.0 ! Reflecting the x1p - 0.0
x2 = 1.0
write(*,*) "x2/x1=", x2/x1
x2 = -1.0
write(*,*) "x2/x1=", x2/x1
end program main
This evaluates to:
x2/x1= Infinity
x2/x1= -Infinity
My fortran version:
$ ifort --version
ifort (IFORT) 19.0.1.144 20181018
Copyright (C) 1985-2018 Intel Corporation. All rights reserved.
And I have three questions:
infinity take the signs shown in outputs for both thetalib::comp_theta and the test program?That there are signed infinite values follows from the compiler supporting IEEE arithmetic with the real type.
For motivation, consider real non-zero numerator and denominator. If these are both of the same sign then the quotient is a real (finite) positive number. If they are of opposite sign the quotient is a real (finite) negative number.
Consider the limit 1/x as x tends to zero from below. For any strictly negative value of x the value is negative. For continuity considerations the limit can be taken to be negative infinity.
So, when the numerator is non-zero, the quotient will be positive infinity if the numerator and denominator are of the same sign, and negative if of opposite sign. Recall also, that the zero denominator may be signed.
If you wish to examine the number, to see whether it is finite you can use the procedure IEEE_IS_FINITE of the intrinsic module ieee_arithmetic. Further, that module has the procedure IEEE_CLASS which provides useful information about its argument. Among other things: