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