| @@ -485,9 +485,6 c between k and the magnetic field (radians) | |||||
| 485 | sum2=sum2+fi(i)*yi(i) |
|
485 | sum2=sum2+fi(i)*yi(i) | |
| 486 | end do |
|
486 | end do | |
| 487 | dl=ak**2*dlf*te/densmks |
|
487 | dl=ak**2*dlf*te/densmks | |
| 488 | c write(*,fmt='("Before YE")') |
|
|||
| 489 | c write(*,*) imode |
|
|||
| 490 | c call exit |
|
|||
| 491 |
|
488 | |||
| 492 | if(imode.eq.1.or.imode.eq.3) then |
|
489 | if(imode.eq.1.or.imode.eq.3) then | |
| 493 |
|
490 | |||
| @@ -496,8 +493,6 c call exit | |||||
| 496 | phi=(omegae/ak)/(sqrt(2.0)*vte) |
|
493 | phi=(omegae/ak)/(sqrt(2.0)*vte) | |
| 497 | psie=(ven/ak)/(sqrt(2.0)*vte) |
|
494 | psie=(ven/ak)/(sqrt(2.0)*vte) | |
| 498 | ye=y_electron(thetae,phi,psie,alpha) |
|
495 | ye=y_electron(thetae,phi,psie,alpha) | |
| 499 | write(*,fmt='("AFTER YE")') |
|
|||
| 500 | c call exit |
|
|||
| 501 |
|
496 | |||
| 502 | else if(imode.eq.2) then |
|
497 | else if(imode.eq.2) then | |
| 503 | c |
|
498 | c | |
| @@ -508,8 +503,7 c | |||||
| 508 | alpha2=abs(pi/2.0-alpha)*180.0/pi |
|
503 | alpha2=abs(pi/2.0-alpha)*180.0/pi | |
| 509 | c write(*,*) "ye: ", ye |
|
504 | c write(*,*) "ye: ", ye | |
| 510 | call collision(densmks, te, freq, alpha2, ye) |
|
505 | call collision(densmks, te, freq, alpha2, ye) | |
| 511 | c write(*,fmt='("AFTER COLLISION")') |
|
506 | c write(*,fmt='(" geobfield: time is before earliest model.")') | |
| 512 | c call exit |
|
|||
| 513 | ye=ye*omega+cmplx(0.0,1.0) |
|
507 | ye=ye*omega+cmplx(0.0,1.0) | |
| 514 |
|
508 | |||
| 515 | end if |
|
509 | end if | |
| @@ -518,12 +512,12 c call exit | |||||
| 518 | p=(cabs(ye))**2*real(sum2)+cabs(sum1+cmplx(0.0,dl))**2*real(ye) |
|
512 | p=(cabs(ye))**2*real(sum2)+cabs(sum1+cmplx(0.0,dl))**2*real(ye) | |
| 519 | p=p/(cabs(yed+sum1))**2 |
|
513 | p=p/(cabs(yed+sum1))**2 | |
| 520 | spect1=p*2.0e0/(omega*pi) |
|
514 | spect1=p*2.0e0/(omega*pi) | |
| 521 | write(*,*) "spect1:",spect1 |
|
515 | ||
| 522 | return |
|
516 | return | |
| 523 | end |
|
517 | end | |
| 524 |
|
518 | |||
| 525 | subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1, |
|
519 | subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1, nion1, | |
| 526 |
& alpha1, dens1, bfld1, acf |
|
520 | & alpha1, dens1, bfld1, acf) | |
| 527 | c |
|
521 | c | |
| 528 | c computes autocorrelation function for given plasma parameters |
|
522 | c computes autocorrelation function for given plasma parameters | |
| 529 | c by integrating real spectrum |
|
523 | c by integrating real spectrum | |
| @@ -537,44 +531,33 c | |||||
| 537 | real pi |
|
531 | real pi | |
| 538 | integer nion1 |
|
532 | integer nion1 | |
| 539 | integer wi1(nion1) |
|
533 | integer wi1(nion1) | |
| 540 |
integer i,j,k |
|
534 | integer i,j,k | |
| 541 | common /mode/imode |
|
|||
| 542 | c |
|
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 | pi=4.0*atan(1.0) |
|
536 | pi=4.0*atan(1.0) | |
| 547 | c |
|
537 | c | |
| 548 | c copy arguments to common block |
|
538 | c copy arguments to common block | |
| 549 | c |
|
539 | c | |
| 550 | ak=2.0*pi/wl |
|
540 | ak=2.0*pi/wl | |
| 551 | imode=2 |
|
541 | imode=2 | |
| 552 | write(*,*) "imode:",imode |
|
|||
| 553 |
|
542 | |||
| 554 | nion=nion1 |
|
543 | nion=nion1 | |
| 555 | alpha=alpha1 |
|
544 | alpha=alpha1 | |
| 556 | te=te1 |
|
545 | te=te1 | |
| 557 | ven=ven1 |
|
546 | ven=ven1 | |
| 558 | c write(*,fmt='("INIT2")') |
|
|||
| 559 | do i=1,nion |
|
547 | do i=1,nion | |
| 560 | c write(*,fmt='("INIT2.5")') |
|
|||
| 561 | ti(i)=ti1(i) |
|
548 | ti(i)=ti1(i) | |
| 562 | fi(i)=fi1(i) |
|
549 | fi(i)=fi1(i) | |
| 563 | vin(i)=vin1(i) |
|
550 | vin(i)=vin1(i) | |
| 564 | wi(i)=wi1(i) |
|
551 | wi(i)=wi1(i) | |
| 565 | end do |
|
552 | end do | |
| 566 | c write(*,fmt='("INIT3")') |
|
|||
| 567 | dens=dens1 |
|
553 | dens=dens1 | |
| 568 | bfld=bfld1 |
|
554 | bfld=bfld1 | |
| 569 |
|
555 | |||
| 570 | c write(*,*) wl,alpha1,bfld1,dens |
|
556 | c write(*,*) wl,alpha1,bfld1,dens | |
| 571 | c call exit |
|
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 | return |
|
561 | return | |
| 579 | end |
|
562 | end | |
| 580 |
|
563 | |||
| @@ -607,7 +590,7 c write(*,*) leniw,lenw | |||||
| 607 | nrmom=0 |
|
590 | nrmom=0 | |
| 608 | ksave=0 |
|
591 | ksave=0 | |
| 609 | momcom=0 |
|
592 | momcom=0 | |
| 610 | write(*,*) "Before qc25f:",acf,imode |
|
593 | ||
| 611 | c much faster, more robust |
|
594 | c much faster, more robust | |
| 612 | c write(*,*) "acf_in: ",acf |
|
595 | c write(*,*) "acf_in: ",acf | |
| 613 | call qc25f(spect1,a,b,tau,integr,nrmom,maxp1,ksave,acf, |
|
596 | call qc25f(spect1,a,b,tau,integr,nrmom,maxp1,ksave,acf, | |
| @@ -615,7 +598,7 c write(*,*) "acf_in: ",acf | |||||
| 615 | c write(*,*) "acf_out: ",acf |
|
598 | c write(*,*) "acf_out: ",acf | |
| 616 | c call qawf(spect1,a,tau,integr,epsabs,acf,abserr,neval, |
|
599 | c call qawf(spect1,a,tau,integr,epsabs,acf,abserr,neval, | |
| 617 | c & ier,limlst,lst,leniw,maxp1,lenw,iwork,work) |
|
600 | c & ier,limlst,lst,leniw,maxp1,lenw,iwork,work) | |
| 618 | write(*,*) "After qc25f:",acf |
|
601 | ||
| 619 | c call exit |
|
602 | ||
| 620 | return |
|
603 | return | |
| 621 | end |
|
604 | end | |
General Comments 0
You need to be logged in to leave comments.
Login now