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
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
|
! Copyright (C) 2005 Barbara Ercolano
!
! Version 2.02
module interpolation_mod
contains
subroutine sortUp(arr)
implicit none
real, dimension(:), intent(inout) :: arr
real, allocatable :: tmp(:)
real :: min
integer :: i, j
integer :: n ! size
n = size(arr)
allocate(tmp(n))
min = 1.e30
do i = 1, n
if (arr(i) < min) min = arr(i)
end do
tmp(1) = min
do j = 2, n
min = 1.e30
do i = 1, n
if (arr(i) < min .and. arr(i) > tmp(j-1)) min = arr(i)
end do
tmp(j) = min
end do
arr = tmp
deallocate(tmp)
end subroutine sortUp
! given an array xa of length n, and given a value x, this
! routine returns a value ns, such that x is located between
! xa(ns) and xa(ns+1)
! ns = 0 or ns=n is returned to indicate that x is out
! of range.
subroutine locate(xa,x,ns)
implicit none
integer, intent(out) :: ns
real, dimension(:), intent(in) :: xa ! input array
real, intent(in) :: x ! variable to be located
! local variables
integer :: n ! size of array xa
n = size(xa)
! first check if x is out of range
if ( x > xa(n) ) then
ns = n
return
end if
if ( x < xa(1) ) then
ns = 0
return
end if
! if not, then locate
! the command finds the location of the smallest positive value of xa-x
! x lies between this location and the previous one
! so subtract one, and then x lies between xa(ns) and xa(ns+1)
ns=max(minloc((xa-x),1,(xa-x).gt.0)-1,1)
end subroutine locate
! this routine will map y(x) onto x_new and return y_new(x_new)
! mapping is carried out by means of linear interpolation
subroutine linearMap(y, x, nx, y_new, x_new, nx_new)
implicit none
real, intent(out) :: y_new(*)
real, intent(in) :: y(*), x(*), x_new(*)
integer, intent(in) :: nx, nx_new
integer :: i, ii
do i = 1, nx_new
call locate(x(1:nx), x_new(i), ii)
if (ii==0) then
y_new(i) = y(1)
else if (ii==nx) then
y_new(i) = y(nx)
else
y_new(i) = y(ii) +(y(ii+1)-y(ii))*(x_new(i)-x(ii))/(x(ii+1)-x(ii))
end if
end do
end subroutine linearMap
! subroutine polint carries out polynomial interpoation
! or extrapolation. given arrays xa nd ya, each of length
! n, and given a value x, it returns a value y and an
! error estimate dy.
subroutine polint(xa, ya, x, y, dy)
implicit none
real, dimension(:), intent(in) :: xa, ya
real, intent(in) :: x
real, intent(out) :: y, dy
! local variables
integer :: i, m ! counters
integer :: n ! size of the arrays
integer :: ns = 1
real :: den, ho, hp, w
real, dimension(size(xa)) :: c, d
! locate the index ns closest to the table entry
call locate(xa, x, ns)
c = ya
d = ya
y = ya(ns)
ns = ns-1
n = size(xa)
! for each entry in the table, loop over the current
! c's and d's and update them
do m=1, n-1
do i = 1, n-m
ho = xa(i) - x
hp = xa(i+m) - x
w = c(i+1) - d(i)
den = ho - hp
if (den == 0) then
print*, "! polint: two input xa's &
& identical (to within roundoff)"
stop
end if
den = w/den
! update c and d
d(i) = hp*den
c(i) = ho*den
end do
! after each column in the table, decide the
! correction c or d to be added to the
! accumulating value of y.
if ( (2*ns) < (n-m) ) then
dy = c(ns+1)
else
dy = d(ns)
ns = ns-1
end if
y = y+dy
end do
end subroutine polint
! subroutine spline given arrays wa and ya each of length
! n and given yp1 and ypn for the first derivatives of
! the interpolating function at points 1 and n
! respectively, returns the array y2a of length n
! containing the second derivatives of the interpolating
! funtion at the tabulated points. to use the condition
! of natural spline (2nd derivatives = 0 at the end
! points) set yp1 and ypn to a value > or = to 1.e30.
subroutine spline(xa, ya, yp1,ypn, y2a)
implicit none
real, dimension(:), intent(inout) :: xa, ya
real, intent(in) :: yp1, ypn
real, dimension(:), intent(out) :: y2a
! local variables
integer :: i, k ! counters
integer :: n ! size of the arrays
integer, parameter :: nmax=500! safety limit
real :: p, qn, sig, un
real, dimension(nmax) :: u
n = size(xa)
if (yp1 > .99e30) then
y2a(1) = 0.
u(1) = 0.
else
y2a(1) = -.5
if ((xa(2)-xa(1))==0.) then
print*, "! spline: bad xa input or x out of range [2,1] "
stop
end if
u(1) = (3./(xa(2)-xa(1)))* &
& ((ya(2)-ya(1))/(xa(2)-xa(1))-yp1)
end if
do i = 2, n-1
if ((xa(i+1)-xa(i-1))==0.) then
print*, "! spline: bad xa input or x out of range [i+1, i-1] ", (i+1), (i-1)
stop
end if
sig = (xa(i)-xa(i-1))/(xa(i+1)-xa(i-1))
p=sig*y2a(i-1)+2.
y2a(i)=(sig-1.)/p
if ((xa(i+1)-xa(i))==0.) then
print*, "! spline: bad xa input or x out of range [i+1, i] ", (i+1), (i)
stop
end if
if ((xa(i)-xa(i-1))==0.) then
print*, "! spline: bad xa input or x out of range [i, i-1] ", (i), (i-1)
stop
end if
u(i)= ya(i+1)
u(i)=u(i) - ya(i)
u(i) = u(i) / (xa(i+1)-xa(i))
u(i)= u(i) - (ya(i)-ya(i-1))/(xa(i)-xa(i-1))
u(i)=(6.*u(i)/(xa(i+1)-xa(i-1)) - sig*u(i-1))/p
end do
if (ypn > .99e30) then
qn=0.
un=0.
else
qn=0.5
if ((xa(n)-xa(n-1))==0.) then
print*, "! spline: bad xa input or x out of range [n, n-1] ", n, (n-1)
stop
end if
un=(3./(xa(n)-xa(n-1))) * &
& (ypn-(ya(n)-ya(n-1))/(xa(n)-xa(n-1)))
end if
y2a(n)=(un-qn*u(n-1))/(qn*y2a(n-1)+1.)
do k = (n-1), 1, -1
y2a(k) = y2a(k)*y2a(k+1)+u(k)
end do
end subroutine spline
! subroutine splint, given the arrays xa and ya each of
! length n and given the array y2a, which is the output
! from the subroutine spline, it returns a cubic spline
! interpolated value y for the given x
! NOTE: if x is out of range splint will return the limit
! values
subroutine splint(xa, ya, y2a, x, y)
real, dimension(:), intent(in) :: xa, ya, y2a
real, intent(in) :: x
real, intent(out) :: y
! local variables
integer :: klo, khi, i, n
integer, parameter :: imax = 1000
real :: a, b, h
n = size(xa)
! check if x is within range
if ( x <= xa(1) ) then
y = ya(1)
return
else if ( x >= xa(n) ) then
y = ya(n)
return
end if
klo = 1
khi = n
do i = 1, imax
if ( (khi-klo) <= 1) exit
k = (khi+klo)/2
if (xa(k) > x) then
khi=k
else
klo=k
end if
end do
h=xa(khi)-xa(klo)
if (h==0.) then
print*, "! splint: bad xa input or x out of range [khi, klow] ", khi, klo
stop
end if
a = (xa(khi) - x)/h
b = (x - xa(klo))/h
y = a*ya(klo)+b*ya(khi) + &
& ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi)) * &
& (h**2)/6.
end subroutine splint
! this function performs a linear interpolation of a scalar quantity in a
! frequency dependent 3D grid. the interpolation coefficients t1, t2, t3 are
! are left to be calculated in the calling program, so that only a small
! number or arguments need to be passed to the procedure
function interpGrid(s3Dgrid, s4DGrid, s5Dgrid, xP, yP, zP, freqP, elem, ion, t1, t2, t3)
implicit none
integer, intent(in) :: xP, yP, zP ! x, y and z axes indeces
integer, intent(in), optional :: freqP ! frequency index
integer, intent(in), optional :: elem ! element index
integer, intent(in), optional :: ion ! ion index
real, intent(in), dimension(:,:,:),&
& optional :: s3Dgrid ! 3D scalar grid
real, intent(in), dimension(:,:,:,:),&
& optional :: s4Dgrid ! 4D scalar grid
real, intent(in), dimension(:,:,:,:,:),&
& optional :: s5Dgrid ! 5D scalar grid
real, intent(in) :: t1, t2, t3 ! interpolation coefficients
real :: interpGrid ! interpolated value
if ( present(freqP) .and. (.not.present(elem)) .and. (.not.present(ion)) ) then
if(.not.present(s4Dgrid)) then
print*, "! interpGrid: insanity occurred - arguments incompatible"
stop
end if
interpGrid = &
& ((1.-t1) * (1.-t2) * (1.-t3))* s4Dgrid(xP , yP , zP , freqP) + &
& ((t1 ) * (1.-t2) * (1.-t3))* s4Dgrid(xP+1, yP , zP , freqP) + &
& ((t1 ) * (t2 ) * (1.-t3))* s4Dgrid(xP+1, yP+1 , zP , freqP) + &
& ((1.-t1) * (t2 ) * (t3 ))* s4Dgrid(xP , yP+1 , zP+1 , freqP) + &
& ((1.-t1) * (t2 ) * (1.-t3))* s4Dgrid(xP , yP+1 , zP , freqP) + &
& ((t1 ) * (1.-t2) * (t3 ))* s4Dgrid(xP+1, yP , zP+1 , freqP) + &
& ((1.-t1) * (1.-t2) * (t3 ))* s4Dgrid(xP , yP , zP+1 , freqP) + &
& ((t1 ) * (t2 ) * (t3 ))* s4Dgrid(xP+1, yP+1 , zP+1 , freqP)
else if ((.not.present(freqP)) .and. present(elem) .and. present(ion) ) then
if(.not.present(s5Dgrid)) then
print*, "! interpGrid: insanity occurred - arguments incompatible"
stop
end if
interpGrid = &
& ((1.-t1) * (1.-t2) * (1.-t3))* s5Dgrid(xP , yP , zP , elem, ion) + &
& ((t1 ) * (1.-t2) * (1.-t3))* s5Dgrid(xP+1, yP , zP , elem, ion) + &
& ((t1 ) * (t2 ) * (1.-t3))* s5Dgrid(xP+1, yP+1 , zP , elem, ion) + &
& ((1.-t1) * (t2 ) * (t3 ))* s5Dgrid(xP , yP+1 , zP+1 , elem, ion) + &
& ((1.-t1) * (t2 ) * (1.-t3))* s5Dgrid(xP , yP+1 , zP , elem, ion) + &
& ((t1 ) * (1.-t2) * (t3 ))* s5Dgrid(xP+1, yP , zP+1 , elem, ion) + &
& ((1.-t1) * (1.-t2) * (t3 ))* s5Dgrid(xP , yP , zP+1 , elem, ion) + &
& ((t1 ) * (t2 ) * (t3 ))* s5Dgrid(xP+1, yP+1 , zP+1 , elem, ion)
else if ((.not.present(freqP)) .and. (.not.present(elem)) .and. (.not.present(ion)) ) then
if(.not.present(s3Dgrid)) then
print*, "! interpGrid: insanity occurred - arguments incompatible"
stop
end if
interpGrid = &
& ((1.-t1) * (1.-t2) * (1.-t3))* s3Dgrid(xP , yP , zP ) + &
& ((t1 ) * (1.-t2) * (1.-t3))* s3Dgrid(xP+1, yP , zP ) + &
& ((t1 ) * (t2 ) * (1.-t3))* s3Dgrid(xP+1, yP+1 , zP ) + &
& ((1.-t1) * (t2 ) * (t3 ))* s3Dgrid(xP , yP+1 , zP+1) + &
& ((1.-t1) * (t2 ) * (1.-t3))* s3Dgrid(xP , yP+1 , zP ) + &
& ((t1 ) * (1.-t2) * (t3 ))* s3Dgrid(xP+1, yP , zP+1) + &
& ((1.-t1) * (1.-t2) * (t3 ))* s3Dgrid(xP , yP , zP+1) + &
& ((t1 ) * (t2 ) * (t3 ))* s3Dgrid(xP+1, yP+1 , zP+1)
end if
end function interpGrid
end module interpolation_mod
|
