I use the same compiler (gfortran) and computer every time I compile it and I don't change the source code, but on different compilations it gives a completely different outcome. Sometimes a reasonable outcome and sometimes a wrong outcome.
For instance I compiled it and got the next very unreasonable outputs.
0.0000000000000000 , 3.0902604013843218E+049 , 3.0902604013843218E+049 , 3.0902604013843218E+049 , 5.3524880238158376E+049
2.0000000000000000 , -6.1610238730665058E+049 , -6.1610238730665058E+049 , -6.1610238730665058E+049 , 1.0671206374795975E+050
4.0000000000000000 , 5.5751160679618236E+049 , 5.5751160679618236E+049 , 5.5751160679618236E+049 , 9.6563842878035016E+049
6.0000000000000000 , -1.0179425493222214E+049 , -1.0179425493222214E+049 , -1.0179425493222214E+049 , 1.7631282146122754E+049
8.0000000000000000 , 2.4002992709421553E+049 , 2.4002992709421553E+049 , 2.4002992709421553E+049 , 4.1574402906423475E+049
10.000000000000000 , 3.5499818567499908E+049 , 3.5499818567499908E+049 , 3.5499818567499908E+049 , 6.1487489418386840E+049
12.000000000000000 , -3.5465339877133604E+049 , -3.5465339877133604E+049 , -3.5465339877133604E+049 , 6.1427770574893967E+049
14.000000000000000 , 3.7523505062483277E+049 , 3.7523505062483277E+049 , 3.7523505062483277E+049 , 6.4992617246289011E+049
Then without changing anything I recompiled the same code and ran it again and got the more reasonable first result:
0.0000000000000000 , -0.20075827532679802 , 5.7540609466025759E-003 , 0.33972754855979093 , 0.39465402770022856
Why does this happen? Is it a problem in my code or the compiler. I leave the code just in case its useful to answer the question. I'm just starting to learn fortran, and sorry for the comments in spanish.
EDIT: In this program I use fgsl, a fortran interface to the gnu scientific library, it can be found here: https://github.com/reinh-bader/fgsl
program trace_of_product
use fgsl
implicit none
real(fgsl_double) :: p1, p2, p3, t, r, omega, tmax, rad, zav, yav, xav, integ1, integ2, integ3
integer(fgsl_int) :: i, j, k, nin, nt, jmax
complex(fgsl_double) :: wigner_func, z_sym, y_sym, x_sym
real(fgsl_double), parameter :: pi = 3.14159265359
character(len=3) :: filenumber
omega = 0.01
nin = 500
nt = 50
jmax = 9
tmax = 100
r = sqrt(1.0d0*jmax)
write(filenumber, '(I0.3)') jmax
open(unit = 1, file = 'data'//trim(filenumber)//'.csv')
do k = 0,nt
t = k*tmax/nt
zav = 0.0d0
yav = 0.0d0
xav = 0.0d0
do j = 1,jmax,2
integ1 = 0.0d0
integ2 = 0.0d0
integ3 = 0.0d0
!
! Esta parte del código calcula la integral del producto de la
! j-función de Wigner y el j-simbolo de los operadores de posición en la región
! [0,2pi)x[0,pi)x[0,4pi) y dV = sin(theta) d(phi) d(theta) d(psi)
!
do i = 0,nin
p1 = 2*pi*rand()
p2 = pi*rand()
p3 = 4*pi*rand()
integ1 = integ1 + realpart(wigner_func(j,r,t,omega,p1,p2,p3))*realpart(z_sym(j,t,omega,p1,p2,p3))&
*sin(p2)
integ2 = integ2 + realpart(wigner_func(j,r,t,omega,p1,p2,p3))*realpart(y_sym(j,t,omega,p1,p2,p3))&
*sin(p2)
integ3 = integ3 + realpart(wigner_func(j,r,t,omega,p1,p2,p3))*realpart(x_sym(j,t,omega,p1,p2,p3))&
*sin(p2)
end do
integ1 = integ1*pi*(j+1)/nin
integ2 = integ2*pi*(j+1)/nin
integ3 = integ3*pi*(j+1)/nin
zav = zav +integ1
yav = yav +integ2
xav = xav +integ3
end do
rad = sqrt(xav**2+yav**2+zav**2)
write(1,*) t,',',xav,',',yav,',',zav,',',rad
write(*,*) t,',',xav,',',yav,',',zav,',',rad
end do
end program
!
! Esta función calcula la j-función de Wigner para un estado coherente de dos
! modos |ab>=|a>|b>, rho = |ab><ab|. N y J están restingidos a que N+J=j
!
function wigner_func(jc,r,t,omega,phi,theta,psi) result(Wp)
use fgsl
implicit none
real(fgsl_double), intent(in) :: r, t, omega, phi, theta, psi
integer(fgsl_int), intent(in) :: jc
complex(fgsl_double) :: wigner_D, Wp
real(fgsl_double) :: wigner_small_d, cg
integer(fgsl_int) :: m, mp, k, l, q, N, J
real(fgsl_double), parameter :: pi = 3.14159265359
do N = 0,jc
J = jc-N
do k = abs(N-J),jc
do m = -J,J
do mp = -N,N
do l = -k,k
do q = -k,k
Wp = Wp + exp(-r**2) * r**(2*jc) * (2*k+1) * sin(pi/8)**(jc-m-mp) * cos(pi/8)**(jc+m+mp)* &
cg(N,mp,J,m,k,l)*wigner_D(k,l,q,phi+omega*t,theta,psi)*wigner_small_d(k,q,N-J,pi/2) / &
sqrt(Gamma(N-mp+1.0d0)*Gamma(J-m+1.0d0)*Gamma(N+mp+1.0d0)*Gamma(J+m+1.0d0)*(jc+1)*(2*N+1))
end do
end do
end do
end do
end do
end do
end function wigner_func
!
! Esta función calcula el j-simbolo del operador z, j = 2l + 1
!
function z_sym(j,t,omega,phi,theta,psi) result(Wz)
use fgsl
implicit none
real(fgsl_double), intent(in) :: t, omega, phi, theta, psi
integer(fgsl_int), intent(in) :: j
complex(fgsl_double) :: Wigner_D, Wz
real(fgsl_double) :: wigner_small_d, cg
integer(fgsl_int) :: l,m,k,qp,q
real(fgsl_double), parameter :: pi = 3.14159265359
if ( mod(j,2) == 1 ) then
l = (j-1)/2
do m = -l,l
do k = 1,j
do qp = -k,k
do q = -k,k
Wz = Wz + sqrt(((l+1.0d0)**2-m**2)/((4*(l+1)**2-1)*(j+1)))*(2*k+1)*&
Wigner_D(k,qp,q,phi+omega*t,theta,psi)* &
(cg(l+1,m,l,m,k,qp)*wigner_small_d(k,q,1,pi/2)/sqrt(j+1.0d0)+&
cg(l,m,l+1,m,k,qp)*wigner_small_d(k,q,-1,pi/2)/sqrt(1.0d0*j))
end do
end do
end do
end do
end if
end function z_sym
!
! Esta función calcula el j-simbolo del operador y, j = 2l + 1
!
function y_sym(j,t,omega,phi,theta,psi) result(Wy)
use fgsl
implicit none
real(fgsl_double), intent(in) :: t, omega, phi, theta, psi
integer(fgsl_int), intent(in) :: j
complex(fgsl_double) :: Wigner_D, Wy
real(fgsl_double) :: wigner_small_d, cg
integer(fgsl_int) :: l, m, k, s, q
real(fgsl_double), parameter :: pi = 3.141596265359
if ( mod(j,2) == 1 ) then
l = (j-1)/2
do m = -l,l
do k = 1,j
do q = -k,k
do s = -k,k
Wy = Wy+cmplx(0,0.5)*sqrt((2*k+1)/((4*(l+1)**2-1)*(j+1.0d0)))*wigner_D(k,s,q,phi+omega*t,theta,psi)*&
(sqrt((l+m+1)*(l+m+2.0d0))*(wigner_small_d(k,q,-1,pi/2)*cg(l,m,l+1,m+1,k,s)/sqrt(1.0d0*j)- &
wigner_small_d(k,q,1,pi/2)*cg(l+1,m+1,l,m,k,s)/sqrt(j+2.0d0)) + sqrt((l-m)*(l-m+1.0d0))*(&
wigner_small_d(k,q,-1,pi/2)*cg(l,m+1,l+1,m,k,s)/sqrt(1.0d0*j)-wigner_small_d(k,q,1,pi/2)*&
cg(l+1,m,l,m+1,k,s)/sqrt(j+2.0d0)))
end do
end do
end do
end do
end if
end function
!
! Esta función calcula el j-simbolo del operador x, j = 2l + 1
!
function x_sym(j,t,omega,phi,theta,psi) result(Wx)
use fgsl
implicit none
real(fgsl_double), intent(in) :: t, omega, phi, theta, psi
integer(fgsl_int), intent(in) :: j
complex(fgsl_double) :: Wigner_D, Wx
real(fgsl_double) :: wigner_small_d, cg
integer(fgsl_int) :: l, m, k, s, q
real(fgsl_double), parameter :: pi = 3.141596265359
if ( mod(j,2) == 1 ) then
l = (j-1)/2
do m = -l,l
do k = 1,j
do q = -k,k
do s = -k,k
Wx = Wx+0.5*sqrt((2*k+1)/((4*(l+1)**2-1)*(j+1.0d0)))*wigner_D(k,s,q,phi+omega*t,theta,psi)*&
(-sqrt((l+m+1)*(l+m+2.0d0))*(wigner_small_d(k,q,-1,pi/2)*cg(l,m,l+1,m+1,k,s)/sqrt(1.0d0*j)+ &
wigner_small_d(k,q,1,pi/2)*cg(l+1,m+1,l,m,k,s)/sqrt(j+2.0d0)) + sqrt((l-m)*(l-m+1.0d0))*(&
wigner_small_d(k,q,-1,pi/2)*cg(l,m+1,l+1,m,k,s)/sqrt(1.0d0*j)+wigner_small_d(k,q,1,pi/2)*&
cg(l+1,m,l,m+1,k,s)/sqrt(j+2.0d0)))
end do
end do
end do
end do
end if
end function
!
! Esta función calcula los coeficientes de Clebsch-Gordan
!
function cg(j,m,j1,m1,j2,m2) result(cgn)
use fgsl
implicit none
integer(fgsl_int), intent(in) :: j, m, j1, m1, j2, m2
real(fgsl_double) :: cgn
cgn = (-1)**(-j1+j2-m)*sqrt(2.0*j+1)*fgsl_sf_coupling_3j(2*j1,2*j2,2*j,2*m1,2*m2,-2*m)
end function cg
!
! Esta función es la d-función de Wigner
!
function wigner_small_d(j,m,mp,theta) result(d)
use fgsl
implicit none
integer(fgsl_int), intent(in) :: j, m, mp
integer(fgsl_int) :: k
real(fgsl_double), intent(in) :: theta
real(fgsl_double) :: d
do k = max(0,m-mp),min(j-mp,j+m)
d = d + (-1)**k * cos(theta/2)**(2*j+m-mp-2*k) * sin(theta/2)**(mp-m+2*k) / &
(Gamma(j-mp-k+1.0d0)*Gamma(k+1.0d0)*Gamma(mp-m+k+1.0d0)*Gamma(j+m-k+1.0d0))
end do
d = sqrt(Gamma(j+m+1.0d0)*Gamma(j-m+1.0d0)*Gamma(j+mp+1.0d0)*Gamma(j-mp+1.0d0))*d
end function wigner_small_d
!
! Esta función es la D-función de Wigner
!
function wigner_D(j,m,mp,phi,theta,psi) result(D)
use fgsl
implicit none
integer(fgsl_int), intent(in) :: j, m, mp
real(fgsl_double), intent(in) :: phi, theta, psi
real(fgsl_double) :: wigner_small_d
complex(fgsl_double) :: D
D = exp(-cmplx(0,1)*phi*m)*wigner_small_d(j,m,mp,theta)*exp(-cmplx(0,1)*psi*mp)
end function wigner_D
