SHARE
TWEET

Untitled

a guest Jan 18th, 2019 56 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. program statistics
  2. implicit none
  3. integer :: NPoints = 10000
  4. integer :: NBlocks = 10000
  5. integer :: istat, xsum, i
  6. integer, parameter :: dp = selected_real_kind(15,300)
  7. real(kind=dp) :: Xavg  !Average
  8. real(kind=dp) :: Std   !Standard Deviation  
  9. real(kind=dp) :: Sign_Level = 0.95
  10. real(kind=dp) :: xbar=0, Nxvar=0, xvar=0, xstd=0
  11. real(kind=dp), dimension(10000) :: X
  12. real(kind=dp), dimension(10000) :: XBavg
  13.  
  14. open(unit=100, file="KK.txt", status="old", action="read",iostat=istat)
  15. if (istat .ne. 0) stop "error opening file"
  16.  
  17. open(unit=200, file = 'err.dat')
  18.  
  19. read(100,*) X
  20.  
  21. call blockav( NPoints, NBlocks, X, XBavg, Xavg, Std)
  22.  
  23.  
  24. print *, Std
  25.  
  26. call meanvar(x,xbar,xvar,Nxvar, xstd)
  27.  
  28.  
  29. contains
  30.  
  31.  subroutine meanvar(x,xbar,xvar,Nxvar, xstd)
  32.   integer, parameter :: dp = selected_real_kind(15,300)
  33.   real(kind=dp), intent(out) :: xbar, Nxvar, xvar, xstd
  34.   integer :: i, j
  35.   real(kind=dp) :: xsum
  36.   real(kind=dp), dimension(10000), intent(in) :: x
  37.  
  38.   do i = 1,10000
  39.  
  40.       xsum = xsum + X(i)
  41.    
  42.       xvar = Nxvar/10000
  43.  
  44.   end do
  45.  
  46.       xbar = xsum/10000
  47.  
  48.   do i = 1,10000
  49.  
  50.       Nxvar = Nxvar + (x(i) - xbar)**2
  51.   end do
  52.  
  53.       xvar = Nxvar/10000
  54.       xstd = xvar**1/2
  55. end subroutine
  56.  
  57. subroutine blockav( NPoints, NBlocks, X, XBavg, Xavg, Std)
  58.  
  59. integer,intent(in) :: NPoints   !Number of values to be averaged
  60. integer,intent(in) :: NBlocks       !Number of blocks to be used
  61. integer, parameter :: dp = selected_real_kind(15,300)
  62. real(kind=dp), dimension(NPoints), intent(in) :: X     !Array for values to be block averaged
  63. real(kind=dp), dimension(NBlocks), intent(out) :: XBavg !Block Averages
  64. real(kind=dp), intent(out) :: Xavg  !Average
  65. real(kind=dp), intent(out) :: Std   !Standard Deviation
  66. real(kind=dp):: Sign_Level = 0.95
  67.  
  68.  
  69. !Local variables.
  70. integer::  iblock, ipoint
  71. integer::  FirstPoint, Pointsperblock
  72.  
  73.  
  74.  
  75.  
  76. Pointsperblock = NPoints / NBlocks   !Integer division
  77.  
  78. if (Pointsperblock.eq.0) then
  79.     write (*,*) 'The number of Points per block is too large for the specified number of blocks'
  80.     write (*,*) 'Please specify a smaller number of Points per block, or increase the number of Points'
  81.     stop
  82. endif
  83.  
  84. FirstPoint = MOD( NPoints, NBlocks )
  85.  
  86.    
  87. !Initializing
  88. XBavg = 0.0
  89. Xavg = 0.0
  90. Std = 0.0
  91.  
  92.  
  93. !Calculate the block averages.
  94. do iblock = 1, NBlocks
  95.     do ipoint = FirstPoint + (iblock - 1) * Pointsperblock + 1, &
  96.                 FirstPoint + iblock * Pointsperblock
  97.         XBavg(iblock) = XBavg(iblock) + X(ipoint)
  98.        
  99.     end do
  100.     XBavg(iblock) = XBavg(iblock) / real(Pointsperblock)
  101.     Std = sqrt((Std + ( XBavg(iblock) - Xavg ) ** 2)/(Nblocks -1)) 
  102.     write(200,*) Std
  103.    
  104. end do                       
  105.  
  106. !Average over the simulation.  
  107. Xavg = SUM( XBavg( 1:NBlocks ) ) / real( NBlocks )
  108.  
  109.  
  110.  
  111.  
  112.  
  113. return
  114.  
  115. end subroutine
  116.  
  117. end program
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top