I hope there is something simple here that I am missing. I have a large project (in high energy physics) that I am retooling to work for a different process and I hit an unexpected roadblock. The details are unimportant, but it seems that when I pass two four-vectors to a DOUBLE COMPLEX function, the second four-vector gets corrupted.
I have simplified the problem to the following sample code:
program complex
double precision p1(4),p2(4)
double complex d0
external sp
p1(1) = -426.6404106d0
p1(2) = -109.0232008d0
p1(3) = -323.4625454d0
p1(4) = +572.4883637d0
p2(1) = +402.2554234d0
p2(2) = -2528.7117650d0
p2(3) = +1250.9838540d0
p2(4) = +2854.8835400d0
write(6,100) p1(1), p1(2), p1(3), p1(4)
write(6,100) p2(1), p2(2), p2(3), p2(4)
write(6,*)
d0 = sp( p1, p2 )
100 format(' p(xyzE)', 4g18.10)
stop
end program complex
double complex function sp( p1, p2 )
double precision p1(4), p2(4)
data x/1/,
& y/2/,
& z/3/,
& E/4/
write(6,*)
write(6,101) p1(x), p1(y), p1(z), p1(E)
write(6,101) p2(x), p2(y), p2(z), p2(E)
write(6,*)
101 format(' p(xyzE)', 4g18.10)
sp = p2(z)*( p1(E) - p1(x) )
& - p1(z)*( p2(E) - p2(x) )
& - (0.d0,1.d0)*( p1(y)*p2(x)+p1(E)*p2(y)
& - p1(y)*p2(E)-p1(x)*p2(y) )
sp = sp / sqrt( ( p1(E)-p1(x) )*( p2(E)-p2(x) ) )
write(6,*) sp
write(6,*)
return
end function sp
This example simply tried to pass the two four-vectors to the FUNCTION SP and then shows that the p1 and p2 seem to be switched around inside and the program crashes.
I have a feeling that this is something simple yet subtle but I would appreciate any advice you could give. I know when I make SP a DOUBLE PRECISION FUNCTION the values are passed correctly, but I need the function to return a complex number.
