-
Notifications
You must be signed in to change notification settings - Fork 14
/
curba_GA.f
2179 lines (1897 loc) · 54.8 KB
/
curba_GA.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! YuP[2020-08-20] Renamed all b into b_GA, to avoid conflict with function b()
subroutine TorGA_curba(rjpd, rjpd0, ratjpd, denom, aspct, enpar,
&tc, thtc, theta, elomom, lh, zeff, model, tol, n0, ig)
cProlog
implicit none
c ----------------------------------------------------------------------
c CURBA version 7.2 12/11/97, edited by Y.R. Lin-Liu
c CURBA version 7.1, 11/21/90, R. H. Cohen, LLNL.
c calculates J/Pd for ECRH and lower hybrid
c from numerical solution for Spitzer-Harm
c problem, including trapped particles, electron-ion collisions,
c and poloidal variation of collision operator. J is flux-
c surface-averaged current density. Reference:
c R.H. Cohen, Effect of Trapped Electrons on Current Drive,
c Phys Fl. 30, 2442 (1987) and erratum in press.
c Usage notes:
c J/Pd normalized to e/(m nu vt), where nu =
c 4 pi n e**4 log lambda / m**2 vt**3, and vt = sqrt (2T/m). Note
c this normalization differs from that in Reference; J/Pd in
c reference is factor of 2 larger.
c When thtc .ne. 1, normalization uses thermal speed and collision
c frequency for bulk electrons.
c When thtc < 0, calculates maximum efficiency for rising buckets.
c Also provides values of Green's function g (Spitzer-Harm
c distribution function including trapped-particle effects)
c --see notes for versions 2.1 and 4.0. Define angular and radial
c parts H and F by:
c g=<rt/R> EXP (-energy/T) c**4 F H/(4 nu vt**3 zhat)
c where rt is mirror ratio from field minimum,
c R is major radius, and zhat=(1+charge)/2.
c For numerical reasons H is normalized differently than in reference;
c it is larger by zhat.
c User may wish to alter function FPP which calculates energy-dependent
c part of Green's function and its derivitive. See defining comments
c in routine.
c OUTPUTS:
c rjpd is current drive efficiency calculated with trapped particles
c rjpd0 is non-relativistic current-drive efficiency in limit of no
c trapped particles, at vparallel=min on relativistic resonance curve.
c ratjpd is the ratio of the above.
c denom is normalized absorbed power (denominator in rjpd).
c INPUTS:
c enpar is k_parallel*c/wave freq. Note enpar**2 < 1-elomom**2
c implies no resonant particles; rjpd set to zero and rjpd0
c determined for vparallel given by nonrelativistic
c resonance condition.
c tc is bulk electron temperature;
c thtc is ratio of temperature along characteristic to that of bulk,
c except for thtc < 0, -thtc is energy (in keV) of bucket rise;
c aspct is inverse aspct ratio;
c theta is poloidal angle at which power is absorbed (measured
c from outside); for thtc<0, theta is poloidal angle
c at which electrons become trapped in bucket; used to
c calculate resonant energy given elomom and enpar. Note
c if calling program fixes minimum resonant energy and
c calculates enpar, then theta has no significance to the
c physics; rjpd depends only on Bmin*Y/B(theta), not on
c theta or Y alone. But it is useful to be able to specify
c theta and Y separately in order to keep track of what higher
c harmonics are doing.
c ABS (lh) is power of e_perp in diffusion coeff; sign governs
c power of p_parallel in diffusion coef: + gives p_par^0;
c - gives p_par^2. For ECRH, |lh| is
c harmonic number (lh=1 for fundamental); + for E_-
c contribution, - for E_parallel; + with lh-->lh+2 for E_+
c component of electric field E. For now, there is no provision
c for the p_par^2 option with lh=0, as the compiler doesn't know
c the difference betweeen +0 and -0. if there is any interest,
c send a message to user 313 and a revision including this option
c will be created.
c elomom is ABS (lh)*cyclotron freq over wave freq for Y<0,
c interpret as evaluated at poloidal angle where electrons are first
c trapped in bucket.
c tol is absolute tolerence for DO2GBF integrator; starting in
c version 1.1 the variables integrated by D02GBE are normalized
c to be O(1), so tol can be viewed as relative error tolerence.
c zeff is ion effective charge.
c model: Absolute value of model selects
c collisional model: 1 for full bounce av, 2 for square well
c (numerical solution); 3 for analytic solution to square well.
c negative model does parallel heating (lower hybrid, fast wave,..).
c n0 is number of mesh points for gaussian integration of
c Green's function; n0=2,4,6,8,10,12,16,20,24,32,48 or 64. n0=128
c gives 64 points over resonance if resonance has only single
c passing-particle segment; 64 points on each if two segments.
c 20,24,32,48, or 64.
c ig selects form for relativistic correction to green's function.
c ig=1 uses model which agrees with Fisch thru first relativistic
c correction to nonrelativistic limit and in relativistic limit;
c ig=2 and 3 are older models; see comments in subroutine fupr.
c version 1.1, 7/17/86, normalizes integrator variables better.
c version 2.0, 7/23/86, adds lower hybrid
c version 2.1 9/15/86, attempted to allow access to the h array on a
c specified mesh of values of eta = mu*B0/energy. Here h(1,) is
c angular part H of generalized Spitzer-Harm function;
c h(2,)=b*dh(1,)/d eta, where b = <sqrt (1-mirror
c ratio * eta)>, and h(3,) was running integral of normalized current.
c This never worked properly; there seems to be no way to force the
c NAG routine d02gbE to maintain a fixed mesh while maintaining
c accuracy (except for the case where the number is the maximum number
c set in the dimension statements). However, as of version 4.0
c the interpolated h is directly available to users.
c version 3.0, 10/12/86, allows different temperature for bulk electrons
c and along resonance curve.
c version 4.0, 10/13/86, integrates using gaussian quadrature in
c EXP (-eperp) to get current, even for model=1,2. As a consequence,
c h is automatically available from calling program using h=hfun(eta,
c hpr), where hpr=d h/d eta (output variable). Also, h matrix is not
c recomputed unless needed.
c To calculate h's only, first call curba from your source program
c with your values for aspct, zeff, model, and tol (a good
c value is usually .002); all other input variables are irrelevant,
c except tc should be nonzero. then h is available as h=hfun(eta,hpr).
c version 5.0, 10/15/86, is relativistic; enpar replaces epart
c as input variable, and tc added.
c version 5.1, 11/11/86, fixes bug so rjpd=0 when no resonant
c particles, and adds denom to calling sequence.
c version 5.2, 11/19/86, provides an improved relativistic correction
c to Green's function and adds switch ig to select Green's function
c models.
c version 5.3, 11/21/86, adds a further-improved relativistic correction
c to Green's function. ig=1 selects new model; 2 and 3 for older models.
c version 5.4, 12/2/86, improves integration along resonance,
c eliminating trapped-particle region from domain.
c version 5.5, 2/16/87, eliminates bug in upper limit for power
c integration
c Version 6.0, 4/2/87, adds bucket-lift current-drive calculation,
c for particles trapped at poloidal angle theta with specified elomom,
c at zero perpendicular energy, and ending on the same characteristic
c with energy increased by -thtc (keV).
c Version 6.1, 12/6/87, corrects comments on version 2.1 and adds
c capability of using model=3 (analytic Legendre function form
c of Green's function) in accessing h via hfun from calling program.
c Also removes possible crash in legfn1 when z**2=1 by substituting
c limiting form for derivative of Legendre function.
c Version 7.0, 3/7/88, adds capability of higher-order flr to
c parallel heating calculation. lh must now be set; diffusion
c coefficient varies as eperp**|lh|. Sign of lh determines
c whether to weight diffusion coefficient with p_parallel^2;
c negative lh adds the weighting. Note that physically we
c expect weighting with v_parallel^2, which is constant on
c a Landau type resonance and hence irrelevant. Hence one
c will almost certainly want to set lh .ge. 0.
c Version 7.1, 11/21/90, corrects bug in ffp, for ig=1. In
c expression for fupr, earlier version had -.5/bfac instead
c of -.5*agrel/afac. Also puts in sign of bucket lift current,
c not tracked previously.
c------------------------------------------------------------------------
c
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
& aspect, sgnj, enpar, charge, zeff, theta, rabs,
&rmax, etam, etamh, et0, thtcu, thtc, theth, tc, aspct, costh, rmin
&, elom, elomom, expeu, enpar2, epst1, epst2, etmax, rjpd, rjpd0,
&ratjps, denom, expel, expe1, expe2, rint, qgauleg_GA, expeld,zhat,
&erise, etal, sgnj0, gammin, etau, dfh, h3fac, sigma, eparnr,
&fjpd0_GA
&, fjpd0h_GA, tol, epar, rl, asp1, ap0, alfsq, elooro, r3, eparsq,
&emin, gammax, ratjpd, ra0, elom2, agrel, bgrel, qgrel, tor
integer kg, ig, model, iheat, modela, l, modl2, ifail, n, n0, ixo,
& lh, lm1, lp1
c
common /TorGA_fcncom/ aspect,epar,rl,rmin,rmax,rabs,etam,charge,
&asp1,zhat,ap0,alfsq,elooro,r3, eparsq,l,iheat,lm1,lp1,modl2
common /TorGA_rela/ enpar2,elom,theth,emin,ra0,elom2,tor, gammin,
&agrel,bgrel,qgrel,ixo,kg
c
cSAP080617
real*8 fcn_GA,fcn
external fcn_GA
external fcd
c
c
cSAP0890617
c write(*,*)'TorGA_curba aspct,enpar',
c & aspct,enpar
c write(*,*)'tc, thtc, theta, elomom',tc, thtc, theta, elomom
c write(*,*)'lh,zeff,model,tol,n0,ig',lh,zeff,model,tol,n0,ig
kg = ig-2
aspect = aspct
sgnj = 1.0d0
if (enpar .lt. 0) sgnj = -1.0d0
charge = zeff
costh = COS (theta)
rabs = 1.0d0+aspect*costh
rmax = 1.0d0+aspect
rmin = 1.0d0-aspect
etam = rmin/rmax
etamh = etam
et0 = 0.0d0
thtcu = thtc
cSAP080617
c write(*,*)'kg,aspect,sgnj,charge',kg,aspect,sgnj,charge
c write(*,*)'costh,rabs,rmax,rmin',costh,rabs,rmax,rmin
c write(*,*)'etam,etamh,et0,thtcu',etam,etamh,et0,thtcu
c
c use th = tc if doing buckets
c
if (thtc .le. 0.0d0) thtcu = 1.0d0
theth = tc*thtcu/511.0d0
cSAP080617
c write(*,*)' theth', theth
if (model .lt. 0) then
iheat = -1
elom = 0
modela = -3
else
iheat = 1
elom = elomom
modela = 3
endif
ixo = 0
l = IABS (lh)
!write(*,*)' TorGA_curba: lh=',lh
c
c set flag to get ppar**(2*ixo) factor in diffusion coef.;
c for ecrh, ixo=0 for E-, E+ and =1 for E_par contributions
c
if (lh .lt. 0) ixo = 1
modl2 = modela-2
call TorGA_setpar
c
c do current-drive integral:
c
ifail = 0
c
c --- set limits of integration for expe = EXP (-epst), where epst=e-emin,
c normalized to th.
c
expeu = 1.0d0
enpar2 = enpar*enpar
call TorGA_limits (epst1, epst2, etmax)
cSAP080617
! write(*,*)'after TorGA_limits epst1, epst2, etmax',
! & epst1, epst2, etmax
if (etmax .le. -1.0d0 .or. ABS(enpar) .lt. 1.d-6) then
rjpd=0.d0
rjpd0=0.d0
ratjps=1.d0
denom=1.d0
return
endif
c
c write(*,*)'etmax,thtc', etmax,thtc
if (etmax .gt. 0.d0 .and. thtc .gt. 0.d0) then
if (epst1 .gt. 0.d0 .and. epst2 .lt. etmax) then
c
c ---calculate integral for resonance divided by trapped region:
c
expel = EXP (-etmax)
expe1 = EXP (-epst1)
expe2 = EXP (-epst2)
n = n0 / 2
rint = qgauleg_GA(fcn_GA,expel,expe2,n,ifail)
! write(*,*)'1,expel,expe2,n,ifail,rint',expel,expe2,n,ifail,rint
rint = rint+qgauleg_GA(fcn_GA,expe1,expeu,n,ifail)
! write(*,*)'2 expe1,expeu,n,ifail,rint', expe1,expeu,n,ifail,rint
else
c
c ---calculate integral for resonance which doesn't enter trapped region
c
n = MIN0 (64, n0)
if (epst2 .le. 0.0d0) expel = EXP (-etmax)
if (epst2 .gt. 0.0d0 .and. epst1 .le. 0.0d0) expel = EXP (-MIN (
& epst2, etmax))
if (epst2 .gt. 0.0d0 .and. epst1 .gt. 0.0d0) expel = EXP (-MIN (
& epst1, etmax))
rint = qgauleg_GA(fcn_GA,expel,expeu,n,ifail)
! write(*,*)'3 expel,expeu,n,ifail,rint', expel,expeu,n,ifail,rint
endif ! (epst1 .gt. 0.d0 .and. epst2 .lt. etmax)
cSAP080617
c write(*,*)'rint',rint
c
c ---calculate denominator
c
expeld = EXP (-etmax)
c%LL call denomf(denom,expeld,etmax)
!write(*,*)'expeld,expeu=',expeld,expeu
denom=theth*qgauleg_GA(fcd,expeld,expeu,n,ifail)
!write(*,*)'denom,zhat,ifail', denom,zhat,ifail
c%LL PRINT *, denom, zdenom
rjpd = -0.125d0*sgnj*thtcu*rmax*rint/(denom*zhat)
c write(*,*)'1 rjpd', rjpd
c
c bucket lift calculation:
else if (etmax .gt. 0.d0 .and. thtc .le. 0) then
erise = -thtc/511.0d0
etal = 0.0d0
sgnj0 = SIGN (1.0d0,(gammin-elom)*enpar)
gammax = gammin+erise
etau = MIN (1.0d0,2.0d0*(rabs/rmax)*elom*erise/ (gammax*gammax-1.0
&d0))
call TorGA_getdfh (dfh, gammin, gammax, etal, etau)
rjpd = -(511.0d0/(8.d0*zhat*tc))*sgnj0*rmax*dfh/erise
c write(*,*)'2 rjpd', rjpd
c
c set j/pd to zero if no resonance
else
rjpd = 0.0d0
gammin = 1.0d0
endif ! (etmax .gt. 0.d0 .and. thtc .gt. 0.d0)
h3fac = 4.0d0 * zhat / (5.0d0 + charge)
rjpd = rjpd*h3fac
c write(*,*)'3 rjpd', rjpd
c
c --- calculate non-relativistic, zero-trapped-particle limits
c
sigma = 1.0d0 - elomom
eparnr = 0.5d0*(1.0d0-elom/gammin)**2/(theth*enpar2)
sgnj0 = SIGN (1.0d0,(gammin-elom)*enpar)
if (iheat .eq. +1) rjpd0 = sgnj0*fjpd0_GA(eparnr,l,sigma,zeff)
& *thtcu
if (iheat .eq. -1) rjpd0 = sgnj0*fjpd0h_GA(eparnr,zeff)*thtcu
ratjpd = rjpd / rjpd0
c write(*,*)' ratjpd,rjpd,rjpd0', ratjpd,rjpd,rjpd0
return
c
end subroutine TorGA_curba
subroutine TorGA_setpar
cProlog
implicit none
c
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
& rl, ra0, rmax, rabs, r3, asp1, aspect, zhat,
&charge, agrel, bgrel, ap0, rmin, eparsq, epar, alfsq, elooro, elom
&, elom2, rlams, etam, alffac, pas, tor, enpar2, theth, emin,
&gammin, qgrel
integer l, lm1, lp1, modl2, nc, nf, iheat, ixo, kg
c
c set eta-independent parameters
ccSAP080617
complex*16
&alpha,chalf,cwun,ci,zarg,cpas,q
common /TorGA_fcncom/ aspect,epar,rl,rmin,rmax,rabs,etam,charge,
&asp1,zhat,ap0,alfsq,elooro,r3, eparsq,l,iheat,lm1,lp1,modl2
common /TorGA_rela/ enpar2,elom,theth,emin,ra0,elom2,tor, gammin,
&agrel,bgrel,qgrel,ixo,kg
common /TorGA_mdl3com/alpha,rlams,pas
data chalf/(-.5d0,0.0d0)/,cwun/(1.0d0,0.0d0)/,ci/(0.0d0,1.0d0)/
c
rl=l
ra0=rmax/rabs
r3=3.0d0*ra0
asp1=1.0d0-aspect
zhat=(1.0d0+charge)/2.d0
agrel=(1.0d0/(5.d0+charge))**2
bgrel=1.0d0-agrel
qgrel=1.0d0/bgrel
ap0 = -0.5d0*sqrt (rmax/rmin)*rmax
eparsq=epar*epar
alfsq=-2.0d0*aspect/asp1
lm1=l-1
lp1=l+1
elooro=elom/ra0
elom2=elom*elom
c parameters for model=3
rlams = sqrt (1.0d0-etam)
alffac=(4.d0/zhat)-.25d0
if (alffac) 10,10,20
10 alpha=chalf+cwun*sqrt (-alffac)
go to 30
20 alpha=chalf+ci*sqrt (alffac)
30 zarg = dcmplx(rlams, 0.0d0)
cSAP080617
c write(*,*)'TorGA_setpar modl2,nf ',modl2,nf
cSAP080618
nf=0
if (modl2 .eq. 1) call TorGA_legfn(alpha,zarg,cpas,q,nc,nf)
cSAP080617
c write(*,*)'after TorGA_legfn alpha,zarg,cpas,q,nc,nf',
c &alpha,zarg,cpas,q,nc,nf
pas = dble(cpas)
cSAP080617
c write(*,*)'pas',pas
tor=2.0d0/ra0
return
c
end
subroutine TorGA_limits (epst1,epst2,etmax)
cProlog
implicit none
c
c calculate range of relativistic range of integration
c etmax is emax-emin, normalized to thot.
c
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
&etmax0, gam1, gam2, gammin, etmax, etcrit, theth,
&gam3, gam4, etam, epst1, epst2, emin, enpar2, aspect, epar, rl,
&rmin, rmax, rabs, charge, asp1, zhat, ap0, alfsq, qgrel, r3,
&eparsq, elom, ra0, elom2, tor, agrel, bgrel, elooro
integer iheat, l, modl2, lm1, lp1, ixo, kg
c
common /TorGA_fcncom/ aspect,epar,rl,rmin,rmax,rabs,etam,charge,
&asp1,zhat,ap0,alfsq,elooro,r3, eparsq,l,iheat,lm1,lp1,modl2
common /TorGA_rela/ enpar2,elom,theth,emin,ra0,elom2,tor, gammin,
&agrel,bgrel,qgrel,ixo,kg
data etmax0/20.d0/
c
if (iheat) 60,60,5
c ECH limits:
c find intersections of resonance with ppar axis:
5 call TorGA_gammas(gam1,gam2,0.0d0)
if (gam1-1.0d0) 10,10,20
10 if (gam2 .le. 1.0d0) go to 40
gammin=gam2
etmax=etmax0
go to 30
20 gammin=gam1
etcrit=(gam2-gam1)/theth
etmax=MIN (etmax0,etcrit)
c find intersections of resonance with passing-trapped separatrix:
30 call TorGA_gammas(gam3,gam4,etam)
epst1=(gam3-gammin)/theth
epst2=(gam4-gammin)/theth
emin=(gammin-1.0d0)/theth
return
c
c --- no resonance if elom+enpar2-1 < 0; set etmax negative as a flag.
c
40 etmax=-1.0d0
return
c Limits for Landau resonance:
60 etmax=etmax0
epst2=-1.0d0
gammin=1.0d0/sqrt (1.0d0-1.0d0/enpar2)
emin=(gammin-1.0d0)/theth
return
c
end
subroutine TorGA_getdfh (dfh, gammin, gammax, etal, etau)
cProlog
implicit none
c
c calculates difference of final and initial greens functions for
c bucket lift from etal,gammin to etau,gammax
c
c
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
& gammin, fl, fupr, gammax, fu, hl, hpr, etal, etau,
& etam, hu, dfh, gam, eta, psq, gamsq, aspect, epar, rl, rmin, rmax
&, rabs, charge, asp1, zhat, ap0, alfsq, elooro,r3, eparsq
integer l, modl2, lm1, lp1, iheat
c
common /TorGA_energy/gam,eta,psq,gamsq
common /TorGA_fcncom/ aspect,epar,rl,rmin,rmax,rabs,etam,charge,
&asp1,zhat,ap0,alfsq,elooro,r3, eparsq,l,iheat,lm1,lp1,modl2
c
call TorGA_setrel(gammin)
call TorGA_ffp(fl,fupr)
call TorGA_setrel(gammax)
call TorGA_ffp(fu,fupr)
c if (modl2) 10,10,20
c 10 hl=hfun(etal,hpr)
c if (etau .ge. etam) go to 30
c hu=hfun(etau,hpr)
c dfh=fu*hu-fl*hl
c return
20 call TorGA_geth3(hl,hpr,etal)
if (etau .ge. etam) go to 30
call TorGA_geth3(hu,hpr,etau)
dfh=fu*hu-fl*hl
return
30 dfh=-fl*hl
return
c
end
subroutine TorGA_legfn (v, z, p, q, nc, nf)
cProlog
implicit none
c
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
& ac, acc, pi, sv, cv, eim, shv, chv, vr, vi, r1, r2
&, r3, r4, rr, vri, eip, pisr
cSAP080617
real*8
& z_imag, zn24
integer nvv, ninter, n23, nff, nf, n24, nfrig, ncvg, nc
c
cSAP080617
real*8
& qinf
cSAP080617
complex*16
& vv,zz,pp,qq,z1,z2,u,a,b_GA,c,gr,cvv,svv,zz1,zz2,srz,
&v,z,p,q, zzs,cisp,cism,pt,vvp,zzz, csqrtk_GA
c....&|--1---------2---------3---------4---------5---------6---------7-|
common /TorGA_legbl/ vv,zz,pp,qq,z1,z2,u,a,b_GA,c,gr,cvv,svv, zz1,
&zz2,srz,acc,pisr,pi,r1,r2,r3,r4,ncvg,nfrig
c
data ac /0.0000001d0/
data qinf /1.79769313486231d+308/
c uses c311bd, legv, leg1, legz, leror
c
cSAP080629: nf initialization
nf=0
cSDAP080617
c write(*,*)'TorGA_legfn v, z, p, q, nc, nf',v, z, p, q, nc, nf
c write(*,*)'pp',pp
call TorGA_c311bd
c write(*,*)'v,ac',v,ac
acc = ac**2
if (dble(v)+0.5d0) 16,17,17
16 vv = -v - 1.0d0
go to 18
17 vv = v
18 nvv = ninter (dble(vv))
cSAP080617
c write(*,*)'18 vv,nvv',vv,nvv
30 vvp = vv * pi
cSAP080617
c write(*,*)'30 vvp',vvp
go to (21, 22, 23,2 4, 25), nvv
21 sv = 1.0d0
go to 26
23 sv = -1.0d0
26 cv = 0.0d0
go to 28
22 cv = -1.0d0
go to 27
24 cv = 1.0d0
27 sv = 0.0d0
go to 28
25 sv = SIN (dble(vvp))
cv = COS (dble(vvp))
28 eip = EXP ( dimag(vvp))
eim = EXP (-dimag(vvp))
shv = 0.5d0 * (eip-eim)
chv = 0.5d0 * (eip+eim)
svv = sv*chv+u*cv*shv
cvv = cv*chv-u*sv*shv
cisp = cvv+u*svv
cism = cvv-u*svv
cSAP080617
c write(*,*)'28 sv,cv',sv,cv
c write(*,*)'eip,eim,vvp',eip,eim,vvp
c write(*,*)'shv,chv',shv,chv
c write(*,*)'svv,cvv,u',svv,cvv,u
c write(*,*)'cisp,cism,z',cisp,cism,z
c write(*,*)'pp',pp
if (dble(z)) 9, 10, 11
9 zz = -z
n23 = 3
nff = -nf
c write(*,*)'9 zz,nff',zz,nff
go to 12
10 n24 = 4
zz = z
if (dimag(zz)) 13, 77, 13
11 zz = z
n23 = 2
nff = nf
12 n24 = 2
13 zzs = zz**2
z1 = (1.0d0-zz) / 2.0d0
c write(*,*)'13 zz,zzs,z1',zz,zzs,z1
c write(*,*)'n23,nf',n23,nf
if (dimag(zz)) 7, 6, 7
6 if ( dble(z1)) 4, 70, 5
4 nfrig = ISIGN (n23, nff)
c write(*,*)'4 n23,nff, nfrig', n23,nff, nfrig
go to 8
5 nfrig = nff
go to 8
7 continue
z_imag = dimag(zz)
c write(*,*)'7 zz,z_imag',zz,z_imag
zn24 = FLOAT (n24)
c write(*,*)'n24,zn24',n24,zn24
nfrig = SIGN (zn24, z_imag)
c write(*,*)'nfrig',nfrig
8 z2 = 1.0d0/zzs
srz = csqrtk_GA (zzs-1.0d0, nfrig, 1)
c write(*,*)'8 zzs,nfrig,srz',zzs,nfrig,srz
zz1 = ( zz+srz) / (2.0d0*srz)
zz2 = (-zz+srz) / (2.0d0*srz)
c write(*,*)'zz,srz,zz1,zz2',zz,srz,zz1,zz2
c write(*,*)' dble(z)',dble(z)
if (dble(z)) 1, 2, 2
1 srz = -srz
2 zzz = zz
zz = z
vr = dble(vv)**2
vi = dimag(vv)**2
r1 = ABS (z1)
r2 = ABS (z2)
r3 = ABS (zz1)
r4 = ABS (zz2)
c
if (r1-r2)61,61,62
61 rr = MAX (r3,r4)/0.045d0-19.5d0
if (rr)66,66,65
65 vri = vr+vi
if (vri)165,67,165
c write(*,*)'before 165 pp',pp
165 if (rr**2-vri+vi**2/(2.0d0*vri)) 66, 67, 67
66 call TorGA_legv
c write(*,*)'66 after TorGA_legv pp',pp
if (ncvg)166,97,166
166 zz = zzz
call TorGA_leg1
c write(*,*)'166 after TorGA_leg1 pp',pp
go to 90
67 zz = zzz
call TorGA_leg1
c write(*,*)'67 after TorGA_leg1 pp',pp
if (ncvg)167,90,167
167 zz = z
call TorGA_legv
c write(*,*)'167 after TorGA_legv pp',pp
go to 97
62 if (vr+vi-16.0d0*r2**2)63,64,64
64 call TorGA_legv
c write(*,*)'64 after TorGA_legv pp',pp
if (ncvg)164,97,164
164 call TorGA_legz
c write(*,*)'164 after TorGA_legz pp',pp
go to 97
63 call TorGA_legz
c write(*,*)'63 after TorGA_legz pp',pp
if (ncvg)163,97,163
163 call TorGA_legv
c write(*,*)'163 after TorGA_legv pp',pp
go to 97
90 if (dble(z)) 93,97,97
93 pt = -2.0d0/pi*svv*qq
if (nfrig) 94,95,96
94 pp = cisp*pp+pt
c write(*,*)'94 pp,cisp,pt',pp,cisp,pt
qq = -cism*qq
go to 97
95 qq = -cvv*qq-pi/2.0d0*svv*pp
pp = cvv*pp+pt
c write(*,*)'95 pp,cvv,pt', pp,cvv,pt
go to 97
96 pp = cism*pp+pt
c write(*,*)'96 pp,cism,pt',pp,cism,pt
qq = -cisp*qq
go to 97
c
97 if (dble(v)+0.5d0) 91,92,92
91 if (ABS (svv) .ne. 0.0d0) go to 98
qq = qinf
go to 92
98 qq = (qq*svv-pi*cvv*pp)/svv
92 p = pp
c write(*,*)'92 p,pp',p,pp
q = qq
nc = ncvg
return
70 qq = qinf
c
if ( dble(z)) 173, 173, 74
173 if (dimag(v)) 71, 73, 71
73 go to (71,72,71,74,71),nvv
71 pp = qq
c write(*,*)'71 pp,qq',pp,qq
go to 82
72 pp = -1.0d0
go to 82
74 pp = 1.0d0
go to 82
c
77 nfrig = nf
call TorGA_legor
c write(*,*)'77 after TorGA_legor pp',pp
82 ncvg = 0
go to 92
c
end
subroutine TorGA_setrel (gamma)
cProlog
implicit none
c
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
& gam, gamma, gamsq, psq, eta
c
common /TorGA_energy/ gam, eta, psq, gamsq
c
gam = gamma
gamsq = gam * gam
psq = gamsq - 1.0d0
return
c
end
subroutine TorGA_ffp (fu, fupr)
cProlog
implicit none
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
& cgrel, bfac, bgrel, psq, bfacq, qgrel, afac, agrel
&, rutaf, fu, fupr, gam, gamsq, denom, eta, enpar2, elom, theth,
&emin, ra0, elom2, tor, gammin
integer kg, ixo
c
c energy-dependent part of Green's function and derivative.
c the Green's function is defined as:
c g=<B_t/R>c^4/4/nuv_t^3B_0)Fu(v/c)h(eta)
c for kg+2=ig=1, fu=p^2(1-1/(1+bgrel*p^2)**qgrel)/sqrt (1+agrel*p^2)
c with agrel=1/(5+Z)^2, bgrel=1-agrel, qgrel=1/bgrel. This form
c gives Green's function which is correct in nonrelativistic
c limit; in limit of no trapped particles reproduces leading-order
c (v^6/c^6) relativistic correction as well as ultrarelativistic
c limit as given by Fisch.
c for kg+2=ig=2, fu=gamma*v^4/c^4, (correct nonrelativistically; correct
c energy dependence at large energy but wrong coefficient)
c for kg+2=ig=3, fu=p^4/(1+cgrel*p^2)^(2/3) with
c p=momentum/mc and cgrel=2/3 (correct to order v^6/c^6 and correct
c energy dependence at large energy but wrong coefficient)
c
common /TorGA_energy/gam,eta,psq,gamsq
common /TorGA_rela/ enpar2,elom,theth,emin,ra0,elom2,tor, gammin,
&agrel,bgrel,qgrel,ixo,kg
data cgrel/.66666666666667d0/
c
if (kg) 10,20,30
c
c --- ig=1 ------------------
c
10 bfac=1.0d0+bgrel*psq
bfacq=bfac**qgrel
afac=1.0d0+agrel*psq
rutaf = sqrt (afac)
fu=psq*(1.0d0-1.0d0/bfacq)/rutaf
fupr=2.0d0*gam*(fu*(1.0d0/psq-.5d0*agrel/afac)+psq/(rutaf*bfac*
&bfacq))
return
c
c --- ig=2 ------------------
c
20 fu=psq*psq/(gamsq*gam)
fupr=psq*(1.0d0+3.0d0/gamsq)/gamsq
return
c
c --- ig=3 ------------------
c
30 denom=1.0d0+cgrel*psq
fu=psq*psq/denom**1.5d0
fupr=fu*gam*(4.d0+cgrel*psq)/(psq*denom)
return
c
end
subroutine TorGA_geth3 (hh, hpr, eta)
cProlog
implicit none
c
c gets angle part of green's function in Legendre-function (square-well)
c approx. Note that h calculated here needs to be multiplied by
c 4/(Z+5) to get h of UCRL-95813
c
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
& rlam, eta, rpa, rppr, hh, rlams, pas, hpr
integer nc, nf
c
cSAP080617
complex*16
&alpha, zarg, pa, ppr, q
common /TorGA_mdl3com/ alpha, rlams, pas
c
c write(*,*)'TorGA_geth3 eta',eta
rlam = sqrt (1.0d0-eta)
zarg = dcmplx(rlam,0.0d0)
cSAP080617
c write(*,*)' rlam,zarg',rlam,zarg
call TorGA_legfn1(alpha,zarg,pa,ppr,q,nc,nf)
c write(*,*)'after TorGA_legfn1 alpha,zarg,pa,ppr,q,nc,nf',
c & alpha,zarg,pa,ppr,q,nc,nf
rpa = dble(pa)
rppr = dble(ppr)
c write(*,*)'pa,rpa',pa,rpa
c write(*,*)'ppr,rppr',ppr,rppr
hh=-rlam+rlams*rpa/pas
c write(*,*)'hh,rlam,rlams,rpa,pas',hh,rlam,rlams,rpa,pas
hpr=(1.0d0-rlams*rppr/pas)*.5d0/rlam
c write(*,*)'hpr,rlams,rppr,pas,rlam',hpr,rlams,rppr,pas,rlam
return
c
end
subroutine TorGA_c311bd
cProlog
implicit none
c
cSAP080617
complex*16
&vv,zz,pp,qq,z1,z2,u,a,b_GA,c,gr,cvv,svv,zz1,zz2,srz
c....&|--1---------2---------3---------4---------5---------6---------7-|
common /TorGA_legbl/ vv,zz,pp,qq,z1,z2,u,a,b_GA,c,gr,cvv,svv, zz1,
&zz2,srz,acc,pisr,pi,r1,r2,r3,r4,ncvg,nfrig
c
c explicit type declaration added 7/12/01 RAJ
cSAP080617
real*8
& pisr, pi, acc,r1, r2, r3, r4
integer ibd, ncvg, nfrig
c
data ibd/0/
c
if (ibd .ne. 0) return
ibd = 1
pisr = 1.77245385090552d0
pi = 3.141592653589793d0
u = (0.0d0,1.0d0)
return
c
end
subroutine TorGA_legv
cProlog
implicit none
c
c explicit type declaration 7/12/01 (RAJ)
cSAP080617
real*8
&acc, sgn, pisr, pi, r1, r2, r3, r4
integer ncvg, nfrig, ncv
cSAP080617
complex*16
& vv,zz,pp,qq,z1,z2,u,a,b_GA,c,gr,cvv,svv,zz1,zz2,srz
c
c....&|--1---------2---------3---------4---------5---------6---------7-|
common /TorGA_legbl/ vv,zz,pp,qq,z1,z2,u,a,b_GA,c,gr,cvv,svv, zz1,
&zz2,srz,acc,pisr,pi,r1,r2,r3,r4,ncvg,nfrig
c
cSAP080617
complex*16
& f1, f2, ssrz, clogok_GA, csqrtk_GA, rgam_GA
c uses hypgm, clogok_GA, csqrtk_GA, rgam_GA
c
a=0.5d0
c=vv+1.5d0
c write(*,*)'in TorGA_legv a,c,zz1,zz2',a,c,zz1,zz2
call TorGA_hypgm(a,a,c,zz1,f1,acc,ncvg)
cyup write(*,*)'f1,ncvg,acc',f1,ncvg,acc
call TorGA_hypgm(a,a,c,zz2,f2,acc,ncv )
c write(*,*)'f2,ncv,acc',f2,ncv,acc
ncvg=ncvg+2*ncv
6 f1=f1*EXP ((vv+0.5d0)*clogok_GA(zz +srz,nfrig,3))
c write(*,*)'zz,srz,nfrig,clogok_GA(zz +srz,nfrig,3)',
c &zz,srz,nfrig,clogok_GA(zz +srz,nfrig,3)
c write(*,*)'f1,zz,srz,nfrig',f1,zz,srz,nfrig