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