V3FIT
splines.f
1  SUBROUTINE spline2d(fun,x,y,nx,ny,kx,coef)
2 c-----------------------------------------------------------------------
3 c setup routine for bicubic spline of fun[x,y]
4 c output of this routine is nx*ny spline coefficients
5 c stored in coef(kx,ny)
6 c-----------------------------------------------------------------------
7  USE precision
8  IMPLICIT NONE
9  REAL(rprec) :: fun(*),x(*),y(*),coef(*)
10  INTEGER :: kx,nx,ny
11  INTEGER :: ind,j
12  DO j=1,ny
13  ind=(j-1)*kx+1
14  CALL spline(x,fun(ind),nx,-1.e30_dbl,-1.e30_dbl,coef(ind))
15  END DO
16  RETURN
17  END
18 c=======================================================================
19  SUBROUTINE spline2dt(fun_new,x_new,y_new,kx_new,nx_new,ny_new,
20  $ fun_old,x_old,y_old,kx_old,nx_old,ny_old,coef)
21 c-----------------------------------------------------------------------
22 c evaluates bicubic spline: determines fun_new[x_new,y_new]
23 c after spline2d has been CALLed to determine coef
24 c in the CALLing PROGRAM fun_new,run_old, and coef are 2d arrays:
25 c fun_old(kx_old,ky_old),fun_new(kx_new,ky_new),coef(kx_old,ky_old)
26 c-----------------------------------------------------------------------
27  USE precision
28  IMPLICIT NONE
29  REAL(rprec) :: fun_new(*),x_new(*),y_new(*)
30  REAL(rprec) :: fun_old(*),x_old(*),y_old(*)
31  REAL(rprec) :: coef(*)
32  INTEGER :: kx_new,nx_new,ny_new
33  INTEGER :: kx_old,nx_old,ny_old
34  INTEGER :: i,j,ind,nlow
35  INTEGER, PARAMETER :: ky_old=2000
36  REAL(rprec) :: ftemp(ky_old),ctemp(ky_old)
37  IF(ny_old.gt.ky_old) THEN
38  WRITE(6,'("dimensioning problem in spline2dt")')
39  stop
40  END IF
41  DO i=1,nx_new
42  DO j=1,ny_old
43  ind=(j-1)*kx_old+1
44  nlow=0
45  CALL splint(x_old,fun_old(ind),coef(ind),nx_old,x_new(i),
46  $ ftemp(j),nlow)
47  END DO
48  CALL spline(y_old,ftemp,ny_old,-1.e30_dbl,-1.e30_dbl,ctemp)
49  DO j=1,ny_new
50  nlow=0
51  CALL splint(y_old,ftemp,ctemp,ny_old,y_new(j),
52  $ fun_new((j-1)*kx_new+i),nlow)
53  END DO
54  END DO
55  RETURN
56  END
57 c=======================================================================
58  SUBROUTINE spline1d(ynew,xnew,nnew,yold,xold,nold,y2old)
59 c-----------------------------------------------------------------------
60 c USE 1d cubic spline on yold[xold] to produce ynew[xnew]
61 c y2old(1:nold) is a work array
62 c ynew(1:nnew) is the output
63 c-----------------------------------------------------------------------
64  USE precision
65  IMPLICIT NONE
66  REAL(rprec) :: ynew(*),yold(*),xnew(*),xold(*)
67  REAL(rprec) :: y2old(*)
68  REAL(rprec) :: yp1,ypn
69  INTEGER :: nnew,nold,i,nlow
70  yp1=-1.e30_dbl
71  ypn=-1.e30_dbl
72  CALL spline(xold,yold,nold,yp1,ypn,y2old)
73  DO i=1,nnew
74  nlow=0
75  CALL splint(xold,yold,y2old,nold,xnew(i),ynew(i),nlow)
76  END DO
77  RETURN
78  END
79 c=======================================================================
80  SUBROUTINE spline1dp(ynew,xnew,nnew,yold,xold,nold,yp1,ypn,y2old)
81 c-----------------------------------------------------------------------
82 c USE 1d cubic spline on yold[xold] to produce ynew[xnew]
83 c y2old(1:nold) is a work array
84 c ynew(1:nnew) is the output
85 c This routine ALLows different boundary conditions
86 c-----------------------------------------------------------------------
87  USE precision
88  IMPLICIT NONE
89  REAL(rprec) :: ynew(*),yold(*),xnew(*),xold(*)
90  REAL(rprec) :: y2old(*)
91  REAL(rprec) :: yp1,ypn
92  INTEGER :: nnew,nold,i,nlow
93 c yp1=-1.e30_dbl
94 c ypn=-1.e30_dbl
95  CALL spline(xold,yold,nold,yp1,ypn,y2old)
96  DO i=1,nnew
97  nlow=0
98  CALL splint(xold,yold,y2old,nold,xnew(i),ynew(i),nlow)
99  END DO
100  RETURN
101  END
102 c=======================================================================
103  SUBROUTINE spline(x,y,n,yp1,ypn,y2)
104 c-----------------------------------------------------------------------
105 c spline routine based upon numerical recipes
106 c this is the setup routine which needs to be CALLed ONLY once
107 c splines y as a FUNCTION of x--both arrays have n elements
108 c yp1 and ypn are boundary conditions on the spline
109 c yp1=y'[x] at x=x[1]
110 c IF yp1>=1.e30 THEN y''[x]=0 at x=x[1] is USEd
111 c IF yp1<=-1.e30_dbl THEN y'[x[1]] is calculated from first four points
112 c ypn=y'[x] at x=x[n]
113 c IF ypn>=1.e30 THEN y''[x]=0 at x=x[n] is USEd
114 c IF ypn<=-1.e30_dbl THEN y'[x[n]] is calculated from last four points
115 c y2[1:n] is calculated array of the second derivatives of the
116 c INTerpolating FUNCTION at the x[i]
117 c-----------------------------------------------------------------------
118  USE precision
119  IMPLICIT NONE
120  INTEGER :: n,k,i
121  INTEGER, PARAMETER :: nmax=2000
122  REAL(rprec) :: x(*),y(*),y2(*)
123  REAL(rprec) :: u(nmax)
124  REAL(rprec) :: yp1,ypn
125  REAL(rprec) :: yp1t,ypnt,p,sig,qn,un
126  IF (n .gt. nmax)then
127  WRITE(6,'("spline; DIMENSIONal error")')
128  stop
129  ENDIF
130  IF (yp1.gt..99e30) THEN
131  y2(1)=0.
132  u(1)=0.
133  ELSE IF(yp1.lt.-.99e30) THEN
134  yp1t=(3*x(1)**2+x(2)*x(3)+
135  $ x(2)*x(4)+x(3)*x(4)-2*x(1)*(x(2)+x(3)+x(4)))*
136  $ y(1)/((x(1)-x(2))*(x(1)-x(3))*(x(1)-x(4)))+
137  $ (-x(1)**2+x(1)*x(3)+x(1)*x(4)-x(3)*x(4))*y(2)/
138  $ ((x(1)-x(2))*(x(2)-x(3))*(x(2)-x(4)))+
139  $ (x(1)**2-x(1)*x(2)-x(1)*x(4)+x(2)*x(4))*y(3)/
140  $ ((x(1)-x(3))*(x(2)-x(3))*(x(3)-x(4)))+
141  $ (-x(1)**2+x(1)*x(2)+x(1)*x(3)-x(2)*x(3))*y(4)/
142  $ ((x(1)-x(4))*(x(2)-x(4))*(x(3)-x(4)))
143  y2(1)=-0.5
144  u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1t)
145  ELSE
146  y2(1)=-0.5
147  u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
148  ENDIF
149  DO i=2,n-1
150  sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
151  p=sig*y2(i-1)+2.
152  y2(i)=(sig-1.)/p
153  u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
154  $ /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
155  END DO
156  IF (ypn.gt..99e30) THEN
157  qn=0.
158  un=0.
159  ELSE IF(ypn.lt.-.99e30) THEN
160  ypnt=(-(x(-2+n)*x(-1+n))+x(-2+n)*x(n)+x(-1+n)*x(n)-x(n)**2)*
161  $ y(-3+n)/
162  $ ((-x(-3+n)+x(-2+n))*(-x(-3+n)+x(-1+n))*
163  $ (-x(-3+n)+x(n)))+
164  $ (x(-3+n)*x(-1+n)-x(-3+n)*x(n)-x(-1+n)*x(n)+x(n)**2)*
165  $ y(-2+n)/
166  $ ((-x(-3+n)+x(-2+n))*(-x(-2+n)+x(-1+n))*
167  $ (-x(-2+n)+x(n)))+
168  $ (-(x(-3+n)*x(-2+n))+x(-3+n)*x(n)+x(-2+n)*x(n)-x(n)**2)*
169  $ y(-1+n)/
170  $ ((-x(-3+n)+x(-1+n))*(-x(-2+n)+x(-1+n))*
171  $ (-x(-1+n)+x(n)))+
172  $ (x(-3+n)*x(-2+n)+x(-3+n)*x(-1+n)+x(-2+n)*x(-1+n)-
173  $ 2*(x(-3+n)+x(-2+n)+x(-1+n))*x(n)+3*x(n)**2)*y(n)/
174  $ ((-x(-3+n)+x(n))*(-x(-2+n)+x(n))*(-x(-1+n)+x(n)))
175  qn=0.5
176  un=(3./(x(n)-x(n-1)))*(ypnt-(y(n)-y(n-1))/(x(n)-x(n-1)))
177  ELSE
178  qn=0.5
179  un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
180  ENDIF
181  y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
182  DO k=n-1,1,-1
183  y2(k)=y2(k)*y2(k+1)+u(k)
184  END DO
185  RETURN
186  END
187 c=======================================================================
188  SUBROUTINE splint(xa,ya,y2a,n,x,y,nlow)
189 c-----------------------------------------------------------------------
190 c cubic spline evaluator--spline must be CALLed first to evaluate
191 c y2a
192 c ya is a FUNCTION of xa--both are arrays of LENgth n
193 c ya2[1:n] contains spline coefficients calculated in spline
194 c x is the argument of y[x] WHERE y is to be evaluated
195 c y=y[x] is the RETURNed value
196 c-----------------------------------------------------------------------
197  USE precision
198  IMPLICIT NONE
199  INTEGER :: n,nlow,nhigh,nave
200  REAL(rprec) :: xa(*),ya(*),y2a(*)
201  REAL(rprec) :: x,y
202  REAL(rprec) :: a,b,h
203 
204  IF(nlow.eq.0) THEN
205  nlow=1
206  nhigh=n
207  100 IF (nhigh-nlow.gt.1) THEN
208  nave=(nhigh+nlow)/2
209  IF(xa(nave).gt.x)then
210  nhigh=nave
211  ELSE
212  nlow=nave
213  ENDIF
214  go to 100
215  ENDIF
216  ELSE
217  nhigh=nlow+1
218  ENDIF
219  h=xa(nhigh)-xa(nlow)
220  IF (h.eq.0.) THEN
221  WRITE(6,'( "two xa = in splinet")')
222  stop
223  ENDIF
224  a=(xa(nhigh)-x)/h
225  b=1.-a
226  y=a*ya(nlow)+b*ya(nhigh)+
227  $ ((a-2.)*y2a(nlow)-(a+1.)*y2a(nhigh))*a*b*h*h/6.
228  RETURN
229  END
230 c=======================================================================
231  SUBROUTINE zspline(xa,ya,y2a,n,zg,ng,za)
232 c-----------------------------------------------------------------------
233 c zspline INTegrates the cubic spline of ya[xa]
234 c it assumes that spline has already been CALLed to evaluate y2a
235 c xa[1:n], ya[1:n], y2a[1:n]
236 c In Mathematica notation z[i]=zg+Integrate[y[x],{x,x[ng],x[i]}]
237 c WHERE both zg and ng are input quantities,
238 c za is calculated here.
239 c za can THEN be USEd in zsplINT to determine z at a
240 c specific x
241 c-----------------------------------------------------------------------
242  USE precision
243  IMPLICIT NONE
244 c
245  REAL(rprec) :: xa(n),ya(n),y2a(n),zg,za(n)
246  INTEGER :: n,ng
247 c
248  INTEGER :: j
249  REAL(rprec) :: const
250 c
251  IF (ng .lt. 0 .or. ng .gt. n)stop 'zspline: wrong ng'
252 c
253  za(1)=0.
254  DO j=2,n
255  za(j)=za(j-1)+0.5*(xa(j)-xa(j-1))*(ya(j)+ya(j-1))-
256  > (xa(j)-xa(j-1))**3*(y2a(j)+y2a(j-1))/24.0
257  ENDdo
258 c
259  const=zg-za(ng)
260  DO j=1,n
261  za(j)=za(j)+const
262  ENDdo
263 c
264  RETURN
265  END
266 c=======================================================================
267  SUBROUTINE zsplint(xa,ya,y2a,za,n,x,y,yp,z)
268 c-----------------------------------------------------------------------
269 c evaluate cubic spline to determine FUNCTION (y),
270 c derivative (yp) and INTegral (z) at location x.
271 c first spline must be CALLed to obtain y2a and
272 c zspline must be CALLed to obtain za
273 c-----------------------------------------------------------------------
274  dimension xa(n),ya(n),y2a(n),za(n)
275  klo=1
276  khi=n
277 1 IF (khi-klo.gt.1) THEN
278  k=(khi+klo)/2
279  IF(xa(k).gt.x)then
280  khi=k
281  ELSE
282  klo=k
283  ENDIF
284  GOTO 1
285  ENDIF
286  h=xa(khi)-xa(klo)
287  IF (h.eq.0.) THEN
288  WRITE(6,'(a)')"bad xa input"
289  stop
290  ENDIF
291  a=(xa(khi)-x)/h
292  b=(x-xa(klo))/h
293  y=a*ya(klo)+b*ya(khi)+
294  * ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
295  yp=(ya(khi)-ya(klo))/h-(3.*a**2-1.)/6.*h*y2a(klo)
296  > +(3.*b**2-1.)/6.*h*y2a(khi)
297  z=za(klo)+0.5*h*ya(klo)*(1.-a**2)+0.5*ya(khi)*h*b**2+
298  > y2a(klo)*h**3*(2.*a**2-a**4-1.)/24.-
299  > y2a(khi)*h**3*(2.*b**2-b**4)/24.
300  RETURN
301  END
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17