##// END OF EJS Templates
schain
RSS
test 5
rflores -
r1640:5807a75fd6c1
parent child
Show More
Side by Side Unified
Context file:
r1640:5807a75fd6c1 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,604 +1,621
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
488 491
489 492 if(imode.eq.1.or.imode.eq.3) then
490 493
491 494 vte=sqrt(bk*te/em)
492 495 thetae=(omega/ak)/(sqrt(2.0)*vte)
493 496 phi=(omegae/ak)/(sqrt(2.0)*vte)
494 497 psie=(ven/ak)/(sqrt(2.0)*vte)
495 498 ye=y_electron(thetae,phi,psie,alpha)
499 write(*,fmt='("AFTER YE")')
500 c call exit
496 501
497 502 else if(imode.eq.2) then
498 503 c
499 504 c use Mike Sulzer's library here: alpha2 is angle off perp in degrees
500 505 c
501 506
502 507 freq=omega/(2.0*pi)
503 508 alpha2=abs(pi/2.0-alpha)*180.0/pi
504 509 c write(*,*) "ye: ", ye
505 510 call collision(densmks, te, freq, alpha2, ye)
506 c write(*,fmt='(" geobfield: time is before earliest model.")')
511 c write(*,fmt='("AFTER COLLISION")')
512 c call exit
507 513 ye=ye*omega+cmplx(0.0,1.0)
508 514
509 515 end if
510 516 yed=ye+cmplx(0.0,dl)
511 517
512 518 p=(cabs(ye))**2*real(sum2)+cabs(sum1+cmplx(0.0,dl))**2*real(ye)
513 519 p=p/(cabs(yed+sum1))**2
514 520 spect1=p*2.0e0/(omega*pi)
515 write(*,*) "spect1 DONE "
521 write(*,*) "spect1:",spect1
516 522 return
517 523 end
518 524
519 subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1, nion1,
520 & alpha1, dens1, bfld1, acf)
525 subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1,
526 & alpha1, dens1, bfld1, acf, nion1)
521 527 c
522 528 c computes autocorrelation function for given plasma parameters
523 529 c by integrating real spectrum
524 530 c tau in sec., alpha in radians, density in cgs, bfield in cgs
525 531 c scattering wavelength (wl) in meters
526 532 c
527 533 include 'fitter.h'
528 534
529 535 real wl,tau,te1,ti1(nion1),fi1(nion1),ven1,vin1(nion1),alpha1,
530 536 & dens1,bfld1,acf
531 537 real pi
532 538 integer nion1
533 539 integer wi1(nion1)
534 integer i,j,k
540 integer i,j,k,imode
541 common /mode/imode
535 542 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")')
536 546 pi=4.0*atan(1.0)
537 547 c
538 548 c copy arguments to common block
539 549 c
540 550 ak=2.0*pi/wl
541 551 imode=2
552 write(*,*) "imode:",imode
542 553
543 554 nion=nion1
544 555 alpha=alpha1
545 556 te=te1
546 557 ven=ven1
558 c write(*,fmt='("INIT2")')
547 559 do i=1,nion
560 c write(*,fmt='("INIT2.5")')
548 561 ti(i)=ti1(i)
549 562 fi(i)=fi1(i)
550 563 vin(i)=vin1(i)
551 564 wi(i)=wi1(i)
552 565 end do
566 c write(*,fmt='("INIT3")')
553 567 dens=dens1
554 568 bfld=bfld1
555 569
556 570 c write(*,*) wl,alpha1,bfld1,dens
557 571 c call exit
558
572 c write(*,fmt='("Before Gauss")')
559 573 call gaussq(tau,acf)
560 574
575 write(*,*) "FINAL acf:",acf
576
577 c write(*,fmt='("After Gauss")')
561 578 return
562 579 end
563 580
564 581 subroutine gaussq(tau,acf)
565 582 c
566 583 c Computes cosine or sine transform of given real function
567 584 c Uses quadpack, tau in sec.
568 585 c
569 586 real a,abserr,epsabs,acf,tau,work,chebmo
570 587 integer ier,integr,iwork,last,leniw,lenw,limit,limlst,
571 588 * lst,maxp1,neval
572 589 integer knorm,ksave,momcom,nrmom
573 590 dimension iwork(300),work(1625),chebmo(61,25)
574 591 external spect1
575 592 include 'fitter.h'
576 593 c
577 594 c write(*,*) "inside gaussq"
578 595 a = 0.0
579 596 b = 2.5e4*ak ! upper integration limit open to debate (was 1.5, now 2.5)
580 597 integr = 1
581 598 epsabs = 1.0e-4
582 599 limlst = 50
583 600 limit = 100
584 601 leniw = limit*2+limlst
585 602 maxp1 = 61
586 603 lenw = leniw*2+maxp1*25
587 604
588 605 c write(*,*) leniw,lenw
589 606
590 607 nrmom=0
591 608 ksave=0
592 609 momcom=0
593
610 write(*,*) "Before qc25f:",acf,imode
594 611 c much faster, more robust
595 612 c write(*,*) "acf_in: ",acf
596 613 call qc25f(spect1,a,b,tau,integr,nrmom,maxp1,ksave,acf,
597 614 & abserr,neval,resabs,resasc,momcom,chebmo)
598 615 c write(*,*) "acf_out: ",acf
599 616 c call qawf(spect1,a,tau,integr,epsabs,acf,abserr,neval,
600 617 c & ier,limlst,lst,leniw,maxp1,lenw,iwork,work)
601
602
618 write(*,*) "After qc25f:",acf
619 c call exit
603 620 return
604 621 end
General Comments 0
You need to be logged in to leave comments. Login now