Исключения с плавающей точкой для основных тригонометрических уравнений

Question

У меня возникли проблемы с исключениями с плавающей запятой в подпрограмме, которая конвертирует из географических в геоцентрические координаты. Переменная geo(x) подается в подпрограмму в виде пар широты и долготы. Переменная xyz(x) выводится как триплет компонентов (x - гринвичский меридиан и экватор; y - 90 градусов долготы и экватора; z -- Северный полюс).

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

Я понимаю что sind а также cosd являются нестандартными функциями, но они связаны с объектным файлом с соответствующими преобразованиями. Я проверял это раньше, и он работает для других кодов, которые используют градусы вместо радиан:

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

Я попытался пройти программу с помощью GDB, и не смог понять проблему. Кажется, что все переменные в уравнении в порядке, но xyz(1) = 1, что не так, если я делаю этот расчет с помощью калькулятора.

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)

Что мне здесь не хватает, что вызывает исключение с плавающей запятой? Я уверен, что это довольно просто, я очень новичок в этом.

Answer 1

Функции sind и cosd возвращают REAL (одинарная точность) с плавающей точкой.

Но subroutine geo2xyz понятия не имею об этом.
С IMPLICIT REAL*8 предполагается, что эти функции возвращают числа с плавающей запятой двойной точности.

Следовательно, вы должны правильно объявить sind и cosd как реальные, прежде чем использовать их, и / или изменить их тип возвращаемого значения на real*8.

Самое смешное, что эта программа могла работать на старой архитектуре x86, потому что возвращаемое значение с одинарной точностью было увеличено до двойного в регистре ST0. Он не будет работать на 64-битной платформе Intel, которая использует SSE2 и регистрирует XMM0: в такой архитектуре нет возможности удвоить ставку.