Line data Source code
1 0 : subroutine lir(z,zi,y,yi,ii,id,nt,l,inter)
2 : c subroutine lir1(z,zi,y,yi,ii,id,nt,l,inter)
3 : c
4 : c
5 : c
6 : c interpolation/extrapolation routine
7 : c
8 : c
9 : c for a such that z=zi(a), sets y(i)=yi(i,a), i=1,ii
10 : c
11 : c zi(n),yi(i,n) must be supplied for n=1,nt and i=1,ii
12 : c id is first dimension of yi
13 : c
14 : c inter is set to 1 for interpolation and 0 for extrapolation
15 : c inter is returned as -1 in case of errors
16 : c
17 : c if l.le.1, scan to find the zi(n) which immediately bound z
18 : c starts at n=1
19 : c if l.gt.1, scan starts from value of n from previous call of lir
20 : c
21 : c
22 : c lir use cubic interpolation/extrapolation unless nt.lt.4
23 : c lir1 use linear interpolation/extrapolation
24 : c
25 : c
26 : c Double precision version.
27 : c ++++++++++++++++++++++++
28 : c
29 : c Dated 10/3/90.
30 : c
31 : c Note: this double precision version of the routine has same name
32 : c as single precision version, but is distinguished by its file name
33 : c
34 : c
35 : implicit double precision(a-h,o-z)
36 : dimension zi(1),y(1),yi(1),a(4)
37 : c
38 : c common defining standard input and output
39 : c
40 : common/cstdio/ istdin, istdou, istdpr, istder
41 : data n/-1/
42 : c
43 0 : il=0
44 0 : go to 1
45 0 : entry lir1(z,zi,y,yi,ii,id,nt,l,inter)
46 0 : il=1
47 : 1 continue
48 0 : ir=1
49 : c
50 : c check nt and reset il if necessary
51 0 : if(nt.lt.2) go to 101
52 0 : if(nt.lt.4) il=1
53 : c
54 : c addressing constants
55 0 : inter=1
56 0 : ir1=ir-1
57 0 : ird=ir*id
58 0 : iir=(ii-1)*ir+1
59 0 : j=(nt-1)*ir+1
60 0 : diff=zi(j)-zi(1)
61 : c
62 : c set index for start of search
63 0 : n=(n-2)*ir+1
64 0 : if(l.le.1.or.n.lt.1) n=1
65 : c
66 : c determine position of z within zi
67 0 : 2 if(n.gt.j) go to 8
68 0 : if(diff) 4,102,3
69 0 : 3 if(zi(n)-z) 5,6,9
70 0 : 4 if(zi(n)-z) 9,6,5
71 0 : 5 n=n+ir
72 0 : go to 2
73 : c
74 : c set y when z lies on a mesh point
75 0 : 6 j=(n-1)*id
76 0 : do 7 i=1,iir
77 0 : y(i)=yi(i+j)
78 0 : 7 if(y(i).eq.0.d0) y(i+ir1)=0.d0
79 0 : go to 30
80 : c
81 : c control when z does not lie on a mesh point
82 0 : 8 inter=0
83 0 : 9 if(n.le.1) inter=0
84 0 : if(il.eq.1) go to 20
85 : c
86 : c cubic interpolation/extrapolation
87 : c
88 : c pivotal point (m) and point (k) closest to z
89 0 : 10 m=n
90 0 : k=3
91 0 : if(n.gt.1+ir) go to 11
92 0 : m=1+ir+ir
93 0 : k=n
94 0 : 11 if(n.lt.j) go to 12
95 0 : m=j-ir
96 0 : k=4
97 : c
98 : c
99 : c weighting factors
100 0 : 12 y1=zi(m-ir*2)
101 0 : y2=zi(m-ir)
102 0 : y3=zi(m)
103 0 : y4=zi(m+ir)
104 : c
105 0 : z1=z-y1
106 0 : z2=z-y2
107 0 : z3=z-y3
108 0 : z4=z-y4
109 : c
110 0 : 13 z12=z1*z2
111 0 : z34=z3*z4
112 : c
113 0 : 14 a(1)=z2*z34/((y1-y2)*(y1-y3)*(y1-y4))
114 0 : a(2)=z1*z34/((y2-y1)*(y2-y3)*(y2-y4))
115 0 : a(3)=z12*z4/((y3-y1)*(y3-y2)*(y3-y4))
116 0 : a(4)=z12*z3/((y4-y1)*(y4-y2)*(y4-y3))
117 : c
118 : c correct a(k)
119 0 : 15 diff=a(1)+a(2)+a(3)+a(4)
120 0 : a(k)=(1.d0+a(k))-diff
121 : c
122 : c compute y
123 0 : 16 m=(m-1)/ir-3
124 0 : m=m*ird
125 0 : do 18 i=1,iir
126 0 : k=i+m
127 0 : yy=0.d0
128 0 : do 17 j=1,4
129 0 : k=k+ird
130 0 : diff=yi(k)
131 0 : 17 yy=yy+a(j)*diff
132 0 : y(i)=yy
133 0 : 18 if(y(i).eq.0.d0) y(i+ir1)=0.d0
134 0 : go to 30
135 : c
136 : c linear interpolation/extrapolation
137 0 : 20 if(n.eq.1) n=1+ir
138 0 : if(n.gt.j) n=j
139 0 : z1=zi(n)
140 0 : y1=(z1-z)/(z1-zi(n-ir))
141 0 : y2=1.0-y1
142 0 : j=(n-1)*id
143 0 : m=j-ird
144 0 : do 21 i=1,iir,ir
145 0 : y(i)=y1*yi(i+m)+y2*yi(i+j)
146 0 : 21 if(y(i).eq.0.d0) y(i+ir1)=0.d0
147 : c
148 : c reset n
149 0 : 30 n=(n+ir-1)/ir
150 0 : return
151 : c
152 : c
153 : c diagnostics
154 0 : 101 if(istdpr.gt.0) write(istdpr,1001) nt
155 0 : inter=-1
156 0 : return
157 0 : 102 if(istdpr.gt.0) write(istdpr,1002) zi(1),nt,zi(j)
158 0 : inter=-1
159 0 : return
160 : c
161 : 1001 format(/1x,10('*'),5x,'there are fewer than two data points in',
162 : * ' lir nt =',i4,5x,10('*')/)
163 : 1002 format(/1x,10('*'),5x,'extreme values of independent variable',
164 : * ' equal in lir',5x,10('*')/16x,'zi(1) =',1pe13.5,', ',
165 : * 'zi(',i4,') =',1pe13.5/)
166 : c
167 : end
|
