V3FIT
spline2d.f
1  SUBROUTINE splnyr(x,y,n,yp,wk)
2  USE precision
3  IMPLICIT NONE
4  INTEGER :: n
5  REAL(rprec) :: x(n),y(n),yp(n),wk(n,5)
6 c...........................................................................
7  REAL(rprec) :: t,v,btmp
8  REAL(rprec) :: a2,b2,c2,an2,bn2,cn2
9  INTEGER :: j
10 c
11  t=(x(2)-x(1))/(x(3)-x(2))
12  v=(x(n)-x(n-2))/(x(n-1)-x(n-2))
13 c
14  DO j=1,n-1
15  yp(j)=1./(x(j+1)-x(j))
16  ENDdo
17 c
18 c set up the maxtrix elements for y'
19  wk(2,2)=(2.+t+1./t)*yp(2)
20  wk(2,3)=(1.+t)*yp(2)
21  wk(2,4)=(y(2)-y(1))*yp(1)**2+(y(3)-y(2))*(2.+3.*t)*yp(2)**2
22  wk(n-1,1)=(v+1./v-2.)*yp(n-2)
23  wk(n-1,2)=(v-1.)*yp(n-2)
24  wk(n-1,4)=(1./v**2+2.*v-3.)*yp(n-2)**2*(y(n-1)-y(n-2))+
25  $ (y(n)-y(n-1))*yp(n-2)**2/v**2
26  IF (n .gt. 4)then
27  DO j=3,n-2
28  wk(j,1)=yp(j-1)
29  wk(j,2)=2.*(yp(j-1)+yp(j))
30  wk(j,3)=yp(j)
31  wk(j,4)=3.*((y(j)-y(j-1))*yp(j-1)**2+(y(j+1)-y(j))*yp(j)**2)
32  ENDdo
33  ENDIF
34 c
35  btmp=wk(2,2)
36  yp(2)=wk(2,4)/btmp
37  DO j=3,n-1
38  wk(j,5)=wk(j-1,3)/btmp
39  btmp=wk(j,2)-wk(j,1)*wk(j,5)
40  yp(j)=(wk(j,4)-wk(j,1)*yp(j-1))/btmp
41  ENDdo
42  DO j=n-2,2,-1
43  yp(j)=yp(j)-yp(j+1)*wk(j+1,5)
44  ENDdo
45 c
46  a2=yp(2)*(x(3)-x(2))
47  b2=3.*(y(3)-y(2))-(x(3)-x(2))*(2.*yp(2)+yp(3))
48  c2=-2.*(y(3)-y(2))+(x(3)-x(2))*(yp(2)+yp(3))
49  yp(1)=(a2+t*(-2.*b2+t*3.*c2))/(x(3)-x(2))
50 c
51  an2=yp(n-2)*(x(n-1)-x(n-2))
52  bn2=3.*(y(n-1)-y(n-2))-(x(n-1)-x(n-2))*(2.*yp(n-2)+yp(n-1))
53  cn2=-2.*(y(n-1)-y(n-2))+(x(n-1)-x(n-2))*(yp(n-2)+yp(n-1))
54  yp(n)=(an2+v*(2.*bn2+v*3.*cn2))/(x(n-1)-x(n-2))
55 c
56  END SUBROUTINE splnyr
57  SUBROUTINE bcspline(x,y,f,nx,ny,nxp,c)
58  USE precision
59  IMPLICIT NONE
60  INTEGER :: nx,ny,nxp
61  REAL(rprec) :: f(nxp,ny),x(nx),y(ny),c(4,nxp,ny)
62 c..........................................................................
63  INTEGER :: kx
64  parameter(kx=515)
65  REAL(rprec) :: g(kx),gp(kx),wk(5*kx)
66  INTEGER :: i,j
67 c
68  IF (nx .lt. 4 .or. nx .gt. kx)then
69  WRITE (6,'(a)').or.'bcspline: nx < 4 nx > kx = 515'
70  stop
71  ENDIF
72  IF (ny .lt. 4 .or. ny .gt. kx)then
73  WRITE (6,'(a)').or.'bcspline: ny < 4 ny > kx = 515'
74  stop
75  ENDIF
76  IF (nx .gt. nxp)then
77  WRITE (6,'(a)')'bcspline: nx > nxp (storage error)'
78  stop
79  ENDIF
80 c
81  IF (x(nx)-x(1) .lt. 0.)then
82  WRITE (6,'(a)')'bcspline: x not monotonically increasing'
83  stop
84  ENDIF
85  IF (y(ny)-y(1) .lt. 0.)then
86  WRITE (6,'(a)')'bcspline: y not monotonically increasing'
87  stop
88  ENDIF
89 c
90  DO j=1,ny
91  DO i=1,nx
92  c(1,i,j)=f(i,j)
93  ENDdo
94  ENDdo
95 c
96  DO j=1,ny
97  DO i=1,nx
98  g(i)=f(i,j)
99  ENDdo
100  CALL splnyr(x,g,nx,gp,wk)
101  DO i=1,nx
102  c(2,i,j)=gp(i)
103  ENDdo
104  ENDdo
105 c
106  DO i=1,nx
107  DO j=1,ny
108  g(j)=f(i,j)
109  ENDdo
110  CALL splnyr(y,g,ny,gp,wk)
111  DO j=1,ny
112  c(3,i,j)=gp(j)
113  ENDdo
114  ENDdo
115 c
116  DO i=1,nx
117  DO j=1,ny
118  g(j)=c(2,i,j)
119  ENDdo
120  CALL splnyr(y,g,ny,gp,wk)
121  DO j=1,ny
122  c(4,i,j)=gp(j)
123  ENDdo
124  ENDdo
125 c
126  END SUBROUTINE bcspline
127  SUBROUTINE bcsplint(x,y,c,nx,ny,nxp,xl,yl,pds,icalc,ier)
128  USE precision
129  IMPLICIT NONE
130  INTEGER :: nx,ny,nxp,icalc,ier
131  REAL(rprec) :: x(nx),y(ny),c(4,nxp,ny),xl,yl,pds(6)
132 c...........................................................................
133  INTEGER :: i,j
134  REAL(rprec) :: hx, hy, sux(2), suy(2),
135  $ su(2), svx(2), sv(2), sxy(2),
136  $ u, v, spln0, spln1, spln2, s0, sh, sp0, sph, h, d
137 c
138  spln0(s0,sh,sp0,sph,h,d) = s0+d*(h*sp0+d*(3.*(sh-s0)-
139  $ (sph+2.*sp0)*h+d*(2.*(s0-sh)+(sph+sp0)*h)))
140  spln1(s0,sh,sp0,sph,h,d) = sp0+d*(6.*(sh-s0)/h-2.0*
141  $ (sph+2.*sp0)+3.*d*(2.*(s0-sh)/h+(sph+sp0)))
142  spln2(s0,sh,sp0,sph,h,d) = 6.*(sh-s0)/h**2-2.0*
143  $ (sph+2.*sp0)/h+d*(2.*(s0-sh)/h**2+(sph+sp0)/h)*6.0
144 c
145  SAVE i,j
146  data i,j/0,0/
147 c
148  ier = 0
149 c --- correlated table search for xl
150  CALL tableintrp (x, nx, xl, i)
151  IF ( i .eq. 0 )then
152  IF (abs(xl-x(1)) .le. abs(x(2)-x(1)))then
153  i=1
154  ier = 3
155  WRITE (6,'(a,i2)')'bcsplint: x < x(1), ier = ',ier
156  go to 10
157  ELSE
158  ier=33
159  WRITE (6,'(a,i2)')'bscplint: x < x(1), ier = ',ier
160  RETURN
161  ENDIF
162  ENDIF
163  IF ( i .eq. nx)then
164  IF (abs(xl-x(nx)) .le. abs(x(nx-1)-x(nx)))then
165  i=nx-1
166  ier=5
167  WRITE (6,'(a,i2)')'bcsplint: x > x(nx), ier = ',ier
168  ELSE
169  ier=35
170  WRITE (6,'(a,i2)')'bcsplint: x > x(nx), ier = ',ier
171  RETURN
172  ENDIF
173  ENDIF
174  10 CONTINUE
175 c --- correlated table search for yl
176  CALL tableintrp (y, ny, yl, j)
177  IF ( j .eq. 0 )then
178  IF (abs(yl-y(1)) .le. abs(y(2)-y(1)))then
179  j=1
180  ier=4
181  WRITE (6,'(a,i2)')'bcsplint: y < y(1), ier = ',ier
182  go to 20
183  ELSE
184  ier=34
185  WRITE (6,'(a,i2)')'bcsplint: y < y(1), ier = ',ier
186  RETURN
187  ENDIF
188  ENDIF
189  IF ( j .eq. ny)then
190  IF (abs(yl-y(ny)) .le. abs(y(ny-1)-y(ny)))then
191  j=ny-1
192  ier=6
193  WRITE (6,'(a,i2)')'bcsplint: y > y(ny), ier = ',ier
194  ELSE
195  ier=36
196  WRITE (6,'(a,i2)')'bcsplint: y > y(ny), ier = ',ier
197  RETURN
198  ENDIF
199  ENDIF
200  20 CONTINUE
201 c
202  hx = x(i+1)-x(i)
203  hy = y(j+1)-y(j)
204  u = (xl-x(i))/hx
205  v = (yl-y(j))/hy
206 c
207  sv(1)= spln0(c(1,i,j),c(1,i,j+1),c(3,i,j),c(3,i,j+1),hy,v)
208  svx(1)=spln0(c(2,i,j),c(2,i,j+1),c(4,i,j),c(4,i,j+1),hy,v)
209  sv(2)= spln0(c(1,i+1,j),c(1,i+1,j+1),c(3,i+1,j),c(3,i+1,j+1),hy,v)
210  svx(2)=spln0(c(2,i+1,j),c(2,i+1,j+1),c(4,i+1,j),c(4,i+1,j+1),hy,v)
211  pds(1)=spln0(sv(1),sv(2),svx(1),svx(2),hx,u)
212  IF (icalc .le. 1)return
213 c
214  su(1)= spln0(c(1,i,j),c(1,i+1,j),c(2,i,j),c(2,i+1,j),hx,u)
215  suy(1)=spln0(c(3,i,j),c(3,i+1,j),c(4,i,j),c(4,i+1,j),hx,u)
216  su(2)= spln0(c(1,i,j+1),c(1,i+1,j+1),c(2,i,j+1),c(2,i+1,j+1),hx,u)
217  suy(2)=spln0(c(3,i,j+1),c(3,i+1,j+1),c(4,i,j+1),c(4,i+1,j+1),hx,u)
218  pds(2)=spln1(sv(1),sv(2),svx(1),svx(2),hx,u)
219  IF (icalc .eq. 2)return
220  pds(3)=spln1(su(1),su(2),suy(1),suy(2),hy,v)
221  IF (icalc .eq. 3)return
222 c
223  sux(1)=spln1(c(1,i,j),c(1,i+1,j),c(2,i,j),c(2,i+1,j),hx,u)
224  sxy(1)=spln1(c(3,i,j),c(3,i+1,j),c(4,i,j),c(4,i+1,j),hx,u)
225  sux(2)=spln1(c(1,i,j+1),c(1,i+1,j+1),c(2,i,j+1),c(2,i+1,j+1),hx,u)
226  sxy(2)=spln1(c(3,i,j+1),c(3,i+1,j+1),c(4,i,j+1),c(4,i+1,j+1),hx,u)
227  pds(4) = spln1(sux(1),sux(2),sxy(1),sxy(2),hy,v)
228  IF (icalc .eq. 4) RETURN
229  pds(5) = spln2(sv(1),sv(2),svx(1),svx(2),hx,u)
230  IF (icalc .eq. 5) RETURN
231  pds(6) = spln2(su(1),su(2),suy(1),suy(2),hy,v)
232 c
233  END SUBROUTINE bcsplint
234 
235 
236  SUBROUTINE tableintrp (xx, n, x, jlo)
237 c
238 c ----------------------------------------------------------------------
239 c --- correlated table search routine. USE jlo from previous CALL to
240 c --- get jlow for current value of x. IF jlo from previous CALL is
241 c --- no good(jlo=0 or jlo=n) THEN DO a binary search.
242 c --- this routine RETURNs jlo so that
243 c xx(jlo) .le. x .le. xx(jlo+1) IF jlo=1
244 c xx(jlo) .lt. x .le. xx(jlo+1) IF jlo =2,3...n-1
245 c it is assumed that xx(j),j=1,2..n is monotonically
246 c increasing,this is NOT checked for.
247 c this is a MODified version of SUBROUTINE HUNT (Numerical Recipes)
248 c ------------------------------------------------------------------ HSJ
249 c
250  USE precision
251  IMPLICIT NONE
252 c
253  INTEGER :: n,jlo,jhi,jmid,inc
254  REAL(rprec) :: x,xx(n)
255 c
256  IF (x .lt. xx(1)) THEN
257  jlo=0 ! indicates out of range below xx(1)
258  ELSE IF (x .le. xx(2)) THEN
259  jlo=1 ! xx(1) .le. x .le. xx(2)
260  ELSE IF (x .le. xx(n)) THEN ! xx(2) .lt. x .le. xx(n)
261 c
262 c check IF jlo from previous CALL is usable
263 c
264  IF (jlo .eq. 0 .or. jlo .eq. n) THEN
265  jlo=2
266  jhi=n
267  go to 15 ! no correlation go directly to bisection
268  ELSE ! 1 .le. jlo .lt. n
269 c
270 c bracket x, THEN USE bisection
271 c start with jlo from previous CALL
272 c
273  inc=1
274  IF (x .gt. xx(jlo)) THEN ! search up
275  4 jhi=jlo+inc
276  IF (jhi .ge. n) THEN
277  jhi=n
278  ELSE IF (x .gt. xx(jhi)) THEN
279  inc=inc+inc
280  jlo=jhi
281  go to 4
282  END IF
283  ELSE
284  5 jhi=jlo
285  jlo=jlo-inc
286  IF (jlo .le. 1) THEN
287  jlo=1
288  ELSE IF (x .le. xx(jlo)) THEN
289  inc=inc+inc
290  go to 5
291  END IF
292  END IF
293 c
294  END IF
295 c
296 c bisection
297 c
298  10 IF (jhi-jlo .eq. 1) RETURN
299  15 jmid=(jhi+jlo)/2
300  IF (x .gt. xx(jmid)) THEN
301  jlo=jmid
302  ELSE ! x .le. xx(jmid)
303  jhi=jmid
304  END IF
305  go to 10
306 c
307  ELSE ! x .gt. xx(n)
308  jlo=n
309  END IF
310 c
311  END SUBROUTINE tableintrp
312 
313 
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17