I am having trouble with floating point exceptions in a subroutine that converts from geographic to geocentric coordinates. The variable geo(x) is fed into the subroutine as pairs of latitude and longitude. The variable xyz(x) is output as a triplet of components (x -- Greenwich meridian and equator; y -- 90 degree longitude and equator; z -- north pole).
subroutine geo2xyz(geo,xyz)
IMPLICIT REAL*8 (A-H,O-Z)
dimension geo(2),xyz(3)
rr=6367443.5
xyz(1)=rr*sind(90.-geo(1))*cosd(geo(2))
xyz(2)=rr*sind(90.-geo(1))*sind(geo(2))
xyz(3)=rr*cosd(90.-geo(1))
return
end
I realize that sind and cosd are non-standard functions, but they are linked to an object file with the appropriate conversions. I have tested this before and it works for other codes that use degrees instead of radians:
real function sind(x)
IMPLICIT REAL*8 (A-H,O-Z)
sind=sin(x*3.141592653589793d0/180.0d0)
return
end
real function cosd(x)
IMPLICIT REAL*8 (A-H,O-Z)
cosd=cos(x*3.141592653589793d0/180.0d0)
return
end
I tried going through the program with GDB, and couldn't figure out the problem. It seems that all the variables in the equation are OK, but xyz(1) = 1, which is not the case if I do that calculation with a calculator.
Program received signal SIGFPE, Arithmetic exception.
0x0000000100001cb8 in geo2xyz (geo=..., xyz=...) at disslip.f:758
758 xyz(1)=rr*sind(90.-geo(1))*cosd(geo(2))
(gdb) print xyz(1)
$1 = 1
(gdb) print rr
$2 = 6367443.5
(gdb) print geo(1)
$3 = 50.350000000000001
(gdb) print geo(2)
$4 = -127.55800000000001
(gdb)
What am I missing here that is causing the floating point exception? I'm sure it's pretty simple, I'm very new to this.
The sind and cosd functions do return a REAL (single precision) floating point.
But subroutine geo2xyz has no clue about that.
With the IMPLICIT REAL*8 it will assume that those functions return a double precision floating point.
Therefore, you must properly declare sind and cosd as real before using them, and/or change their return type to real*8.
The funny thing is that this program might have worked on old x86 architecture because the single precision value returned were promoted to double in the ST0 register. It won't work on 64 bit intel which is using SSE2 and register XMM0: there is no promotion to double in such architecture.