##// END OF EJS Templates
schain
RSS
lagp.f update 2
rflores -
r1621:c024a8a81a90
parent child
Show More
Side by Side Unified
Context file:
r1621:c024a8a81a90 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,566 +1,566
1 subroutine lagp(plag,wl,r0,dr,nl,nrange)
1 subroutine lagp(plag,wl,r0,dr,nl,nrange)
2 c
2 c
3 c predict lagged product matrix using state parameters stored
3 c predict lagged product matrix using state parameters stored
4 c as cubic b-splines
4 c as cubic b-splines
5 c read weights for numeric integration from a file
5 c read weights for numeric integration from a file
6 c
6 c
7 c faster version - assumes samples taken at half-integer lags, ranges
7 c faster version - assumes samples taken at half-integer lags, ranges
8 c
8 c
9 include 'fitter.h'
9 include 'fitter.h'
10 include 'spline.h'
10 include 'spline.h'
11 include 'bfield.h'
11 include 'bfield.h'
12 integer nl,nrange
12 integer nl,nrange
13 real plag(nl,nrange),taus(nl)
13 real plag(nl,nrange),taus(nl)
14 parameter(nw=1000)
14 parameter(nw=1000)
15 integer ilag(nw)
15 integer ilag(nw)
16 real weight(nw),dtau(nw),drange(nw)
16 real weight(nw),dtau(nw),drange(nw)
17 logical first
17 logical first
18 real store(2*nl,2*nrange,2)
18 real store(2*nl,2*nrange,2)
19
19
20 character(1024) :: fqual_temp
20 character(1024) :: fqual_temp
21 character(:), allocatable :: fqual
21 character(:), allocatable :: fqual
22
22
23 data first/.true./
23 data first/.true./
24 c
24 c
25
25
26 c write(*,*) "dens_before: ",dens
26 c write(*,*) "dens_before: ",dens
27 c call exit
27 c call exit
28 c zero out array
28 c zero out array
29 c write(*,*) "Starting lagp"
29 c write(*,*) "Starting lagp"
30 do i=1,nrange
30 do i=1,nrange
31 do j=1,nl
31 do j=1,nl
32 plag(j,i)=0.0
32 plag(j,i)=0.0
33 end do
33 end do
34 end do
34 end do
35
35
36 c
36 c
37 c compute all required acfs
37 c compute all required acfs
38 c
38 c
39 do i=1,2*nrange
39 do i=1,2*nrange
40 ii=(i-1)/2+1
40 ii=(i-1)/2+1
41 alt=r0+float(i-1)*dr*0.5 ! half integer ranges
41 alt=r0+float(i-1)*dr*0.5 ! half integer ranges
42 c write(*,*) "dens_before: ",dens
42 c write(*,*) "dens_before: ",dens
43 c call exit
43 c call exit
44 c write(*,*) "alt: ",alt
44 c write(*,*) "alt: ",alt
45 c write(*,*) "dens: ",dens
45 c write(*,*) "dens: ",dens
46 c write(*,*) "te: ",te
46 c write(*,*) "te: ",te
47 c write(*,*) "ti1: ",ti1
47 c write(*,*) "ti1: ",ti1
48 c write(*,*) "hf: ",hf
48 c write(*,*) "hf: ",hf
49 c write(*,*) "hef: ",hef
49 c write(*,*) "hef: ",hef
50 c call exit
50 c call exit
51 call get_spline(alt,dens,te,ti1,hf,hef)
51 call get_spline(alt,dens,te,ti1,hf,hef)
52 c write(*,*) "dens_after: ",dens
52 c write(*,*) "dens_after: ",dens
53 c call exit
53 c call exit
54 do k=1,nion
54 do k=1,nion
55 ti(k)=ti1
55 ti(k)=ti1
56 end do
56 end do
57 fi(2)=hf
57 fi(2)=hf
58 fi(3)=hef
58 fi(3)=hef
59 fi(1)=1.0-hf-hef
59 fi(1)=1.0-hf-hef
60
60
61 do j=1,2*nl
61 do j=1,2*nl
62 tau=float(j-1)*(dr/1.5e5)*0.5 ! half integer lags
62 tau=float(j-1)*(dr/1.5e5)*0.5 ! half integer lags
63
63
64 c just consider a single azimuth angle for now
64 c just consider a single azimuth angle for now
65
65
66 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
66 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
67 c & alpha_prof(ii)+1.74e-2,dens,bfld_prof(ii),rho)
67 c & alpha_prof(ii)+1.74e-2,dens,bfld_prof(ii),rho)
68 c store(j,i,1)=rho*dens*(100.0/alt)**2
68 c store(j,i,1)=rho*dens*(100.0/alt)**2
69
69
70 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
70 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
71 c & alpha_prof(ii)-1.74e-2,dens,bfld_prof(ii),rho)
71 c & alpha_prof(ii)-1.74e-2,dens,bfld_prof(ii),rho)
72 c store(j,i,2)=rho*dens*(100.0/alt)**2
72 c store(j,i,2)=rho*dens*(100.0/alt)**2
73
73
74 call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
74 call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
75 & alpha_prof(ii),dens,bfld_prof(ii),rho)
75 & alpha_prof(ii),dens,bfld_prof(ii),rho)
76 store(j,i,1)=rho*dens*(100.0/alt)**2
76 store(j,i,1)=rho*dens*(100.0/alt)**2
77 c write(*,*) "store", store(j,i,1)
77 c write(*,*) "store", store(j,i,1)
78 c call exit
78 c call exit
79 end do
79 end do
80 end do
80 end do
81
81
82 c construct lagged products
82 c construct lagged products
83 c
83 c
84 c write(*,*) "before weights"
84 c write(*,*) "before weights"
85 c write(*,*) first
85 c write(*,*) first
86 c write(*,*) "weight before: ",weight
86 c write(*,*) "weight before: ",weight
87 c write(*,*) "BEFORE GET_PATH"
87 c write(*,*) "BEFORE GET_PATH"
88 call get_path(fqual_temp)
88 call get_path(fqual_temp)
89 c write(*,*) "L_BEF: ", fqual_temp, "L_BEF_end"
89 c write(*,*) "L_BEF: ", fqual_temp, "L_BEF_end"
90 fqual = TRIM(fqual_temp)
90 fqual = TRIM(fqual_temp)
91 fqual = fqual//'/weights.dat'
91 fqual = fqual//'weights.dat'
92 c write(*,*) "Final fqual",fqual,"FINAL"
92 c write(*,*) "Final fqual",fqual,"FINAL"
93
93
94 if(first) then
94 if(first) then
95 open(unit=25,file=fqual,
95 open(unit=25,file=fqual,
96 & status='old')
96 & status='old')
97 read (25,*) nwt
97 read (25,*) nwt
98 do i=1,nwt
98 do i=1,nwt
99 read(25,*) weight(i),dtau(i),drange(i),ilag(i)
99 read(25,*) weight(i),dtau(i),drange(i),ilag(i)
100 c write(*,*) "ilag: ",ilag(i)
100 c write(*,*) "ilag: ",ilag(i)
101 end do
101 end do
102 c flush(25)
102 c flush(25)
103 close(25)
103 close(25)
104 c first=.false.
104 c first=.false.
105 end if
105 end if
106
106
107 c write(*,*) "after weights"
107 c write(*,*) "after weights"
108
108
109 c write(*,*) "weight(3): ",weight(3)
109 c write(*,*) "weight(3): ",weight(3)
110 c call exit
110 c call exit
111 c write(*,*) "nl: ",nl,"nrange: ",nrange,"nwt: ",nwt
111 c write(*,*) "nl: ",nl,"nrange: ",nrange,"nwt: ",nwt
112
112
113 do i=1+nl,nrange ! don't worry about redundancy here
113 do i=1+nl,nrange ! don't worry about redundancy here
114
114
115 do j=1,nwt
115 do j=1,nwt
116 c write(*,*) "j: ",j
116 c write(*,*) "j: ",j
117 alt=r0+(float(i-1)+drange(j))*dr
117 alt=r0+(float(i-1)+drange(j))*dr
118 tau=abs(dtau(j)*dr/1.5e5)
118 tau=abs(dtau(j)*dr/1.5e5)
119 lag=ilag(j)+1
119 lag=ilag(j)+1
120 c write(*,*) tau
120 c write(*,*) tau
121
121
122 k=1+2.0*(tau/(dr/1.5e5))
122 k=1+2.0*(tau/(dr/1.5e5))
123 l=1+2.0*(alt-r0)/dr
123 l=1+2.0*(alt-r0)/dr
124
124
125 do m=1,1 ! 2
125 do m=1,1 ! 2
126
126
127 c call exit
127 c call exit
128 c write(*,*) i
128 c write(*,*) i
129 c call exit
129 c call exit
130 c if (i .eq. nrange) then
130 c if (i .eq. nrange) then
131 c if (j .eq. 227) then
131 c if (j .eq. 227) then
132
132
133 c write(*,*) "lag: ",lag
133 c write(*,*) "lag: ",lag
134 c write(*,*) "ilag: ",ilag(j)
134 c write(*,*) "ilag: ",ilag(j)
135 c call exit
135 c call exit
136 c write(*,*) "plag: ",plag(lag,i)
136 c write(*,*) "plag: ",plag(lag,i)
137 c write(*,*) "weight: ",weight
137 c write(*,*) "weight: ",weight
138 c
138 c
139 c write(*,*) "weight: ",weight(j)
139 c write(*,*) "weight: ",weight(j)
140 c write(*,*) "store: ",store(k,l,m)
140 c write(*,*) "store: ",store(k,l,m)
141 c call exit
141 c call exit
142 c end if
142 c end if
143 c end if
143 c end if
144 plag(lag,i)=plag(lag,i)+weight(j)*store(k,l,m)
144 plag(lag,i)=plag(lag,i)+weight(j)*store(k,l,m)
145 c if (i .eq. nrange) then
145 c if (i .eq. nrange) then
146 c write(*,*) "lag: ",lag
146 c write(*,*) "lag: ",lag
147 c write(*,*) "plag: ",plag(lag,i)
147 c write(*,*) "plag: ",plag(lag,i)
148 c write(*,*) "plag: ",plag(lag,i)
148 c write(*,*) "plag: ",plag(lag,i)
149 c call exit
149 c call exit
150 c end if
150 c end if
151 end do
151 end do
152 c write(*,*) "plag: ",plag(lag,i)
152 c write(*,*) "plag: ",plag(lag,i)
153 c write(*,*) "plag: ",plag(12,72)
153 c write(*,*) "plag: ",plag(12,72)
154 c call exit
154 c call exit
155 end do
155 end do
156 c call exit
156 c call exit
157 c write(*,*) "plag: ",plag(12,72)
157 c write(*,*) "plag: ",plag(12,72)
158 end do
158 end do
159 c write(*,*) "plag: ",plag(12,72)
159 c write(*,*) "plag: ",plag(12,72)
160 c call exit
160 c call exit
161 c write(*,*) "End LAGP"
161 c write(*,*) "End LAGP"
162 return
162 return
163 end
163 end
164
164
165 subroutine lagp_old(plag,wl,r0,dr,nl,nrange)
165 subroutine lagp_old(plag,wl,r0,dr,nl,nrange)
166 c
166 c
167 c predict lagged product matrix using state parameters stored
167 c predict lagged product matrix using state parameters stored
168 c as cubic b-splines
168 c as cubic b-splines
169 c read weights for numeric integration from a file
169 c read weights for numeric integration from a file
170 c
170 c
171 c general version - samples can be anywhere
171 c general version - samples can be anywhere
172 c
172 c
173 include 'fitter.h'
173 include 'fitter.h'
174 include 'spline.h'
174 include 'spline.h'
175 include 'bfield.h'
175 include 'bfield.h'
176 integer nl,nrange
176 integer nl,nrange
177 real plag(nl,nrange),taus(nl)
177 real plag(nl,nrange),taus(nl)
178 parameter(nw=1000)
178 parameter(nw=1000)
179 integer ilag(nw)
179 integer ilag(nw)
180 real weight(nw),dtau(nw),drange(nw)
180 real weight(nw),dtau(nw),drange(nw)
181 logical first
181 logical first
182 data first/.true./
182 data first/.true./
183 c
183 c
184
184
185 c zero out array
185 c zero out array
186
186
187 do i=1,nrange
187 do i=1,nrange
188 do j=1,nl
188 do j=1,nl
189 plag(j,i)=0.0
189 plag(j,i)=0.0
190 end do
190 end do
191 end do
191 end do
192
192
193 c
193 c
194 c construct lagged products
194 c construct lagged products
195 c
195 c
196
196
197 if(first) then
197 if(first) then
198 open(unit=25,file='weights.dat',status='old')
198 open(unit=25,file='weights.dat',status='old')
199 read (25,*) nwt
199 read (25,*) nwt
200 do i=1,nwt
200 do i=1,nwt
201 read(25,*) weight(i),dtau(i),drange(i),ilag(i)
201 read(25,*) weight(i),dtau(i),drange(i),ilag(i)
202 end do
202 end do
203 close(25)
203 close(25)
204 first=.false.
204 first=.false.
205 end if
205 end if
206
206
207 c can avoid recalculating splines and/or ACFs when parameters repeat
207 c can avoid recalculating splines and/or ACFs when parameters repeat
208
208
209 altp=-1.0
209 altp=-1.0
210 taup=-1.0
210 taup=-1.0
211
211
212 do i=1+nl,nrange
212 do i=1+nl,nrange
213 do j=1,nwt
213 do j=1,nwt
214 alt=r0+(float(i-1)+drange(j))*dr
214 alt=r0+(float(i-1)+drange(j))*dr
215 ii=(alt-r0)/dr+1
215 ii=(alt-r0)/dr+1
216
216
217 if(alt.ne.altp) then
217 if(alt.ne.altp) then
218 call get_spline(alt,dens,te,ti1,hf,hef)
218 call get_spline(alt,dens,te,ti1,hf,hef)
219
219
220 c write(*,*) alt,dens,te,ti1,hf,hef
220 c write(*,*) alt,dens,te,ti1,hf,hef
221
221
222 do k=1,nion
222 do k=1,nion
223 ti(k)=ti1
223 ti(k)=ti1
224 end do
224 end do
225 fi(2)=hf
225 fi(2)=hf
226 fi(3)=hef
226 fi(3)=hef
227 fi(1)=1.0-hf-hef
227 fi(1)=1.0-hf-hef
228 end if
228 end if
229
229
230 tau=abs(dtau(j)*dr/1.5e5)
230 tau=abs(dtau(j)*dr/1.5e5)
231 lag=ilag(j)+1
231 lag=ilag(j)+1
232
232
233 c write(*,*) i,alpha_prof(i),bfld_prof(i)
233 c write(*,*) i,alpha_prof(i),bfld_prof(i)
234
234
235 if(tau.ne.taup.or.alt.ne.altp) then
235 if(tau.ne.taup.or.alt.ne.altp) then
236
236
237 c two- or three- point Gauss Hermite quadrature rule
237 c two- or three- point Gauss Hermite quadrature rule
238
238
239 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
239 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
240 c & alpha_prof(ii)+1.74e-2,dens,bfld_prof(ii),rho1)
240 c & alpha_prof(ii)+1.74e-2,dens,bfld_prof(ii),rho1)
241
241
242 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
242 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
243 c & alpha_prof(ii)-1.74e-2,dens,bfld_prof(ii),rho2)
243 c & alpha_prof(ii)-1.74e-2,dens,bfld_prof(ii),rho2)
244
244
245 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
245 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
246 c & alpha_prof(ii),dens,bfld_prof(ii),rho1)
246 c & alpha_prof(ii),dens,bfld_prof(ii),rho1)
247
247
248 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
248 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
249 c & alpha_prof(ii)+1.9e-2,dens,bfld_prof(ii),rho2)
249 c & alpha_prof(ii)+1.9e-2,dens,bfld_prof(ii),rho2)
250
250
251 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
251 c call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
252 c & alpha_prof(ii)-1.9e-2,dens,bfld_prof(ii),rho3)
252 c & alpha_prof(ii)-1.9e-2,dens,bfld_prof(ii),rho3)
253
253
254 call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
254 call acf2(wl,tau,te,ti,fi,ven,vin,wi,nion,
255 & alpha_prof(ii),dens,bfld_prof(ii),rho)
255 & alpha_prof(ii),dens,bfld_prof(ii),rho)
256
256
257 c rho=(rho1+rho2)/2.0
257 c rho=(rho1+rho2)/2.0
258 c rho=(rho2+rho3)*0.29541+rho1*1.18164)/1.77246
258 c rho=(rho2+rho3)*0.29541+rho1*1.18164)/1.77246
259
259
260 end if
260 end if
261
261
262 plag(lag,i)=plag(lag,i)+weight(j)*rho*dens*(100.0/alt)**2
262 plag(lag,i)=plag(lag,i)+weight(j)*rho*dens*(100.0/alt)**2
263
263
264 altp=alt
264 altp=alt
265 taup=tau
265 taup=tau
266
266
267 end do
267 end do
268 end do
268 end do
269
269
270 return
270 return
271 end
271 end
272
272
273 function atanh(x)
273 function atanh(x)
274 c
274 c
275 real atanh,x
275 real atanh,x
276 c
276 c
277 atanh=log(sqrt((1.+x)/(1.-x)))
277 atanh=log(sqrt((1.+x)/(1.-x)))
278 return
278 return
279 end
279 end
280
280
281 subroutine get_spline(alt,dens,te,ti,hf,hef)
281 subroutine get_spline(alt,dens,te,ti,hf,hef)
282 c
282 c
283 c routines for handling cubic b-spline interpolation
283 c routines for handling cubic b-spline interpolation
284 c gets values consistent with stored coefficients
284 c gets values consistent with stored coefficients
285 c
285 c
286
286
287 include 'spline.h'
287 include 'spline.h'
288
288
289 c
289 c
290 c bspline values from 200-1500 km in 15-km intervals
290 c bspline values from 200-1500 km in 15-km intervals
291 c must specify five knots above ceiling
291 c must specify five knots above ceiling
292 c note offset ... start accessing splines above bottom two knots
292 c note offset ... start accessing splines above bottom two knots
293 c
293 c
294 c five banks of splines, for Ne, Te, Ti, H+, and He+ versus altitude
294 c five banks of splines, for Ne, Te, Ti, H+, and He+ versus altitude
295 c need to initialize ta somewhere
295 c need to initialize ta somewhere
296 c
296 c
297 c write(*,*) "ta: ",ta
297 c write(*,*) "ta: ",ta
298 c write(*,*) "bcoef(1,1): ",bcoef(1,1)
298 c write(*,*) "bcoef(1,1): ",bcoef(1,1)
299 c write(*,*) "nspline: ",nspline
299 c write(*,*) "nspline: ",nspline
300 c write(*,*) "norder: ",norder
300 c write(*,*) "norder: ",norder
301 c write(*,*) "alt: ",alt
301 c write(*,*) "alt: ",alt
302 dens=bvalue(ta,bcoef(1,1),nspline,norder,alt,0)
302 dens=bvalue(ta,bcoef(1,1),nspline,norder,alt,0)
303 c write(*,*) "dens_inside_get_spline: ",dens
303 c write(*,*) "dens_inside_get_spline: ",dens
304 c write(*,*) dens
304 c write(*,*) dens
305 c call exit
305 c call exit
306 dens=10.0**MAX(dens,2.0)
306 dens=10.0**MAX(dens,2.0)
307 c write(*,*) dens
307 c write(*,*) dens
308 c call exit
308 c call exit
309 te=bvalue(ta,bcoef(1,2),nspline,norder,alt,0)
309 te=bvalue(ta,bcoef(1,2),nspline,norder,alt,0)
310 te=t0+t1*(1.0+tanh(te))/2.0 ! DLH 10/14 was 3500
310 te=t0+t1*(1.0+tanh(te))/2.0 ! DLH 10/14 was 3500
311 c tr=bvalue(ta,bcoef(1,3),nspline,norder,alt,0)
311 c tr=bvalue(ta,bcoef(1,3),nspline,norder,alt,0)
312 c tr=exp(tr)
312 c tr=exp(tr)
313 c ti=te/tr
313 c ti=te/tr
314 ti=bvalue(ta,bcoef(1,3),nspline,norder,alt,0)
314 ti=bvalue(ta,bcoef(1,3),nspline,norder,alt,0)
315 ti=t0+t1*(1.0+tanh(ti))/2.0 ! DLH 10/14 was 3500
315 ti=t0+t1*(1.0+tanh(ti))/2.0 ! DLH 10/14 was 3500
316 hf=bvalue(ta,bcoef(1,4),nspline,norder,alt,0)
316 hf=bvalue(ta,bcoef(1,4),nspline,norder,alt,0)
317 hf=(1.0+tanh(hf))/2.0
317 hf=(1.0+tanh(hf))/2.0
318 hef=bvalue(ta,bcoef(1,5),nspline,norder,alt,0)
318 hef=bvalue(ta,bcoef(1,5),nspline,norder,alt,0)
319 hef=(1.0+tanh(hef))/2.0
319 hef=(1.0+tanh(hef))/2.0
320 return
320 return
321 end
321 end
322
322
323
323
324 real function bvalue ( t, bcoef, n, k, x, jderiv )
324 real function bvalue ( t, bcoef, n, k, x, jderiv )
325 c from * a practical guide to splines * by c. de boor
325 c from * a practical guide to splines * by c. de boor
326 calls interv
326 calls interv
327 c
327 c
328 calculates value at x of jderiv-th derivative of spline from b-repr.
328 calculates value at x of jderiv-th derivative of spline from b-repr.
329 c the spline is taken to be continuous from the right, EXCEPT at the
329 c the spline is taken to be continuous from the right, EXCEPT at the
330 c rightmost knot, where it is taken to be continuous from the left.
330 c rightmost knot, where it is taken to be continuous from the left.
331 c
331 c
332 c****** i n p u t ******
332 c****** i n p u t ******
333 c t, bcoef, n, k......forms the b-representation of the spline f to
333 c t, bcoef, n, k......forms the b-representation of the spline f to
334 c be evaluated. specifically,
334 c be evaluated. specifically,
335 c t.....knot sequence, of length n+k, assumed nondecreasing.
335 c t.....knot sequence, of length n+k, assumed nondecreasing.
336 c bcoef.....b-coefficient sequence, of length n .
336 c bcoef.....b-coefficient sequence, of length n .
337 c n.....length of bcoef and dimension of spline(k,t),
337 c n.....length of bcoef and dimension of spline(k,t),
338 c a s s u m e d positive .
338 c a s s u m e d positive .
339 c k.....order of the spline .
339 c k.....order of the spline .
340 c
340 c
341 c w a r n i n g . . . the restriction k .le. kmax (=20) is imposed
341 c w a r n i n g . . . the restriction k .le. kmax (=20) is imposed
342 c arbitrarily by the dimension statement for aj, dl, dr below,
342 c arbitrarily by the dimension statement for aj, dl, dr below,
343 c but is n o w h e r e c h e c k e d for.
343 c but is n o w h e r e c h e c k e d for.
344 c
344 c
345 c x.....the point at which to evaluate .
345 c x.....the point at which to evaluate .
346 c jderiv.....integer giving the order of the derivative to be evaluated
346 c jderiv.....integer giving the order of the derivative to be evaluated
347 c a s s u m e d to be zero or positive.
347 c a s s u m e d to be zero or positive.
348 c
348 c
349 c****** o u t p u t ******
349 c****** o u t p u t ******
350 c bvalue.....the value of the (jderiv)-th derivative of f at x .
350 c bvalue.....the value of the (jderiv)-th derivative of f at x .
351 c
351 c
352 c****** m e t h o d ******
352 c****** m e t h o d ******
353 c The nontrivial knot interval (t(i),t(i+1)) containing x is lo-
353 c The nontrivial knot interval (t(i),t(i+1)) containing x is lo-
354 c cated with the aid of interv . The k b-coeffs of f relevant for
354 c cated with the aid of interv . The k b-coeffs of f relevant for
355 c this interval are then obtained from bcoef (or taken to be zero if
355 c this interval are then obtained from bcoef (or taken to be zero if
356 c not explicitly available) and are then differenced jderiv times to
356 c not explicitly available) and are then differenced jderiv times to
357 c obtain the b-coeffs of (d**jderiv)f relevant for that interval.
357 c obtain the b-coeffs of (d**jderiv)f relevant for that interval.
358 c Precisely, with j = jderiv, we have from x.(12) of the text that
358 c Precisely, with j = jderiv, we have from x.(12) of the text that
359 c
359 c
360 c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) )
360 c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) )
361 c
361 c
362 c where
362 c where
363 c / bcoef(.), , j .eq. 0
363 c / bcoef(.), , j .eq. 0
364 c /
364 c /
365 c bcoef(.,j) = / bcoef(.,j-1) - bcoef(.-1,j-1)
365 c bcoef(.,j) = / bcoef(.,j-1) - bcoef(.-1,j-1)
366 c / ----------------------------- , j .gt. 0
366 c / ----------------------------- , j .gt. 0
367 c / (t(.+k-j) - t(.))/(k-j)
367 c / (t(.+k-j) - t(.))/(k-j)
368 c
368 c
369 c Then, we use repeatedly the fact that
369 c Then, we use repeatedly the fact that
370 c
370 c
371 c sum ( a(.)*b(.,m,t)(x) ) = sum ( a(.,x)*b(.,m-1,t)(x) )
371 c sum ( a(.)*b(.,m,t)(x) ) = sum ( a(.,x)*b(.,m-1,t)(x) )
372 c with
372 c with
373 c (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
373 c (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
374 c a(.,x) = ---------------------------------------
374 c a(.,x) = ---------------------------------------
375 c (x - t(.)) + (t(.+m-1) - x)
375 c (x - t(.)) + (t(.+m-1) - x)
376 c
376 c
377 c to write (d**j)f(x) eventually as a linear combination of b-splines
377 c to write (d**j)f(x) eventually as a linear combination of b-splines
378 c of order 1 , and the coefficient for b(i,1,t)(x) must then be the
378 c of order 1 , and the coefficient for b(i,1,t)(x) must then be the
379 c desired number (d**j)f(x). (see x.(17)-(19) of text).
379 c desired number (d**j)f(x). (see x.(17)-(19) of text).
380 c
380 c
381 parameter (kmax = 40)
381 parameter (kmax = 40)
382 integer jderiv,k,n, i,ilo,imk,j,jc,jcmin,jcmax,jj,kmj,km1
382 integer jderiv,k,n, i,ilo,imk,j,jc,jcmin,jcmax,jj,kmj,km1
383 * ,mflag,nmi,jdrvp1
383 * ,mflag,nmi,jdrvp1
384 c integer kmax
384 c integer kmax
385 C real bcoef(n),t(1),x, aj(20),dl(20),dr(20),fkmj
385 C real bcoef(n),t(1),x, aj(20),dl(20),dr(20),fkmj
386 real bcoef(n),x, aj(kmax),dl(kmax),dr(kmax),fkmj
386 real bcoef(n),x, aj(kmax),dl(kmax),dr(kmax),fkmj
387 dimension t(n+k)
387 dimension t(n+k)
388 c former fortran standard made it impossible to specify the length of t
388 c former fortran standard made it impossible to specify the length of t
389 c precisely without the introduction of otherwise superfluous addition-
389 c precisely without the introduction of otherwise superfluous addition-
390 c al arguments.
390 c al arguments.
391 bvalue = 0.
391 bvalue = 0.
392 if (jderiv .ge. k) go to 99
392 if (jderiv .ge. k) go to 99
393 c
393 c
394 c *** Find i s.t. 1 .le. i .lt. n+k and t(i) .lt. t(i+1) and
394 c *** Find i s.t. 1 .le. i .lt. n+k and t(i) .lt. t(i+1) and
395 c t(i) .le. x .lt. t(i+1) . If no such i can be found, x lies
395 c t(i) .le. x .lt. t(i+1) . If no such i can be found, x lies
396 c outside the support of the spline f , hence bvalue = 0.
396 c outside the support of the spline f , hence bvalue = 0.
397 c (The asymmetry in this choice of i makes f rightcontinuous, except
397 c (The asymmetry in this choice of i makes f rightcontinuous, except
398 c at t(n+k) where it is leftcontinuous.)
398 c at t(n+k) where it is leftcontinuous.)
399 call interv ( t, n+k, x, i, mflag )
399 call interv ( t, n+k, x, i, mflag )
400 if (mflag .ne. 0) go to 99
400 if (mflag .ne. 0) go to 99
401 c *** if k = 1 (and jderiv = 0), bvalue = bcoef(i).
401 c *** if k = 1 (and jderiv = 0), bvalue = bcoef(i).
402 km1 = k - 1
402 km1 = k - 1
403 if (km1 .gt. 0) go to 1
403 if (km1 .gt. 0) go to 1
404 bvalue = bcoef(i)
404 bvalue = bcoef(i)
405 go to 99
405 go to 99
406 c
406 c
407 c *** store the k b-spline coefficients relevant for the knot interval
407 c *** store the k b-spline coefficients relevant for the knot interval
408 c (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j),
408 c (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j),
409 c dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
409 c dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
410 c from input to zero. set any t.s not obtainable equal to t(1) or
410 c from input to zero. set any t.s not obtainable equal to t(1) or
411 c to t(n+k) appropriately.
411 c to t(n+k) appropriately.
412 1 jcmin = 1
412 1 jcmin = 1
413 imk = i - k
413 imk = i - k
414 if (imk .ge. 0) go to 8
414 if (imk .ge. 0) go to 8
415 jcmin = 1 - imk
415 jcmin = 1 - imk
416 do 5 j=1,i
416 do 5 j=1,i
417 5 dl(j) = x - t(i+1-j)
417 5 dl(j) = x - t(i+1-j)
418 do 6 j=i,km1
418 do 6 j=i,km1
419 aj(k-j) = 0.
419 aj(k-j) = 0.
420 6 dl(j) = dl(i)
420 6 dl(j) = dl(i)
421 go to 10
421 go to 10
422 8 do 9 j=1,km1
422 8 do 9 j=1,km1
423 9 dl(j) = x - t(i+1-j)
423 9 dl(j) = x - t(i+1-j)
424 c
424 c
425 10 jcmax = k
425 10 jcmax = k
426 nmi = n - i
426 nmi = n - i
427 if (nmi .ge. 0) go to 18
427 if (nmi .ge. 0) go to 18
428 jcmax = k + nmi
428 jcmax = k + nmi
429 do 15 j=1,jcmax
429 do 15 j=1,jcmax
430 15 dr(j) = t(i+j) - x
430 15 dr(j) = t(i+j) - x
431 do 16 j=jcmax,km1
431 do 16 j=jcmax,km1
432 aj(j+1) = 0.
432 aj(j+1) = 0.
433 16 dr(j) = dr(jcmax)
433 16 dr(j) = dr(jcmax)
434 go to 20
434 go to 20
435 18 do 19 j=1,km1
435 18 do 19 j=1,km1
436 19 dr(j) = t(i+j) - x
436 19 dr(j) = t(i+j) - x
437 c
437 c
438 20 do 21 jc=jcmin,jcmax
438 20 do 21 jc=jcmin,jcmax
439 21 aj(jc) = bcoef(imk + jc)
439 21 aj(jc) = bcoef(imk + jc)
440 c
440 c
441 c *** difference the coefficients jderiv times.
441 c *** difference the coefficients jderiv times.
442 if (jderiv .eq. 0) go to 30
442 if (jderiv .eq. 0) go to 30
443 do 23 j=1,jderiv
443 do 23 j=1,jderiv
444 kmj = k-j
444 kmj = k-j
445 fkmj = float(kmj)
445 fkmj = float(kmj)
446 ilo = kmj
446 ilo = kmj
447 do 23 jj=1,kmj
447 do 23 jj=1,kmj
448 aj(jj) = ((aj(jj+1) - aj(jj))/(dl(ilo) + dr(jj)))*fkmj
448 aj(jj) = ((aj(jj+1) - aj(jj))/(dl(ilo) + dr(jj)))*fkmj
449 23 ilo = ilo - 1
449 23 ilo = ilo - 1
450 c
450 c
451 c *** compute value at x in (t(i),t(i+1)) of jderiv-th derivative,
451 c *** compute value at x in (t(i),t(i+1)) of jderiv-th derivative,
452 c given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).
452 c given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).
453 30 if (jderiv .eq. km1) go to 39
453 30 if (jderiv .eq. km1) go to 39
454 jdrvp1 = jderiv + 1
454 jdrvp1 = jderiv + 1
455 do 33 j=jdrvp1,km1
455 do 33 j=jdrvp1,km1
456 kmj = k-j
456 kmj = k-j
457 ilo = kmj
457 ilo = kmj
458 do 33 jj=1,kmj
458 do 33 jj=1,kmj
459 aj(jj) = (aj(jj+1)*dl(ilo) + aj(jj)*dr(jj))/(dl(ilo)+dr(jj))
459 aj(jj) = (aj(jj+1)*dl(ilo) + aj(jj)*dr(jj))/(dl(ilo)+dr(jj))
460 33 ilo = ilo - 1
460 33 ilo = ilo - 1
461 39 bvalue = aj(1)
461 39 bvalue = aj(1)
462 c
462 c
463 99 return
463 99 return
464 end
464 end
465 subroutine interv ( xt, lxt, x, left, mflag )
465 subroutine interv ( xt, lxt, x, left, mflag )
466 c from * a practical guide to splines * by C. de Boor
466 c from * a practical guide to splines * by C. de Boor
467 computes left = max( i : xt(i) .lt. xt(lxt) .and. xt(i) .le. x ) .
467 computes left = max( i : xt(i) .lt. xt(lxt) .and. xt(i) .le. x ) .
468 c
468 c
469 c****** i n p u t ******
469 c****** i n p u t ******
470 c xt.....a real sequence, of length lxt , assumed to be nondecreasing
470 c xt.....a real sequence, of length lxt , assumed to be nondecreasing
471 c lxt.....number of terms in the sequence xt .
471 c lxt.....number of terms in the sequence xt .
472 c x.....the point whose location with respect to the sequence xt is
472 c x.....the point whose location with respect to the sequence xt is
473 c to be determined.
473 c to be determined.
474 c
474 c
475 c****** o u t p u t ******
475 c****** o u t p u t ******
476 c left, mflag.....both integers, whose value is
476 c left, mflag.....both integers, whose value is
477 c
477 c
478 c 1 -1 if x .lt. xt(1)
478 c 1 -1 if x .lt. xt(1)
479 c i 0 if xt(i) .le. x .lt. xt(i+1)
479 c i 0 if xt(i) .le. x .lt. xt(i+1)
480 c i 0 if xt(i) .lt. x .eq. xt(i+1) .eq. xt(lxt)
480 c i 0 if xt(i) .lt. x .eq. xt(i+1) .eq. xt(lxt)
481 c i 1 if xt(i) .lt. xt(i+1) .eq. xt(lxt) .lt. x
481 c i 1 if xt(i) .lt. xt(i+1) .eq. xt(lxt) .lt. x
482 c
482 c
483 c In particular, mflag = 0 is the 'usual' case. mflag .ne. 0
483 c In particular, mflag = 0 is the 'usual' case. mflag .ne. 0
484 c indicates that x lies outside the CLOSED interval
484 c indicates that x lies outside the CLOSED interval
485 c xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the
485 c xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the
486 c intervals is due to the decision to make all pp functions cont-
486 c intervals is due to the decision to make all pp functions cont-
487 c inuous from the right, but, by returning mflag = 0 even if
487 c inuous from the right, but, by returning mflag = 0 even if
488 C x = xt(lxt), there is the option of having the computed pp function
488 C x = xt(lxt), there is the option of having the computed pp function
489 c continuous from the left at xt(lxt) .
489 c continuous from the left at xt(lxt) .
490 c
490 c
491 c****** m e t h o d ******
491 c****** m e t h o d ******
492 c The program is designed to be efficient in the common situation that
492 c The program is designed to be efficient in the common situation that
493 c it is called repeatedly, with x taken from an increasing or decrea-
493 c it is called repeatedly, with x taken from an increasing or decrea-
494 c sing sequence. This will happen, e.g., when a pp function is to be
494 c sing sequence. This will happen, e.g., when a pp function is to be
495 c graphed. The first guess for left is therefore taken to be the val-
495 c graphed. The first guess for left is therefore taken to be the val-
496 c ue returned at the previous call and stored in the l o c a l varia-
496 c ue returned at the previous call and stored in the l o c a l varia-
497 c ble ilo . A first check ascertains that ilo .lt. lxt (this is nec-
497 c ble ilo . A first check ascertains that ilo .lt. lxt (this is nec-
498 c essary since the present call may have nothing to do with the previ-
498 c essary since the present call may have nothing to do with the previ-
499 c ous call). Then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left =
499 c ous call). Then, if xt(ilo) .le. x .lt. xt(ilo+1), we set left =
500 c ilo and are done after just three comparisons.
500 c ilo and are done after just three comparisons.
501 c Otherwise, we repeatedly double the difference istep = ihi - ilo
501 c Otherwise, we repeatedly double the difference istep = ihi - ilo
502 c while also moving ilo and ihi in the direction of x , until
502 c while also moving ilo and ihi in the direction of x , until
503 c xt(ilo) .le. x .lt. xt(ihi) ,
503 c xt(ilo) .le. x .lt. xt(ihi) ,
504 c after which we use bisection to get, in addition, ilo+1 = ihi .
504 c after which we use bisection to get, in addition, ilo+1 = ihi .
505 c left = ilo is then returned.
505 c left = ilo is then returned.
506 c
506 c
507 integer left,lxt,mflag, ihi,ilo,istep,middle
507 integer left,lxt,mflag, ihi,ilo,istep,middle
508 real x,xt(lxt)
508 real x,xt(lxt)
509 data ilo /1/
509 data ilo /1/
510 save ilo
510 save ilo
511 ihi = ilo + 1
511 ihi = ilo + 1
512 if (ihi .lt. lxt) go to 20
512 if (ihi .lt. lxt) go to 20
513 if (x .ge. xt(lxt)) go to 110
513 if (x .ge. xt(lxt)) go to 110
514 if (lxt .le. 1) go to 90
514 if (lxt .le. 1) go to 90
515 ilo = lxt - 1
515 ilo = lxt - 1
516 ihi = lxt
516 ihi = lxt
517 c
517 c
518 20 if (x .ge. xt(ihi)) go to 40
518 20 if (x .ge. xt(ihi)) go to 40
519 if (x .ge. xt(ilo)) go to 100
519 if (x .ge. xt(ilo)) go to 100
520 c
520 c
521 c **** now x .lt. xt(ilo) . decrease ilo to capture x .
521 c **** now x .lt. xt(ilo) . decrease ilo to capture x .
522 istep = 1
522 istep = 1
523 31 ihi = ilo
523 31 ihi = ilo
524 ilo = ihi - istep
524 ilo = ihi - istep
525 if (ilo .le. 1) go to 35
525 if (ilo .le. 1) go to 35
526 if (x .ge. xt(ilo)) go to 50
526 if (x .ge. xt(ilo)) go to 50
527 istep = istep*2
527 istep = istep*2
528 go to 31
528 go to 31
529 35 ilo = 1
529 35 ilo = 1
530 if (x .lt. xt(1)) go to 90
530 if (x .lt. xt(1)) go to 90
531 go to 50
531 go to 50
532 c **** now x .ge. xt(ihi) . increase ihi to capture x .
532 c **** now x .ge. xt(ihi) . increase ihi to capture x .
533 40 istep = 1
533 40 istep = 1
534 41 ilo = ihi
534 41 ilo = ihi
535 ihi = ilo + istep
535 ihi = ilo + istep
536 if (ihi .ge. lxt) go to 45
536 if (ihi .ge. lxt) go to 45
537 if (x .lt. xt(ihi)) go to 50
537 if (x .lt. xt(ihi)) go to 50
538 istep = istep*2
538 istep = istep*2
539 go to 41
539 go to 41
540 45 if (x .ge. xt(lxt)) go to 110
540 45 if (x .ge. xt(lxt)) go to 110
541 ihi = lxt
541 ihi = lxt
542 c
542 c
543 c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval.
543 c **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval.
544 50 middle = (ilo + ihi)/2
544 50 middle = (ilo + ihi)/2
545 if (middle .eq. ilo) go to 100
545 if (middle .eq. ilo) go to 100
546 c note. it is assumed that middle = ilo in case ihi = ilo+1 .
546 c note. it is assumed that middle = ilo in case ihi = ilo+1 .
547 if (x .lt. xt(middle)) go to 53
547 if (x .lt. xt(middle)) go to 53
548 ilo = middle
548 ilo = middle
549 go to 50
549 go to 50
550 53 ihi = middle
550 53 ihi = middle
551 go to 50
551 go to 50
552 c**** set output and return.
552 c**** set output and return.
553 90 mflag = -1
553 90 mflag = -1
554 left = 1
554 left = 1
555 return
555 return
556 100 mflag = 0
556 100 mflag = 0
557 left = ilo
557 left = ilo
558 return
558 return
559 110 mflag = 1
559 110 mflag = 1
560 if (x .eq. xt(lxt)) mflag = 0
560 if (x .eq. xt(lxt)) mflag = 0
561 left = lxt
561 left = lxt
562 111 if (left .eq. 1) return
562 111 if (left .eq. 1) return
563 left = left - 1
563 left = left - 1
564 if (xt(left) .lt. xt(lxt)) return
564 if (xt(left) .lt. xt(lxt)) return
565 go to 111
565 go to 111
566 end
566 end
General Comments 0
You need to be logged in to leave comments. Login now