V3FIT
doakima.f
1  SUBROUTINE doakima(s, f, nin, ier, case)
2  USE precision
3  USE knots
4  IMPLICIT NONE
5  REAL(rprec), DIMENSION(nin), INTENT(in) :: s, f
6  INTEGER, INTENT(in) :: nin
7  CHARACTER(len=*), INTENT(in) :: case
8  INTEGER :: id, m
9  INTEGER, PARAMETER :: nset1 = 11, nset2 = 20
10  REAL(rprec) :: ds
11  REAL(rprec), DIMENSION(nin) :: s3
12  REAL(rprec), DIMENSION(:), ALLOCATABLE :: d1f, d2f, br
13  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: co
14  REAL, DIMENSION(:), ALLOCATABLE :: xplt,yplt
15  REAL siz
16  INTEGER, DIMENSION(:), ALLOCATABLE :: slctn, olli, span
17  INTEGER :: n, i, lastspan, tot, j, count, ier
18  LOGICAL :: match
19  REAL(rprec) :: fvn_spline_eval, fvn_spline_devaldx
20  REAL(rprec) :: fvn_spline_d2evaldx2
21  REAL(rprec) :: c
22  CHARACTER*99 xlabel, ylabel, slabel*60, mlabel*50
23  REAL(rprec) :: rand
24  print*," Akima Spline for ", case
25  CALL pgslw(2)
26  CALL pgscf(2)
27  n = nin
28  ALLOCATE( br(n),stat = ier)
29  ALLOCATE( co(4,n),stat = ier)
30  ALLOCATE( slctn(n),stat = ier)
31  ALLOCATE(f2(n), s2(n),stat = ier)
32  slctn = 0
33  CALL fvn_akima(n,s,f,br,co)
34  ds = s(2)-s(1)
35  DO i = 1,n
36  s2(i) = s(i) !+rand(i)*ds
37  f2(i) = fvn_spline_eval(s(i),n-1,br,co)
38  ENDDO
39  xlabel = 's'
40  ylabel = case
41  mlabel = "Test that spline of ALL points would be successful."
42  slabel = ' '
43  CALL graf2(real(s),real(f),real(f2),n,'s',
44  & trim(ylabel),trim(mlabel)," ")
45  IF (ALLOCATED(xplt)) DEALLOCATE(xplt,yplt)
46  ALLOCATE(xplt(2),yplt(2))
47  CALL pgsave
48  CALL pgsci(2)
49  CALL pgslw(3)
50  yplt(1) = maxval(f2)/6
51  yplt(2) = yplt(1)
52  CALL pgline(2, (/.02,.08/), yplt)
53  CALL pgtext (0.1, yplt(1) ,"Spline Fit")
54  CALL pgsci(3)
55  CALL pgslw(3)
56  CALL pgsls(5)
57  xplt(2) = 0.05
58  yplt(1) = maxval(f2)/5
59  CALL pgpt1 (xplt(2),yplt(1),2)
60  CALL pgtext (0.1, yplt(1), "knot location/values")
61  call pgsch(3)
62  CALL pgpt1 (xplt(2),yplt(1),2)
63  CALL pgunsa
64  IF (ALLOCATED(xplt)) DEALLOCATE(xplt,yplt)
65 ! Now start thinning INTo a new s2.
66  slctn(1) = 1; slctn(2) = SIZE(f)
67  lastspan = 2; tot = 2
68  ALLOCATE(span(nin/nset1))
69  span = 0
70  span = (/(i, i = 1,nin,nset1)/)
71  lastspan = lastspan+nin/nset1-1
72  DO i = 1,nin
73  IF(slctn(i) == 0) EXIT
74  ENDDO
75  slctn(i:i+SIZE(span)-1) = span
76 ! df SIGN change
77  DO i = 1,nin
78  IF(slctn(i) > nin)slctn(i) = 0
79  IF(slctn(i) == 0) cycle
80  tot = i
81  ENDDO
82  lastspan = tot
83  ALLOCATE(olli(n))
84  olli = 0
85  DO i = 2,nin-1
86  IF((f(i)-f(i-1)) == 0.)cycle
87  c = (f(i+1)-f(i))/(f(i)-f(i-1))
88  IF(c < 0 ) olli(i) = -1
89  ENDDO
90  IF (sum(olli) /= 0) THEN
91  DO i = 1,nin
92  IF(olli(i) == 0) cycle
93  slctn(lastspan+1) = i
94  lastspan = lastspan+1
95  ENDDO
96  ENDIF
97  DO i = 1,nin
98  IF(slctn(i) == 0) cycle
99  tot = i
100  ENDDO
101 
102 ! first derivative magnitudes
103  s2 = s
104  f2 = s
105  DO i = 1,nin
106  f2(i) = fvn_spline_devaldx(s(i),n-1,br,co)
107  ENDDO
108  s2 = log10(abs(f2))
109 ! SIGN changes again
110  olli = 0
111  DO i = 2,nin-1
112  IF((f(i)-f(i-1)) == 0.)cycle
113  c = (f2(i+1)-f2(i))/(f2(i)-f2(i-1))
114  IF(c < 0 ) olli(i) = -1
115  ENDDO
116  IF (sum(olli) /= 0) THEN
117  DO i = 1,nin
118  IF(olli(i) == 0) cycle
119  slctn(lastspan+1) = i
120  lastspan = lastspan+1
121  ENDDO
122  ENDIF
123  DO i = 1,nin
124  IF(slctn(i) == 0) cycle
125  tot = i
126  ENDDO
127  DEALLOCATE(span)
128  ALLOCATE(span(SIZE(f2)))
129  span = (/(i,i = 1,nin)/)
130  s3 = s2
131  CALL ssort( s2, span, nin)
132 ! WRITE(*,2000)
133  2000 FORMAT(' Original ',5x,' Sorted ',/,' Array ',
134  $ 5x,' Array')
135 ! WRITE(*,2001) (i, s3(i), span(i), s2(i), i = 1,nin)
136  2001 FORMAT(i3,2x,f5.3,5x, i3,2x, f5.3)
137  IF(ALLOCATED(olli))DEALLOCATE(olli)
138 ! Now start thinning INTo a new s2.
139  ALLOCATE(olli(nin/nset1))
140  olli = (/(i, i = 1,nin,nset1)/)
141  lastspan = lastspan+nin/nset1
142  DO i = 1,nin
143  IF(slctn(i) == 0) EXIT
144  ENDDO
145  slctn(i:i+SIZE(olli)-1) = olli
146  DO i = 1,nin
147  IF(slctn(i) == 0) cycle
148  tot = i
149  ENDDO
150 ! Back to largest d/ds values
151  slctn(lastspan+1:lastspan+nset1) = span(1:nset1)
152  lastspan = lastspan+nset1
153  DO i = 1,nin
154  IF(slctn(i) == 0) cycle
155  tot = i
156  ENDDO
157  IF(ALLOCATED(olli))DEALLOCATE(olli)
158  ALLOCATE(olli(lastspan))
159  olli = (/(i,i = 1,lastspan)/)
160  CALL ssort(dble(slctn(1:lastspan)),olli,lastspan)
161  count = 0
162  DO i = 1,lastspan
163  DO j = 1, lastspan
164  match = (slctn(olli(i)) == slctn(olli(j)) .AND. i /= j)
165  IF(.not. match) cycle
166  slctn(olli(j)) = 0
167  count = count+1
168  ENDDO
169  ENDDO
170  tot = 0
171  DO i = 1,nin
172  IF(slctn(i) == 0) cycle
173  tot = tot+1
174  ENDDO
175  IF(ALLOCATED(span))DEALLOCATE(span)
176  ALLOCATE(span(lastspan))
177  span = slctn(olli(lastspan:1:-1))
178 ! final sort of INDEX
179  IF(ALLOCATED(olli)) DEALLOCATE(olli)
180  ALLOCATE(olli(lastspan))
181  olli = (/(i, i = 1,lastspan)/)
182  CALL ssort(dble(span),olli,lastspan)
183  DO i = 1,lastspan
184  IF(span(olli(i)) == 0) EXIT
185  ENDDO
186  lastspan = i-1
187  IF(ALLOCATED(slctn))DEALLOCATE(slctn)
188  ALLOCATE(slctn(lastspan))
189  slctn = span(olli(lastspan:1:-1))
190  ALLOCATE(xplt(lastspan),yplt(lastspan))
191  xplt = s(slctn); yplt = f(slctn)
192  WRITE(xlabel,fmt = '(i3.3)')lastspan
193  CALL pgqch(siz)
194  CALL pgsch(3*siz/2)
195  IF (lastspan > 100) slabel = "Thinning FAILED!!"
196  CALL graf1pt(xplt,yplt,lastspan,'s',trim(ylabel),
197  5 "Akima spline -- knots thinned to "//trim(xlabel),trim(slabel))
198  slabel = ' '
199  CALL pgsch(siz)
200 ! now spline thinned data set
201  IF(ALLOCATED(co))DEALLOCATE(co,br,s2,f2)
202  n = lastspan
203  ALLOCATE(s2(n),f2(n),br(n),co(4,n))
204  s2 = s(slctn)
205  f2 = f(slctn)
206  CALL fvn_akima(n,s2,f2,br,co)
207  DEALLOCATE(s2,f2)
208  ALLOCATE(s2(nin),f2(nin))
209  s2 = s+ds/2
210  s2(nin) = s(nin)
211  s3=0
212  call random_number(s3)
213  DO i = 1,nin
214  s2(i) = s(i)+s3(i)*ds
215  f2(i) = fvn_spline_eval(s(i),n-1,br,co)
216  ENDDO
217  CALL pgsave
218  CALL pgsci(4)
219  CALL pgslw(7)
220  IF( ALLOCATED(xplt) )DEALLOCATE(xplt,yplt)
221  ALLOCATE(xplt(nin),yplt(nin))
222  xplt=s2; yplt=f2
223  CALL pgline(nin,xplt,yplt)
224  CALL pgsci(3)
225  CALL pgslw(4)
226  CALL pgsls(5)
227  CALL pgline (nin,real(s),real(f))
228  CALL pgunsa
229  CALL pgsci(2)
230  CALL pgpt1(0.1,real(maxval(f2))/5,2)
231  CALL pgtext (0.18, real(maxval(f2))/5, "knot location/values")
232  CALL pgsci(4)
233  IF( ALLOCATED(xplt) )DEALLOCATE(xplt,yplt)
234  ALLOCATE(xplt(2),yplt(2))
235  yplt(1) = real(maxval(f2)/6)
236  yplt(2) = yplt(1)
237  xplt(1) = 0.02
238  xplt(2) = 0.12
239  CALL pgline(2,xplt,yplt)
240  CALL pgtext (0.18, real(maxval(f2))/6, "Akima spline result")
241  DEALLOCATE(xplt,yplt)
242  DEALLOCATE(s2,f2)
243  ALLOCATE(s2(SIZE(slctn)),f2(SIZE(slctn)))
244  s2 = s(slctn); f2 = f(slctn); n2 = SIZE(slctn)
245  IF (lastspan .GT. 100) ier = -lastspan
246  END SUBROUTINE doakima
247 
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17