• API
• FAQ
• Tools
• Archive
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.
15. if (istat .ne. 0) stop "error opening file"
16.
17. open(unit=200, file = 'err.dat')
18.
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.

Top