1
0
mirror of https://github.com/TREX-CoE/qmc-lttc.git synced 2024-10-14 20:12:12 +02:00
qmc-lttc/qmc_stats.f90

58 lines
1.2 KiB
Fortran
Raw Normal View History

2021-01-26 00:22:37 +01:00
subroutine ave_error(x,n,ave,err)
implicit none
2021-01-07 10:01:55 +01:00
integer, intent(in) :: n
double precision, intent(in) :: x(n)
double precision, intent(out) :: ave, err
2021-01-29 13:23:00 +01:00
double precision :: variance
2021-01-26 00:22:37 +01:00
if (n < 1) then
stop 'n<1 in ave_error'
2021-01-29 13:23:00 +01:00
2021-01-26 00:22:37 +01:00
else if (n == 1) then
2021-01-07 10:01:55 +01:00
ave = x(1)
err = 0.d0
2021-01-29 13:23:00 +01:00
2021-01-07 10:01:55 +01:00
else
2021-01-29 13:23:00 +01:00
ave = sum(x(:)) / dble(n)
variance = sum((x(:) - ave)**2) / dble(n-1)
err = dsqrt(variance/dble(n))
2021-01-07 10:01:55 +01:00
endif
end subroutine ave_error
2021-01-07 11:07:18 +01:00
subroutine random_gauss(z,n)
implicit none
integer, intent(in) :: n
double precision, intent(out) :: z(n)
double precision :: u(n+1)
double precision, parameter :: two_pi = 2.d0*dacos(-1.d0)
integer :: i
call random_number(u)
2021-01-29 13:23:00 +01:00
2021-01-07 11:07:18 +01:00
if (iand(n,1) == 0) then
! n is even
do i=1,n,2
z(i) = dsqrt(-2.d0*dlog(u(i)))
2021-01-12 10:55:00 +01:00
z(i+1) = z(i) * dsin( two_pi*u(i+1) )
z(i) = z(i) * dcos( two_pi*u(i+1) )
2021-01-07 11:07:18 +01:00
end do
2021-01-29 13:23:00 +01:00
2021-01-07 11:07:18 +01:00
else
! n is odd
do i=1,n-1,2
z(i) = dsqrt(-2.d0*dlog(u(i)))
2021-01-12 10:55:00 +01:00
z(i+1) = z(i) * dsin( two_pi*u(i+1) )
z(i) = z(i) * dcos( two_pi*u(i+1) )
2021-01-07 11:07:18 +01:00
end do
2021-01-29 13:23:00 +01:00
2021-01-07 11:07:18 +01:00
z(n) = dsqrt(-2.d0*dlog(u(n)))
2021-01-12 10:55:00 +01:00
z(n) = z(n) * dcos( two_pi*u(n+1) )
2021-01-29 13:23:00 +01:00
2021-01-07 11:07:18 +01:00
end if
2021-01-29 13:23:00 +01:00
2021-01-07 11:07:18 +01:00
end subroutine random_gauss