##// END OF EJS Templates
schain
RSS
fitacf.f update
rflores -
r1636:43055b8aca7f
parent child
Show More
Side by Side Unified
Context file:
r1636:43055b8aca7f Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,621 +1,604
1 1 subroutine guess(acf,tau,npts,zero,amin,te,tr)
2 2 c
3 3 c find zero crossing (zero), depth of minimum (amin), height of maximum
4 4 c
5 5 real acf(npts),tau(npts)
6 6
7 7 zero=0.0
8 8 amin=1.0
9 9 tmin=0.0
10 10 jmin=0
11 11
12 12 do i=npts,2,-1
13 13 if(acf(i)*acf(i-1).lt.0.0) then
14 14 zero=(tau(i-1)*acf(i)-tau(i)*acf(i-1))/(acf(i)-acf(i-1))
15 15 end if
16 16 if(acf(i).lt.amin) then
17 17 amin=acf(i)
18 18 jmin=i
19 19 end if
20 20 end do
21 21
22 22 if(jmin.gt.0) then
23 23 call parab1(tau(jmin-1),acf(jmin-1),a,b,c)
24 24 tmin=-b/(2.0*a)
25 25 amin=c+tmin*(b+tmin*a)
26 26 end if
27 27
28 28 tr=cdtr1(-amin)
29 29 te=czte1(zero*1000.0,tr)
30 30 return
31 31 end
32 32
33 33 subroutine parab1(x,y,a,b,c)
34 34 C-----
35 35 dimension x(3),y(3)
36 36 delta=x(1)-x(2)
37 37 a=(y(1)-2.*y(2)+y(3))/(2.*delta*delta)
38 38 b=(y(1)-y(2))/delta - a*(x(1)+x(2))
39 39 c=y(1)-a*x(1)*x(1)-b*x(1)
40 40 return
41 41 end
42 42
43 43 real function cdtr1(depth)
44 44 C-----convert depth to te/ti ratio
45 45 dimension tr(4)
46 46 C according to the curve published in farley et al 1967
47 47 c modified for 2004 conditions on axis
48 48 c data tr/7.31081,3.53286,5.92271,.174/
49 49 data tr/9.5,4.0,8.5,.3/
50 50 data nt/4/
51 51 cdtr1=tr(1)
52 52 do i=2,nt
53 53 cdtr1=cdtr1*depth + tr(i)
54 54 end do
55 55 return
56 56 end
57 57
58 58 real function czte1(zlag,tr)
59 59 C-----convert zero crossing point to te
60 60 C according to the curve published in farley et al 1967
61 61 c modified for 2004 conditions on axis
62 62 dimension dt(4)
63 63 c data dt/0.00945025,-0.0774338,.203626,.812397/,nd/4/
64 64 data dt/0.00945025,-0.0774338,.2,0.9/,nd/4/
65 65 data t0/1000./
66 66 tr1=min(abs(tr),5.)
67 67 if(zlag .eq. 0)then
68 68 czte1=1000000.
69 69 else
70 70 dt0=dt(1)
71 71 do i=2,nd
72 72 dt0=dt0*tr1 + dt(i)
73 73 end do
74 74 czte1=t0*(dt0/zlag)**2
75 75 end if
76 76 return
77 77 end
78 78
79 79 subroutine fit(wl,taup,rhop,covar,cinv,sigma2p,paramp,ebp,
80 80 & bfldp,alphap,densp,alt,time,nl,ifitp,ist)
81 81 c
82 82 c subroutine to fit measured ACF
83 83 c wavelength wl (m),lags taup (s), normalized acf rhop,
84 84 c experimental variances sigma2p, covariances covar,
85 85 c fit parameters params, error bars ebp, magnetic field bfldp (gauss),
86 86 c alphap B-field angle (radians), density densp (gcs),
87 87 c altitude alt (km), time (LT hours), nl lags
88 88 c ifitp determines which parameters are fit (see below)
89 89 c
90 90 real tol,pi,wl
91 91 parameter(nlmax=100,npmax=10,lwa=2000,tol=1.0e-5)
92 92 external fcn
93 93
94 94 real covar(nl,nl),cinv(nl,nl)
95 95 real ev(nlmax*nlmax),ap(nlmax*(nlmax+1)/2)
96 96 real fv1(nlmax),fv2(nlmax),det(2),w(nlmax)
97 97
98 98 real taup(nl),rhop(nl),sigma2p(nl),paramp(npmax)
99 99 real p2(npmax),ebp(npmax),alphap,bfldp,densp
100 100 real tau(nlmax),rho(nlmax),sigma2(nlmax),params(npmax)
101 101 real fvec(nlmax),wa(lwa),wa2(lwa),eb(npmax)
102 102
103 103 integer iwa(npmax),ifitp(npmax),ifit(npmax)
104 104
105 105 include 'fitter.h'
106 106 common /mode/imode
107 107
108 108 common/fitter/tau,rho,sigma2,params,ifit
109 109 common/trans/ev
110 110 C-----
111 111
112 112 c set ifit(1-5) to unity to fit:
113 113 c 1 normalization
114 114 c 2 Te
115 115 c 3 Ti
116 116 c 4 H+
117 117 c 5 He+
118 118
119 119 imode=2
120 120 c
121 121 pi=4.0*atan(1.0)
122 122 ak=2.0*pi/wl
123 123
124 124 wi(1)=16
125 125 wi(2)=1
126 126 wi(3)=4
127 127 nion=3
128 128
129 129 c
130 130 c invert covariances and find eigenvalues and eigenvectors
131 131 c
132 132 l = 0
133 133 do 20 j = 1, nl
134 134 do 10 i = 1, j
135 135 l = l + 1
136 136 ap(l) = covar(i,j)
137 137 10 continue
138 138 20 continue
139 139
140 140 call sppfa(ap,nl,info)
141 141 call sppdi(ap,nl,det,01)
142 142
143 143 l = 0
144 144 do 40 j = 1, nl
145 145 do 30 i = 1, j
146 146 l = l + 1
147 147 cinv(i,j)=ap(l)
148 148 cinv(j,i)=ap(l)
149 149 30 continue
150 150 40 continue
151 151
152 152 c
153 153 c transformation matrix is inverse (transpose) of eigenvectors
154 154 c
155 155
156 156 matz=1
157 157 call rs(nl,nl,cinv,w,matz,ev,fv1,fv2,ierr)
158 158
159 159 do i=1,nl
160 160 sigma2(i)=1.0/w(i)
161 161 end do
162 162
163 163 c
164 164 c extract desired parameters
165 165 c
166 166 iparm=0
167 167 do i=1,5
168 168 eb(i)=0.0
169 169 if(ifitp(i).eq.1) then
170 170 iparm=iparm+1
171 171 p2(iparm)=paramp(i)
172 172 end if
173 173 ifit(i)=ifitp(i)
174 174 params(i)=paramp(i)
175 175 end do
176 176 np=iparm
177 177
178 178 alpha=alphap
179 179 dens=densp
180 180 do i=1,nl
181 181 tau(i)=taup(i)
182 182 rho(i)=rhop(i)
183 183 sigma2p(i)=sigma2(i)
184 184 end do
185 185
186 186 bfld=bfldp
187 187
188 188 c
189 189 c no. equations is no. lags - do nlls fit
190 190 c
191 191
192 192 call lmdif1(fcn,nl,np,p2,fvec,tol,info,iwa,wa,lwa)
193 193 ist=info
194 194
195 195 c
196 196 c generate error bars here
197 197 c
198 198 call fdjac2(fcn,nl,np,p2,fvec,wa,nl,iflag,0.0e0,wa2)
199 199
200 200 do i=1,np
201 201 err=0.0
202 202 do j=1,nl
203 203 err=err+(wa(j+(i-1)*nl))**2
204 204 end do
205 205 if(err.gt.0.0) eb(i)=sqrt(err**-1)
206 206 end do
207 207
208 208 c reorder results
209 209 iparm=0
210 210 do i=1,5
211 211 if(ifit(i).eq.1) then
212 212 iparm=iparm+1
213 213 paramp(i)=p2(iparm)
214 214 ebp(i)=eb(iparm)
215 215 else
216 216 ebp(i)=0.0
217 217 paramp(i)=params(i)
218 218 end if
219 219 end do
220 220
221 221 9 continue
222 222 c write(*,*) dens
223 223 c write(*,*) te
224 224
225 225 return
226 226 end
227 227
228 228
229 229 subroutine fcn(m,n,x,fvec,iflag)
230 230 c
231 231 c provides m functions (fvec) in n variables (x) to minimize
232 232 c in least-squares sense
233 233 c
234 234
235 235 parameter(nlmax=100,npmax=10)
236 236 integer m,n,iflag
237 237 real fvec(m)
238 238
239 239 integer ifit(npmax)
240 240 real x(n),fv(nlmax),anorm,params(npmax),ev(nlmax*nlmax)
241 241 real tau(nlmax),rho(nlmax),sigma2(nlmax)
242 242 real chisq,acf
243 243
244 244 include 'fitter.h'
245 245
246 246 common/fitter/tau,rho,sigma2,params,ifit
247 247 common /errs/ chisq
248 248 common /trans/ev
249 249
250 250 c 1 0 0 0 0 fit zero lag (otherwise set to default)
251 251 c 0 1 0 0 0 fit Te (otherwise default)
252 252 c 0 0 1 0 0 fit Ti (otherwise set to Te)
253 253 c 0 0 0 1 0 fit H+ (otherwize default)
254 254 c 0 0 0 0 1 fit He+ (otherwise default)
255 255
256 256 8 continue
257 257
258 258 c
259 259 c ignore collisions
260 260 c
261 261 c write(*,*) "starting fcn"
262 262 ven=0.0
263 263 do i=1,3
264 264 vin(i)=0.0
265 265 end do
266 266
267 267 iparm=0
268 268
269 269 c
270 270 c zero lag normalization constant
271 271 c
272 272 if(ifit(1).eq.1) then
273 273 iparm=iparm+1
274 274 c=x(iparm)
275 275 else
276 276 c=params(1)
277 277 end if
278 278 c
279 279 c Te
280 280 c
281 281 if(ifit(2).eq.1) then
282 282 iparm=iparm+1
283 283 te=x(iparm)
284 284 else
285 285 te=params(2)
286 286 end if
287 287 c
288 288 c Ti - default is Te rather than initial value
289 289 c
290 290 if(ifit(3).eq.1) then
291 291 iparm=iparm+1
292 292 do i=1,3
293 293 ti(i)=x(iparm)
294 294 end do
295 295 else
296 296 do i=1,3
297 297 ti(i)=te
298 298 end do
299 299 params(3)=te
300 300 end if
301 301
302 302 c three-ion plasma
303 303
304 304 nion=3
305 305 c
306 306 c composition H+ first
307 307 c
308 308 fi(1)=1.0
309 309 if(ifit(4).eq.1) then
310 310 iparm=iparm+1
311 311 fi(2)=(x(iparm))
312 312 fi(1)=fi(1)-(x(iparm))
313 313 else
314 314 fi(2)=params(4)
315 315 fi(1)=fi(1)-params(4)
316 316 end if
317 317 c
318 318 c He+
319 319 c
320 320 if(ifit(5).eq.1) then
321 321 iparm=iparm+1
322 322 fi(3)=(x(iparm))
323 323 fi(1)=fi(1)-(x(iparm))
324 324 else
325 325 fi(3)=params(5)
326 326 fi(1)=fi(1)-params(5)
327 327 end if
328 328
329 329 c write(*,*) x,bfld,alpha,dens,ak,m
330 330
331 331 call gaussq(0.0,anorm)
332 332 anorm=anorm/c
333 333
334 334 if(anorm.eq.0.0.or.m.eq.0)return
335 335 chisq=0.0
336 336 do i=1,m
337 337 call gaussq(tau(i),acf)
338 338 fv(i)=(acf/anorm-rho(i))
339 339 end do
340 340
341 341 c
342 342 c transform to space where s2 is diagonal
343 343 c (are i and j transposed below?)
344 344 c
345 345 do i=1,m
346 346 fvec(i)=0.0
347 347 do j=1,m
348 348 fvec(i)=fvec(i)+fv(j)*ev(j+(i-1)*m)/sqrt(sigma2(i))
349 349 end do
350 350 chisq=chisq+fvec(i)**2
351 351 end do
352 352
353 353 chisq=chisq/float(m)
354 354
355 355 c write(*,*) fvec
356 356 c write(*,*) chisq
357 357 c stop
358 358
359 359 return
360 360 end
361 361
362 362 c
363 363 c not really fond of above code
364 364 c
365 365
366 366
367 367 complex function cj_ion(theta,psi)
368 368 c
369 369 real theta,phi,psi,alpha
370 370 complex z
371 371 complex*16 zz
372 372 z=zz(dcmplx(-theta,psi))
373 373 cj_ion=z*cmplx(0.0,-1.0)
374 374 return
375 375 end
376 376
377 377 complex function cj_electron(theta,phi,psi,alpha)
378 378 c
379 379 c Theta, phi, and psi are the normalized frequency, gyrofrequency,
380 380 c and collision frequency. Alpha is angle between wavevector and
381 381 c magnetic field in radians.
382 382 c
383 383 parameter(nterms=10)
384 384 real theta,phi,psi,alpha
385 385 real arg
386 386 complex cj,cy(0:nterms)
387 387 complex z
388 388 complex*16 zz
389 389 integer m,nz,ierr
390 390 integer imode
391 391 common/mode/imode
392 392 c
393 393 if(imode.eq.3) then
394 394
395 395 arg=0.5*(sin(alpha)/phi)**2
396 396
397 397 c calculate modified Bessel functions using amos library
398 398
399 399 call cbesi(cmplx(arg,0.0),0.0,2,nterms+1,cy,nz,ierr)
400 400
401 401 cj=cmplx(0.0,0.0)
402 402
403 403 do m=-nterms,nterms
404 404 z=zz(dcmplx(-(theta-float(m)*phi)/cos(alpha),
405 405 & psi/cos(alpha)))
406 406 cj=cj+z*cy(iabs(m))
407 407 end do
408 408
409 409 cj_electron=cj*cmplx(0.0,-1.0/cos(alpha))
410 410
411 411 else
412 412
413 413 z=zz(dcmplx(-theta/cos(alpha),psi/cos(alpha)))
414 414 cj_electron=z*cmplx(0.0,-1.0/cos(alpha))
415 415 cj_electron=cj_electron*(1.0-(sin(alpha)/phi)**2/2.0)
416 416
417 417 end if
418 418
419 419 return
420 420 end
421 421
422 422 complex function y_ion(theta,psi)
423 423 c
424 424 real theta,phi,psi,alpha
425 425 complex cj, cj_ion
426 426 cj=cj_ion(theta,psi)
427 427 y_ion=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)
428 428 y_ion=y_ion/(1-psi*cj)
429 429 return
430 430 end
431 431
432 432 complex function y_electron(theta,phi,psi,alpha)
433 433 c
434 434 real theta,phi,psi,alpha
435 435 complex cj, cj_electron
436 436 cj=cj_electron(theta,phi,psi,alpha)
437 437 y_electron=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)
438 438 y_electron=y_electron/(1.0-psi*cj)
439 439 return
440 440 end
441 441
442 442 real function spect1(omega)
443 443 c
444 444 c Function to generate IS ion-line spectrum
445 445 c Provisions for nimax ions
446 446 c No provisions for drifts
447 447 c
448 448 c imode=1 Farley B-field treatment for electrons
449 449 c 2 use Mike Sulzer's model for ye
450 450 c 3 calculate ye using sums of Bessel functions
451 451 c
452 452 include 'fitter.h'
453 453
454 454 real omega,thetae,thetai(nimax),psie,psii(nimax),phi,p
455 455 real tr(nimax),vti(nimax),vte,omegae
456 456 real bk,em,e,dlf,dl,pi
457 457 real alpha2,densmks,freq
458 458 complex ye,yed,yet,yi(nimax),sum1,sum2,sumdl
459 459 complex y_ion, y_electron, y_esum
460 460 integer i,j,k,imode
461 461 common/mode/ imode
462 462 data bk,em,e,dlf/1.38e-23,9.1e-31,1.6e-19,4772.9/
463 463
464 464 c fi is ion fraction, wi is ion atomic weight
465 465 c dens is electron density (cgs), alpha is angle
466 466 c between k and the magnetic field (radians)
467 467
468 468 pi=4.0*atan(1.0)
469 469
470 470 if(omega.eq.0.0) omega=1.0
471 471
472 472 omegae=e*bfld*1.0e-4/em
473 473 densmks=dens*1.0e6
474 474
475 475 sum1=0.0
476 476 sum2=0.0
477 477
478 478 do i=1,nion
479 479 tr(i)=te/ti(i)
480 480 vti(i)=sqrt(bk*ti(i)/(1.67e-27*wi(i)))
481 481 thetai(i)=(omega/ak)/(sqrt(2.0)*vti(i))
482 482 psii(i)=(vin(i)/ak)/(sqrt(2.0)*vti(i))
483 483 yi(i)=y_ion(thetai(i),psii(i))
484 484 sum1=sum1+fi(i)*tr(i)*yi(i)
485 485 sum2=sum2+fi(i)*yi(i)
486 486 end do
487 487 dl=ak**2*dlf*te/densmks
488 c write(*,fmt='("Before YE")')
489 c write(*,*) imode
490 c call exit
491 488
492 489 if(imode.eq.1.or.imode.eq.3) then
493 490
494 491 vte=sqrt(bk*te/em)
495 492 thetae=(omega/ak)/(sqrt(2.0)*vte)
496 493 phi=(omegae/ak)/(sqrt(2.0)*vte)
497 494 psie=(ven/ak)/(sqrt(2.0)*vte)
498 495 ye=y_electron(thetae,phi,psie,alpha)
499 write(*,fmt='("AFTER YE")')
500 c call exit
501 496
502 497 else if(imode.eq.2) then
503 498 c
504 499 c use Mike Sulzer's library here: alpha2 is angle off perp in degrees
505 500 c
506 501
507 502 freq=omega/(2.0*pi)
508 503 alpha2=abs(pi/2.0-alpha)*180.0/pi
509 504 c write(*,*) "ye: ", ye
510 505 call collision(densmks, te, freq, alpha2, ye)
511 c write(*,fmt='("AFTER COLLISION")')
512 c call exit
506 c write(*,fmt='(" geobfield: time is before earliest model.")')
513 507 ye=ye*omega+cmplx(0.0,1.0)
514 508
515 509 end if
516 510 yed=ye+cmplx(0.0,dl)
517 511
518 512 p=(cabs(ye))**2*real(sum2)+cabs(sum1+cmplx(0.0,dl))**2*real(ye)
519 513 p=p/(cabs(yed+sum1))**2
520 514 spect1=p*2.0e0/(omega*pi)
521 write(*,*) "spect1:",spect1
515
522 516 return
523 517 end
524 518
525 subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1,
526 & alpha1, dens1, bfld1, acf, nion1)
519 subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1, nion1,
520 & alpha1, dens1, bfld1, acf)
527 521 c
528 522 c computes autocorrelation function for given plasma parameters
529 523 c by integrating real spectrum
530 524 c tau in sec., alpha in radians, density in cgs, bfield in cgs
531 525 c scattering wavelength (wl) in meters
532 526 c
533 527 include 'fitter.h'
534 528
535 529 real wl,tau,te1,ti1(nion1),fi1(nion1),ven1,vin1(nion1),alpha1,
536 530 & dens1,bfld1,acf
537 531 real pi
538 532 integer nion1
539 533 integer wi1(nion1)
540 integer i,j,k,imode
541 common /mode/imode
534 integer i,j,k
542 535 c
543 c write(*,*) "INITIAL acf:",wl,tau,te1,ti1,fi1,ven1,vin1,wi1,alpha1
544 write(*,*) "INITIAL acf:",dens1, bfld1, acf, nion1
545 c write(*,fmt='("INIT")')
546 536 pi=4.0*atan(1.0)
547 537 c
548 538 c copy arguments to common block
549 539 c
550 540 ak=2.0*pi/wl
551 541 imode=2
552 write(*,*) "imode:",imode
553 542
554 543 nion=nion1
555 544 alpha=alpha1
556 545 te=te1
557 546 ven=ven1
558 c write(*,fmt='("INIT2")')
559 547 do i=1,nion
560 c write(*,fmt='("INIT2.5")')
561 548 ti(i)=ti1(i)
562 549 fi(i)=fi1(i)
563 550 vin(i)=vin1(i)
564 551 wi(i)=wi1(i)
565 552 end do
566 c write(*,fmt='("INIT3")')
567 553 dens=dens1
568 554 bfld=bfld1
569 555
570 556 c write(*,*) wl,alpha1,bfld1,dens
571 557 c call exit
572 c write(*,fmt='("Before Gauss")')
573 call gaussq(tau,acf)
574 558
575 write(*,*) "FINAL acf:",acf
559 call gaussq(tau,acf)
576 560
577 c write(*,fmt='("After Gauss")')
578 561 return
579 562 end
580 563
581 564 subroutine gaussq(tau,acf)
582 565 c
583 566 c Computes cosine or sine transform of given real function
584 567 c Uses quadpack, tau in sec.
585 568 c
586 569 real a,abserr,epsabs,acf,tau,work,chebmo
587 570 integer ier,integr,iwork,last,leniw,lenw,limit,limlst,
588 571 * lst,maxp1,neval
589 572 integer knorm,ksave,momcom,nrmom
590 573 dimension iwork(300),work(1625),chebmo(61,25)
591 574 external spect1
592 575 include 'fitter.h'
593 576 c
594 577 c write(*,*) "inside gaussq"
595 578 a = 0.0
596 579 b = 2.5e4*ak ! upper integration limit open to debate (was 1.5, now 2.5)
597 580 integr = 1
598 581 epsabs = 1.0e-4
599 582 limlst = 50
600 583 limit = 100
601 584 leniw = limit*2+limlst
602 585 maxp1 = 61
603 586 lenw = leniw*2+maxp1*25
604 587
605 588 c write(*,*) leniw,lenw
606 589
607 590 nrmom=0
608 591 ksave=0
609 592 momcom=0
610 write(*,*) "Before qc25f:",acf,imode
593
611 594 c much faster, more robust
612 595 c write(*,*) "acf_in: ",acf
613 596 call qc25f(spect1,a,b,tau,integr,nrmom,maxp1,ksave,acf,
614 597 & abserr,neval,resabs,resasc,momcom,chebmo)
615 598 c write(*,*) "acf_out: ",acf
616 599 c call qawf(spect1,a,tau,integr,epsabs,acf,abserr,neval,
617 600 c & ier,limlst,lst,leniw,maxp1,lenw,iwork,work)
618 write(*,*) "After qc25f:",acf
619 c call exit
601
602
620 603 return
621 604 end
General Comments 0
You need to be logged in to leave comments. Login now