Line data Source code
1 24 : subroutine derive(x,y,dydx,nn,iy,idy,ipy,ipdy)
2 : c subroutine deriv2(x,y,dydx,nn,iy,idy,ipy,ipdy)
3 : c
4 : c
5 : c derive sets dydx(1,n) = first derivative of y(1,n) w.r.t. x(n)
6 : c deriv2 sets dydx(1,n) = second derivative of y(1,n) w.r.t. x(n)
7 : c n=1,nn
8 : c
9 : c iy is the first dimension of y in the calling programme
10 : c idy is the first dimension of dydx in the calling programme
11 : c
12 : c second order accuracy differences are used at interior points
13 : c at end points third order differences are used for first
14 : c derivative. second derivative is obtained by quadratic
15 : c interpolation from the interior
16 : c
17 : c
18 : c revised on 6/9 1984 to include scaling by interval length,
19 : c to avoid underflow and consequent divide errors.
20 : c
21 : c ****************************************
22 : c
23 : c notes on precision parameters:
24 : c
25 : c The arguments ipy and ipdy are kept for consistency with previous
26 : c versions of routine. However they have no effect.
27 : c
28 : c .............................................................................
29 : c
30 : c Double precision version.
31 : c ++++++++++++++++++++++++
32 : c
33 : c Dated 10/3/90.
34 : c
35 : c Note: this double precision version of the routine has same name
36 : c as single precision version, but is distinguished by its file name
37 : c
38 : implicit double precision(a-h,o-z)
39 : logical second
40 : dimension x(*), y(*),dydx(*)
41 : c
42 : c common defining standard input and output
43 : c
44 : common/cstdio/ istdin, istdou, istdpr, istder
45 : data epsil/1.d-8/
46 : c
47 : c
48 12 : second=.false.
49 12 : go to 10
50 : c
51 0 : entry deriv2(x,y,dydx,nn,iy,idy,ipy,ipdy)
52 0 : second=.true.
53 : c
54 : c
55 12 : 10 ky=iy*ipy
56 12 : kdy=idy*ipdy
57 12 : if(ipdy.eq.1) go to 30
58 : c
59 : c if dydx is real*8, set insignificant half to zero
60 0 : j=1-kdy
61 0 : 20 do 21 n=1,nn
62 0 : j=j+kdy
63 0 : 21 dydx(j)=0.d0
64 : c
65 12 : 30 n1=nn-1
66 12 : n=0
67 12 : i=1
68 12 : k=1+ky
69 12 : hn =x(2)-x(1)
70 12 : xn= abs(x(1))
71 12 : hna=abs(hn)
72 12 : e=y(k)-y(1)
73 12 : if(second) go to 37
74 : c
75 : c first derivative at end points
76 12 : hn1=x(3)-x(2)
77 12 : h3=x(4)-x(3)
78 12 : nx1=0
79 12 : xn1= abs(x(4))
80 : c
81 12 : hna1=abs(hn1)
82 12 : ha3=abs(h3)
83 : c
84 : c rescale differences
85 : c
86 12 : if(hn.eq.0) go to 361
87 12 : hnin=1.d0/hn
88 12 : hn1=hnin*hn1
89 12 : h3=hnin*h3
90 : c
91 24 : 31 hn2=1.d0+hn1
92 24 : hn3=hn2+h3
93 : c test size of intervals
94 24 : xxa=(xn+xn1)*epsil
95 24 : if(hna.lt.xxa.or.hna1.lt.xxa.or.ha3.lt.xxa) go to 361
96 : c=hn2*hn3*(hn2*((hn1+h3)*(hn3+1.d0)-2.d0*h3-hn1*hn3)-
97 24 : . (h3*h3+hn1*hn3))
98 24 : if(second) go to 34
99 24 : 32 d=hn2*hn2
100 24 : dysc=hnin
101 24 : f=1.d0
102 24 : b=f*h3
103 24 : f=f*hn1
104 24 : hn1=hn3*hn3
105 24 : a=d*hn1*h3
106 24 : b=-hn1*(b+f)
107 24 : d=d*f
108 24 : if(n) 35,35,33
109 12 : 33 c=-c
110 12 : go to 36
111 0 : 34 a=-hn2*hn3*h3*(hn2+hn3)
112 0 : b=hn3*(hn1+h3)*(hn3+1.d0)
113 0 : d=-hn1*hn2*(1.d0+hn2)
114 0 : c=0.5d0*c
115 0 : dysc=hnin*hnin
116 12 : if(n) 35,35,36
117 12 : 35 j=1
118 12 : l=k+ky*2
119 24 : 36 dydx(i)=dysc*(a*e+b*(y(k+ky)-y(j))+d*(y(l)-y(j)))/c
120 24 : go to 362
121 : c zero interval. write diagnostics
122 0 : 361 nj1=nx1+1
123 0 : nj2=nx1+4
124 0 : if(istdpr.gt.0)
125 0 : * write(istdpr,1000) nj1,nj2,(x(nx1+j),j=1,4)
126 0 : dydx(i)=0.d0
127 : c
128 24 : 362 if(n.ne.0) return
129 : c
130 : c derivative at interior points
131 17604 : 37 do 42 n=2,n1
132 17592 : i=i+kdy
133 17592 : j=k
134 17592 : k=k+ky
135 17592 : xn1=xn
136 17592 : xn= abs(x(n))
137 17592 : xxa=(xn1+xn)*epsil
138 17592 : d=e
139 17592 : e=y(k)-y(j)
140 17592 : hn1=hn
141 17592 : hn=x(n+1)-x(n)
142 17592 : hna1=hna
143 17592 : hna=abs(hn)
144 : c test size of intervals
145 17592 : if(hna.ge.xxa.and.hna1.ge.xxa) go to 371
146 0 : dydx(i)=dydx(i-kdy)
147 0 : nj1=n-1
148 0 : nj2=n+1
149 0 : if(istdpr.gt.0)
150 0 : * write(istdpr,1000) nj1,nj2,x(nj1),x(n),x(nj2)
151 0 : go to 42
152 : c
153 : c rescale differences
154 : c
155 17592 : 371 hnin=1.d0/hn
156 17592 : hn1=hnin*hn1
157 : c
158 17592 : c=hn1*(1.d0+hn1)
159 17592 : if(second) go to 39
160 17592 : 38 a=hn1*hn1
161 17592 : b=1.d0
162 17592 : dysc=hnin
163 17592 : go to 40
164 0 : 39 a=hn1
165 0 : b=-1.d0
166 0 : c=0.5d0*c
167 0 : dysc=hnin*hnin
168 17592 : 40 dydx(i)=dysc*(a*e+b*d)/c
169 12 : 42 continue
170 : c
171 12 : h3=x(nn-2)-x(nn-3)
172 : c
173 12 : ha3=abs(h3)
174 12 : h3=hnin*h3
175 : c
176 12 : xn= abs(x(nn))
177 12 : xn1= abs(x(nn-3))
178 12 : nx1=nn-4
179 12 : if(second) go to 50
180 : c
181 : c storage indices for first derivative at last point
182 12 : i=i+kdy
183 12 : j=k
184 12 : k=k-ky*3
185 12 : l=k
186 12 : e=-e
187 12 : go to 31
188 : c
189 : c second derivative at end points
190 0 : 50 j=i
191 0 : k=i-kdy
192 0 : l=k-kdy
193 0 : i=i+kdy
194 0 : 51 a=1.d0+hn1
195 0 : b=a+h3
196 0 : c=hn1+h3
197 0 : xxa=(xn1+xn)*epsil
198 : c test size of intervals
199 0 : if(hna.ge.xxa.and.hna1.ge.xxa.and.ha3.ge.xxa) go to 52
200 0 : 51100 nj1=nx1+1
201 0 : nj2=nx1+4
202 0 : if(istdpr.gt.0)
203 0 : * write(istdpr,1000) nj1,nj2,(x(nx1+j),j=1,4)
204 0 : dydx(i)=0.d0
205 0 : go to 53
206 : 52 dydx(i)=(a*b/(hn1*c))*dydx(j)-(b/(hn1*h3))*dydx(k)
207 0 : . +(a/(c*h3))*dydx(l)
208 0 : 53 if(i.eq.1) return
209 0 : i=1
210 0 : j=i+kdy
211 0 : k=j+kdy
212 0 : l=k+kdy
213 0 : hn=x(2)-x(1)
214 0 : hn1=x(3)-x(2)
215 0 : h3=x(4)-x(3)
216 0 : xn= abs(x(1))
217 0 : xn1= abs(x(4))
218 0 : nx1=0
219 0 : hna=abs(hn)
220 0 : hna1=abs(hn1)
221 0 : ha3=abs(h3)
222 : c
223 : c rescale differences
224 : c
225 0 : if(hn.eq.0) go to 51100
226 0 : hnin=1.d0/hn
227 0 : hn1=hnin*hn1
228 0 : h3=hnin*h3
229 : c
230 0 : go to 51
231 : c
232 : 1000 format(' **** from derive: degeneracy among x(',i5,') - x(',
233 : . i5,') = ',1p4e16.8)
234 : end
|
