PROGRAM var

! This program will illustrate some pathologies of floating point
!  arithmetic.

IMPLICIT NONE

INTEGER, PARAMETER :: real_kind=4

REAL (KIND=real_kind) :: p  ! this is some varying 'data' variable
REAL (KIND=real_kind) :: sum1, sum2 ! sums of first and second powers of p
REAL (KIND=real_kind) :: mean, std_dev, variance
REAL (KIND=real_kind) :: pi
INTEGER :: count   
INTEGER :: n

pi = 4._real_kind*atan(1.0_real_kind)

sum1 = 0.0
sum2 = 0.0
count = 0
DO n = 1,1280
  p = 1000._real_kind + sin(pi*float(n)/64._real_kind)
!  PRINT *,'n, p ',n,p
  sum1 = sum1 + p
  sum2 = sum2 + p**2
  count = count + 1
ENDDO

!analytically, variance = 0.5

mean = sum1 / real(count)
PRINT *,'mean ', mean
variance = sum2 / real(count) - mean**2
PRINT *,'variance ', variance
std_dev = sqrt(variance)
PRINT *,'std. deviation ', std_dev

END PROGRAM var
