-
Notifications
You must be signed in to change notification settings - Fork 0
/
statistics.f90
191 lines (171 loc) · 4.91 KB
/
statistics.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
subroutine getsstatistics(coord,nsteps,skip,nevalfluc,dt,Z,H0,minsegmentlenght,goodrav,gooddevav)
implicit none
real (kind=4), DIMENSION(nsteps) :: coord
integer :: skip, nsteps, nskip, nmax, nevalfluc, nsegment, segmentlenght, minsegmentlenght
double precision, allocatable, dimension(:) :: segmentedcoord, segmenteddev, coordav
double precision, dimension(38) :: test
double precision :: dt, S, m, Z, xi, xj, U, V, gooddevav, goodrav
integer :: i,j,k,start,end,count,N
logical :: toosmall, nsegmentok
logical, intent(inout) :: H0
! nmax=0
H0=.False.
toosmall=.False.
nskip=skip
do while((.not.H0) .and. (.not.toosmall))
if (allocated(coordav)) deallocate(coordav)
allocate(coordav(nsteps-nskip))
coordav=0.d0
do k=nskip+1,nsteps
coordav(k-nskip)=(coordav(k-nskip-1)*(dble(k-1-nskip))+coord(k))/dble(k-nskip)
end do
if ((nsteps-nskip).le.nevalfluc) toosmall=.True.
if ((nsteps-nskip).le.nevalfluc) write(*,*) "Not enough data 1"
nmax=0
do i=2, nevalfluc-1
if ((coordav(i) .gt. coordav(i-1)) .and. (coordav(i) .gt. coordav(i+1))) nmax=nmax+1
end do
if (nmax.eq.0) then
segmentlenght=minsegmentlenght
else
segmentlenght=nevalfluc/nmax
end if
write(9999,*) "Fluctuation time is:", dble(segmentlenght*dt)
if (segmentlenght.lt.minsegmentlenght) segmentlenght=minsegmentlenght
nsegment=(nsteps-nskip)/segmentlenght
if (nsegment.gt.100) then
nsegment=100
segmentlenght=(nsteps-nskip)/nsegment
endif
write(9999,*) "Number of segments:", nsegment
write(9999,*) "Number of points in each segment:", segmentlenght
if (nsegment.lt.24) toosmall=.True.
if (nsegment.lt.24) write(*,*) "Not enough data 2"
if (allocated(segmentedcoord)) deallocate(segmentedcoord)
if (allocated(segmenteddev)) deallocate(segmenteddev)
allocate(segmentedcoord(nsegment),segmenteddev(nsegment))
!computes mean
do i=1,nsegment
segmentedcoord(i)=0.d0
start=(i-1)*segmentlenght+1
end=start+segmentlenght-1
do j=start,end
count=1
segmentedcoord(i)=(segmentedcoord(i)*(dble(count-1))+coordav(j))/dble(count)
count=count+1
end do
end do
!computes standard deviation
do i=1,nsegment
segmenteddev(i)=0.d0
start=(i-1)*segmentlenght+1
end=start+segmentlenght-1
do j=start,end
segmenteddev(i)=segmenteddev(i)+(coordav(j)-segmentedcoord(i))**2
end do
segmenteddev(i)=dsqrt(segmenteddev(i)/(segmentlenght-1))
end do
!do Mann-Kendall test for trend
! nsegmentok=.False.
! N=24
! do while ((.not. nsegmentok) .and. (N.le.nsegment))
N=nsegment
S=0.d0
do i=1,N-1
do j=i+1,N
xi=segmentedcoord(i)
xj=segmentedcoord(j)
if ((xj-xi).gt.0.d0) S = S+1.d0
if ((xj-xi).lt.0.d0) S = S-1.d0
if ((xj-xi).eq.0.d0) S = S+0.d0
end do
end do
if (S.gt.0.d0) m=-1.d0
if (S.lt.0.d0) m=1.d0
if (S.eq.0.d0) m=0.d0
V=dsqrt((N*(N-1))/18.d0)
V=V*dsqrt(dble(2*N+5))
Z=(S-m)/V
U=1.96d0
H0=.False.
if (abs(Z).lt.U) then
H0=.True.
! nsegmentok=.True.
! write(*,*) N,nsegmentok
! write(*,*) "Gotta use:",N*segmentlenght, " steps (",N, "segments) Drop ", nskip, " steps"
! else
! nskip=nskip+1000
endif
! N=N+1
! end do
if (H0) then
write(*,*) "Steps used for analysis:",N*segmentlenght
write(*,*) "Segments:", N
write(*,*) "Droped:", nskip
write(*,*) "Total:", nskip+N*segmentlenght
else
nskip=nskip+1000
endif
end do
N=N-1
! do i=1,N
! write(8888,*) i, segmentedcoord(i)
! end do
! write(8888,*)
if(toosmall) then
write(*,*) "Computing straight average"
if (allocated(coordav)) deallocate(coordav)
allocate(coordav(nsteps-skip))
do k=skip+1,nsteps
coordav(k-skip)=(coordav(k-skip-1)*(dble(k-1-skip))+coord(k))/dble(k-skip)
end do
goodrav=coordav(nsteps-skip)
else
goodrav=0.d0
do i=1,N
goodrav=goodrav+segmentedcoord(i)
end do
goodrav=goodrav/dble(N)
endif
gooddevav=0.d0
do j=1,N
gooddevav=gooddevav+(segmentedcoord(j)-goodrav)**2
end do
gooddevav=dsqrt(gooddevav/(N-1))
end subroutine getsstatistics
subroutine getmaxstd(nrestr,nrep,rep,fav,devav,maxstd,maxstdat,forcedev)
implicit none
double precision, dimension(3,nrestr,nrep), intent(in) :: fav, devav
integer, intent(in) :: nrestr, nrep, rep
double precision, intent(out) :: maxstd
double precision :: fmod,std
integer :: i,j,maxstdat
logical :: forcedev
maxstd=0.d0
maxstdat=1
do i=1,nrestr
if(forcedev) then! propagates std on |F|; in that case devav should be kref*devav
std=(fav(1,i,rep)*devav(1,i,rep))**2+&
(fav(2,i,rep)*devav(2,i,rep))**2+&
(fav(3,i,rep)*devav(3,i,rep))**2
fmod=fav(1,i,rep)**2+&
fav(2,i,rep)**2+&
fav(3,i,rep)**2
fmod=dsqrt(fmod)
std=dsqrt(std)
std=std/fmod
if (std.gt.maxstd) then
maxstd=std
maxstdat=i
end if
else
do j=1,3
std=devav(j,i,rep)
if (std.gt.maxstd) then
maxstd=std
maxstdat=i
end if
end do
end if
end do
end subroutine getmaxstd