R2D2f

Fortran source

0001 C
0002 C   CODE R2D2
0003 C
0004 C#####################################################################
0005 C          evaluating the non-relativistic and relativistic          #
0006 C        dielectric tensor elements (two different formalisms)       #
0007 C   (only for waves in the electron cyclotron range of frequencies)  #
0008 C#####################################################################
0009 C
0010 C---------------------------------------------------------------------
0011 C
0012 C VERSION 5/3J:September 23, 2003
0013 C             May 2003
0014 C              (modified the integration routine to check for
0015 C               convergence -- seems to be doing fine)
0016 C             June 5, 2003
0017 C              (incorporated the weiss' way of doing things)
0018 C             June 9, 2003
0019 C              (incorporated solving of dispersion
0020 C               relation at different positions in a plasma)
0021 C             June 10, 2003
0022 C              (seems to be working)
0023 C              (Note: need to reduce accuracy of root finder for
0024 C                     some relativistic cases to 1.0e-6;
0025 C                     true specially near doppler-shifted harmonics)
0026 C              (Note: need to keep at most 20 terms in the 
0027 C                     non-relativistic dielectric tensor elements;
0028 C                     accuracy issues for higher harmonic terms)
0029 C             June 24, 2003
0030 C              (added polarization subroutine for diffusion
0031 C               coefficient calculations)
0032 C             July 27, 2003
0033 C             (extended to arbitrary number of roots -- use the
0034 C              n_rts to change the number of roots in various
0035 C              routines)
0036 C             September 19, 2003
0037 C              (extended evaluation of polarizations to the
0038 C               relativistic dispersion tensor)
0039 C             September 23, 2003
0040 C              (the number of terms to be kept in nonrelativistic
0041 C               tensor element sums is now an input parameter)
0042 C              (Note: sometimes to follow the EBW near the resonance
0043 C                     may need 150 or so elements instead of the 20
0044 C                     mentioned above.)
0045 C              (zeroed small elements in the polarization vector)
0046 C             October 10, 2003
0047 C              (added non-relativistic power flow/ power dissipated
0048 C               routines including solving the hermitian dispersion
0049 C               tensor for roots)
0050 C              (note: in the routine for solving the hermitian 
0051 C                     tensor one can change the number of intervals
0052 C                     in which to search for roots, and the range 
0053 C                     over which to solve for the roots)
0054 C                     (finished and tested October 14, 2003)
0055 C             December 2, 2004
0056 C              (fixed the modified Bessel function identity in
0057 C               besk2 in three places -- courtesy of Eric)
0058 C             April 5, 2005
0059 C              (added parallel and y-direction power flows in the 
0060 C               non-relativistic power flow package -- pwr_nr)
0061 C             April to October 2005
0062 C              (added relativistic power flow package)
0063 C
0064 C---------------------------------------------------------------------
0065 C
0066 C
0067       subroutine R2D2(rout,pini,pinr,gues,pos)
0068 C
0069       implicit double precision (a-h, o-z)
0070 C
0071       parameter (n_rts=10)
0072 C
0073       double precision guesr(n_rts),guesi(n_rts)
0074       double complex croots(4),roots(n_rts),frt(n_rts)
0075 C
0076       character*20 tokamak,tokparam
0077       integer l_tokamak
0078 C
0079       real*8 rout(2*(4+2*n_rts))
0080       real*8 pini(21),pinr(26),gues(2*n_rts),pos(32768*5)
0081 C
0082       common /para1/ w0
0083       common /para2/ npos
0084       common /para4/ den0,bt0,te0
0085       common /para5/ wpe0,wce0,vte0
0086 C
0087       common /nrels/ nr_nterms
0088 C
0089       common /trub0/ ktrub
0090       common /trub1/ xi_lo,dxi,nxi
0091       common /trub2/ diff_err,navg
0092       common /trerr/ errabs0,errrel0,irule0
0093 C
0094       common /ppara/ p_para_min,p_para_max
0095       common /pperp/ p_perp_min,p_perp_max
0096       common /weerr/ errabs1,errrel1,irule1
0097 C
0098       common /nrts1/ tol_ht,knr,nrt_nr,nrt_fr,nd_ht
0099 C
0100       common /chck1/ kparam
0101       common /chck2/ cnpa_in,cnpa_fi,n_cnpa
0102       common /chck3/ omega_in,omega_fi,n_omega
0103       common /chck4/ alpha_in,alpha_fi,n_alpha
0104 C
0105       common /rterr/ errt_abs,errt_rel,it_max
0106 C
0107       common /wave0/ cnpa
0108 C
0109       common /guinp/ guesr,guesi,kgues
0110 C
0111       common /prnt/  kpr
0112 C
0113 C   Simulation parameters
0114 C
0115 C   open(20,file='simparam.m')'simparam.m')
0116 C
0117 C   Code parameters
0118 C
0119 C      open(22,file='codeparam.m')'codeparam.m')
0120 C
0121 C....  output non-relativistic and relativistic roots .....
0122 C....... (data output is the local position and all the roots) .......
0123 C.... (kparam = 1: local cnpa, omega, or alpha and all the roots) ....
0124 C
0125 C      open(5,file='nr1.m',status='replace')'replace')
0126 C      open(8,file='r1.m',status='replace')'replace')
0127 C
0128 C........ the non-relativistic and relativistic polarizations ........
0129 C..... polarizations are w.r.t. E_x and output as E_+, E_-, E_z ......
0130 C.. the first element is either the position, cnpa, omega, or alpha ..
0131 C...... the second/third elements are cnpe -- the complex root .......
0132 C......... the other elements are the complex polarizations ..........
0133 C
0134       open(4,file='pn.m',status='replace')
0135       open(2,file='pr.m',status='replace')
0136 C
0137 C..... output file containing information on elements and roots ......
0138 C
0139       open(6,file='d1.m',status='replace')
0140 C
0141 C... output file containing local plasma parameters with position ....
0142 C......... position, wpe, wce, wuh, 1/wce, vte, bt, den, te ..........
0143 C
0144       open(7,file='p1.m',status='replace')
0145 C
0146 C...... output file containing all the run-time error messages .......
0147 C
0148       open(9,file='errmsg1',status='replace')
0149 C
0150 C....... output file for checking various parts of the program .......
0151 C.................... (presently not being used) .....................
0152 C
0153       open(10,file='out1',status='replace')
0154 C
0155 C
0156 C........ controlling whether solving the dispersion relation ........
0157 C.......... in the conventional form (i.e. along x) or as a ..........
0158 C.............. function of n_parallel, omega, or alpha ..............
0159 C
0160 C........ if knp = 0 ==> cnpa is actually k_para in cm**(-1) .........
0161 C...... if knp = 1 ==> cnpa is normalized; cnpa = (c*k_para/w) .......
0162 C...... (NOTE: knp is used only for the conventional d.r. form) ......
0163 C....... (NOTE: otherwise, it is assumed that cnpa is n_para) ........
0164 C.............. (omega = w/wce,  alpha = (wpe/wce)**2) ...............
0165 C
0166 C...... kparam = 0 solves the conventional dispersion relation .......
0167 C..... (here the input n_para or k_para is as given in cnpa_in) ......
0168 C..... kparam = 1 solves the d.r. for a range in n_para (param1) .....
0169 C....... (omega and alpha are fixed at omega_in and alpha_in) ........
0170 C..... kparam = 2 solves the d.r. for a range in omega (param2) ......
0171 C....... (n_para and alpha are fixed at cnpa_in and alpha_in) ........
0172 C..... kparam = 3 solves the d.r. for a range in alpha (param3) ......
0173 C....... (n_para and omega are fixed at cnpa_in and omega_in) ........
0174 C...... initial n_para, final n_para, number of steps in n_para ......
0175 C....... initial omega, final omega, number of steps in omega ........
0176 C....... initial alpha, final alpha, number of steps in alpha ........
0177 C.. (NOTE: for kparam.ge.1 the temperature used is te0 given below) ..
0178 C
0179 C...... output file containing power flow and power dissipated .......
0180 C...... (output: input parameter e.g. position, cnpe, sx, pdis) ......
0181 C
0182       open(11,file='pf.m',status='replace')
0183       open(13,file='pfr.m',status='replace')
0184 C
0185 C      namelist /param0/ kparam
0186 C      namelist /param1/ knp,cnpa_in,cnpa_fi,n_cnpa
0187 C      namelist /param2/ omega_in,omega_fi,n_omega
0188 C      namelist /param3/ alpha_in,alpha_fi,n_alpha
0189 C
0190 C............ if kw0 = 0 ==> w0 is the frequency in Hertz ............
0191 C............ if kw0 = 1 ==> w0 is the angular frequency  ............
0192 C
0193 C      namelist /para/ kw0
0194 C
0195 C... error limits for root finder and maximum number of iterations ...
0196 C
0197 C      namelist /rtfnd/ errt_abs,errt_rel,it_max
0198 C
0199 C............... for the trubnikov integration package ...............
0200 C........... ktrub = 1 uses the trubnikov form of the d.r. ...........
0201 C............. (this is a single time-history integral) ..............
0202 C........... ktrub = 0 uses the momentum  form of the d.r. ...........
0203 C.......... (this is the double integral in momenta space) ...........
0204 C.......... ktrub = 2 gives results for both representation ..........
0205 C.... (NOTE: used only by the rel_chk for checking d.r. elements) ....
0206 C......... the upper and lower limit for the xi-integration ..........
0207 C..... nxi = number of intervals in which this range is divided ......
0208 C.... absolute and relative errors and the rule for integration ......
0209 C.. navg = number of intervals for averaging the integration result ..
0210 C.... diff_err = error tolerance between three of these averages .....
0211 C
0212 C      namelist /trint1/ ktrub,xi_lo,xi_hi,nxi
0213 C      namelist /trint2/ errabs0,errrel0,irule0
0214 C      namelist /trint3/ navg,diff_err
0215 C
0216 C... these are parameters when doing the integration over momenta ....
0217 C.......... minimum and maximum values of parallel momentum ..........
0218 C....... minimum and maximum values of perpendicular momentum ........
0219 C.... (NOTE: p_perp_min has to be equal to or greater than zero) .....
0220 C.... absolute and relative errors and the rule for integration ......
0221 C
0222 C      namelist /pint1/ p_para_min,p_para_max,p_perp_min,p_perp_max
0223 C      namelist /pint2/ errabs1,errrel1,irule1
0224 C
0225 C......... kpr  = 0 suppresses printing out of various data ..........
0226 C... kcheck = 1 checks the dispersion elements for a given n_perp ....
0227 C... cnpe_r and cnpe_i are the real and imaginary parts of n_perp ....
0228 C
0229 C      namelist /pwrit/ kpr,kcheck,cnpe_r,cnpe_i
0230 C
0231 C...... number of normalized radial positions ......
0232 C.............. (NOTE: these are used when kparam = 0) ...............
0233 C
0234 C      read (20,*) npos
0235        npos=int(pini(1))
0236 C
0237 C... knr = 1 solves the non-relativistic disperison relation only ....
0238 C.. knr = 0 solves the fully relativistic dispersion relation also ... 
0239 C..... number of roots to be found in the non-relativistic case ......
0240 C....... number of roots to be found in the relativistic case ........
0241 C
0242 C       read (20,*) tol_ht 
0243        knr=int(pini(2)) 
0244        nrt_nr=int(pini(3)) 
0245        nrt_fr=int(pini(4))
0246        nd_ht=int(pini(5))
0247        tol_ht=pinr(1) 
0248 C
0249 C.......... number of terms to be used in the sums for the ...........
0250 C............ non-relativistic dielectric tensor elements ............
0251 C
0252 C   read (20,*) nr_nterms
0253        nr_nterms=int(pini(6))
0254 C
0255 C............ kflag = 1 : these initial guesses are used .............
0256 C....... when solving the non-relativistic dispersion relation .......
0257 C......... (real part in guesr and imaginary part in guesi) ..........
0258 C....... kflag = 0 : these guesses are ignored and the guesses .......
0259 C.... obtained from the approximate dispersion relation are used .....
0260 C
0261 C   read (20,*) kgues
0262       kgues = int(pini(7))
0263 C
0264 C............... input guesses of roots if so desired ................
0265 C
0266 C      read (20,*) guesr(1:nrt_nr)
0267 C      read (20,*) guesi(1:nrt_nr)
0268         do i=1,nrt_nr
0269           guesr(i) = gues(2*i-1)
0270           guesi(i) = gues(2*i)
0271         end do
0272 C
0273 C.............. reference magnetic field,  ..............
0274 C.............. peak electron temperature, peak density ..............
0275 C............... frequency ...........................................
0276 C........ (NOTE: this te0 is used for all kparam given above) ........
0277 C.... (NOTE: the other parameters are used only when kparam = 0) .....
0278 C
0279 C
0280 C   read (20,*) bt0,te0,den0,w0
0281       bt0 = pinr(2)
0282       te0 = pinr(3)
0283       den0 = pinr(4)
0284       w0 = pinr(5)
0285 C
0286 C      read (22,param0)
0287 C      read (22,param1)
0288 C      read (22,param2)
0289 C      read (22,param3)
0290 C      read (22,para)
0291 C      read (22,rtfnd)
0292 C      read (22,trint1)
0293 C      read (22,trint2)
0294 C      read (22,trint3)
0295 C      read (22,pint1)
0296 C      read (22,pint2)
0297 C      read (22,pwrit)
0298 C      
0299       kparam = int(pini(8))
0300       knp = int(pini(9))
0301       n_cnpa = int(pini(10))
0302       n_omega = int(pini(11))
0303       n_alpha = int(pini(12))
0304       kw0 = int(pini(13))
0305       it_max = int(pini(14))
0306       ktrub = int(pini(15))
0307       nxi = int(pini(16))
0308       irule0 = int(pini(17))
0309       navg = int(pini(18))
0310       irule1 = int(pini(19))
0311       kpr = int(pini(20))
0312       kcheck = int(pini(21))
0313 C
0314       cnpa_in = pinr(6)
0315       cnpa_fi = pinr(7)
0316       omega_in = pinr(8)
0317       omega_fi = pinr(9)
0318       alpha_in = pinr(10)
0319       alpha_fi = pinr(11)
0320       errt_abs = pinr(12)
0321       errt_rel = pinr(13)
0322       xi_lo = pinr(14)
0323       xi_hi = pinr(15)
0324       errabs0 = pinr(16)
0325       errrel0 = pinr(17)
0326       diff_err = pinr(18)
0327       p_para_min = pinr(19)
0328       p_para_max = pinr(20)
0329       p_perp_min = pinr(21)
0330       p_perp_max = pinr(22)
0331       errabs1 = pinr(23)
0332       errrel1 = pinr(24)
0333       cnpe_r = pinr(25)
0334       cnpe_i = pinr(26)
0335 C
0336       twopi=6.283185307179586476925287
0337       clight=2.99792458d10
0338 C
0339 C..... getting the angular frequency (if kw0 = 0) in radians/sec .....
0340 C
0341       if (kparam.eq.0.and.kw0.eq.0) then
0342         w0=twopi*w0
0343       end if
0344 C
0345 C...... the integration grid for the trubnikov form (if needed) ......
0346 C
0347       dxi=(xi_hi-xi_lo)/dfloat(nxi)
0348       ci=dcmplx(0.d0,1.d0)
0349 C
0350       if (kcheck.eq.1) then
0351         call rel_chk(cnpa_in,omega_in,alpha_in,te0,cnpe_r,cnpe_i)
0352         go to 99
0353       end if
0354 C
0355 C... solving the dispersion relation for two different conditions ....
0356 C
0357       if (kparam.eq.0) then
0358         call rel_disp(croots,roots,frt,pos)
0359       else
0360         call rel_eff(croots,roots,frt,te0)
0361       end if
0362 C
0363    99 continue
0364 C
0365       do i=1,4
0366         rout(2*i-1) = dreal(croots(i))
0367         rout(2*i) = dimag(croots(i))
0368       end do
0369 C      
0370       do i=1,n_rts
0371         rout(8+2*i-1) = dreal(roots(i))
0372         rout(8+2*i) = dimag(roots(i))
0373         rout(28+2*i-1) = dreal(frt(i))
0374         rout(28+2*i) = dimag(frt(i))
0375       end do
0376 C
0377       return
0378       end
0379 C
0380 C
0381 C=====================================================================
0382 C
0383       subroutine rel_chk(cnpa_in,omega_in,alpha_in,te0,cnpe_r,cnpe_i)
0384 C
0385       implicit double precision (a-h, o-z)
0386 C
0387       double complex dcmplx,cnpe,ep(6),ci
0388 C
0389       common /para9/ alpha,omega,vt,besk2
0390       common /wave0/ cnpa
0391       common /wave1/ cnpe
0392       common /nrels/ nr_nterms
0393       common /trub0/ ktrub
0394       common /prnt/  kpr
0395 C
0396       external dbsk0e,dbsk1e,besk_a,umach,erset
0397 C
0398       cnpa=cnpa_in
0399       omega=omega_in
0400       alpha=alpha_in
0401       vt=4.19396657d7/2.99792458d10*dsqrt(te0)
0402 C
0403       cnpe=dcmplx(cnpe_r,cnpe_i)
0404 C      
0405       arg=1./vt**2
0406       if (arg.lt.50.d0) then
0407         call umach(-3,9)
0408         besk2=dbsk0e(arg)+2.*dbsk1e(arg)/arg
0409         call erset(0,1,1)
0410       else
0411         besk2=besk_a(2,arg)
0412       end if
0413 C
0414       call d_nr(nr_nterms,cnpe)
0415 C
0416       if (ktrub.eq.1.or.ktrub.eq.2) then
0417         call trub(ep)
0418         write (6,701) ep(1:2)
0419         write (6,701) ep(3:4)
0420         write (6,701) ep(5:6)
0421   701   format(1x,"trub x's = ",4g18.10)
0422       end if
0423 C
0424       if (ktrub.eq.0.or.ktrub.eq.2) then
0425         call weiss(ep)
0426         write (6,801) ep(1:2)
0427         write (6,801) ep(3:4)
0428         write (6,801) ep(5:6)
0429   801   format(1x,"weiss x's = ",4g18.10)
0430       end if
0431 C
0432 C
0433       return
0434       end
0435 C
0436 C=====================================================================
0437 C
0438       subroutine rel_disp(croots,roots,frt,pos)
0439 C
0440       implicit double precision (a-h, o-z)
0441 C
0442       parameter (n_rts=10)
0443 C
0444       double precision guesr(n_rts),guesi(n_rts)
0445 C
0446       real*8 pos(32768*5)
0447 C     
0448       double complex dcmplx,ci
0449 C
0450       double complex croots(4)
0451       double complex rguess(n_rts),roots(n_rts)
0452       double complex xrguess(n_rts,1+npos),xroots(n_rts,1+npos)
0453       double complex frtgues(n_rts),frt(n_rts)
0454       double complex frt_herm_guess(n_rts),frt_herm(n_rts)
0455       double complex frt_herm_guess_in(n_rts)
0456       double complex rguess_in(n_rts),frtgues_in(n_rts)
0457 C
0458       common /para1/ w0
0459       common /para2/ npos
0460       common /para8/ wp,wc,vte
0461       common /para9/ alpha,omega,vt,besk2
0462       common /nrts1/ tol_ht,knr,nrt_nr,nrt_fr,nd_ht
0463       common /guinp/ guesr,guesi,kgues
0464       common /prnt/  kpr
0465 C
0466       ci=dcmplx(0.d0,1.d0)
0467 C
0468       call prof_in
0469 C
0470       call pp(r,1,pos,0)
0471 C
0472 C..... solving for the cold plasma and approximate plasma roots ......
0473 C
0474       call coldrt(croots)
0475       call zero_roots(croots,4)
0476       call app_nr(rguess)
0477       call zero_roots(rguess,6)
0478 C      write (6,401) croots(1:2)
0479 C      write (6,401) croots(3:4)
0480       write (6,402) rguess(1:2)
0481       write (6,402) rguess(3:4)
0482       write (6,402) rguess(5:6)
0483   401 format("cold rt = ",4g18.10)
0484   402 format("warm rt = ",4g18.10)
0485 C
0486 C........... setting up the initial guesses if so desired ............
0487 C
0488       if (kgues.eq.1) then
0489         do 1 i=1,nrt_nr
0490           rguess(i)=dcmplx(guesr(i),guesi(i))
0491     1   continue
0492       end if
0493 C
0494 C......... solving the non-relativistic dispersion relation ..........
0495 C
0496       call nr_disp(rguess,roots,nrt_nr)
0497       call zero_roots(roots,nrt_nr)
0498       xroots(1:nrt_nr,1) = roots(1:nrt_nr)
0499 C      write (5,107) r,roots(1:nrt_nr)
0500       do 21 j=1,nrt_nr
0501         rguess(j)=roots(j)
0502         rguess_in(j)=rguess(j)
0503         write (6,501) r,roots(j)
0504   501   format(1x,"r, nr-root = ",3g18.10)
0505         call wave_polar(r,roots(j),0)
0506    21 continue
0507       call nr_disp_herm(r,roots,nrt_nr)
0508 C
0509 C........ solving the fully-relativistic dispersion relation .........
0510 C
0511       if (knr.eq.0) then
0512        do 31 j=1,nrt_fr
0513           frtgues(j)=rguess(j)
0514    31   continue
0515         call fr_disp(frtgues,frt,nrt_fr)
0516         call zero_roots(frt,nrt_fr)
0517 C        write (8,108) r,frt(1:nrt_fr)
0518         do 41 j=1,nrt_fr
0519           frtgues(j)=frt(j)
0520           frtgues_in(j)=frtgues(j)
0521           write (6,601) r,frt(j)
0522   601     format(1x,"r, rel-root = ",3g18.10)
0523           call wave_polar(r,frt(j),1)
0524           frt_herm_guess(j)=dcmplx(dreal(frt(j)),0.d0)
0525    41   continue
0526         call fr_herm_disp(r,frt_herm_guess,frt_herm,nrt_fr)
0527         do 51 j=1,nrt_fr
0528           write (6,701) r,frt_herm(j)
0529   701     format(1x,"r, hermitian rel-root = ",3g18.10)
0530           frt_herm_guess_in(j)=frt_herm(j)
0531    51   continue
0532       end if
0533 C
0534 C...... do-loop for calculating the roots at various positions .......
0535 C
0536       do 3 ipos=1,npos
0537 C
0538 C.......... reading next position r and the corresponding plasma parameters ..........
0539 C
0540         call pp(r,1,pos,ipos)
0541 C
0542 C......... solving the non-relativistic dispersion relation ..........
0543 C
0544         call nr_disp(rguess,roots,nrt_nr)
0545         call zero_roots(roots,nrt_nr)
0546 C        write (5,107) r,roots(1:nrt_nr)
0547         do 22 j=1,nrt_nr
0548           rguess(j)=roots(j)
0549           if (ipos.eq.1.or.kpr.eq.1) write (6,501) r,roots(j)
0550           call wave_polar(r,roots(j),0)
0551    22   continue
0552         call nr_disp_herm(r,roots,nrt_nr)
0553 C
0554 C........ solving the fully-relativistic dispersion relation .........
0555 C
0556         if (knr.eq.0) then
0557           call fr_disp(frtgues,frt,nrt_fr)
0558           call zero_roots(frt,nrt_fr)
0559 C          write (8,108) r,frt(1:nrt_fr)
0560           do 42 j=1,nrt_fr
0561             frtgues(j)=frt(j)
0562             if (ipos.eq.1.or.kpr.eq.1) write (6,601) r,frt(j)
0563             call wave_polar(r,frt(j),1)
0564             frt_herm_guess(j)=frt_herm(j)
0565    42     continue
0566           call fr_herm_disp(r,frt_herm_guess,frt_herm,nrt_fr)
0567         end if
0568 C
0569     3 continue
0570 C
0571     4 continue
0572 C
0573   107 format(<2*nrt_nr+1>g18.10)
0574   108 format(<2*nrt_fr+1>g18.10)
0575 C
0576 C
0577       return
0578       end
0579 C
0580 C=====================================================================
0581 C
0582       subroutine rel_eff(croots,roots,frt,te0)
0583 C
0584       implicit double precision (a-h, o-z)
0585 C
0586       parameter (n_rts=10)
0587 C
0588       double precision guesr(n_rts),guesi(n_rts)
0589 C
0590       double complex dcmplx,ci
0591 C
0592       double complex croots(4)
0593       double complex rguess(n_rts),roots(n_rts)
0594       double complex frtgues(n_rts),frt(n_rts)
0595       double complex frt_herm_guess(n_rts),frt_herm(n_rts)
0596 C
0597       common /para9/ alpha,omega,vt,besk2
0598       common /nrts1/ tol_ht,knr,nrt_nr,nrt_fr,nd_ht
0599       common /chck1/ kparam
0600       common /chck2/ cnpa_in,cnpa_fi,n_cnpa
0601       common /chck3/ omega_in,omega_fi,n_omega
0602       common /chck4/ alpha_in,alpha_fi,n_alpha
0603       common /wave0/ cnpa
0604       common /guinp/ guesr,guesi,kgues
0605       common /prnt/  kpr
0606 C
0607       external dbsk0e,dbsk1e,besk_a,umach,erset
0608 C
0609       clight=2.99792458d10
0610       ci=dcmplx(0.d0,1.d0)
0611 C
0612 C................ vt is the electron thermal velocity ................
0613 C
0614       vt=4.19388628d7*dsqrt(te0)/clight
0615 C
0616 C... the modified bessel function in the denominator of integrals ....
0617 C...... (this is multiplied by exp(x) where x is the argument) .......
0618 C
0619       arg=1./vt**2
0620       if (arg.lt.50.d0) then
0621         call umach(-3,9)
0622         besk2=dbsk0e(arg)+2.*dbsk1e(arg)/arg
0623         call erset(0,1,1)
0624       else
0625         besk2=besk_a(2,arg)
0626       end if
0627 C
0628 C.... kparam = 1 ==> dispersion solved for different n_para       ....
0629 C....        = 2 ==> dispersion solved for different w0/wce       ....
0630 C....        = 3 ==> dispersion solved for different (wpe/wce)**2 ....
0631 C
0632       if (kparam.eq.1) then
0633 C
0634          omega=omega_in
0635          alpha=alpha_in
0636 C
0637         if (n_cnpa.eq.1) then
0638           d_cnpa=0.d0
0639         else
0640           d_cnpa=(cnpa_fi-cnpa_in)/dfloat(n_cnpa-1)
0641         end if
0642 C
0643         cnpa=cnpa_in-d_cnpa
0644 C
0645         do 1 i=1,n_cnpa
0646 C
0647           cnpa=cnpa+d_cnpa
0648 C
0649           write (6,301) cnpa
0650   301     format(1x,"cnpa = ",g18.10)
0651 C
0652 C..... solving for the cold plasma and approximate plasma roots ......
0653 C
0654           if (i.eq.1) then
0655             call coldrt(croots)
0656             call zero_roots(croots,4)
0657             call app_nr(rguess)
0658             call zero_roots(rguess,6)
0659             write (6,401) croots(1:2)
0660             write (6,401) croots(3:4)
0661             write (6,402) rguess(1:2)
0662             write (6,402) rguess(3:4)
0663             write (6,402) rguess(5:6)
0664   401       format("cold rt = ",4g18.10)
0665   402       format("warm rt = ",4g18.10)
0666           end if
0667 C
0668 C........... setting up the initial guesses if so desired ............
0669 C
0670       if (i.eq.1.and.kgues.eq.1) then
0671         do 61 ig=1,nrt_nr
0672           rguess(ig)=dcmplx(guesr(ig),guesi(ig))
0673    61   continue
0674       end if
0675 C
0676 C......... solving the non-relativistic dispersion relation ..........
0677 C
0678           call nr_disp(rguess,roots,nrt_nr)
0679           call zero_roots(roots,nrt_nr)
0680 C          write (5,107) cnpa,roots(1:nrt_nr)
0681           do 21 j=1,nrt_nr
0682             rguess(j)=roots(j)
0683             if (i.eq.1.or.kpr.eq.1) write (6,501) cnpa,roots(j)
0684   501       format(1x,"cnpa, nr-root = ",3g18.10)
0685             call wave_polar(cnpa,roots(j),0)
0686    21     continue
0687           call nr_disp_herm(cnpa,roots,nrt_nr)
0688 C
0689           if (knr.eq.1) go to 1
0690 C
0691 C........ solving the fully relativistic dispersion relation .........
0692 C
0693           if (i.eq.1) then
0694             do 31 j=1,nrt_fr
0695               frtgues(j)=rguess(j)
0696    31        continue
0697           end if
0698 CC          frtgues(1)=dcmplx(93.76,-0.56e-2)
0699           call fr_disp(frtgues,frt,nrt_fr)
0700           call zero_roots(frt,nrt_fr)
0701 C          write (8,108) cnpa,frt(1:nrt_fr)
0702           do 41 j=1,nrt_fr
0703             frtgues(j)=frt(j)
0704             if (i.eq.1.or.kpr.eq.1) write (6,601) cnpa,frt(j)
0705   601       format(1x,"cnpa, rel-root = ",3g18.10)
0706             call wave_polar(cnpa,frt(j),1)
0707           if (i.eq.1) then
0708             frt_herm_guess(j)=dcmplx(dreal(frt(j)),0.d0)
0709           else
0710             frt_herm_guess(j)=frt_herm(j)
0711           end if
0712    41     continue
0713           call fr_herm_disp(cnpa,frt_herm_guess,frt_herm,nrt_fr)
0714 C
0715     1   continue
0716 C
0717       end if
0718 C
0719 C
0720       if (kparam.eq.2) then
0721 C
0722          cnpa=cnpa_in
0723          alpha=alpha_in
0724 C
0725         if (n_omega.eq.1) then
0726           d_omega=0.d0
0727         else
0728           d_omega=(omega_fi-omega_in)/dfloat(n_omega-1)
0729         end if
0730 C
0731         omega=omega_in-d_omega
0732 C
0733         do 2 i=1,n_omega
0734 C
0735           omega=omega+d_omega
0736 C
0737           write (6,302) omega
0738   302     format(1x,"omega/omega_c = ",g18.10)
0739 C
0740 C..... solving for the cold plasma and approximate plasma roots ......
0741 C
0742           if (i.eq.1) then
0743             call coldrt(croots)
0744             call zero_roots(croots,4)
0745             call app_nr(rguess)
0746             call zero_roots(rguess,6)
0747             write (6,401) croots(1:2)
0748             write (6,401) croots(3:4)
0749             write (6,402) rguess(1:2)
0750             write (6,402) rguess(3:4)
0751             write (6,402) rguess(5:6)
0752           end if
0753 C
0754 C........... setting up the initial guesses if so desired ............
0755 C
0756       if (i.eq.1.and.kgues.eq.1) then
0757         do 62 ig=1,nrt_nr
0758           rguess(ig)=dcmplx(guesr(ig),guesi(ig))
0759    62   continue
0760       end if
0761 C
0762 C......... solving the non-relativistic dispersion relation ..........
0763 C
0764           call nr_disp(rguess,roots,nrt_nr)
0765           call zero_roots(roots,nrt_nr)
0766 C          write (5,107) omega,roots(1:nrt_nr)
0767           do 22 j=1,nrt_nr
0768             rguess(j)=roots(j)
0769             if (i.eq.1.or.kpr.eq.1) write (6,502) omega,roots(j)
0770   502       format(1x,"omega, nr-root = ",3g18.10)
0771             call wave_polar(omega,roots(j),0)
0772    22     continue
0773           call nr_disp_herm(omega,roots,nrt_nr)
0774 C
0775           if (knr.eq.1) go to 2
0776 C
0777 C........ solving the fully relativistic dispersion relation .........
0778 C
0779           if (i.eq.1) then
0780             do 32 j=1,nrt_fr
0781               frtgues(j)=rguess(j)
0782    32       continue
0783           end if
0784 CC          frtgues(1)=dcmplx(93.76,-0.56e-2)
0785           call fr_disp(frtgues,frt,nrt_fr)
0786           call zero_roots(frt,nrt_fr)
0787 C          write (8,108) omega,frt(1:nrt_fr)
0788           do 42 j=1,nrt_fr
0789             frtgues(j)=frt(j)
0790             if (i.eq.1.or.kpr.eq.1) write (6,602) omega,frt(j)
0791   602       format(1x,"cnpa, rel-root = ",3g18.10)
0792             call wave_polar(omega,frt(j),1)
0793           if (i.eq.1) then
0794             frt_herm_guess(j)=dcmplx(dreal(frt(j)),0.d0)
0795           else
0796             frt_herm_guess(j)=frt_herm(j)
0797           end if
0798    42     continue
0799           call fr_herm_disp(omega,frt_herm_guess,frt_herm,nrt_fr)
0800 C
0801     2   continue
0802 C
0803       end if
0804 C
0805 C
0806       if (kparam.eq.3) then
0807 C
0808          cnpa=cnpa_in
0809          omega=omega_in
0810 C
0811         if (n_alpha.eq.1) then
0812           d_alpha=0.d0
0813         else
0814           d_alpha=(alpha_fi-alpha_in)/dfloat(n_alpha-1)
0815         end if
0816 C
0817         alpha=alpha_in-d_alpha
0818 C
0819         do 3 i=1,n_alpha
0820 C
0821           alpha=alpha+d_alpha
0822 C
0823           write (6,302) alpha
0824   303     format(1x,"alpha = ",g18.10)
0825 C
0826 C..... solving for the cold plasma and approximate plasma roots ......
0827 C
0828           if (i.eq.1) then
0829             call coldrt(croots)
0830             call zero_roots(croots,4)
0831             call app_nr(rguess)
0832             call zero_roots(rguess,6)
0833             write (6,401) croots(1:2)
0834             write (6,401) croots(3:4)
0835             write (6,402) rguess(1:2)
0836             write (6,402) rguess(3:4)
0837             write (6,402) rguess(5:6)
0838           end if
0839 C
0840 C........... setting up the initial guesses if so desired ............
0841 C
0842       if (i.eq.1.and.kgues.eq.1) then
0843         do 63 ig=1,nrt_nr
0844           rguess(ig)=dcmplx(guesr(ig),guesi(ig))
0845    63   continue
0846       end if
0847 C
0848 C......... solving the non-relativistic dispersion relation ..........
0849 C
0850           call nr_disp(rguess,roots,nrt_nr)
0851           call zero_roots(roots,nrt_nr)
0852 C          write (5,107) alpha,roots(1:nrt_nr)
0853           do 23 j=1,nrt_nr
0854             rguess(j)=roots(j)
0855             if (i.eq.1.or.kpr.eq.1) write (6,503) alpha,roots(j)
0856   503       format(1x,"alpha, nr-root = ",3g18.10)
0857             call wave_polar(alpha,roots(j),0)
0858    23    continue
0859           call nr_disp_herm(alpha,roots,nrt_nr)
0860 C
0861           if (knr.eq.1) go to 3
0862 C
0863 C........ solving the fully relativistic dispersion relation .........
0864 C
0865           if (i.eq.1) then
0866             do 33 j=1,nrt_fr
0867               frtgues(j)=rguess(j)
0868    33      continue
0869           end if
0870 CC          frtgues(1)=dcmplx(93.76,-0.56e-2)
0871           call fr_disp(frtgues,frt,nrt_fr)
0872           call zero_roots(frt,nrt_fr)
0873 C          write (8,108) alpha,frt(1:nrt_fr)
0874           do 43 j=1,nrt_fr
0875             frtgues(j)=frt(j)
0876             if (i.eq.1.or.kpr.eq.1) write (6,603) alpha,frt(j)
0877   603       format(1x,"alpha, rel-root = ",3g18.10)
0878             call wave_polar(alpha,frt(j),1)
0879           if (i.eq.1) then
0880             frt_herm_guess(j)=dcmplx(dreal(frt(j)),0.d0)
0881           else
0882             frt_herm_guess(j)=frt_herm(j)
0883           end if
0884    43     continue
0885           call fr_herm_disp(alpha,frt_herm_guess,frt_herm,nrt_fr)
0886 C
0887     3   continue
0888 C
0889       end if
0890 C
0891 C
0892   101 format("om, rel-eps = ",3g18.10)
0893   103 format("------------------------------------------------------")
0894   104 format("om, nr-eps  = ",3g18.10)
0895   105 format("om, nr-sin  = ",3g18.10)
0896   106 format("om, rel-sin = ",3g18.10)
0897   107 format(<2*nrt_nr+1>g18.10)
0898   108 format(<2*nrt_fr+1>g18.10)
0899   201 format("sum/d11,d12 = ",4g18.10)
0900   202 format("sum/d13,d22 = ",4g18.10)
0901   203 format("sum/d23,d33 = ",4g18.10)
0902   204 format("sum/x11,x12 = ",4g18.10)
0903   205 format("sum/x13,x22 = ",4g18.10)
0904   206 format("sum/x23,x33 = ",4g18.10)
0905 C
0906 C
0907       return
0908       end
0909 C
0910 C=====================================================================
0911 C
0912 C
0913 C---------------------------------------------------------------------
0914 C
0915 C..... initial data for density, magnetic field, and temperature .....
0916 C
0917       subroutine prof_in
0918 C
0919       implicit double precision (a-h, o-z)
0920 C
0921       common /para1/ w0
0922       common /para4/ den0,bt0,te0
0923       common /para5/ wpe0,wce0,vte0
0924 C
0925       pi=3.141592653589793238462643d0
0926 C
0927 C...... electron plasma and cycloton frequencies at the center .......
0928 C........................ (normalised to w0) .........................
0929 C
0930       wpe0=5.6414187d4/w0*dsqrt(den0)
0931       wce0=1.7588103d7/w0*bt0
0932 C
0933 C......... thermal velocity at the center (normalised to c) ..........
0934 C
0935       vte0=4.19396657d7/2.99792458d10*dsqrt(te0)
0936 C
0937 C
0938       return
0939       end
0940 C
0941 C---------------------------------------------------------------------
0942 C
0943 C...... calculating the various plasma parameters at position r ......
0944 C
0945       subroutine pp(r,kout,pos,ipos)
0946 C
0947       implicit double precision (a-h, o-z)
0948 C
0949       real*8 pos(32768*5)
0950 C
0951       common /para1/ w0
0952       common /para4/ den0,bt0,te0
0953       common /para5/ wpe0,wce0,vte0
0954       common /para8/ wpe,wce,vte
0955       common /para9/ alpha,omega,vt,besk2
0956       common /wave0/ cnpa
0957 
0958 C
0959       external dbsk0e,dbsk1e,besk_a,umach,erset
0960 C
0961 C... distances normalised to plasma radius .....
0962 C.............. profiles for density and magnetic field ..............
0963 C.............. (normalised to their central values) .................
0964 C................ extending into the scrape-off layer ................
0965 C... (in this region the density/temperatures are the edge values) ...
0966 C
0967 C   read(20,*) r, te, den, bt, cnpa
0968     r = pos(ipos*5 + 1)
0969     te = pos(ipos*5 + 2)
0970     den = pos(ipos*5 + 3)
0971     bt = pos(ipos*5 + 4)
0972     cnpa = pos(ipos*5 + 5)
0973 C   read(21,*) rcm, tem, denm, qm
0974 C   read(21,*) rcp, tep, denp, qp
0975 
0976 C   te = ((rcp - r)*tem + (r - rcm)*tep)/(rcp - rcm)
0977 C   den = ((rcp - r)*denm + (r - rcm)*denp)/(rcp - rcm)
0978 C   q = ((rcp - r)*qm + (r - rcm)*qp)/(rcp - rcm)
0979 
0980 
0981 
0982 C
0983 C      bpol=r/r0/q/(1.+r/r0*cos(theta))
0984 C.......Modified j.d. 9/16/03. Poloidal dependence added for B field 
0985 C      btor=1./(1.+r/r0*cos(theta))
0986 
0987 C      bt=dsqrt(bpol**2+btor**2)
0988 
0989 C
0990 C....... various frequencies normalised to the wave frequency ........
0991 C
0992       wpe=wpe0*dsqrt(den)
0993       wce=wce0*bt
0994 C
0995 C........... alpha is the overdense parameter (wpe/wce)**2 ...........
0996 C
0997       alpha=(wpe/wce)**2
0998 C
0999       omega=1.d0/wce
1000 C
1001 C............. electron thermal velocity normalized to c .............
1002 C
1003       vte=vte0*dsqrt(te)
1004       vt=vte
1005 C
1006 C.... this is K_2(c**2/vt**2) needed for relativistic Maxwellians ....
1007 C... the modified bessel function in the denominator of integrals ....
1008 C...... (this is multiplied by exp(x) where x is the argument) .......
1009 C
1010       arg=1./vte**2
1011       if (arg.lt.50.d0) then
1012         call umach(-3,9)
1013         besk2=dbsk0e(arg)+2.*dbsk1e(arg)/arg
1014         call erset(0,1,1)
1015       else
1016         besk2=besk_a(2,arg)
1017       end if
1018 C
1019       if (kout.eq.1) then
1020         write (7,101) r,wpe,wce,dsqrt(wpe**2+wce**2),
1021      %                1./wce,vte,bt*bt0,den*den0,te*te0
1022       end if
1023 C
1024 C
1025   101 format(9g18.10)
1026 C
1027       return
1028       end
1029 C
1030 C---------------------------------------------------------------------
1031 C
1032 C... aymptotic form for the modified Bessel function (second kind) ...
1033 C.......... this is for real argument and is exp(x)*K(nu,x) ..........
1034 C
1035       double precision function besk_a(nu,z)
1036 C
1037       implicit double precision (a-h,o-z)
1038 C
1039       double precision mu
1040 C
1041 C.......................... fac=sqrt(pi/2) ...........................
1042 C
1043       fac=1.253314137315500251207883
1044 C
1045       mu=dfloat(4*nu**2)
1046       besk_a=fac/dsqrt(z)*(1.+(mu-1.)/(8.*z)+
1047      %                        (mu-1.)*(mu-9.)/(128.*z**2)+
1048      %                        (mu-1.)*(mu-9.)*(mu-25.)/(3072.*z**3))
1049 C
1050 C
1051       return
1052       end
1053 C
1054 C---------------------------------------------------------------------
1055 C
1056 C... aymptotic form for the modified Bessel function (second kind) ...
1057 C........ this is for complex argument and is exp(x)*K(nu,x) .........
1058 C
1059       double complex function cbesk_a(nu,z)
1060 C
1061       implicit double precision (a-h,o-z)
1062 C
1063       double precision mu
1064 C
1065       double complex z,cdsqrt
1066 C
1067 C.......................... fac=sqrt(pi/2) ...........................
1068 C
1069       fac=1.253314137315500251207883
1070 C
1071       mu=dfloat(4*nu**2)
1072       cbesk_a=fac/cdsqrt(z)*(1.+(mu-1.)/(8.*z)+
1073      %                          (mu-1.)*(mu-9.)/(128.*z**2)+
1074      %                          (mu-1.)*(mu-9.)*(mu-25.)/(3072.*z**3))
1075 C
1076 C
1077       return
1078       end
1079 C
1080 C---------------------------------------------------------------------
1081 C---------------------------------------------------------------------
1082 C
1083 C........ solving the fully relativistic dispersion relation .........
1084 C
1085       subroutine fr_disp(rguess,roots,nrts)
1086 C
1087       implicit double precision (a-h,o-z)
1088 C
1089       parameter (n_rts=10)
1090 C
1091       integer info(n_rts)
1092 C
1093       double complex fr_func
1094       double complex rguess(n_rts),roots(n_rts)
1095 C
1096       common /rterr/ errt_abs,errt_rel,it_max
1097       common /prnt/  kpr
1098 C
1099       external fr_func,dzanly,umach,erset
1100 C
1101       errabs=errt_abs
1102       errrel=errt_rel
1103       itmax=it_max
1104       nknown=0
1105       nnew=nrts
1106       nguess=nrts
1107 C
1108       call umach(-3,9)
1109       call dzanly(fr_func,errabs,errrel,nknown,nnew,nguess,
1110      %            rguess,itmax,roots,info)
1111       call erset(0,1,1)
1112 C
1113       call order_roots(roots,nrts)
1114 C
1115 CC      if (kpr.eq.1) then
1116 CC        write (6,101) roots(1:2)
1117 CC        write (6,101) roots(3:4)
1118 CC        write (6,101) roots(5:6)
1119 CC        write (6,101) roots(7:8)
1120 CC      end if
1121   101 format(1x,"full r-r = ",4g18.10)
1122 C
1123 C
1124       return
1125       end
1126 C
1127 C---------------------------------------------------------------------
1128 C
1129 C............ the fully-relativistic dispersion function .............
1130 C
1131       double complex function fr_func(roots)
1132 C
1133       implicit double precision (a-h, o-z)
1134 C
1135       double complex roots,cnpe
1136       double complex d11,d12,d13,d22,d23,d33
1137       double complex x11,x12,x13,x22,x23,x33
1138 C
1139       double complex ep(6)
1140 C
1141       common /trub0/ ktrub
1142       common /wave0/ cnpa
1143       common /wave1/ cnpe
1144       common /coutd/ d11,d12,d13,d22,d23,d33
1145 C
1146       cnpe=roots
1147 C
1148       if (ktrub.eq.1) then
1149         call trub(ep)
1150       else
1151         call weiss(ep)
1152       end if
1153 C
1154       x11=ep(1)
1155       x12=ep(2)
1156       x13=ep(3)
1157       x22=ep(4)
1158       x23=ep(5)
1159       x33=ep(6)
1160 C
1161 C........ the dielectric tensor elements (satisfying D.E = 0) ........
1162 C
1163       d11=1.-cnpa**2+x11
1164       d12=x12
1165       d13=cnpe*cnpa+x13
1166       d22=1.-cnpe**2-cnpa**2+x22
1167       d23=x23
1168       d33=1.-cnpe**2+x33
1169 C
1170       fr_func=d11*(d22*d33+d23*d23)+
1171      %        d12*(d23*d13+d12*d33)+
1172      %        d13*(d12*d23-d22*d13)
1173 C
1174 C
1175       return
1176       end
1177 C
1178 C---------------------------------------------------------------------
1179 C---------------------------------------------------------------------
1180 C
1181 C......... solving the non-relativistic dispersion relation ..........
1182 C
1183       subroutine nr_disp(rguess,roots,nrts)
1184 C
1185       implicit double precision (a-h, o-z)
1186 C
1187       parameter (n_rts=10)
1188 C
1189       integer info(n_rts)
1190 C
1191       double complex nr_func
1192       double complex rguess(n_rts),roots(n_rts)
1193 C
1194       common /rterr/ errt_abs,errt_rel,it_max
1195       common /prnt/  kpr
1196 C
1197       external nr_func,dzanly,umach,erset
1198 C
1199       errabs=errt_abs
1200       errrel=errt_rel
1201       itmax=it_max
1202       nknown=0
1203       nnew=nrts
1204       nguess=nrts
1205 C
1206       call umach(-3,9)
1207       call dzanly(nr_func,errabs,errrel,nknown,nnew,nguess,
1208      %            rguess,itmax,roots,info)
1209       call erset(0,1,1)
1210 C
1211       call order_roots(roots,nrts)
1212 C
1213 CC      if (kpr.eq.1) then
1214 CC        write (6,101) roots(1:2)
1215 CC        write (6,101) roots(3:4)
1216 CC        write (6,101) roots(5:6)
1217 CC        write (6,101) roots(7:8)
1218 CC      end if
1219 C
1220   101 format(1x,"full nr-rts = ",4g18.10)
1221 C
1222 C
1223       return
1224       end
1225 C
1226 C---------------------------------------------------------------------
1227 C
1228 C............. the non-relativistic dispersion function ..............
1229 C
1230       double complex function nr_func(roots)
1231 C
1232       implicit double precision (a-h, o-z)
1233 C
1234       double complex roots
1235       double complex d11,d12,d13,d22,d23,d33
1236 C
1237       common /nrels/ nr_nterms
1238       common /coutd/ d11,d12,d13,d22,d23,d33
1239 C
1240       call d_nr(nr_nterms,roots)
1241 C
1242       nr_func=d11*(d22*d33+d23*d23)+
1243      %        d12*(d23*d13+d12*d33)+
1244      %        d13*(d12*d23-d22*d13)
1245 C
1246 C
1247       return
1248       end
1249 C
1250 C---------------------------------------------------------------------
1251 C
1252 C................. non-relativistic tensor elements ..................
1253 C
1254       subroutine d_nr(nn,cnpe)
1255 C
1256       implicit double precision (a-h, o-z)
1257 C
1258       double complex dcmplx,cnpe
1259       double complex ci,scmu,cmu
1260       double complex x11,x12,x13,x22,x23,x33
1261       double complex d11,d12,d13,d22,d23,d33
1262 C
1263       double complex cbi(0:1000),dcbi(0:1000)
1264       double complex zm(0:1000),zp(0:1000)
1265       double complex zargm(0:1000),zargp(0:1000)
1266       double complex csum(6)
1267 C
1268       common /para9/ alpha,omega,vt,besk2
1269       common /coutd/ d11,d12,d13,d22,d23,d33
1270       common /coutx/ x11,x12,x13,x22,x23,x33
1271       common /wave0/ cnpa
1272       common /prnt/  kpr
1273 C
1274       ci=dcmplx(0.d0,1.d0)
1275 C
1276 C.. scmu is the sqrt (cmu) (defined to take care of sign of k_perp) ..
1277 C.............. cmu = (k_perp * v_te)**2 / omega_ce**2 ...............
1278 C............ cy0 = omega / ( sqrt(2)*abs(k_para)*v_te ) .............
1279 C
1280       scmu=cnpe*vt*omega
1281       cmu=scmu**2
1282       cy0=1./(dsqrt(2.d0)*dabs(cnpa)*vt)
1283 C
1284 C...... the number of terms to be calculated in the summations .......
1285 C....... this is based on when exp(-mu)*I_n(mu) becomes small ........
1286 C
1287       nn=idint(cdabs(cmu)/dsqrt(2.d0))*5+20
1288 C
1289 C.......... arguments of the z-functions in the summations ...........
1290 C
1291       do 1 n=0,nn
1292         zargm(n)=dcmplx(cy0*(1.-dfloat(n)/omega),0.d0)
1293         zargp(n)=dcmplx(cy0*(1.+dfloat(n)/omega),0.d0)
1294     1 continue
1295 C
1296       call ebes(nn,cmu,cbi,dcbi)
1297 C
1298       call pdisp(nn,zargm,zargp,zm,zp)
1299 C
1300       csum(1)=dcmplx(0.d0,0.d0)
1301       csum(2)=dcmplx(0.d0,0.d0)
1302       csum(3)=dcmplx(0.d0,0.d0)
1303       csum(4)=2.*(cbi(0)-dcbi(0))*zp(0)
1304       csum(5)=zargp(0)*zp(0)*(dcbi(0)-cbi(0))
1305       csum(6)=zargp(0)*(1.d0+cbi(0)*zargp(0)*zp(0))
1306 C
1307       do 2 n=1,nn
1308 C
1309         csum(1)=csum(1)+
1310      %          dfloat(n**2)*cbi(n)*(zp(n)+zm(n))
1311         csum(2)=csum(2)+
1312      %          dfloat(n)*(dcbi(n)-cbi(n))*(zm(n)-zp(n))
1313         csum(3)=csum(3)+
1314      %          dfloat(n)*cbi(n)*(zargm(n)*zm(n)-
1315      %                            zargp(n)*zp(n))
1316         csum(4)=csum(4)+
1317      %          ((2.+dfloat(n**2)/cmu**2)*cbi(n)-2.*dcbi(n))*
1318      %          (zp(n)+zm(n))
1319         csum(5)=csum(5)+
1320      %          (dcbi(n)-cbi(n))*(zargp(n)*zp(n)+zargm(n)*zm(n))
1321         csum(6)=csum(6)+
1322      %          cbi(n)*(zargp(n)**2*zp(n)+zargm(n)**2*zm(n))
1323 C
1324     2 continue
1325 C
1326 C.......................... sign of k_para ...........................
1327 C
1328       sk=cnpa/(dabs(cnpa))
1329 C
1330       fac=alpha/omega**2
1331 C
1332       x11=fac*cy0/cmu*csum(1)
1333       x12=-ci*fac*cy0*csum(2)
1334       x13=sk*fac*dsqrt(2.d0)/scmu*cy0*csum(3)
1335       x22=fac*cmu*cy0*csum(4)
1336       x23=ci*sk*fac*dsqrt(2.d0)*scmu*cy0*csum(5)
1337       x33=2.*fac*cy0*csum(6)
1338 C
1339       if (kpr.eq.1) then
1340         write (6,601) x11,x12
1341         write (6,601) x13,x22
1342         write (6,601) x23,x33
1343   601   format("non-rel x's = ",4g18.10)
1344       end if
1345 C
1346 C........ the dielectric tensor elements (satisfying D.E = 0) ........
1347 C
1348       d11=1.-cnpa**2+x11
1349       d12=x12
1350       d13=cnpe*cnpa+x13
1351       d22=1.-cnpe**2-cnpa**2+x22
1352       d23=x23
1353       d33=1.-cnpe**2+x33
1354 C
1355 C
1356       return
1357       end
1358 C
1359 C---------------------------------------------------------------------
1360 C
1361 C....... evaluating the Bessel functions and their derivatives .......
1362 C
1363       subroutine ebes(n,carg,cbi,dcbi)
1364 C
1365       implicit double precision (a-h, o-z)
1366 C
1367       double complex cdexp
1368       double complex carg,factor
1369 C
1370       double complex cbi(0:1000),dcbi(0:1000)
1371 C
1372       external dcbis
1373 C
1374       xnu=0.
1375       numb=n+2
1376 C
1377 C............... evaluating the bessel function itself ...............
1378 C
1379       call dcbis(xnu,carg,numb,cbi)
1380 C
1381 C......... calculating the derivates of the bessel functions .........
1382 C
1383       dcbi(0)=cbi(1)
1384       do 1 i=1,numb-1
1385         dcbi(i)=0.5*(cbi(i-1)+cbi(i+1))
1386     1 continue
1387 C
1388 C.............. multiplying by exp(-mu_e) for the sums ...............
1389 C
1390       factor=cdexp(-carg)
1391       do 2 i=0,numb-1
1392         cbi(i)=cbi(i)*factor
1393         dcbi(i)=dcbi(i)*factor
1394     2 continue
1395 C
1396 C
1397       return
1398       end
1399 C
1400 C---------------------------------------------------------------------
1401 C
1402 C............ calculating the plasma dispersion functions ............
1403 C
1404       subroutine pdisp(nar,cargm,cargp,zm,zp)
1405 C
1406       implicit double precision (a-h, o-z)
1407 C
1408       double complex dcmplx
1409       double complex carg,csqpi,zerfe
1410 C
1411       double complex cargm(0:1000),cargp(0:1000)
1412       double complex zm(0:1000),zp(0:1000)
1413 C
1414       external zerfe
1415 C
1416       sqpi=1.772453850905516027298167d0
1417         pi=3.141592653589793238462643d0
1418 C
1419       csqpi=dcmplx(0.d0,sqpi)
1420 C     
1421       do 1 i=0,nar
1422         zm(i)=csqpi*zerfe(cargm(i))
1423         zp(i)=csqpi*zerfe(cargp(i))
1424     1 continue
1425 C
1426 C
1427         return
1428         end
1429 C
1430 C---------------------------------------------------------------------
1431 C---------------------------------------------------------------------
1432 C
1433 C................. calculating the cold plasma roots .................
1434 C
1435       subroutine coldrt(croots)
1436 C
1437       implicit double precision (a-h, o-z)
1438 C
1439       double complex dcmplx,cdsqrt,sc,s2,cdummy
1440       double complex cnpe
1441 C
1442       double complex cr(2),croots(4)
1443 C
1444       common /para9/ alpha,omega,vt,besk2
1445       common /wave0/ cnpa
1446       common /wave1/ cnpe
1447 C
1448 C..................... the stix tensor elements ......................
1449 C
1450       sx=1.-alpha/(omega**2-1.)
1451       sy=-alpha/(omega*(omega**2-1.))
1452       sz=1.-alpha/omega**2
1453 C
1454       s1=0.5*(sx+sz-cnpa**2)-0.5*(cnpa**2*sz+sy**2)/sx
1455       sdummy=((cnpa**2-sx)*(sx-sz)+sy**2)**2+
1456      %         4.*cnpa**2*sy**2*sz
1457       sc=dcmplx(sdummy,0.d0)
1458       s2=0.5*cdsqrt(sc)/sx
1459 C
1460       cr(1)=s1+s2
1461       cr(2)=s1-s2
1462 C
1463 C...... this orders the roots in terms of their absolute values ......
1464 C
1465       if (cdabs(cr(1)).gt.cdabs(cr(2))) then
1466         cdummy=cr(1)
1467         cr(1)=cr(2)
1468         cr(2)=cdummy
1469       end if
1470 C
1471 C.................. and now all the roots of n-perp ..................
1472 C... roots(1) and (2) are for fast wave; (3) and (4) for slow wave ...
1473 C
1474       croots(1)= cdsqrt(cr(1))
1475       croots(2)=-croots(1)
1476       croots(3)= cdsqrt(cr(2))
1477       croots(4)=-croots(3)
1478 C
1479       return
1480       end
1481 C
1482 C---------------------------------------------------------------------
1483 C
1484 C......... approximate non-relativistic dispersion equation ..........
1485 C
1486       subroutine app_nr(roots)
1487 C
1488       implicit double precision (a-h, o-z)
1489 C
1490 C......... this is the order of the polynomial to be solved ..........
1491 C....... (should be set to at least 3 and can be 4 if desired) .......
1492 C
1493       parameter (n_order=3)
1494 C
1495       double complex dcmplx,ci,cdsqrt,zero
1496       double complex cnpe
1497 C
1498       double complex zm(0:1000),zp(0:1000)
1499       double complex zargm(0:1000),zargp(0:1000)
1500       double complex xi(0:5),be(0:5)
1501       double complex ga(0:5),de(0:5),ph(0:5)
1502       double complex d11(0:8),d12(0:8),d13(0:8)
1503       double complex d22(0:8),d23(0:8),d33(0:8)
1504       double complex x11(0:8),x12(0:8),x13(0:8)
1505       double complex x22(0:8),x23(0:8),x33(0:8)
1506       double complex co(0:4),coo(0:n_order)
1507       double complex rts(n_order),roots(8)
1508 C
1509       common /para9/ alpha,omega,vt,besk2
1510       common /wave0/ cnpa
1511       common /wave1/ cnpe
1512 C
1513       external dzpocc,umach,erset
1514 C
1515       ci=dcmplx(0.d0,1.d0)
1516 C
1517 C.. scmu is the sqrt (cmu) (defined to take care of sign of k_perp) ..
1518 C............ cmu = (v_te / c)**2 / (omega_ce / omega)**2 ............
1519 C............ cy0 = omega / ( sqrt(2)*abs(k_para)*v_te ) .............
1520 C
1521       scmu=vt*omega
1522       cmu=scmu**2
1523       cy0=1./(dsqrt(2.d0)*dabs(cnpa)*vt)
1524 C
1525 C.......................... sign of k_para ...........................
1526 C
1527       sk=cnpa/(dabs(cnpa))
1528 C
1529 C.......... arguments of the z-functions in the summations ...........
1530 C
1531       do 1 n=0,6
1532         zargm(n)=dcmplx(cy0*(1.-dfloat(n)/omega),0.d0)
1533         zargp(n)=dcmplx(cy0*(1.+dfloat(n)/omega),0.d0)
1534     1 continue
1535 C
1536       nn=7
1537       call pdisp(nn,zargm,zargp,zm,zp)
1538 C
1539       do 2 n=0,5
1540         be(n)=zm(n)-zp(n)
1541         ga(n)=zargm(n)*zm(n)-zargp(n)*zp(n)
1542         if (n.eq.0) then
1543           xi(0)=zm(0)
1544           de(0)=zargm(0)*zm(0)
1545           ph(0)=zargm(0)**2*zm(0)
1546         else
1547           xi(n)=zm(n)+zp(n)
1548           de(n)=zargm(n)*zm(n)+zargp(n)*zp(n)
1549           ph(n)=zargm(n)**2*zm(n)+zargp(n)**2*zp(n)
1550         end if
1551     2 continue
1552 C
1553       zero=dcmplx(0.d0,0.d0)
1554       do 3 i=0,8
1555         d11(i)=zero
1556         d12(i)=zero
1557         d13(i)=zero
1558         d22(i)=zero
1559         d23(i)=zero
1560         d33(i)=zero
1561         x11(i)=zero
1562         x12(i)=zero
1563         x13(i)=zero
1564         x22(i)=zero
1565         x23(i)=zero
1566         x33(i)=zero
1567     3 continue
1568 C
1569       d11(0)=1.-cnpa**2
1570       d13(1)=cnpa
1571       d22(0)=1.-cnpa**2
1572       d22(2)=-1.
1573       d33(0)=1.
1574       d33(2)=-1.
1575 C
1576       x11(0)=xi(1)
1577       x11(2)=cmu*(-xi(1)+xi(2))
1578       x11(4)=cmu**2*(0.625*xi(1)-xi(2)+0.375*xi(3))
1579       x11(6)=cmu**3/24.*(-7.*xi(1)+14.*xi(2)-
1580      %                    9.*xi(3)+2.*xi(4))
1581       x11(8)=cmu**4/384.*(42.*xi(1)-96.*xi(2)+81.*xi(3)-
1582      %                     32.*xi(4)+5.*xi(5))
1583 C
1584 C........ this is i*D_12 (or D_12-bar) as set up in the notes ........
1585 C
1586       x12(0)=be(1)
1587       x12(2)=cmu*(-2.*be(1)+be(2))
1588       x12(4)=cmu**2*(1.875*be(1)-1.5*be(2)+0.375*be(3))
1589       x12(6)=cmu**3/12.*(-14.*be(1)+14.*be(2)-6.*be(3)+be(4))
1590       x12(8)=cmu**4/384.*(210.*be(1)-240.*be(2)+125.*be(3)-
1591      %                        40.*be(4)+5.*be(5))
1592 C
1593       cdum=sk*scmu*dsqrt(2.d0)
1594       x13(1)=cdum*ga(1)
1595       x13(3)=cdum*cmu*(-ga(1)+0.5*ga(2))
1596       x13(5)=cdum*cmu**2*0.125*(5.*ga(1)-4.*ga(2)+ga(3))
1597       x13(7)=cdum*cmu**3/48.*(-14.*ga(1)+14.*ga(2)-
1598      %                         6.*ga(3)+ga(4))
1599 C
1600       x22(0)=xi(1)
1601       x22(2)=cmu*(4.*xi(0)-3.*xi(1)+xi(2))
1602       x22(4)=cmu**2*(-6.*xi(0)+4.625*xi(1)-
1603      %                2.*xi(2)+0.375*xi(3))
1604       x22(6)=cmu**3*(120.*xi(0)-97.*xi(1)+50.*xi(2)-
1605      %               15.*xi(3)+2.*xi(4))
1606       x22(8)=cmu**4/384.*(-560.*xi(0)+938.*xi(1)-544.*xi(2)+
1607      %                     209.*xi(3)-48.*xi(4)+5.*xi(5))
1608 C
1609 C
1610 C........ this is i*D_23 (or D_23-bar) as set up in the notes ........
1611 C
1612       cdum=-sk*dsqrt(2.d0)*scmu
1613       x23(1)=cdum*(-2.*de(0)+de(1))
1614       x23(3)=cdum*cmu*(3.*de(0)-2.*de(1)+0.5*de(2))
1615       x23(5)=cdum*cmu**2*(-2.5*de(0)+1.875*de(1)-
1616      %                    0.75*de(2)+0.125*de(3))
1617       x23(7)=cdum*cmu**3/48.*(70.*de(0)-56.*de(1)+28.*de(2)-
1618      %                        8.*de(3)+de(4))
1619 C
1620       x33(0)=4.*(cy0+ph(0))
1621       x33(2)=cmu*2.*(-2.*ph(0)+ph(1))
1622       x33(4)=cmu**2*(3.*ph(0)-2.*ph(1)+0.5*ph(2))
1623       x33(6)=cmu**3/12.*(-20.*ph(0)+15.*ph(1)-
1624      %                    6.*ph(2)+ph(3))
1625       x33(8)=cmu**4/96.*(70.*ph(0)-56.*ph(1)+28.*ph(2)-
1626      %                    8.*ph(3)+ph(4))
1627 C
1628       fac=0.5*alpha/omega**2*cy0
1629 C
1630       do 4 i=0,8
1631         x11(i)=fac*x11(i)
1632         x12(i)=fac*x12(i)
1633         x13(i)=fac*x13(i)
1634         x22(i)=fac*x22(i)
1635         x23(i)=fac*x23(i)
1636         x33(i)=fac*x33(i)
1637     4 continue
1638 C
1639       do 5 i=0,8
1640         d11(i)=d11(i)+x11(i)
1641         d12(i)=d12(i)+x12(i)
1642         d13(i)=d13(i)+x13(i)
1643         d22(i)=d22(i)+x22(i)
1644         d23(i)=d23(i)+x23(i)
1645         d33(i)=d33(i)+x33(i)
1646     5 continue
1647 C
1648 C.... the coefficients of (n-perp)**2 in the dispersion function .....
1649 C........ disp-fn = co(0)+co(1)*(n-perp)**2+co(2)*(n-perp)**4 ........
1650 C.......             + co(3)*(n-perp)**6+co(4)*(n-perp)**8    ........
1651 C............ the co(i)'s were dertermined using macsyma .............
1652 C
1653       co(0) = d11(0)*d22(0)*d33(0)-d12(0)**2*d33(0)
1654 C
1655       co(1) = d11(0)*d22(0)*d33(2)-d12(0)**2*d33(2)+
1656      %        d11(0)*d33(0)*d22(2)-2.*d12(0)*d33(0)*d12(2)+
1657      %        d22(0)*d33(0)*d11(2)-d11(0)*d23(1)**2-
1658      %        2.*d12(0)*d13(1)*d23(1)-d22(0)*d13(1)**2
1659 C
1660       co(2) = d11(0)*d22(0)*d33(4)-d12(0)**2*d33(4)+
1661      %        d11(0)*d33(0)*d22(4)-2.*d12(0)*d33(0)*d12(4)+
1662      %        d22(0)*d33(0)*d11(4)-2.*d11(0)*d23(1)*d23(3)-
1663      %        2.*d12(0)*d13(1)*d23(3)-2.*d12(0)*d23(1)*d13(3)-
1664      %        2.*d22(0)*d13(1)*d13(3)+d11(0)*d22(2)*d33(2)-
1665      %        2.*d12(0)*d12(2)*d33(2)+d22(0)*d11(2)*d33(2)+
1666      %        d33(0)*d11(2)*d22(2)-d13(1)**2*d22(2)-
1667      %        d33(0)*d12(2)**2-2.*d13(1)*d23(1)*d12(2)-
1668      %        d23(1)**2*d11(2)
1669 C
1670       co(3) = d11(0)*d22(0)*d33(6)-d12(0)**2*d33(6)+
1671      %        d11(0)*d33(0)*d22(6)-2.*d12(0)*d33(0)*d12(6)+
1672      %        d22(0)*d33(0)*d11(6)-2.*d11(0)*d23(1)*d23(5)-
1673      %        2.*d12(0)*d13(1)*d23(5)-2.*d12(0)*d23(1)*d13(5)-
1674      %        2.*d22(0)*d13(1)*d13(5)+d11(0)*d22(2)*d33(4)-
1675      %        2.*d12(0)*d12(2)*d33(4)+d22(0)*d11(2)*d33(4)+
1676      %        d11(0)*d33(2)*d22(4)+d33(0)*d11(2)*d22(4)-
1677      %        d13(1)**2*d22(4)-2.*d12(0)*d33(2)*d12(4)-
1678      %        2.*d33(0)*d12(2)*d12(4)-2.*d13(1)*d23(1)*d12(4)+
1679      %        d22(0)*d33(2)*d11(4)+d33(0)*d22(2)*d11(4)-
1680      %        d23(1)**2*d11(4)-d11(0)*d23(3)**2-
1681      %        2.*d12(0)*d13(3)*d23(3)-2.*d13(1)*d12(2)*d23(3)-
1682      %        2.*d23(1)*d11(2)*d23(3)-d22(0)*d13(3)**2-
1683      %        2.*d13(1)*d22(2)*d13(3)-2.*d23(1)*d12(2)*d13(3)+
1684      %        d11(2)*d22(2)*d33(2)-d12(2)**2*d33(2)
1685 C
1686       co(4) = d11(0)*d22(0)*d33(8)-d12(0)**2*d33(8)+
1687      %        d11(0)*d33(0)*d22(8)-2.*d12(0)*d33(0)*d12(8)+
1688      %        d22(0)*d33(0)*d11(8)-2.*d11(0)*d23(1)*d23(7)-
1689      %        2.*d12(0)*d13(1)*d23(7)-2.*d12(0)*d23(1)*d13(7)-
1690      %        2.*d22(0)*d13(1)*d13(7)+d11(0)*d22(2)*d33(6)-
1691      %        2.*d12(0)*d12(2)*d33(6)+d22(0)*d11(2)*d33(6)+
1692      %        d11(0)*d33(2)*d22(6)+d33(0)*d11(2)*d22(6)-
1693      %        d13(1)**2*d22(6)-2.*d12(0)*d33(2)*d12(6)-
1694      %        2.*d33(0)*d12(2)*d12(6)-2.*d13(1)*d23(1)*d12(6)+
1695      %        d22(0)*d33(2)*d11(6)+d33(0)*d22(2)*d11(6)-
1696      %        d23(1)**2*d11(6)-2.*d11(0)*d23(3)*d23(5)-
1697      %        2.*d12(0)*d13(3)*d23(5)-2.*d13(1)*d12(2)*d23(5)-
1698      %        2.*d23(1)*d11(2)*d23(5)-2.*d12(0)*d23(3)*d13(5)-
1699      %        2.*d22(0)*d13(3)*d13(5)-2.*d13(1)*d22(2)*d13(5)-
1700      %        2.*d23(1)*d12(2)*d13(5)+d11(0)*d22(4)*d33(4)-
1701      %        2.*d12(0)*d12(4)*d33(4)+d22(0)*d11(4)*d33(4)+
1702      %        d11(2)*d22(2)*d33(4)-d12(2)**2*d33(4)+
1703      %        d33(0)*d11(4)*d22(4)-2.*d13(1)*d13(3)*d22(4)+
1704      %        d11(2)*d33(2)*d22(4)-d33(0)*d12(4)**2-
1705      %        2.*d13(1)*d23(3)*d12(4)-2.*d23(1)*d13(3)*d12(4)-
1706      %        2.*d12(2)*d33(2)*d12(4)-2.*d23(1)*d23(3)*d11(4)+
1707      %        d22(2)*d33(2)*d11(4)-d11(2)*d23(3)**2-
1708      %        2.*d12(2)*d13(3)*d23(3)-d22(2)*d13(3)**2
1709 C
1710       do 6 i=0,n_order
1711         coo(i)=co(i)
1712     6 continue
1713 C
1714       ndeg=n_order
1715       call umach(-3,9)
1716       call dzpocc(ndeg,coo,rts)
1717       call erset(0,1,1)
1718 C
1719       do 7 i=1,n_order
1720         roots(2*i-1)=cdsqrt(rts(i))
1721         roots(2*i)=-roots(2*i-1)
1722     7 continue
1723 C
1724 C
1725       return
1726       end
1727 C
1728 C---------------------------------------------------------------------
1729 C---------------------------------------------------------------------
1730 C
1731 C......... this is for zeroing out components that are small .........
1732 C
1733       subroutine zero_roots(roots,nrt0)
1734 C
1735       implicit double precision (a-h, o-z)
1736 C
1737       parameter (n_rts=10)
1738 C
1739       double complex dcmplx
1740 C
1741       double complex roots(n_rts)
1742 C
1743       do 1 i=1,nrt0
1744         if (dabs(dreal(roots(i))).lt.1.0d-15) 
1745      %      roots(i)=dcmplx(0.d0,dimag(roots(i)))
1746         if (dabs(dimag(roots(i))).lt.1.0d-15) 
1747      %      roots(i)=dcmplx(dreal(roots(i)),0.d0)
1748     1 continue
1749 C
1750 C
1751       return
1752       end
1753 C
1754 C---------------------------------------------------------------------
1755 C
1756 C........... calculating polarizations with respect to E_x ...........
1757 C
1758       subroutine wave_polar(data,cnpe,kfl_fr)
1759 C
1760       implicit double precision (a-h, o-z)
1761 C
1762       double complex dcmplx,ci,cnpe,cdummy
1763       double complex d11,d12,d13,d22,d23,d33
1764       double complex ey,ez
1765       double complex pol_p,pol_m,pol_z
1766 C
1767       double complex nr_func,fr_func
1768 C
1769       double complex det(3),pdum(3)
1770 C
1771       common /wave0/ cnpa
1772       common /coutd/ d11,d12,d13,d22,d23,d33
1773 C
1774       external nr_func,fr_func
1775 C
1776       ci=dcmplx(0.d0,1.d0)
1777 C
1778       if (kfl_fr.eq.0) cdummy=nr_func(cnpe)
1779       if (kfl_fr.eq.1) cdummy=fr_func(cnpe)
1780 C
1781       det(1)=-d12*d23-d13*d22
1782       det(2)= d22*d33-d23*d23
1783       det(3)= d12*d33+d13*d23
1784 C
1785       if (cdabs(det(1)).gt.cdabs(det(2))) then
1786         if (cdabs(det(1)).gt.cdabs(det(3))) then
1787           nn=1
1788         else
1789           nn=3
1790         end if
1791       else
1792         if (cdabs(det(2)).gt.cdabs(det(3))) then
1793           nn=2
1794         else
1795           nn=3
1796         end if
1797       end if
1798 C
1799       if (nn.eq.1) then
1800         ey=(d11*d23-d13*d12)/det(1)
1801         ez=(d12*d12+d11*d22)/det(1)
1802       end if
1803       if (nn.eq.2) then
1804         ey=(d12*d33-d13*d23)/det(2) 
1805         ez=(d12*d23-d13*d22)/det(2)
1806       end if
1807       if (nn.eq.3) then
1808         ey= (d13*d13-d11*d33)/det(3)
1809         ez=-(d12*d13+d11*d23)/det(3)
1810       end if
1811 C
1812 C....... zeroing out the small components of the polarizations .......
1813 C
1814       pdum(1)=(1.d0+ci*ey)/dsqrt(2.d0)
1815       pdum(2)=(1.d0-ci*ey)/dsqrt(2.d0)
1816       pdum(3)=ez
1817 C
1818       do 1 i=1,3
1819         if (dabs(dreal(pdum(i))).lt.1.0d-15) 
1820      %      pdum(i)=dcmplx(0.d0,dimag(pdum(i)))
1821         if (dabs(dimag(pdum(i))).lt.1.0d-15) 
1822      %      pdum(i)=dcmplx(dreal(pdum(i)),0.d0)
1823     1 continue
1824 C
1825       pol_p=pdum(1)
1826       pol_m=pdum(2)
1827       pol_z=pdum(3)
1828 C
1829       if (kfl_fr.eq.0) 
1830      %  write (4,101) data,cnpe,pol_p,pol_m,pol_z
1831       if (kfl_fr.eq.1) 
1832      %  write (2,101) data,cnpe,pol_p,pol_m,pol_z
1833 C
1834   101 format(9g18.10)
1835 C
1836 C
1837       return
1838       end
1839 C
1840 C---------------------------------------------------------------------
1841 C---------------------------------------------------------------------
1842 C
1843 C......  ordering the roots found from the dispersion relation .......
1844 C
1845       subroutine order_roots(roots,nrt0)
1846 C
1847       implicit double precision (a-h, o-z)
1848 C
1849       parameter (n_rts=10)
1850 C
1851       double complex cdummy,ctemp
1852 C
1853       double complex roots(n_rts)
1854 C
1855 C... arranging the real and imaginary parts in terms of magnitude ....
1856 C... values decrease from highest positive to the lowest negative ....
1857 C
1858       do 1 i=1,nrt0
1859         cdummy=roots(i)
1860         do 2 j=i+1,nrt0
1861           if (dreal(roots(j)).gt.dreal(cdummy)) then
1862             ctemp=roots(j)
1863             roots(j)=cdummy
1864             cdummy=ctemp
1865           end if
1866           if (dreal(roots(j)).eq.dreal(cdummy)) then
1867             if (dimag(roots(j)).eq.dimag(cdummy)) then
1868               ctemp=roots(j)
1869               roots(j)=cdummy
1870               cdummy=ctemp
1871             end if
1872           end if
1873     2   continue
1874         roots(i)=cdummy
1875     1 continue
1876 C
1877       call zero_roots(roots,nrt0)
1878 C
1879 C
1880       return
1881       end
1882 C
1883 C---------------------------------------------------------------------
1884 C---------------------------------------------------------------------
1885 C
1886 C..... evaluating the trubnikov form of the relativistic tensor ......
1887 C
1888       subroutine trub(epsilon)
1889 C
1890       implicit double precision (a-h, o-z)
1891 C
1892       double precision stot(3),sumr(0:1000)
1893 C
1894       double complex dcmplx,csum
1895 C
1896       double complex epsilon(6),co(0:1)
1897 C
1898       common /trub1/ xi_lo,dxi,nxi
1899       common /trub2/ diff_err,navg
1900       common /trerr/ errabs0,errrel0,irule0
1901       common /trubc/ kel,kreal
1902       common /prnt/  kpr
1903 C
1904       external dqdag,ftrub,umach,erset
1905 C
1906       co(0)=dcmplx(0.d0,1.d0)
1907       co(1)=dcmplx(1.d0,0.d0)
1908 C
1909       do 1 j=1,6
1910         kel=j
1911         csum=dcmplx(0.d0,0.d0)
1912 C
1913         do 2 k=0,1
1914           kreal=k
1915 C
1916 C.. ntr keeps track of the number of averages (3) that are compared ..
1917 C
1918           ntr=0
1919 C
1920           xi_0=xi_lo
1921           xi_1=xi_0+dxi
1922 C
1923           sum=0.d0
1924 C
1925           do 3 i=1,nxi
1926 C
1927             errabs=errabs0
1928             errrel=errrel0
1929             irule=irule0
1930 C
1931             call umach(-3,9)
1932             call dqdag(ftrub,xi_0,xi_1,errabs,errrel,irule,
1933      %                 t_real,errest)
1934             call erset(0,1,1)
1935 C
1936             sum=sum+t_real
1937 C
1938             jj=mod(i-1,navg)
1939             sumr(jj)=sum
1940             if (jj.eq.navg-1) then
1941               ntr=ntr+1
1942               stot(ntr)=0.d0
1943               do 4 ij=0,navg-1
1944                 stot(ntr)=stot(ntr)+sumr(ij)
1945     4         continue
1946               stot(ntr)=stot(ntr)/dfloat(navg)
1947               if (ntr.eq.3) then
1948                 diff1=dabs(dabs(stot(1))-dabs(stot(2)))
1949                 diff2=dabs(dabs(stot(2))-dabs(stot(3)))
1950                 if (diff1.lt.diff_err.and.diff2.lt.diff_err) then
1951                   sum_real=stot(3)
1952                   go to 5
1953                 else
1954                   stot(1)=stot(2)
1955                   stot(2)=stot(3)
1956                   ntr=2
1957                 end if
1958               end if
1959             end if
1960 C
1961             xi_0=xi_1
1962             xi_1=xi_0+dxi
1963 C
1964     3     continue
1965 C
1966     5     continue
1967 C
1968           if (kpr.eq.1) then
1969             if (kreal.eq.0) then
1970               write (6,201) i,sum
1971   201         format("imag integral used ",i9,
1972      %               "steps to get ",g18.10)
1973             else
1974               write (6,202) i,sum
1975   202         format("real integral used ",i9,
1976      %               "steps to get ",g18.10)
1977             end if
1978           end if
1979 C
1980           csum=csum+co(k)*sum
1981     2   continue
1982 C
1983         epsilon(kel)=csum
1984 C
1985     1 continue
1986 C
1987 C
1988       return
1989       end
1990 C
1991 C---------------------------------------------------------------------
1992 C
1993       double precision function ftrub(xi)
1994 C
1995       implicit double precision (a-h, o-z)
1996 C
1997       double complex dcmplx,ci,factor
1998       double complex tx_11,tx_12,tx_13
1999       double complex tx_22,tx_23,tx_33
2000 C
2001       common /para9/ alpha,omega,vt,besk2
2002       common /trubc/ kel,kreal
2003 C
2004       external tx_11,tx_12,tx_13
2005       external tx_22,tx_23,tx_33
2006 C
2007       ci=dcmplx(0.d0,1.d0)
2008 C
2009 C...... factor is the multiplier in the dielectric tensor form .......
2010 C.............. factor = i * {(wpe/wce)**2} * (wce/w0) ...............
2011 C
2012       factor=ci*alpha/omega
2013 C
2014       if (kreal.eq.0) then
2015         if (kel.eq.1) ftrub=dimag(factor*tx_11(xi))
2016         if (kel.eq.2) ftrub=dimag(factor*tx_12(xi))
2017         if (kel.eq.3) ftrub=dimag(factor*tx_13(xi))
2018         if (kel.eq.4) ftrub=dimag(factor*tx_22(xi))
2019         if (kel.eq.5) ftrub=dimag(factor*tx_23(xi))
2020         if (kel.eq.6) ftrub=dimag(factor*tx_33(xi))
2021       else
2022         if (kel.eq.1) ftrub=dreal(factor*tx_11(xi))
2023         if (kel.eq.2) ftrub=dreal(factor*tx_12(xi))
2024         if (kel.eq.3) ftrub=dreal(factor*tx_13(xi))
2025         if (kel.eq.4) ftrub=dreal(factor*tx_22(xi))
2026         if (kel.eq.5) ftrub=dreal(factor*tx_23(xi))
2027         if (kel.eq.6) ftrub=dreal(factor*tx_33(xi))
2028       end if
2029 C
2030 C
2031       return
2032       end
2033 C
2034 C---------------------------------------------------------------------
2035 C
2036       double complex function tx_11(xi)
2037 C
2038       implicit double precision (a-h, o-z)
2039 C
2040       double complex dcmplx,cdexp,cdsqrt
2041       double complex cnpe
2042       double complex ci,r,sr,t2
2043       double complex ck2,ck3,factor
2044       double complex ce,cfac,cbesk_a
2045 C
2046       double complex ckbs(2)
2047 C
2048       common /para9/ alpha,omega,vt,besk2
2049       common /wave0/ cnpa
2050       common /wave1/ cnpe
2051 C
2052       external dcbks,umach,erset
2053 C
2054       ci=dcmplx(0.d0,1.d0)
2055 C
2056       twopi=6.283185307179586476925287
2057       xim=dmod(xi,twopi)
2058 C
2059 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
2060 C.... this takes care of the normalization factors in the tensor .....
2061 C... sr is the square-root of r (argument of the Bessel functions) ...
2062 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
2063 C
2064       dom=omega*vt**2
2065       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
2066      %  (cnpa*dom*xi)**2
2067 CC      r=(1./vt**2-ci*xi*omega)**2+
2068 CC     %  2.*(cnpe*omega)**2*(1.-dcos(xim))+
2069 CC     %  (cnpa*omega*xi)**2
2070       sr=cdsqrt(r)/vt**2
2071       r=r*besk2
2072 C
2073       if (dreal(sr).lt.50.d0) then
2074 C
2075         numb=2
2076         xnu=2.
2077         call umach(-3,9)
2078         call dcbks(xnu,sr,numb,ckbs)
2079         call erset(0,1,1)
2080 C
2081         ck2=ckbs(1)*dexp(1./vt**2)
2082         ck3=ckbs(2)*dexp(1./vt**2)
2083 C
2084       else
2085 C
2086         ce=1./vt**2-sr
2087         cfac=cdexp(ce)
2088         ck2=cbesk_a(2,sr)*cfac
2089         ck3=cbesk_a(3,sr)*cfac
2090 C
2091       end if
2092 C
2093       t1=dcos(xim)
2094 C
2095       co=omega**2
2096 C
2097       t2=co*(cnpe*dsin(xim))**2
2098 C
2099 CC      tx_11=ck2/(r*vt**4*besk2)*t1-ck3/(r*sr*vt**4*besk2)*t2
2100       tx_11=ck2/r*t1-ck3/(r*sr)*t2
2101 C
2102 C
2103       return
2104       end
2105 C
2106 C---------------------------------------------------------------------
2107 C
2108       double complex function tx_12(xi)
2109 C
2110       implicit double precision (a-h, o-z)
2111 C
2112       double complex dcmplx,cdexp,cdsqrt
2113       double complex cnpe
2114       double complex ci,r,sr,t2
2115       double complex ck2,ck3
2116       double complex ce,cfac,cbesk_a
2117 C
2118       double complex ckbs(2)
2119 C
2120       common /para9/ alpha,omega,vt,besk2
2121       common /wave0/ cnpa
2122       common /wave1/ cnpe
2123 C
2124       external dcbks,umach,erset
2125 C
2126       ci=dcmplx(0.d0,1.d0)
2127 C
2128       twopi=6.283185307179586476925287
2129       xim=dmod(xi,twopi)
2130 C
2131 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
2132 C.... this takes care of the normalization factors in the tensor .....
2133 C... sr is the square-root of r (argument of the Bessel functions) ...
2134 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
2135 C
2136       dom=omega*vt**2
2137       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
2138      %  (cnpa*dom*xi)**2
2139 CC      r=(1./vt**2-ci*xi*omega)**2+
2140 CC     %  2.*(cnpe*omega)**2*(1.-dcos(xim))+
2141 CC     %  (cnpa*omega*xi)**2
2142       sr=cdsqrt(r)/vt**2
2143       r=r*besk2
2144 C
2145       if (dreal(sr).lt.50.d0) then
2146 C
2147         numb=2
2148         xnu=2.
2149         call umach(-3,9)
2150         call dcbks(xnu,sr,numb,ckbs)
2151         call erset(0,1,1)
2152 C
2153         ck2=ckbs(1)*dexp(1./vt**2)
2154         ck3=ckbs(2)*dexp(1./vt**2)
2155 C
2156       else
2157 C
2158         ce=1./vt**2-sr
2159         cfac=cdexp(ce)
2160         ck2=cbesk_a(2,sr)*cfac
2161         ck3=cbesk_a(3,sr)*cfac
2162 C
2163       end if
2164 C
2165       t1=-dsin(xim)
2166 C
2167       co=omega**2
2168 C
2169       t2=-co*cnpe**2*dsin(xim)*(1.-dcos(xim))
2170 C
2171 CC      tx_12=ck2/(r*vt**4*besk2)*t1-ck3/(r*sr*vt**4*besk2)*t2
2172       tx_12=ck2/r*t1-ck3/(r*sr)*t2
2173 C
2174 C
2175       return
2176       end
2177 C
2178 C---------------------------------------------------------------------
2179 C
2180       double complex function tx_13(xi)
2181 C
2182       implicit double precision (a-h, o-z)
2183 C
2184       double complex dcmplx,cdexp,cdsqrt
2185       double complex cnpe
2186       double complex ci,r,sr,t2
2187       double complex ck3
2188       double complex ce,cfac,cbesk_a
2189 C
2190       double complex ckbs(2)
2191 C
2192       common /para9/ alpha,omega,vt,besk2
2193       common /wave0/ cnpa
2194       common /wave1/ cnpe
2195 C
2196       external dcbks,umach,erset
2197 C
2198       ci=dcmplx(0.d0,1.d0)
2199 C
2200       twopi=6.283185307179586476925287
2201       xim=dmod(xi,twopi)
2202 C
2203 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
2204 C.... this takes care of the normalization factors in the tensor .....
2205 C... sr is the square-root of r (argument of the Bessel functions) ...
2206 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
2207 C
2208       dom=omega*vt**2
2209       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
2210      %  (cnpa*dom*xi)**2
2211 CC      r=(1./vt**2-ci*xi*omega)**2+
2212 CC     %  2.*(cnpe*omega)**2*(1.-dcos(xim))+
2213 CC     %  (cnpa*omega*xi)**2
2214       sr=cdsqrt(r)/vt**2
2215       r=r*besk2
2216 C
2217       if (dreal(sr).lt.50.d0) then
2218 C
2219         numb=2
2220         xnu=2.
2221         call umach(-3,9)
2222         call dcbks(xnu,sr,numb,ckbs)
2223         call erset(0,1,1)
2224 C
2225         ck2=ckbs(1)*dexp(1./vt**2)
2226         ck3=ckbs(2)*dexp(1./vt**2)
2227 C
2228       else
2229 C
2230         ce=1./vt**2-sr
2231         cfac=cdexp(ce)
2232         ck2=cbesk_a(2,sr)*cfac
2233         ck3=cbesk_a(3,sr)*cfac
2234 C
2235       end if
2236 C
2237       co=omega**2
2238 C
2239       t2=co*cnpe*cnpa*xi*dsin(xim)
2240 C
2241 CC      tx_13=-ck3/(r*sr*vt**4*besk2)*t2
2242       tx_13=-ck3/(r*sr)*t2
2243 C
2244 C
2245       return
2246       end
2247 C
2248 C---------------------------------------------------------------------
2249 C
2250       double complex function tx_22(xi)
2251 C
2252       implicit double precision (a-h, o-z)
2253 C
2254       double complex dcmplx,cdexp,cdsqrt
2255       double complex cnpe
2256       double complex ci,r,sr,t2
2257       double complex ck2,ck3
2258       double complex ce,cfac,cbesk_a
2259 C
2260       double complex ckbs(2)
2261 C
2262       common /para9/ alpha,omega,vt,besk2
2263       common /wave0/ cnpa
2264       common /wave1/ cnpe
2265 C
2266       external dcbks,umach,erset
2267 C
2268       ci=dcmplx(0.d0,1.d0)
2269 C
2270       twopi=6.283185307179586476925287
2271       xim=dmod(xi,twopi)
2272 C
2273 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
2274 C.... this takes care of the normalization factors in the tensor .....
2275 C... sr is the square-root of r (argument of the Bessel functions) ...
2276 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
2277 C
2278       dom=omega*vt**2
2279       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
2280      %  (cnpa*dom*xi)**2
2281 CC      r=(1./vt**2-ci*xi*omega)**2+
2282 CC     %  2.*(cnpe*omega)**2*(1.-dcos(xim))+
2283 CC     %  (cnpa*omega*xi)**2
2284       sr=cdsqrt(r)/vt**2
2285       r=r*besk2
2286 C
2287       if (dreal(sr).lt.50.d0) then
2288 C
2289         numb=2
2290         xnu=2.
2291         call umach(-3,9)
2292         call dcbks(xnu,sr,numb,ckbs)
2293         call erset(0,1,1)
2294 C
2295         ck2=ckbs(1)*dexp(1./vt**2)
2296         ck3=ckbs(2)*dexp(1./vt**2)
2297 C
2298       else
2299 C
2300         ce=1./vt**2-sr
2301         cfac=cdexp(ce)
2302         ck2=cbesk_a(2,sr)*cfac
2303         ck3=cbesk_a(3,sr)*cfac
2304 C
2305       end if
2306 C
2307       t1=dcos(xim)
2308 C
2309       co=omega**2
2310 C
2311       t2=-co*cnpe**2*(1.-dcos(xim))**2
2312 C
2313 CC      tx_22=ck2/(r*vt**4*besk2)*t1-ck3/(r*sr*vt**4*besk2)*t2
2314       tx_22=ck2/r*t1-ck3/(r*sr)*t2
2315 C
2316 C
2317       return
2318       end
2319 C
2320 C---------------------------------------------------------------------
2321 C
2322       double complex function tx_23(xi)
2323 C
2324       implicit double precision (a-h, o-z)
2325 C
2326       double complex dcmplx,cdexp,cdsqrt
2327       double complex cnpe
2328       double complex ci,r,sr,t2
2329       double complex ck3
2330       double complex ce,cfac,cbesk_a
2331 C
2332       double complex ckbs(2)
2333 C
2334       common /para9/ alpha,omega,vt,besk2
2335       common /wave0/ cnpa
2336       common /wave1/ cnpe
2337 C
2338       external dcbks,umach,erset
2339 C
2340       ci=dcmplx(0.d0,1.d0)
2341 C
2342       twopi=6.283185307179586476925287
2343       xim=dmod(xi,twopi)
2344 C
2345 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
2346 C.... this takes care of the normalization factors in the tensor .....
2347 C... sr is the square-root of r (argument of the Bessel functions) ...
2348 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
2349 C
2350       dom=omega*vt**2
2351       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
2352      %  (cnpa*dom*xi)**2
2353 CC      r=(1./vt**2-ci*xi*omega)**2+
2354 CC     %  2.*(cnpe*omega)**2*(1.-dcos(xim))+
2355 CC     %  (cnpa*omega*xi)**2
2356       sr=cdsqrt(r)/vt**2
2357       r=r*besk2
2358 C
2359       if (dreal(sr).lt.50.d0) then
2360 C
2361         numb=2
2362         xnu=2.
2363         call umach(-3,9)
2364         call dcbks(xnu,sr,numb,ckbs)
2365         call erset(0,1,1)
2366 C
2367         ck2=ckbs(1)*dexp(1./vt**2)
2368         ck3=ckbs(2)*dexp(1./vt**2)
2369 C
2370       else
2371 C
2372         ce=1./vt**2-sr
2373         cfac=cdexp(ce)
2374         ck2=cbesk_a(2,sr)*cfac
2375         ck3=cbesk_a(3,sr)*cfac
2376 C
2377       end if
2378 C
2379       co=omega**2
2380 C
2381       t2=co*cnpe*cnpa*xi*(1.-dcos(xim))
2382 C
2383 CC      tx_23=-ck3/(r*sr*vt**4*besk2)*t2
2384       tx_23=-ck3/(r*sr)*t2
2385 C
2386 C
2387       return
2388       end
2389 C
2390 C---------------------------------------------------------------------
2391 C
2392       double complex function tx_33(xi)
2393 C
2394       implicit double precision (a-h, o-z)
2395 C
2396       double complex dcmplx,cdexp,cdsqrt
2397       double complex cnpe
2398       double complex ci,r,sr,t2
2399       double complex ck2,ck3
2400       double complex ce,cfac,cbesk_a
2401 C
2402       double complex ckbs(2)
2403 C
2404       common /para9/ alpha,omega,vt,besk2
2405       common /wave0/ cnpa
2406       common /wave1/ cnpe
2407 C
2408       external dcbks,umach,erset
2409 C
2410       ci=dcmplx(0.d0,1.d0)
2411 C
2412       twopi=6.283185307179586476925287
2413       xim=dmod(xi,twopi)
2414 C
2415 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
2416 C.... this takes care of the normalization factors in the tensor .....
2417 C... sr is the square-root of r (argument of the Bessel functions) ...
2418 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
2419 C
2420       dom=omega*vt**2
2421       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
2422      %  (cnpa*dom*xi)**2
2423 CC      r=(1./vt**2-ci*xi*omega)**2+
2424 CC     %  2.*(cnpe*omega)**2*(1.-dcos(xim))+
2425 CC     %  (cnpa*omega*xi)**2
2426       sr=cdsqrt(r)/vt**2
2427       r=r*besk2
2428 C
2429       if (dreal(sr).lt.50.d0) then
2430 C
2431         numb=2
2432         xnu=2.
2433         call umach(-3,9)
2434         call dcbks(xnu,sr,numb,ckbs)
2435         call erset(0,1,1)
2436 C
2437         ck2=ckbs(1)*dexp(1./vt**2)
2438         ck3=ckbs(2)*dexp(1./vt**2)
2439 C
2440       else
2441 C
2442         ce=1./vt**2-sr
2443         cfac=cdexp(ce)
2444         ck2=cbesk_a(2,sr)*cfac
2445         ck3=cbesk_a(3,sr)*cfac
2446 C
2447       end if
2448 C
2449       co=omega**2
2450 C
2451       t2=co*(cnpa*xi)**2
2452 C
2453 CC      tx_33=ck2/(r*vt**4*besk2)-ck3/(r*sr*vt**4*besk2)*t2
2454       tx_33=ck2/r-ck3/(r*sr)*t2
2455 C
2456 C
2457       return
2458       end
2459 C
2460 C---------------------------------------------------------------------
2461 C
2462 C#####################################################################
2463 C       evaluating the relativistic dielectric tensor elements        
2464 C          this part of the code is for the Weiss' formalism          
2465 C#####################################################################
2466 C
2467 C.....................................................................
2468 C
2469 C.............. doing the two dimensional integrations ...............
2470 C
2471       subroutine weiss(epsilon)
2472 C
2473       implicit double precision (a-h, o-z)
2474 C
2475       double complex cnpe,ckpe,dcmplx,csum,ci
2476 C
2477       double complex co(0:1),ep(6),epsilon(6)
2478 C
2479       common /para9/ alpha,omega,vt,besk2
2480       common /ppara/ p_para_min,p_para_max
2481       common /pperp/ p_perp_min,p_perp_max
2482       common /weerr/ errabs0,errrel0,irule0
2483       common /wave0/ cnpa
2484       common /wave1/ cnpe
2485       common /wkpa/  ckpa
2486       common /wkpe/  ckpe
2487       common /nar2/  ireal
2488       common /nar4/  kel
2489 C
2490       external dqdag,fperp,umach,erset
2491 C
2492       ckpa=cnpa*omega
2493       ckpe=cnpe*omega
2494 C
2495       pi=3.141592653589793238462643d0
2496       factor=0.5*pi*alpha/(vt**4*omega*besk2)
2497 C
2498       co(0)=dcmplx(0.d0,1.d0)
2499       co(1)=dcmplx(1.d0,0.d0)
2500 C
2501       do 1 j=1,6
2502         kel=j
2503         csum=dcmplx(0.d0,0.d0)
2504 C
2505         do 2 i=0,1
2506           ireal=i
2507           errabs=errabs0
2508           errrel=errrel0
2509           irule=irule0
2510           call umach(-3,9)
2511           call dqdag(fperp,p_perp_min,p_perp_max,
2512      %              errabs,errrel,irule,result,errest)
2513           call erset(0,1,1)
2514           csum=csum+co(i)*result
2515     2   continue
2516 C
2517         ep(kel)=csum*factor
2518 C
2519     1 continue
2520 C
2521       ci=dcmplx(0.d0,1.d0)
2522 C
2523       epsilon(1)=0.5d0*(ep(1)+2.d0*ep(2)+ep(4))
2524       epsilon(2)=ci*0.5d0*(ep(1)-ep(4))
2525       epsilon(3)=1.d0/dsqrt(2.d0)*(ep(3)-ep(5))
2526       epsilon(4)=0.5d0*(ep(1)-2.d0*ep(2)+ep(4))
2527       epsilon(5)=-ci/dsqrt(2.d0)*(ep(3)+ep(5))
2528       epsilon(6)=ep(6)
2529 C
2530 C
2531       return
2532       end
2533 C
2534 C=====================================================================
2535 C
2536       double precision function fperp(p_perp)
2537 C
2538       implicit double precision (a-h, o-z)
2539 C
2540       double precision pa(200)
2541 C
2542       double complex ci,dcmplx,sum
2543 C
2544       double complex co(0:1)
2545 C
2546       common /weerr/ errabs0,errrel0,irule0
2547       common /ppara/ p_para_min,p_para_max
2548       common /var1/  ppe,ppa0
2549       common /nar1/  nsing
2550       common /nar2/  ireal
2551       common /nar3/  kreal
2552 C
2553       external dqdawc,fpara,dqdag,umach,erset
2554 C
2555       pi=3.141592653589793238462643d0
2556       ci=dcmplx(0.d0,1.d0)
2557 C
2558       co(0)=ci
2559       co(1)=dcmplx(1.d0,0.d0)
2560 C
2561       ppe=p_perp
2562 C
2563       call sing(ppe,p_para_min,p_para_max,pa,index)
2564 C
2565 C...... nsing keeps track of whether we have singularity or not ......
2566 C.... if index is not zero then all intervals in the integration .....
2567 C.... contain a singular point, otherwise there is no singularity ....
2568 C... for no singularity the integration routine used is different ....
2569 C
2570       if (index.eq.0) then
2571         nsing=0
2572       else
2573         nsing=1
2574       end if
2575 C
2576       if (index.ne.0) then
2577         if (index.eq.1) then
2578           sum=dcmplx(0.d0,0.d0)
2579           do 1 i=0,1
2580             kreal=i
2581             errabs=errabs0
2582             errrel=errrel0
2583             pin=p_para_min
2584             pfi=p_para_max
2585             ppa0=pa(1)
2586             call umach(-3,9)
2587             call dqdawc(fpara,pin,pfi,ppa0,
2588      %                  errabs,errrel,result,errest)
2589             call erset(0,1,1)
2590             sum=sum+co(i)*(result+ci*pi*fpara(ppa0))
2591     1     continue
2592           go to 5
2593         else
2594           sum=dcmplx(0.d0,0.d0)
2595           pin=p_para_min
2596           do 2 i=1,index-1
2597             do 3 j=0,1
2598               kreal=j
2599               pfi=(pa(i+1)+pa(i))/2.
2600               errabs=errabs0
2601               errrel=errrel0
2602               ppa0=pa(i)
2603               call umach(-3,9)
2604 CC              write (10,201) pin,pfi,ppa0
2605   201         format("fperp: pin,pfi,ppa0 =",3g18.10)
2606               call dqdawc(fpara,pin,pfi,ppa0,
2607      %                     errabs,errrel,result,errest)
2608               call erset(0,1,1)
2609               sum=sum+co(j)*(result+ci*pi*fpara(ppa0))
2610 CC              write (10,202) sum
2611   202         format("sum = ",2g18.10)
2612     3       continue
2613             pin=pfi
2614     2     continue
2615           do 4 i=0,1
2616             kreal=i
2617             pfi=p_para_max
2618             errabs=errabs0
2619             errrel=errrel0
2620             ppa0=pa(index)
2621 CC            write (10,201) pin,pfi,ppa0
2622             call umach(-3,9)
2623             call dqdawc(fpara,pin,pfi,ppa0,
2624      %                  errabs,errrel,result,errest)
2625             call erset(0,1,1)
2626             sum=sum+co(i)*(result+ci*pi*fpara(ppa0))
2627 CC            write (10,202) sum
2628     4     continue
2629         end if
2630       end if
2631 C
2632       if (index.eq.0) then
2633         errabs=errabs0
2634         errrel=errrel0
2635         irule=irule0
2636 C
2637         pin=p_para_min
2638         pfi=p_para_max
2639         call umach(-3,9)
2640         call dqdag(fpara,pin,pfi,
2641      %             errabs,errrel,irule,result,errest)
2642         call erset(0,1,1)
2643         sum=dcmplx(result,0.d0)
2644       end if
2645 C
2646     5 continue
2647 C
2648       if (ireal.eq.0) then
2649         fperp=dimag(sum)
2650       else
2651         fperp=dreal(sum)
2652       end if
2653 C
2654 C
2655       return
2656       end
2657 C
2658 C=====================================================================
2659 C
2660       double precision function fpara(p_para)
2661 C
2662       implicit double precision (a-h, o-z)
2663 C
2664       double complex sigma_11,sigma_12,sigma_13
2665       double complex sigma_22,sigma_23,sigma_33
2666 C
2667       common /nar3/ kreal
2668       common /nar4/ kel
2669 C
2670       external sigma_11,sigma_12,sigma_13
2671       external sigma_22,sigma_23,sigma_33
2672 C
2673       if (kreal.eq.0) then
2674         if (kel.eq.1) fpara=dimag(sigma_11(p_para))
2675         if (kel.eq.2) fpara=dimag(sigma_12(p_para))
2676         if (kel.eq.3) fpara=dimag(sigma_13(p_para))
2677         if (kel.eq.4) fpara=dimag(sigma_22(p_para))
2678         if (kel.eq.5) fpara=dimag(sigma_23(p_para))
2679         if (kel.eq.6) fpara=dimag(sigma_33(p_para))
2680       else
2681         if (kel.eq.1) fpara=dreal(sigma_11(p_para))
2682         if (kel.eq.2) fpara=dreal(sigma_12(p_para))
2683         if (kel.eq.3) fpara=dreal(sigma_13(p_para))
2684         if (kel.eq.4) fpara=dreal(sigma_22(p_para))
2685         if (kel.eq.5) fpara=dreal(sigma_23(p_para))
2686         if (kel.eq.6) fpara=dreal(sigma_33(p_para))
2687       end if
2688 C
2689 C
2690       return
2691       end
2692 C
2693 C=====================================================================
2694 C
2695       double complex function sigma_11(p_para)
2696 C
2697       implicit double precision (a-h, o-z)
2698 C
2699       double complex ckpe,carg,besfun
2700 C
2701       double complex cbesj(1),cbesy(1)
2702 C
2703       common /para9/ alpha,omega,vt,besk2
2704       common /wkpa/  ckpa
2705       common /wkpe/  ckpe
2706       common /var1/  p_perp,ppa0
2707       common /nar1/  nsing
2708       common /nar3/  kreal
2709 C
2710       external sinom
2711       external dcbjs,dcbys
2712 C
2713       pi=3.141592653589793238462643d0
2714 C
2715       gamma=dsqrt(1.+p_perp**2+p_para**2)
2716 C
2717       omt=omega*gamma-ckpa*p_para
2718       carg=ckpe*p_perp
2719 C
2720       ord=dabs(omt+1.d0)
2721       numb=1
2722       call dcbjs(ord,carg,numb,cbesj)
2723       numb=1
2724       call dcbys(ord,carg,numb,cbesy)
2725 C
2726       besfun=cbesj(1)*(cbesj(1)*dcos(pi*ord)-
2727      %                 cbesy(1)*dsin(pi*ord))
2728 C
2729       sigma_11=p_perp**3/(2.d0*gamma)*
2730      %         dexp(-(gamma-1.)/vt**2)*
2731      %         sinom(p_perp,p_para,ppa0)*besfun
2732 C
2733 C
2734       return
2735       end
2736 C
2737 C=====================================================================
2738 C
2739       double complex function sigma_12(p_para)
2740 C
2741       implicit double precision (a-h, o-z)
2742 C
2743       double complex ckpe,carg,besfun,cdum
2744 C
2745       double complex cbesj(1),cbesy(1)
2746 C
2747       common /para9/ alpha,omega,vt,besk2
2748       common /wkpa/  ckpa
2749       common /wkpe/  ckpe
2750       common /var1/  p_perp,ppa0
2751       common /nar1/  nsing
2752       common /nar3/  kreal
2753 C
2754       external sinom
2755       external dcbjs,dcbys
2756 C
2757       pi=3.141592653589793238462643d0
2758 C
2759       gamma=dsqrt(1.+p_perp**2+p_para**2)
2760 C
2761       omt=omega*gamma-ckpa*p_para
2762       carg=ckpe*p_perp
2763 C
2764       ord=1.d0+dabs(omt)
2765       numb=1
2766       call dcbjs(ord,carg,numb,cbesj)
2767       cdum=cbesj(1)
2768       if (dabs(omt).le.1.d0) then
2769         ord=1.d0-dabs(omt)
2770         numb=1
2771         call dcbjs(ord,carg,numb,cbesj)
2772         besfun=cdum*cbesj(1)
2773       else
2774         ord=dabs(omt)-1.d0
2775         numb=1
2776         call dcbjs(ord,carg,numb,cbesj)
2777         numb=1
2778         call dcbys(ord,carg,numb,cbesy)
2779         besfun=cdum*(cbesj(1)*dcos(pi*ord)-
2780      %               cbesy(1)*dsin(pi*ord))
2781       end if
2782 C
2783       sigma_12=p_perp**3/(2.d0*gamma)*
2784      %         dexp(-(gamma-1.)/vt**2)*
2785      %         sinom(p_perp,p_para,ppa0)*besfun
2786 C
2787 C
2788       return
2789       end
2790 C
2791 C=====================================================================
2792 C
2793       double complex function sigma_13(p_para)
2794 C
2795       implicit double precision (a-h, o-z)
2796 C
2797       double complex ckpe,carg,besfun,cdum
2798 C
2799       double complex cbesj(1),cbesy(1)
2800 C
2801       common /para9/ alpha,omega,vt,besk2
2802       common /wkpa/  ckpa
2803       common /wkpe/  ckpe
2804       common /var1/  p_perp,ppa0
2805       common /nar1/  nsing
2806       common /nar3/  kreal
2807 C
2808       external sinom
2809       external dcbjs,dcbys
2810 C
2811       pi=3.141592653589793238462643d0
2812 C
2813       gamma=dsqrt(1.+p_perp**2+p_para**2)
2814 C
2815       omt=omega*gamma-ckpa*p_para
2816       carg=ckpe*p_perp
2817 C
2818       if (omt.ge.0.d0) then
2819         ord=omt+1.d0
2820         numb=1
2821         call dcbjs(ord,carg,numb,cbesj)
2822         cdum=cbesj(1)
2823         ord=omt
2824         numb=1
2825         call dcbjs(ord,carg,numb,cbesj)
2826         numb=1
2827         call dcbys(ord,carg,numb,cbesy)
2828         besfun=cdum*(cbesj(1)*dcos(pi*ord)-
2829      %               cbesy(1)*dsin(pi*ord))
2830       end if
2831 C
2832       if (omt.gt.-1.d0.and.omt.lt.0.d0) then
2833         ord=omt+1.d0
2834         numb=1
2835         call dcbjs(ord,carg,numb,cbesj)
2836         cdum=cbesj(1)
2837         ord=-omt
2838         numb=1
2839         call dcbjs(ord,carg,numb,cbesj)
2840         besfun=cdum*cbesj(1)
2841       end if
2842 C
2843       if (omt.le.-1.d0) then
2844         ord=-omt
2845         numb=1
2846         call dcbjs(ord,carg,numb,cbesj)
2847         cdum=cbesj(1)
2848         ord=-1.d0-omt
2849         numb=1
2850         call dcbjs(ord,carg,numb,cbesj)
2851         numb=1
2852         call dcbys(ord,carg,numb,cbesy)
2853         besfun=cdum*(cbesj(1)*dcos(pi*ord)-
2854      %               cbesy(1)*dsin(pi*ord))
2855       end if
2856 C
2857       sigma_13=-p_perp**2*p_para/(dsqrt(2.d0)*gamma)*
2858      %          dexp(-(gamma-1.)/vt**2)*
2859      %          sinom(p_perp,p_para,ppa0)*besfun
2860 C
2861 C
2862       return
2863       end
2864 C
2865 C=====================================================================
2866 C
2867       double complex function sigma_22(p_para)
2868 C
2869       implicit double precision (a-h, o-z)
2870 C
2871       double complex ckpe,carg,besfun,cdum
2872 C
2873       double complex cbesj(1),cbesy(1)
2874 C
2875       common /para9/ alpha,omega,vt,besk2
2876       common /wkpa/  ckpa
2877       common /wkpe/  ckpe
2878       common /var1/  p_perp,ppa0
2879       common /nar1/  nsing
2880       common /nar3/  kreal
2881 C
2882       external sinom
2883       external dcbjs,dcbys
2884 C
2885       pi=3.141592653589793238462643d0
2886 C
2887       gamma=dsqrt(1.+p_perp**2+p_para**2)
2888 C
2889       omt=omega*gamma-ckpa*p_para
2890       carg=ckpe*p_perp
2891 C
2892       ord=dabs(omt-1.d0)
2893       numb=1
2894       call dcbjs(ord,carg,numb,cbesj)
2895       numb=1
2896       call dcbys(ord,carg,numb,cbesy)
2897       besfun=cbesj(1)*(cbesj(1)*dcos(pi*ord)-
2898      %                 cbesy(1)*dsin(pi*ord))
2899 C
2900       sigma_22=p_perp**3/(2.d0*gamma)*
2901      %         dexp(-(gamma-1.)/vt**2)*
2902      %         sinom(p_perp,p_para,ppa0)*besfun
2903 C
2904 C
2905       return
2906       end
2907 C
2908 C=====================================================================
2909 C
2910       double complex function sigma_23(p_para)
2911 C
2912       implicit double precision (a-h, o-z)
2913 C
2914       double complex ckpe,carg,besfun,cdum
2915 C
2916       double complex cbesj(1),cbesy(1)
2917 C
2918       common /para9/ alpha,omega,vt,besk2
2919       common /wkpa/  ckpa
2920       common /wkpe/  ckpe
2921       common /var1/  p_perp,ppa0
2922       common /nar1/  nsing
2923       common /nar3/  kreal
2924 C
2925       external sinom
2926       external dcbjs,dcbys
2927 C
2928       pi=3.141592653589793238462643d0
2929 C
2930       gamma=dsqrt(1.+p_perp**2+p_para**2)
2931 C
2932       omt=omega*gamma-ckpa*p_para
2933       carg=ckpe*p_perp
2934 C
2935       if (omt.le.0.d0) then
2936         ord=1.d0-omt
2937         numb=1
2938         call dcbjs(ord,carg,numb,cbesj)
2939         cdum=cbesj(1)
2940         ord=-omt
2941         numb=1
2942         call dcbjs(ord,carg,numb,cbesj)
2943         numb=1
2944         call dcbys(ord,carg,numb,cbesy)
2945         besfun=cdum*(cbesj(1)*dcos(pi*ord)-
2946      %               cbesy(1)*dsin(pi*ord))
2947       end if
2948 C
2949       if (omt.gt.0.d0.and.omt.lt.1.d0) then
2950         ord=omt
2951         numb=1
2952         call dcbjs(ord,carg,numb,cbesj)
2953         cdum=cbesj(1)
2954         ord=1.d0-omt
2955         numb=1
2956         call dcbjs(ord,carg,numb,cbesj)
2957         besfun=cdum*cbesj(1)
2958       end if
2959 C
2960       if (omt.ge.1.d0) then
2961         ord=omt
2962         numb=1
2963         call dcbjs(ord,carg,numb,cbesj)
2964         cdum=cbesj(1)
2965         ord=omt-1.d0
2966         numb=1
2967         call dcbjs(ord,carg,numb,cbesj)
2968         numb=1
2969         call dcbys(ord,carg,numb,cbesy)
2970         besfun=cdum*(cbesj(1)*dcos(pi*ord)-
2971      %               cbesy(1)*dsin(pi*ord))
2972       end if
2973 C
2974       sigma_23=-p_perp**2*p_para/(dsqrt(2.d0)*gamma)*
2975      %          dexp(-(gamma-1.)/vt**2)*
2976      %          sinom(p_perp,p_para,ppa0)*besfun
2977 C
2978 C
2979       return
2980       end
2981 C
2982 C=====================================================================
2983 C
2984       double complex function sigma_33(p_para)
2985 C
2986       implicit double precision (a-h, o-z)
2987 C
2988       double complex ckpe,carg,besfun
2989 C
2990       double complex cbesj(1),cbesy(1)
2991 C
2992       common /para9/ alpha,omega,vt,besk2
2993       common /wkpa/  ckpa
2994       common /wkpe/  ckpe
2995       common /var1/  p_perp,ppa0
2996       common /nar1/  nsing
2997       common /nar3/  kreal
2998 C
2999       external sinom
3000       external dcbjs,dcbys
3001 C
3002       pi=3.141592653589793238462643d0
3003 C
3004       gamma=dsqrt(1.+p_perp**2+p_para**2)
3005 C
3006       omt=omega*gamma-ckpa*p_para
3007       carg=ckpe*p_perp
3008 C
3009       ord=dabs(omt)
3010       numb=1
3011       call dcbjs(ord,carg,numb,cbesj)
3012       numb=1
3013       call dcbys(ord,carg,numb,cbesy)
3014       besfun=cbesj(1)*(cbesj(1)*dcos(pi*ord)-
3015      %                 cbesy(1)*dsin(pi*ord))
3016 C
3017       sigma_33=-p_perp*p_para**2/gamma*dexp(-(gamma-1.)/vt**2)*
3018      %          sinom(p_perp,p_para,ppa0)*besfun
3019 C
3020 C
3021       return
3022       end
3023 C
3024 C=====================================================================
3025 C
3026 C........ this is the sinusoidal part of the tensor elements .........
3027 C
3028       double precision function sinom(ppe,ppa,ppa0)
3029 C
3030       implicit double precision (a-h,o-z)
3031 C
3032       common /para9/ alpha,omega,vt,besk2
3033       common /wkpa/  ckpa
3034       common /nar1/  nsing
3035 C
3036       pi=3.141592653589793238462643d0
3037 C
3038       gamma=dsqrt(1.+ppe**2+ppa**2)
3039       omt=(omega*gamma-ckpa*ppa)
3040 C
3041       if (nsing.eq.1) then
3042         domt=omega*ppa/gamma-ckpa
3043         if (ppa.eq.ppa0) then
3044           fun=1./(pi*dcos(pi*omt)*domt)
3045         else
3046           fun=(ppa-ppa0)/(dsin(pi*omt))
3047         end if
3048       else
3049         fun=1./(dsin(pi*omt))
3050       end if
3051 C
3052       sinom=fun
3053 C
3054 C
3055       return
3056       end
3057 C
3058 C=====================================================================
3059 C
3060 C...... finding the location of the singularities in p_parallel ......
3061 C.......... this is the case when n_para is not equal to 1 ...........
3062 C
3063       subroutine sing(p_perp,p_para_min,p_para_max,pres,index)
3064 C
3065       implicit double precision (a-h, o-z)
3066 C
3067       double precision pres(200)
3068 C
3069       common /para9/ alpha,omega,vt,besk2
3070       common /wkpa/  ckpa
3071 C
3072 C....... calculating the range in l for a range in p_parallel ........
3073 C
3074       gamma=dsqrt(1.+p_perp**2+p_para_min**2)
3075       l1=int(gamma*omega-ckpa*p_para_min)
3076       gamma=dsqrt(1.+p_perp**2+p_para_max**2)
3077       l2=int(gamma*omega-ckpa*p_para_max)
3078 C
3079       lmax=max(iabs(l1),iabs(l2))
3080 C
3081       index=0
3082 C
3083       do 1 i=0,lmax
3084         fl=dfloat(i)
3085         sqp=(ckpa**2-omega**2)*(1.+p_perp**2)+fl**2
3086         if (sqp.lt.0.) then
3087           go to 2
3088         else
3089           dummy=(ckpa*fl+omega*dsqrt(sqp))/(omega**2-ckpa**2)
3090           if (dummy.gt.p_para_min.and.dummy.lt.p_para_max) then
3091             index=index+1
3092             pres(index)=dummy
3093           end if
3094           dummy=(ckpa*fl-omega*dsqrt(sqp))/(omega**2-ckpa**2)
3095           if (dummy.gt.p_para_min.and.dummy.lt.p_para_max) then
3096             index=index+1
3097             pres(index)=dummy
3098           end if
3099         end if
3100 C
3101     2   continue
3102 C
3103     1 continue
3104 C
3105 C...... arranging the pres array in order of increasing p_para .......
3106 C
3107       nr=0
3108 C
3109       do 3 i=1,index
3110         dum1=pres(i)
3111         do 4 j=i+1,index
3112           if (pres(j).lt.dum1) then
3113             temp1=pres(j)
3114             pres(j)=dum1
3115             dum1=temp1
3116           end if
3117     4   continue
3118         pres(i)=dum1
3119     3 continue
3120 C
3121 C
3122       return
3123       end
3124 C
3125 C=====================================================================
3126 C=====================================================================
3127 C=====================================================================
3128 C=====================================================================
3129 C
3130 C........ for plotting the data for a given tensor component .........
3131 C
3132       subroutine fun_plot
3133 C
3134       implicit double precision (a-h,o-z)
3135 C
3136       common /trub1/ xi_lo,dxi,nxi
3137 C
3138       external rinteg,ainteg
3139 C
3140       xi=xi_lo-dxi
3141 C
3142       do 1 i=1,nxi
3143         xi=xi+dxi
3144         dumr=rinteg(xi)
3145         dumi=ainteg(xi)
3146         if (dabs(dumr).lt.1.d-40) dumr=0.d0
3147         if (dabs(dumi).lt.1.d-40) dumi=0.d0
3148 C        write (5,101) xi,dumr,dumi
3149     1 continue
3150 C
3151   101 format(3g15.10)
3152 C
3153 C
3154       return
3155       end
3156 C
3157 C
3158 C---------------------------------------------------------------------
3159 C
3160 C.......... evaluating the integral for each tensor element ..........
3161 C
3162       subroutine cint(c_int)
3163 C
3164       implicit double precision (a-h, o-z)
3165 C
3166       double precision sumr(0:1000),stot(3)
3167 C
3168       double complex c_int,csum
3169       double complex dcmplx
3170 C
3171       common /trub1/ xi_lo,dxi,nxi
3172       common /trub2/ diff_err,navg
3173       common /trerr/ errabs0,errrel0,irule0
3174 C
3175       external rinteg,ainteg,dqdag,umach,erset
3176 C
3177 C.. ntr keeps track of the number of averages (3) that are compared ..
3178 C
3179       ntr=0
3180 C
3181       xi_0=xi_lo
3182       xi_1=xi_0+dxi
3183 C
3184       sum=0.d0
3185 C
3186       do 1 i=1,nxi
3187 C
3188         errabs=errabs0
3189         errrel=errrel0
3190         irule=irule0
3191 C
3192         call umach(-3,9)
3193         call dqdag(rinteg,xi_0,xi_1,errabs,errrel,irule,
3194      %             t_real,errest)
3195         call erset(0,1,1)
3196 C
3197         sum=sum+t_real
3198 C
3199         jj=mod(i-1,navg)
3200         sumr(jj)=sum
3201         if (jj.eq.navg-1) then
3202           ntr=ntr+1
3203           stot(ntr)=0.d0
3204           do 2 j=0,navg-1
3205             stot(ntr)=stot(ntr)+sumr(j)
3206     2     continue
3207           stot(ntr)=stot(ntr)/dfloat(navg)
3208           if (ntr.eq.3) then
3209             diff1=dabs(dabs(stot(1))-dabs(stot(2)))
3210             diff2=dabs(dabs(stot(2))-dabs(stot(3)))
3211             if (diff1.lt.diff_err.and.diff2.lt.diff_err) then
3212               sum_real=stot(3)
3213               go to 3
3214             else
3215               stot(1)=stot(2)
3216               stot(2)=stot(3)
3217               ntr=2
3218             end if
3219           end if
3220         end if
3221 C
3222         xi_0=xi_1
3223         xi_1=xi_0+dxi
3224 C
3225     1 continue
3226 C
3227     3 continue
3228 C
3229       write (6,201) i,sum
3230   201 format("real integral used ",i9," steps to get ",g18.10)
3231 C
3232       if (i.ge.nxi) sum_real=sum
3233 C
3234       xi_0=xi_lo
3235       xi_1=xi_0+dxi
3236 C
3237       sum=0.d0
3238       ntr=0
3239 C
3240       do 4 i=1,nxi
3241 C
3242         errabs=errabs0
3243         errrel=errrel0
3244         irule=irule0
3245 C
3246         call umach(-3,9)
3247         call dqdag(ainteg,xi_0,xi_1,errabs,errrel,irule,
3248      %             t_imag,errest)
3249         call erset(0,1,1)
3250 C
3251         sum=sum+t_imag
3252 C
3253         jj=mod(i-1,navg)
3254         sumr(jj)=sum
3255         if (jj.eq.navg-1) then
3256           ntr=ntr+1
3257           stot(ntr)=0.d0
3258           do 5 j=0,navg-1
3259             stot(ntr)=stot(ntr)+sumr(j)
3260     5     continue
3261           stot(ntr)=stot(ntr)/dfloat(navg)
3262           if (ntr.eq.3) then
3263             diff1=dabs(dabs(stot(1))-dabs(stot(2)))
3264             diff2=dabs(dabs(stot(2))-dabs(stot(3)))
3265             if (diff1.lt.diff_err.and.diff2.lt.diff_err) then
3266               sum_imag=stot(3)
3267               go to 6
3268             else
3269               stot(1)=stot(2)
3270               stot(2)=stot(3)
3271               ntr=2
3272             end if
3273           end if
3274         end if
3275 C
3276         xi_0=xi_1
3277         xi_1=xi_0+dxi
3278 C
3279     4 continue
3280 C
3281     6 continue
3282 C
3283       if (i.ge.nxi) sum_imag=sum
3284 C
3285       csum=dcmplx(sum_real,sum_imag)
3286 C
3287       write (6,202) i,sum
3288   202 format("imaginary integral used ",i9," steps to get ",g18.10)
3289 C
3290       c_int=csum
3291 C
3292 C
3293       return
3294       end
3295 C
3296 C---------------------------------------------------------------------
3297 C
3298 C... real and imaginary parts of the integrand (given separately) ....
3299 C
3300 C
3301       double precision function rinteg(xi)
3302 C
3303       implicit double precision (a-h, o-z)
3304 C
3305       double precision t1(3,3)
3306 C
3307       double complex dcmplx,cdexp,cdsqrt
3308       double complex cnpe
3309       double complex ci,r,sr
3310       double complex ck2,ck3,factor
3311       double complex ce,cfac,cbesk_a
3312       double complex delta_t,cd_del
3313       double complex rt,srt
3314 C
3315       double complex ckbs(2)
3316       double complex t2(3,3)
3317 C
3318       common /para9/ alpha,omega,vt,besk2
3319       common /wave0/ cnpa
3320       common /wave1/ cnpe
3321       common /nele/  n1,n2
3322 C
3323       external dcbks,umach,erset
3324       external cbesk_a
3325 C
3326       ci=dcmplx(0.d0,1.d0)
3327 C
3328       twopi=6.283185307179586476925287
3329       xim=dmod(xi,twopi)
3330 C
3331 C...... factor is the multiplier in the dielectric tensor form .......
3332 C
3333       factor=ci*alpha/omega
3334 C
3335 C... sr is the square-root of r (argument of the Bessel functions) ...
3336 C
3337       r=(1./vt**2-ci*xi*omega)**2+
3338      %  2.*(cnpe*omega)**2*(1.-dcos(xim))+
3339      %  (cnpa*omega*xi)**2
3340       sr=cdsqrt(r)
3341 C
3342       if (dreal(sr).lt.50.d0) then
3343 C
3344         numb=2
3345         xnu=2.
3346         call umach(-3,9)
3347         call dcbks(xnu,sr,numb,ckbs)
3348         call erset(0,1,1)
3349 C
3350         ck2=ckbs(1)*dexp(1./vt**2)
3351         ck3=ckbs(2)*dexp(1./vt**2)
3352 C
3353       else
3354 C
3355         ce=1./vt**2-sr
3356         cfac=cdexp(ce)
3357         ck2=cbesk_a(2,sr)*cfac
3358         ck3=cbesk_a(3,sr)*cfac
3359 C
3360       end if
3361 C
3362       t1(1,1)=dcos(xim)
3363       t1(1,2)=-dsin(xim)
3364       t1(1,3)=0.d0
3365       t1(2,1)=-t1(1,2)
3366       t1(2,2)=t1(1,1)
3367       t1(2,3)=0.d0
3368       t1(3,1)=0.d0
3369       t1(3,2)=0.d0
3370       t1(3,3)=1.d0
3371 C
3372       co=omega**2
3373 C
3374       t2(1,1)=co*(cnpe*dsin(xim))**2
3375       t2(1,2)=-co*cnpe**2*dsin(xim)*(1.-dcos(xim))
3376       t2(1,3)=co*cnpe*cnpa*xi*dsin(xim)
3377       t2(2,1)=-t2(1,2)
3378       t2(2,2)=-co*cnpe**2*(1.-dcos(xim))**2
3379       t2(2,3)=co*cnpe*cnpa*xi*(1.-dcos(xim))
3380       t2(3,1)=t2(1,3)
3381       t2(3,2)=-t2(2,3)
3382       t2(3,3)=co*(cnpa*xi)**2
3383 C
3384 C.................... the relativistic integrand .....................
3385 C
3386       rinteg_r=dreal(factor*(ck2/(r*vt**4*besk2)*t1(n1,n2)-
3387      %                       ck3/(r*sr*vt**4*besk2)*t2(n1,n2)))
3388 C
3389       rinteg=rinteg_r
3390 C
3391 C
3392       return
3393       end
3394 C
3395 C---------------------------------------------------------------------
3396 C
3397 C... real and imaginary parts of the integrand (given separately) ....
3398 C
3399 C
3400       double precision function ainteg(xi)
3401 C
3402       implicit double precision (a-h, o-z)
3403 C
3404       double precision t1(3,3)
3405 C
3406       double complex dcmplx,cdexp,cdsqrt
3407       double complex cnpe
3408       double complex ci,r,sr
3409       double complex ck2,ck3,factor
3410       double complex ce,cfac,cbesk_a
3411       double complex delta_t,cd_del
3412       double complex rt,srt
3413 C
3414       double complex ckbs(2)
3415       double complex t2(3,3)
3416 C
3417       common /para9/ alpha,omega,vt,besk2
3418       common /wave0/ cnpa
3419       common /wave1/ cnpe
3420       common /nele/  n1,n2
3421 C
3422       external dcbks,umach,erset
3423       external cbesk_a
3424 C
3425       ci=dcmplx(0.d0,1.d0)
3426 C
3427 C...... factor is the multiplier in the dielectric tensor form .......
3428 C
3429       factor=ci*alpha/omega
3430 C
3431 C... sr is the square-root of r (argument of the Bessel functions) ...
3432 C
3433       r=(1./vt**2-ci*xi*omega)**2+
3434      %  2.*(cnpe*omega)**2*(1.-dcos(xi))+
3435      %  (cnpa*omega*xi)**2
3436       sr=cdsqrt(r)
3437 C
3438       if (dreal(sr).lt.50.d0) then
3439 C
3440         numb=2
3441         xnu=2.
3442         call umach(-3,9)
3443         call dcbks(xnu,sr,numb,ckbs)
3444         call erset(0,1,1)
3445 C
3446         ck2=ckbs(1)*dexp(1./vt**2)
3447         ck3=ckbs(2)*dexp(1./vt**2)
3448 C
3449       else
3450 C
3451         ce=1./vt**2-sr
3452         cfac=cdexp(ce)
3453         ck2=cbesk_a(2,sr)*cfac
3454         ck3=cbesk_a(3,sr)*cfac
3455 C
3456       end if
3457 C
3458       t1(1,1)=dcos(xi)
3459       t1(1,2)=-dsin(xi)
3460       t1(1,3)=0.d0
3461       t1(2,1)=-t1(1,2)
3462       t1(2,2)=t1(1,1)
3463       t1(2,3)=0.d0
3464       t1(3,1)=0.d0
3465       t1(3,2)=0.d0
3466       t1(3,3)=1.d0
3467 C
3468       co=omega**2
3469 C
3470       t2(1,1)=co*(cnpe*dsin(xi))**2
3471       t2(1,2)=-co*cnpe**2*dsin(xi)*(1.-dcos(xi))
3472       t2(1,3)=co*cnpe*cnpa*xi*dsin(xi)
3473       t2(2,1)=-t2(1,2)
3474       t2(2,2)=-co*cnpe**2*(1.-dcos(xi))**2
3475       t2(2,3)=co*cnpe*cnpa*xi*(1.-dcos(xi))
3476       t2(3,1)=t2(1,3)
3477       t2(3,2)=-t2(2,3)
3478       t2(3,3)=co*(cnpa*xi)**2
3479 C
3480 C.................... the relativistic integrand .....................
3481 C
3482       ainteg_r=dimag(factor*(ck2/(r*vt**4*besk2)*t1(n1,n2)-
3483      %                       ck3/(r*sr*vt**4*besk2)*t2(n1,n2)))
3484 C
3485       ainteg=ainteg_r
3486 C
3487 C
3488       return
3489       end
3490 C
3491 C---------------------------------------------------------------------
3492 C
3493 C---------------------------------------------------------------------
3494 C.......... non-relativistic power flow and energy density ...........
3495 C---------------------------------------------------------------------
3496 C
3497 C.... solving the non-relativistic hermitian dispersion relation .....
3498 C
3499       subroutine nr_disp_herm(data,rguess,nrts)
3500 C
3501       implicit double precision (a-h, o-z)
3502 C
3503       parameter (n_rts=10)
3504 C
3505       integer info(n_rts)
3506 C
3507       double precision nr_func_herm
3508 C
3509       double precision re_rt(n_rts)
3510 C
3511       double complex rguess(n_rts),roots(n_rts)
3512 C
3513       common /rterr/ errt_abs,errt_rel,it_max
3514       common /nrels/ nr_nterms
3515       common /prnt/  kpr
3516 C
3517       external nr_func_herm,dzbren,umach,erset
3518 C
3519 C.......... taking the real part of the complex input roots ..........
3520 C
3521       do 1 i=1,nrts
3522         re_rt(i)=dreal(rguess(i))
3523     1 continue
3524 C
3525 C... finding the maximum and minimum limits for finding the roots ....
3526 C
3527       rtmin=re_rt(1)
3528       rtmax=re_rt(1)
3529       do 2 i=2,nrts
3530         rtmin=dmin1(rtmin,re_rt(i))
3531         rtmax=dmax1(rtmax,re_rt(i))
3532     2 continue
3533 C
3534       rtmin=(1.d0-0.25*rtmin/dabs(rtmin))*rtmin
3535       rtmax=(1.d0+0.25*rtmax/dabs(rtmax))*rtmax
3536 C
3537 C....... finding the real n_perp roots for the hermitian tensor ......
3538 C
3539       newrt=0
3540       nu_int=11
3541       drt=(rtmax-rtmin)/dfloat(nu_int)
3542       rti=rtmin-drt
3543       do 3 i=1,nu_int
3544         rti=rti+drt
3545         rtf=rti+drt
3546         if (nr_func_herm(rti)*nr_func_herm(rtf).lt.0.d0) then
3547           errabs=0.d0
3548           errrel=errt_rel
3549           maxfn=it_max
3550           xin=rti
3551           xfi=rtf
3552 C
3553           call umach(-3,9)
3554           call dzbren(nr_func_herm,errabs,errrel,xin,xfi,maxfn)
3555           call erset(0,1,1)
3556 C
3557           newrt=newrt+1
3558           roots(newrt)=dcmplx(xfi,0.d0)
3559 C
3560         end if
3561 C
3562     3 continue
3563 C
3564       nn=nr_nterms
3565       do 4 i=1,newrt
3566         call pwr_nr(nn,roots(i),data)
3567     4 continue
3568 C
3569 C
3570       return
3571       end
3572 C
3573 C---------------------------------------------------------------------
3574 C
3575 C........ the non-relativistic hermitian dispersion function .........
3576 C
3577       double precision function nr_func_herm(re_roots)
3578 C
3579       implicit double precision (a-h, o-z)
3580 C
3581       double complex dcmplx,roots
3582       double complex x11,x12,x13,x22,x23,x33
3583 C
3584       common /nrels/ nr_nterms
3585       common /wave0/ cnpa
3586       common /coutx/ x11,x12,x13,x22,x23,x33
3587 C
3588       cnpe=re_roots
3589       roots=dcmplx(cnpe,0.0d0)
3590 C
3591       call d_nr(nr_nterms,roots)
3592 C
3593 C........... hermitian and anti-hermitian tensor elements ............
3594 C
3595       d11h=1.-cnpa**2+dreal(x11)
3596       d12h=dimag(x12)
3597       d13h=cnpe*cnpa+dreal(x13)
3598       d22h=1.-cnpe**2-cnpa**2+dreal(x22)
3599       d23h=dimag(x23)
3600       d33h=1.-cnpe**2+dreal(x33)
3601 C
3602       nr_func_herm=d11h*(d22h*d33h+d23h*d23h)-
3603      %             d12h*(d23h*d13h+d12h*d33h)-
3604      %             d13h*(d12h*d23h+d22h*d13h)
3605 C
3606 C
3607       return
3608       end
3609 C
3610 C---------------------------------------------------------------------
3611 C
3612 C................. non-relativistic tensor elements ..................
3613 C
3614       subroutine pwr_nr(nn,cnpe,data)
3615 C
3616       implicit double precision (a-h, o-z)
3617 C
3618       double complex dcmplx,dconjg,cnpe
3619       double complex ci,scmu,cmu
3620       double complex d11,d12,d13,d22,d23,d33
3621       double complex d11h,d12h,d13h,d22h,d23h,d33h
3622       double complex up,um,dup,dum,uzp,uzm,duzp,duzm
3623       double complex uz2p,duz2p,duz2m,duz3p
3624       double complex ex,ey,ez
3625       double complex sx,sy,sz,pdis
3626 C
3627       double complex x11(4),x12(4),x13(4)
3628       double complex x22(4),x23(4),x33(4)
3629       double complex x11h(4),x12h(4),x13h(4)
3630       double complex x22h(4),x23h(4),x33h(4)
3631       double complex x11a(4),x12a(4),x13a(4)
3632       double complex x22a(4),x23a(4),x33a(4)
3633       double complex cbi(0:1000),dcbi(0:1000),dcbi2(0:1000)
3634       double complex zm(0:1000),zp(0:1000)
3635       double complex zargm(0:1000),zargp(0:1000)
3636       double complex dzm(0:1000),dzp(0:1000)
3637       double complex csum(6,4)
3638 C
3639       common /para1/ w0,rp,r0
3640       common /para9/ alpha,omega,vt,besk2
3641       common /wave0/ cnpa
3642       common /prnt/  kpr
3643       common /pwrdi/ d11h,d12h,d13h,d22h,d23h,d33h
3644 C
3645       ci=dcmplx(0.d0,1.d0)
3646 C
3647 C.. scmu is the sqrt (cmu) (defined to take care of sign of k_perp) ..
3648 C.............. cmu = (k_perp * v_te)**2 / omega_ce**2 ...............
3649 C............ cy0 = omega / ( sqrt(2)*abs(k_para)*v_te ) .............
3650 C
3651       scmu=cnpe*vt*omega
3652       cmu=scmu**2
3653       cy0=1./(dsqrt(2.d0)*dabs(cnpa)*vt)
3654 C
3655 C...... the number of terms to be calculated in the summations .......
3656 C....... this is based on when exp(-mu)*I_n(mu) becomes small ........
3657 C
3658       nn=idint(cdabs(cmu)/dsqrt(2.d0))*5+20
3659 C
3660 C.......... arguments of the z-functions in the summations ...........
3661 C
3662       do 1 n=0,nn
3663         zargm(n)=dcmplx(cy0*(1.-dfloat(n)/omega),0.d0)
3664         zargp(n)=dcmplx(cy0*(1.+dfloat(n)/omega),0.d0)
3665     1 continue
3666 C
3667       call ebes2(nn,cmu,cbi,dcbi,dcbi2)
3668 C
3669       call pdisp2(nn,zargm,zargp,zm,zp,dzm,dzp)
3670 C
3671 C..... the first element denotes the component of the 3x3 tensor .....
3672 C.. (1 is for xx, 2 for xy, 3 for xz, 4 for yy, 5 for yz, 6 for zz) ..
3673 C............. the second element is defined as follows ..............
3674 C.............. 1 gives the dielectric tensor elements ...............
3675 C............. 2 gives the derivative with respect to mu .............
3676 C......... 3 gives omega * derivative with respect to omega ..........
3677 C..... 4 gives k_parallel*derivative with respect to k_parallel ......
3678 C
3679       csum(1,1)=dcmplx(0.d0,0.d0)
3680       csum(1,2)=dcmplx(0.d0,0.d0)
3681       csum(1,3)=dcmplx(0.d0,0.d0)
3682       csum(1,4)=dcmplx(0.d0,0.d0)
3683 C
3684       csum(2,1)=dcmplx(0.d0,0.d0)
3685       csum(2,2)=dcmplx(0.d0,0.d0)
3686       csum(2,3)=dcmplx(0.d0,0.d0)
3687       csum(2,4)=dcmplx(0.d0,0.d0)
3688 C
3689       csum(3,1)=dcmplx(0.d0,0.d0)
3690       csum(3,2)=dcmplx(0.d0,0.d0)
3691       csum(3,3)=dcmplx(0.d0,0.d0)
3692       csum(3,4)=dcmplx(0.d0,0.d0)
3693 C
3694       csum(4,1)=2.*(cbi(0)-dcbi(0))*zp(0)
3695       csum(4,2)=2.*(dcbi(0)-dcbi2(0))*zp(0)
3696       csum(4,3)=2.*(cbi(0)-dcbi(0))*dzp(0)
3697       csum(4,4)=2.*(cbi(0)-dcbi(0))*cy0*dzp(0)
3698 C
3699       csum(5,1)=(dcbi(0)-cbi(0))*cy0*zp(0)
3700       csum(5,2)=(dcbi2(0)-dcbi(0))*cy0*zp(0)
3701       csum(5,3)=(dcbi(0)-cbi(0))*(zp(0)+cy0*dzp(0))
3702       csum(5,4)=(dcbi(0)-cbi(0))*cy0*(zp(0)+cy0*dzp(0))
3703 C
3704       csum(6,1)=cy0*(1.d0+cbi(0)*cy0*zp(0))
3705       csum(6,2)=cy0*(1.d0+dcbi(0)*cy0*zp(0))
3706       csum(6,3)=1.d0+cbi(0)*cy0*(2.d0*zp(0)+cy0*dzp(0))
3707       csum(6,4)=cy0*(-1.d0+cbi(0)*cy0**2*dzp(0))
3708 C
3709       do 2 n=1,nn
3710 C
3711         up=zm(n)+zp(n)
3712         um=zm(n)-zp(n)
3713         dup=dzm(n)+dzp(n)
3714         dum=dzm(n)-dzp(n)
3715         uzp=zargm(n)*zm(n)+zargp(n)*zp(n)
3716         uzm=zargm(n)*zm(n)-zargp(n)*zp(n)
3717         uz2p=zargm(n)**2*zm(n)+zargp(n)**2*zp(n)
3718         duzp=zargm(n)*dzm(n)+zargp(n)*dzp(n)
3719         duzm=zargm(n)*dzm(n)-zargp(n)*dzp(n)
3720         duz2p=zargm(n)**2*dzm(n)+zargp(n)**2*dzp(n)
3721         duz2m=zargm(n)**2*dzm(n)-zargp(n)**2*dzp(n)
3722         duz3p=zargm(n)**3*dzm(n)+zargp(n)**3*dzp(n)
3723 C
3724         csum(1,1)=csum(1,1)+dfloat(n**2)*cbi(n)*up
3725         csum(1,2)=csum(1,2)+dfloat(n**2)*dcbi(n)*up
3726         csum(1,3)=csum(1,3)+dfloat(n**2)*cbi(n)*dup
3727         csum(1,4)=csum(1,4)+dfloat(n**2)*cbi(n)*duzp
3728 C
3729         csum(2,1)=csum(2,1)+dfloat(n)*(dcbi(n)-cbi(n))*um
3730         csum(2,2)=csum(2,2)+dfloat(n)*(dcbi2(n)-dcbi(n))*um
3731         csum(2,3)=csum(2,3)+dfloat(n)*(dcbi(n)-cbi(n))*dum
3732         csum(2,4)=csum(2,4)+dfloat(n)*(dcbi(n)-cbi(n))*duzm
3733 C
3734         csum(3,1)=csum(3,1)+dfloat(n)*cbi(n)*uzm
3735         csum(3,2)=csum(3,2)+dfloat(n)*dcbi(n)*uzm
3736         csum(3,3)=csum(3,3)+dfloat(n)*cbi(n)*(um+duzm)
3737         csum(3,4)=csum(3,4)+dfloat(n)*cbi(n)*duz2m
3738 C
3739         csum(4,1)=csum(4,1)+
3740      %            ((2.+dfloat(n**2)/cmu**2)*cbi(n)-2.*dcbi(n))*up
3741         csum(4,2)=csum(4,2)+
3742      %            (-2.*dfloat(n**2)/cmu**3*cbi(n)+
3743      %             (2.+dfloat(n**2)/cmu**2)*dcbi(n)-2.*dcbi2(n))*up
3744         csum(4,3)=csum(4,3)+
3745      %            ((2.+dfloat(n**2)/cmu**2)*cbi(n)-2.*dcbi(n))*dup
3746         csum(4,4)=csum(4,4)+
3747      %            ((2.+dfloat(n**2)/cmu**2)*cbi(n)-2.*dcbi(n))*duzp
3748 C
3749         csum(5,1)=csum(5,1)+(dcbi(n)-cbi(n))*uzp
3750         csum(5,2)=csum(5,2)+(dcbi2(n)-dcbi(n))*uzp
3751         csum(5,3)=csum(5,3)+(dcbi(n)-cbi(n))*(up+duzp)
3752         csum(5,4)=csum(5,4)+(dcbi(n)-cbi(n))*(uzp+duz2p)
3753 C
3754         csum(6,1)=csum(6,1)+cbi(n)*uz2p
3755         csum(6,2)=csum(6,2)+dcbi(n)*uz2p
3756         csum(6,3)=csum(6,3)+cbi(n)*(2.*uzp+duz2p)
3757         csum(6,4)=csum(6,4)+cbi(n)*duz3p
3758 C
3759     2 continue
3760 C
3761 C.......................... sign of k_para ...........................
3762 C
3763       sk=cnpa/(dabs(cnpa))
3764 C
3765       fac=alpha/omega**2
3766 C
3767 C........................ the tensor elements ........................
3768 C
3769       x11(1)=fac*cy0/cmu*csum(1,1)
3770       x12(1)=-ci*fac*cy0*csum(2,1)
3771       x13(1)=sk*fac*dsqrt(2.d0)/scmu*cy0*csum(3,1)
3772       x22(1)=fac*cmu*cy0*csum(4,1)
3773       x23(1)=ci*sk*fac*dsqrt(2.d0)*scmu*cy0*csum(5,1)
3774       x33(1)=2.*fac*cy0*csum(6,1)
3775 C
3776 C....... derivatives of the tensor elements with respect to mu .......
3777 C
3778       x11(2)=-(1.+1./cmu)*x11(1)
3779      %       +fac*cy0/cmu*csum(1,2)
3780       x12(2)=-x12(1)
3781      %       -ci*fac*cy0*csum(2,2)
3782       x13(2)=-(1.+0.5/cmu)*x13(1)
3783      %       +sk*fac*dsqrt(2.d0)/scmu*cy0*csum(3,2)
3784       x22(2)= (1./cmu-1.)*x22(1)
3785      %       +fac*cmu*cy0*csum(4,2)
3786       x23(2)= (0.5/cmu-1.)*x23(1)
3787      %       +ci*sk*fac*dsqrt(2.d0)*scmu*cy0*csum(5,2)
3788       x33(2)=-x33(1)
3789      %       +2.*fac*cy0*csum(6,2)
3790 C
3791 C... omega * derivatives of tensor elements with respect to omega ....
3792 C
3793       x11(3)=-x11(1)
3794      %       +fac*cy0**2/cmu*csum(1,3)
3795       x12(3)=-x12(1)
3796      %       -ci*fac*cy0**2*csum(2,3)
3797       x13(3)=-x13(1)
3798      %       +sk*fac*dsqrt(2.d0)/scmu*cy0**2*csum(3,3)
3799       x22(3)=-x22(1)
3800      %       +fac*cmu*cy0**2*csum(4,3)
3801       x23(3)=-x23(1)
3802      %       +ci*sk*fac*dsqrt(2.d0)*scmu*cy0**2*csum(5,3)
3803       x33(3)=-x33(1)
3804      %       +2.*fac*cy0**2*csum(6,3)
3805 C
3806 C.. k_para * derivatives of tensor elements with respect to k_para ...
3807 C
3808       x11(4)=-x11(1)
3809      %       -fac*cy0/cmu*csum(1,4)
3810       x12(4)=-x12(1)
3811      %       +ci*fac*cy0*csum(2,4)
3812       x13(4)=-2.*x13(1)
3813      %       -sk*fac*dsqrt(2.d0)/scmu*cy0*csum(3,4)
3814       x22(4)=-x22(1)
3815      %       -fac*cmu*cy0*csum(4,4)
3816       x23(4)=-x23(1)
3817      %       -ci*sk*fac*dsqrt(2.d0)*scmu*cy0*csum(5,4)
3818       x33(4)=-3.*x33(1)
3819      %       -2.*fac*cy0*csum(6,4)
3820 C
3821 C.... evaluating the hermitian and anti-hermitian tensor elements ....
3822 C
3823       do 3 i=1,4
3824         x11h(i)=dcmplx(dreal(x11(i)),0.d0)
3825         x11a(i)=dcmplx(0.d0,dimag(x11(i)))
3826         x12h(i)=dcmplx(0.d0,dimag(x12(i)))
3827         x12a(i)=dcmplx(dreal(x12(i)),0.d0)
3828         x13h(i)=dcmplx(dreal(x13(i)),0.d0)
3829         x13a(i)=dcmplx(0.d0,dimag(x13(i)))
3830         x22h(i)=dcmplx(dreal(x22(i)),0.d0)
3831         x22a(i)=dcmplx(0.d0,dimag(x22(i)))
3832         x23h(i)=dcmplx(0.d0,dimag(x23(i)))
3833         x23a(i)=dcmplx(dreal(x23(i)),0.d0)
3834         x33h(i)=dcmplx(dreal(x33(i)),0.d0)
3835         x33a(i)=dcmplx(0.d0,dimag(x33(i)))
3836     3 continue
3837 C
3838 C........ the dielectric tensor elements (satisfying D.E = 0) ........
3839 C
3840       d11=1.-cnpa**2+x11(1)
3841       d12=x12(1)
3842       d13=cnpe*cnpa+x13(1)
3843       d22=1.-cnpe**2-cnpa**2+x22(1)
3844       d23=x23(1)
3845       d33=1.-cnpe**2+x33(1)
3846 C
3847 C...... dispersion tensor for the hermitian permittivity tensor ......
3848 C
3849       d11h=1.-cnpa**2+dcmplx(dreal(x11(1)),0.d0)
3850       d12h=dcmplx(0.d0,dimag(x12(1)))
3851       d13h=cnpe*cnpa+dcmplx(dreal(x13(1)),0.d0)
3852       d22h=1.-cnpe**2-cnpa**2+dcmplx(dreal(x22(1)),0.d0)
3853       d23h=dcmplx(0.d0,dimag(x23(1)))
3854       d33h=1.-cnpe**2+dcmplx(dreal(x33(1)),0.d0)
3855 C
3856 C............. evaluating the electric field components ..............
3857 C
3858       call wave_polar_pwr(cnpe,ex,ey,ez)
3859 C
3860 C.......................... speed of light ...........................
3861       c=2.99792458d10
3862 C
3863 C........... power flow in the x-direction (in cgs units) ............
3864 C
3865       sx=0.5*c*dreal(cnpe*(cdabs(ey)**2+cdabs(ez)**2)-
3866      %               cnpa*dconjg(ex)*ez)-
3867      %   0.5*c*cmu/cnpe*(cdabs(ex)**2*x11h(2)+cdabs(ey)**2*x22h(2)+
3868      %                   cdabs(ez)**2*x33h(2)+
3869      %                   2.*(ci*dimag(dconjg(ex)*ey)*x12h(2)+
3870      %                          dreal(dconjg(ex)*ez)*x13h(2)+
3871      %                       ci*dimag(dconjg(ey)*ez)*x23h(2)))
3872 C
3873 C........... power flow in the y-direction (in cgs units) ............
3874 C
3875       sy=-0.5*c*dreal(cnpa*dconjg(ey)*ez+cnpe*ex*dconjg(ey))
3876 C
3877 C........... power flow in the y-direction (in cgs units) ............
3878 C
3879       sz=0.5*c*dreal(cnpa*(cdabs(ex)**2+cdabs(ey)**2)-
3880      %               cnpe*ex*dconjg(ez))-
3881      %   0.25*c/cnpa*(cdabs(ex)**2*x11h(4)+cdabs(ey)**2*x22h(4)+
3882      %                cdabs(ez)**2*x33h(4)+
3883      %                2.*(ci*dimag(dconjg(ex)*ey)*x12h(4)+
3884      %                       dreal(dconjg(ex)*ez)*x13h(4)+
3885      %                    ci*dimag(dconjg(ey)*ez)*x23h(4)))
3886 C
3887 C.............. converting the power flow to mks units ...............
3888 C..... this gives the power flow density (in watts/meter**2) per .....
3889 C.... unit electric field amplitude squared (in {volts/meter}**2) ....
3890 C...... the conversion factor is = 1.e-11/(4*pi*2.99792458**2) .......
3891 C
3892       fact=8.854187817620391d-14
3893 C
3894       sxr=dreal(sx)
3895       sxi=dimag(sx)
3896       if (dabs(sxr).lt.1.d-20) sxr=0.d0
3897       if (dabs(sxi).lt.1.d-20) sxi=0.d0
3898       syr=dreal(sy)
3899       syi=dimag(sy)
3900       if (dabs(syr).lt.1.d-20) syr=0.d0
3901       if (dabs(syi).lt.1.d-20) syi=0.d0
3902       szr=dreal(sz)
3903       szi=dimag(sz)
3904       if (dabs(szr).lt.1.d-20) szr=0.d0
3905       if (dabs(szi).lt.1.d-20) szi=0.d0
3906 C
3907       sx=dcmplx(sxr,sxi)
3908       sy=dcmplx(syr,syi)
3909       sz=dcmplx(szr,szi)
3910       sx=sx*fact
3911       sy=sy*fact
3912       sz=sz*fact
3913 C
3914 C
3915 C......................... power disspiated ..........................
3916 C
3917       pdis=-0.5*ci*w0*(x11a(1)*cdabs(ex)**2+x22a(1)*cdabs(ey)**2+
3918      %                 x33a(1)*cdabs(ez)**2+
3919      %                 2.*ci*x12a(1)*dimag(dconjg(ex)*ey)+
3920      %                 2.*x13a(1)*dreal(dconjg(ex)*ez)+
3921      %                 2.*ci*x23a(1)*dimag(dconjg(ey)*ez))
3922 C...................... converting to mks units ......................
3923 C.... this gives the power dissipated density (in watts/meter**3) ....
3924 C.... per electric field amplitude squared (in {volts/meter}**2) .....
3925 C........ the conversion factor = 1.e-9/(4.*pi*2.99792458**2) ........
3926 C
3927       fact=8.854187817620391d-12
3928       pdisr=dreal(pdis)
3929       pdisi=dimag(pdis)
3930       if (dabs(pdisr).lt.1.d-20) pdisr=0.d0
3931       if (dabs(pdisi).lt.1.d-20) pdisi=0.d0
3932       pdis=dcmplx(pdisr,pdisi)
3933       pdis=fact*pdis
3934 C
3935       write (11,101) data,cnpe,sx,sy,sz,pdis,ex,ey,ez
3936   101 format(17g18.10)
3937 C
3938 C
3939       return
3940       end
3941 C
3942 C---------------------------------------------------------------------
3943 C
3944 C...... Bessel functions and their first and second derivatives ......
3945 C
3946       subroutine ebes2(n,carg,cbi,dcbi,dcbi2)
3947 C
3948       implicit double precision (a-h, o-z)
3949 C
3950       double complex cdexp
3951       double complex carg,factor
3952 C
3953       double complex cbi(0:1000),dcbi(0:1000),dcbi2(0:1000)
3954 C
3955       external dcbis
3956 C
3957       xnu=0.
3958       numb=n+3
3959 C
3960 C............... evaluating the bessel function itself ...............
3961 C
3962       call dcbis(xnu,carg,numb,cbi)
3963 C
3964 C......... calculating the derivates of the bessel functions .........
3965 C
3966       dcbi(0)=cbi(1)
3967       do 1 i=1,numb-1
3968         dcbi(i)=0.5*(cbi(i-1)+cbi(i+1))
3969     1 continue
3970 C
3971 C..... calculating the second derivates of the bessel functions ......
3972 C
3973       dcbi2(0)=dcbi(1)
3974       do 2 i=1,numb-2
3975         dcbi2(i)=0.5*(dcbi(i-1)+dcbi(i+1))
3976     2 continue
3977 C
3978 C.............. multiplying by exp(-mu_e) for the sums ...............
3979 C
3980       factor=cdexp(-carg)
3981       do 3 i=0,numb-1
3982         cbi(i)=cbi(i)*factor
3983         dcbi(i)=dcbi(i)*factor
3984         dcbi2(i)=dcbi2(i)*factor
3985     3 continue
3986 C
3987 C
3988       return
3989       end
3990 C
3991 C---------------------------------------------------------------------
3992 C
3993 C....... the plasma dispersion functions and their derivatives .......
3994 C
3995       subroutine pdisp2(nar,cargm,cargp,zm,zp,dzm,dzp)
3996 C
3997       implicit double precision (a-h, o-z)
3998 C
3999       double complex dcmplx
4000       double complex carg,csqpi,zerfe
4001 C
4002       double complex cargm(0:1000),cargp(0:1000)
4003       double complex zm(0:1000),zp(0:1000)
4004       double complex dzm(0:1000),dzp(0:1000)
4005 C
4006       external zerfe
4007 C
4008       sqpi=1.772453850905516027298167d0
4009         pi=3.141592653589793238462643d0
4010 C
4011       csqpi=dcmplx(0.d0,sqpi)
4012 C     
4013       do 1 i=0,nar
4014         zm(i)=csqpi*zerfe(cargm(i))
4015         zp(i)=csqpi*zerfe(cargp(i))
4016     1 continue
4017 C
4018 C......... this is d/dy[Z(y)] for y being y_plus and y_minus .........
4019 C
4020       do 2 i=0,nar
4021         dzm(i)=-2.*(1.+cargm(i)*zm(i))
4022         dzp(i)=-2.*(1.+cargp(i)*zp(i))
4023     2 continue
4024 C
4025 C
4026         return
4027         end
4028 C
4029 C---------------------------------------------------------------------
4030 C
4031 C........... calculating polarizations with respect to E_x ...........
4032 C
4033       subroutine wave_polar_pwr(cnpe,ex,ey,ez)
4034 C
4035       implicit double precision (a-h, o-z)
4036 C
4037       double complex dcmplx,ci,cnpe,cdummy
4038       double complex d11,d12,d13,d22,d23,d33
4039       double complex ex,ey,ez
4040       double complex pol_p,pol_m,pol_z
4041 C
4042       double complex det(3),pdum(3)
4043 C
4044       common /wave0/ cnpa
4045       common /pwrdi/ d11,d12,d13,d22,d23,d33
4046 C
4047       ci=dcmplx(0.d0,1.d0)
4048 C
4049       det(1)=-d12*d23-d13*d22
4050       det(2)= d22*d33-d23*d23
4051       det(3)= d12*d33+d13*d23
4052 C
4053       if (cdabs(det(1)).gt.cdabs(det(2))) then
4054         if (cdabs(det(1)).gt.cdabs(det(3))) then
4055           nn=1
4056         else
4057           nn=3
4058         end if
4059       else
4060         if (cdabs(det(2)).gt.cdabs(det(3))) then
4061           nn=2
4062         else
4063           nn=3
4064         end if
4065       end if
4066 C
4067       if (nn.eq.1) then
4068         ey=(d11*d23-d13*d12)/det(1)
4069         ez=(d12*d12+d11*d22)/det(1)
4070       end if
4071       if (nn.eq.2) then
4072         ey=(d12*d33-d13*d23)/det(2) 
4073         ez=(d12*d23-d13*d22)/det(2)
4074       end if
4075       if (nn.eq.3) then
4076         ey= (d13*d13-d11*d33)/det(3)
4077         ez=-(d12*d13+d11*d23)/det(3)
4078       end if
4079 C
4080 C....... zeroing out the small components of the polarizations .......
4081 C
4082       pdum(1)=(1.d0+ci*ey)/dsqrt(2.d0)
4083       pdum(2)=(1.d0-ci*ey)/dsqrt(2.d0)
4084       pdum(3)=ez
4085 C
4086       do 1 i=1,3
4087         if (dabs(dreal(pdum(i))).lt.1.0d-15) 
4088      %      pdum(i)=dcmplx(0.d0,dimag(pdum(i)))
4089         if (dabs(dimag(pdum(i))).lt.1.0d-15) 
4090      %      pdum(i)=dcmplx(dreal(pdum(i)),0.d0)
4091     1 continue
4092 C
4093       ex=dcmplx(1.0d0,0.d0)
4094       ey=ci*(dsqrt(2.d0)*pdum(2)-1.d0)
4095       ez=pdum(3)
4096 C
4097 C
4098       return
4099       end
4100 C
4101 C---------------------------------------------------------------------
4102 C
4103 C
4104 C---------------------------------------------------------------------
4105 C=====================================================================
4106 C=====================================================================
4107 C
4108 C...... solving the hermitian relativistic dispersion relation .......
4109 C
4110       subroutine fr_herm_disp(data,rguess,roots,nrts)
4111 C
4112       implicit double precision (a-h,o-z)
4113 C
4114       parameter (n_rts=10)
4115 C
4116       integer info(n_rts)
4117 C
4118       double complex fr_herm_func
4119       double complex rguess(n_rts),roots(n_rts)
4120 C
4121       common /rterr/ errt_abs,errt_rel,it_max
4122       common /prnt/  kpr
4123 C
4124       external fr_herm_func,dzanly,umach,erset
4125 C
4126       errabs=errt_abs
4127       errrel=errt_rel
4128       itmax=it_max
4129       nknown=0
4130       nnew=nrts
4131       nguess=nrts
4132 C
4133       call umach(-3,9)
4134       call dzanly(fr_herm_func,errabs,errrel,nknown,nnew,nguess,
4135      %            rguess,itmax,roots,info)
4136       call erset(0,1,1)
4137 C
4138       call order_roots(roots,nrts)
4139 C
4140 CCCCCCCCCCCCCCCCCCC
4141 CC      do 421 j=1,nrts
4142 CC        write (10,422) roots(j)
4143 CC  422   format("hermitian root = ",2g18.10)
4144 CC  421 continue
4145 CCCCCCCCCCCCCCCCCCC
4146       do 1 i=1,nrts
4147         call pwr_rel(roots(i),data)
4148     1 continue
4149 C
4150 C
4151       return
4152       end
4153 C
4154 C---------------------------------------------------------------------
4155 C
4156 C............ the fully-relativistic dispersion function .............
4157 C
4158       double complex function fr_herm_func(roots)
4159 C
4160       implicit double precision (a-h, o-z)
4161 C
4162       double complex roots,cnpe
4163       double complex dh11,dh12,dh13,dh22,dh23,dh33
4164       double complex xh11,xh12,xh13,xh22,xh23,xh33
4165 C
4166       double complex ep(6)
4167 C
4168       common /wave0/ cnpa
4169       common /wave1/ cnpe
4170 C
4171       cnpe=roots
4172 C
4173       kherm=1
4174       call trub_herm(kherm,ep)
4175 C
4176       xh11=ep(1)
4177       xh12=ep(2)
4178       xh13=ep(3)
4179       xh22=ep(4)
4180       xh23=ep(5)
4181       xh33=ep(6)
4182 C
4183 C........ the dielectric tensor elements (satisfying D.E = 0) ........
4184 C
4185       dh11=1.-cnpa**2+xh11
4186       dh12=xh12
4187       dh13=cnpe*cnpa+xh13
4188       dh22=1.-cnpe**2-cnpa**2+xh22
4189       dh23=xh23
4190       dh33=1.-cnpe**2+xh33
4191 C
4192       fr_herm_func=dh11*(dh22*dh33+dh23*dh23)+
4193      %             dh12*(dh23*dh13+dh12*dh33)+
4194      %             dh13*(dh12*dh23-dh22*dh13)
4195 C
4196 C
4197       return
4198       end
4199 C
4200 C
4201 C---------------------------------------------------------------------
4202 C
4203 C..... evaluating the trubnikov form of the relativistic tensor ......
4204 C
4205       subroutine trub_herm(kherm,epsilon)
4206 C
4207       implicit double precision (a-h, o-z)
4208 C
4209       double precision stot(3),sumr(0:1000)
4210 C
4211       double complex dcmplx
4212 C
4213       double complex epsilon(6),co(0:1)
4214 C
4215       common /trub1/ xi_lo,dxi,nxi
4216       common /trub2/ diff_err,navg
4217       common /trerr/ errabs0,errrel0,irule0
4218       common /trubh/ kel,kreal
4219       common /prnt/  kpr
4220 C
4221       external dqdag,ftrub_herm,umach,erset
4222 C
4223       co(0)=dcmplx(0.d0,1.d0)
4224       co(1)=dcmplx(1.d0,0.d0)
4225 C
4226       do 1 j=1,6
4227         kel=j
4228 C
4229         if (kherm.eq.1) then
4230           if (j.eq.1.or.j.eq.3.or.j.eq.4.or.j.eq.6) kreal=1
4231           if (j.eq.2.or.j.eq.5) kreal=0
4232         else
4233           if (j.eq.1.or.j.eq.3.or.j.eq.4.or.j.eq.6) kreal=0
4234           if (j.eq.2.or.j.eq.5) kreal=1
4235         end if
4236 C
4237         k=kreal
4238 C
4239 C.. ntr keeps track of the number of averages (3) that are compared ..
4240 C
4241         ntr=0
4242 C
4243         xi_0=xi_lo
4244         xi_1=xi_0+dxi
4245 C
4246         sum=0.d0
4247 CCCCCCCCCCCCCCCCCCCCCCC
4248 CC        ncount=0
4249 CCCCCCCCCCCCCCCCCCCCCCC
4250 C
4251         do 2 i=1,nxi
4252 C
4253 CCCCCCCCCCCCCCCCCCCCCCC
4254 CC          ncount=ncount+1
4255 CCCCCCCCCCCCCCCCCCCCCCC
4256           errabs=errabs0
4257           errrel=errrel0
4258           irule=irule0
4259 C
4260           call umach(-3,9)
4261           call dqdag(ftrub_herm,xi_0,xi_1,errabs,errrel,irule,
4262      %               t_real,errest)
4263           call erset(0,1,1)
4264 C
4265           sum=sum+t_real
4266 C
4267           jj=mod(i-1,navg)
4268           sumr(jj)=sum
4269           if (jj.eq.navg-1) then
4270             ntr=ntr+1
4271             stot(ntr)=0.d0
4272             do 3 ij=0,navg-1
4273               stot(ntr)=stot(ntr)+sumr(ij)
4274     3       continue
4275             stot(ntr)=stot(ntr)/dfloat(navg)
4276             if (ntr.eq.3) then
4277               diff1=dabs(dabs(stot(1))-dabs(stot(2)))
4278               diff2=dabs(dabs(stot(2))-dabs(stot(3)))
4279               if (diff1.lt.diff_err.and.diff2.lt.diff_err) then
4280                 sum_real=stot(3)
4281                 go to 4
4282               else
4283                 stot(1)=stot(2)
4284                 stot(2)=stot(3)
4285                 ntr=2
4286               end if
4287             end if
4288           end if
4289 C
4290           xi_0=xi_1
4291           xi_1=xi_0+dxi
4292 C
4293     2   continue
4294 C
4295     4   continue
4296 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4297 CC        write (10,389) kel,ncount
4298 CC  389   format("inside trub_herm: element, int_count = ",2i8)
4299 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4300 C
4301         epsilon(kel)=co(k)*sum
4302 C
4303     1 continue
4304 C
4305 C
4306       return
4307       end
4308 C
4309 C---------------------------------------------------------------------
4310 C
4311       double precision function ftrub_herm(xi)
4312 C
4313       implicit double precision (a-h, o-z)
4314 C
4315       double complex dcmplx,ci,factor
4316       double complex tx_11,tx_12,tx_13
4317       double complex tx_22,tx_23,tx_33
4318 C
4319       common /para9/ alpha,omega,vt,besk2
4320       common /trubh/ kel,kreal
4321 C
4322       external tx_11,tx_12,tx_13
4323       external tx_22,tx_23,tx_33
4324 C
4325       ci=dcmplx(0.d0,1.d0)
4326 C
4327 C...... factor is the multiplier in the dielectric tensor form .......
4328 C.............. factor = i * {(wpe/wce)**2} * (wce/w0) ...............
4329 C
4330       factor=ci*alpha/omega
4331 C
4332       if (kreal.eq.0) then
4333         if (kel.eq.1) ftrub_herm=dimag(factor*tx_11(xi))
4334         if (kel.eq.2) ftrub_herm=dimag(factor*tx_12(xi))
4335         if (kel.eq.3) ftrub_herm=dimag(factor*tx_13(xi))
4336         if (kel.eq.4) ftrub_herm=dimag(factor*tx_22(xi))
4337         if (kel.eq.5) ftrub_herm=dimag(factor*tx_23(xi))
4338         if (kel.eq.6) ftrub_herm=dimag(factor*tx_33(xi))
4339       else
4340         if (kel.eq.1) ftrub_herm=dreal(factor*tx_11(xi))
4341         if (kel.eq.2) ftrub_herm=dreal(factor*tx_12(xi))
4342         if (kel.eq.3) ftrub_herm=dreal(factor*tx_13(xi))
4343         if (kel.eq.4) ftrub_herm=dreal(factor*tx_22(xi))
4344         if (kel.eq.5) ftrub_herm=dreal(factor*tx_23(xi))
4345         if (kel.eq.6) ftrub_herm=dreal(factor*tx_33(xi))
4346       end if
4347 C
4348 C
4349       return
4350       end
4351 C
4352 C
4353 C=====================================================================
4354 C
4355 C...... subroutines for calculating various derivatives of the .......
4356 C... hermitian (and anti-hermitian) relativistic dielectric tensor ...
4357 C
4358 C..... the derivatives are x(i)*[d/dx(i)] where i is as follows ......
4359 C..... i = 1 is w (frequency), i = 2 is k_perp, i = 3 is k_para ......
4360 C...... i = 4 is wc (cyclotron), i = 5 is vt (thermal velocity) ......
4361 C.................. i = 6 is wp (plasma frequency) ...................
4362 C
4363 C
4364 C.... working with the trubnikov form of the relativistic tensor .....
4365 C
4366 C..... here kdiff plays the same role as i in x(i) listed above ......
4367 C
4368       subroutine d_drel(kherm,kdiff,eps,deps)
4369 C
4370       implicit double precision (a-h, o-z)
4371 C
4372       double complex dcmplx,cz
4373 C
4374       double complex eps(6),deps(6),d_eps(6),dg(6)
4375 C
4376       common /para9/ alpha,omega,vt,besk2
4377 C
4378       external dbsk0e,dbsk1e,besk_a,umach,erset
4379 C
4380 C..... dg is the derivative of the term multiplying the integral .....
4381 C
4382       cz=dcmplx(0.d0,0.d0)
4383       do 1 i=1,6
4384         dg(i)=cz
4385     1 continue
4386 C
4387       if (kdiff.eq.6) then
4388         do 2 i=1,6
4389           d_eps(i)=cz
4390     2  continue
4391       else
4392         call d_trub(kherm,kdiff,d_eps)
4393       end if
4394 C
4395 C.......... g has derivatives w.r.t. w, wc, vt, and wp only ..........
4396 C
4397       if (kdiff.eq.1.or.kdiff.eq.4) then
4398         do 3 i=1,6
4399           dg(i)=-eps(i)
4400     3   continue
4401       end if
4402 C
4403       if (kdiff.eq.5) then
4404         arg=1./vt**2
4405         if (arg.lt.50.d0) then
4406           call umach(-3,9)
4407           bk2=dbsk0e(arg)+2.*dbsk1e(arg)/arg
4408           bk1=dbsk1e(arg)
4409           call erset(0,1,1)
4410         else
4411           bk1=besk_a(1,arg)
4412           bk2=besk_a(2,arg)
4413         end if
4414         fac=-2.*(4.+arg*bk1/bk2)
4415         do 4 i=1,6
4416           dg(i)=fac*eps(i)
4417     4   continue
4418       end if
4419 C
4420       if (kdiff.eq.6) then
4421         do 5 i=1,6
4422           dg(i)=2.*eps(i)
4423     5   continue
4424       end if
4425 C
4426 C.................. the derivatives of all the xi's ..................
4427 C
4428       do 6 i=1,6
4429         deps(i)=dg(i)+d_eps(i)
4430     6 continue
4431 C
4432 C
4433       return
4434       end
4435 C
4436 C---------------------------------------------------------------------
4437 C
4438       subroutine d_trub(kherm,kdiff,d_eps)
4439 C
4440       implicit double precision (a-h, o-z)
4441 C
4442       double precision stot(3),sumr(0:1000)
4443 C
4444       double complex dcmplx
4445 C
4446       double complex d_eps(6),co(0:1)
4447 C
4448       common /trub1/ xi_lo,dxi,nxi
4449       common /trub2/ diff_err,navg
4450       common /trerr/ errabs0,errrel0,irule0
4451       common /trubp/ kel,kreal
4452       common /trubd/ nder
4453       common /prnt/  kpr
4454 C
4455       external dqdag,d_ftrub,umach,erset
4456 C
4457 C............. kdiff is the particular derivative needed .............
4458 C.... kdiff = i calculates x(i)*[d/dx(i)] of the tensor elements .....
4459 C
4460       nder=kdiff
4461 C
4462       co(0)=dcmplx(0.d0,1.d0)
4463       co(1)=dcmplx(1.d0,0.d0)
4464 C
4465       do 1 j=1,6
4466         kel=j
4467 C
4468         if (kherm.eq.1) then
4469           if (j.eq.1.or.j.eq.3.or.j.eq.4.or.j.eq.6) kreal=1
4470           if (j.eq.2.or.j.eq.5) kreal=0
4471         else
4472           if (j.eq.1.or.j.eq.3.or.j.eq.4.or.j.eq.6) kreal=0
4473           if (j.eq.2.or.j.eq.5) kreal=1
4474         end if
4475 C
4476         k=kreal
4477 C
4478 C.. ntr keeps track of the number of averages (3) that are compared ..
4479 C
4480         ntr=0
4481 C
4482         xi_0=xi_lo
4483         xi_1=xi_0+dxi
4484 C
4485         sum=0.d0
4486 C
4487 CCCCCCCCCCCCCCCCCCCCC
4488 CC        write (10,433)
4489 CC  433   format(1x,"inside d_trub")
4490 CC        write (10,444) xi_0,xi_1,nder
4491 CC  444   format("xi_0, xi_1, nder = ",2g18.10,5x,i4)
4492 CC        ncount=0
4493 CCCCCCCCCCCCCCCCCCCCC
4494         do 2 i=1,nxi
4495 C
4496 CCCCCCCCCCCCCCCCCCCCC
4497 CC          ncount=ncount+1
4498 CCCCCCCCCCCCCCCCCCCCC
4499           errabs=errabs0
4500           errrel=errrel0
4501           irule=irule0
4502 C
4503           call umach(-3,9)
4504           call dqdag(d_ftrub,xi_0,xi_1,errabs,errrel,irule,
4505      %               t_real,errest)
4506           call erset(0,1,1)
4507 C
4508           sum=sum+t_real
4509 C
4510           jj=mod(i-1,navg)
4511           sumr(jj)=sum
4512           if (jj.eq.navg-1) then
4513             ntr=ntr+1
4514             stot(ntr)=0.d0
4515             do 3 ij=0,navg-1
4516               stot(ntr)=stot(ntr)+sumr(ij)
4517     3       continue
4518             stot(ntr)=stot(ntr)/dfloat(navg)
4519             if (ntr.eq.3) then
4520               diff1=dabs(dabs(stot(1))-dabs(stot(2)))
4521               diff2=dabs(dabs(stot(2))-dabs(stot(3)))
4522 CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4523 CC              write (10,520) kel
4524 CC  520         format("d_trub element = ",i2)
4525 CC              write (10,519) stot(1),stot(2),stot(3)
4526 CC  519         format("s1,s2,s3 = ",3g18.10)
4527 CC              write (10,521) (stot(1)+stot(2))/2.,
4528 CC     %         (stot(2)+stot(3))/2.
4529 CC  521         format("s12,s23 = ",2g18.10)
4530 CC              s12=(stot(1)+stot(2))/2.
4531 CC              s23=(stot(2)+stot(3))/2.
4532 CC              write (10,522) diff1,diff2,dabs(s12-s23)
4533 CC  522         format("diff1,diff2,diff12=",3g18.10)
4534 CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4535               if (diff1.lt.diff_err.and.diff2.lt.diff_err) then
4536                 sum_real=stot(3)
4537                 go to 4
4538               else
4539                 stot(1)=stot(2)
4540                 stot(2)=stot(3)
4541                 ntr=2
4542               end if
4543             end if
4544           end if
4545 C
4546           xi_0=xi_1
4547           xi_1=xi_0+dxi
4548 CCCCCCCCCCCCCCCCCCC
4549 CC          if (kel.eq.6) then
4550 CC            write (10,589) i,xi_0-dxi,xi_0,co(k)*sum
4551 CC  589       format("d_trub: i,xi0,xi1,sum ",i8,4g18.10)
4552 CC          end if
4553 CCCCCCCCCCCCCCCCCCC
4554 C
4555     2   continue
4556 C
4557     4   continue
4558 C
4559         d_eps(kel)=co(k)*sum
4560 CCCCCCCCCCCCCCCCCCC
4561 CC        write (10,389) kel,ncount
4562 CC  389   format("inside d_trub: element, int_count = ",2i8)
4563 CC        write (10,434) j,kel,d_eps(kel)
4564 CC  434   format("j, kel, d_eps = ",2i4,2g18.10)
4565 CCCCCCCCCCCCCCCCCCC
4566 C
4567     1 continue
4568 C
4569 C
4570       return
4571       end
4572 C
4573 C---------------------------------------------------------------------
4574 C
4575       double precision function d_ftrub(xi)
4576 C
4577       implicit double precision (a-h, o-z)
4578 C
4579       double complex dcmplx,ci,factor
4580       double complex dtx_11,dtx_12,dtx_13
4581       double complex dtx_22,dtx_23,dtx_33
4582 C
4583       common /para9/ alpha,omega,vt,besk2
4584       common /trubp/ kel,kreal
4585 C
4586       external dtx_11,dtx_12,dtx_13
4587       external dtx_22,dtx_23,dtx_33
4588 C
4589       ci=dcmplx(0.d0,1.d0)
4590 C
4591 C...... factor is the multiplier in the dielectric tensor form .......
4592 C.............. factor = i * {(wpe/wce)**2} * (wce/w0) ...............
4593 C
4594       factor=ci*alpha/omega
4595 C
4596       if (kreal.eq.0) then
4597         if (kel.eq.1) d_ftrub=dimag(factor*dtx_11(xi))
4598         if (kel.eq.2) d_ftrub=dimag(factor*dtx_12(xi))
4599         if (kel.eq.3) d_ftrub=dimag(factor*dtx_13(xi))
4600         if (kel.eq.4) d_ftrub=dimag(factor*dtx_22(xi))
4601         if (kel.eq.5) d_ftrub=dimag(factor*dtx_23(xi))
4602         if (kel.eq.6) d_ftrub=dimag(factor*dtx_33(xi))
4603       else
4604         if (kel.eq.1) d_ftrub=dreal(factor*dtx_11(xi))
4605         if (kel.eq.2) d_ftrub=dreal(factor*dtx_12(xi))
4606         if (kel.eq.3) d_ftrub=dreal(factor*dtx_13(xi))
4607         if (kel.eq.4) d_ftrub=dreal(factor*dtx_22(xi))
4608         if (kel.eq.5) d_ftrub=dreal(factor*dtx_23(xi))
4609         if (kel.eq.6) d_ftrub=dreal(factor*dtx_33(xi))
4610       end if
4611 C
4612 C
4613       return
4614       end
4615 C
4616 C---------------------------------------------------------------------
4617 C
4618       double complex function dtx_11(xi)
4619 C
4620       implicit double precision (a-h, o-z)
4621 C
4622       double complex dcmplx,cdexp,cdsqrt
4623       double complex cnpe
4624       double complex ci,cz,r,sr,t2
4625       double complex ck3,ck4,dr,dt2
4626       double complex cm1,cm2
4627       double complex ce,cfac,cbesk_a
4628 C
4629       double complex ckbs(2)
4630 C
4631       common /para9/ alpha,omega,vt,besk2
4632       common /wave0/ cnpa
4633       common /wave1/ cnpe
4634       common /trubd/ nder
4635 C
4636       external dcbks,umach,erset
4637 C
4638       ci=dcmplx(0.d0,1.d0)
4639       cz=dcmplx(0.d0,0.d0)
4640 C
4641       twopi=6.283185307179586476925287
4642 C
4643 C..... xim is xi modulo two*pi (needed in the sines and cosines) .....
4644 C
4645       xim=dmod(xi,twopi)
4646 C
4647 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
4648 C.... this takes care of the normalization factors in the tensor .....
4649 C... sr is the square-root of r (argument of the Bessel functions) ...
4650 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
4651 C
4652       dom=omega*vt**2
4653       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
4654      %  (cnpa*dom*xi)**2
4655       sr=cdsqrt(r)/vt**2
4656       r=r*besk2
4657 C
4658       if (dreal(sr).lt.50.d0) then
4659 C
4660         numb=2
4661         xnu=3.
4662         call umach(-3,9)
4663         call dcbks(xnu,sr,numb,ckbs)
4664         call erset(0,1,1)
4665 C
4666         ck3=ckbs(1)*dexp(1./vt**2)
4667         ck4=ckbs(2)*dexp(1./vt**2)
4668 C
4669       else
4670 C
4671         ce=1./vt**2-sr
4672         cfac=cdexp(ce)
4673         ck3=cbesk_a(3,sr)*cfac
4674         ck4=cbesk_a(4,sr)*cfac
4675 C
4676       end if
4677 C
4678       cm1=-ck3/(2.*sr*r)
4679       cm2= ck4/(2.*sr**2*r)
4680 C
4681       t1=dcos(xim)
4682       t2=(cnpe*omega*dsin(xim))**2
4683 C
4684 C.... derivatives of r (un-normalized) w.r.t. various parameters .....
4685 C.. derivatives of t2 and tensor element w.r.t. various parameters ...
4686 C
4687       if (nder.eq.1) then
4688         dr=-2.*ci*xi*omega*(1./vt**2-ci*xi*omega)
4689         dtx_11=(cm1*t1+cm2*t2)*dr
4690       end if
4691       if (nder.eq.2) then
4692         dr=4.*(cnpe*omega)**2*(1.-dcos(xim))
4693         dt2=2.*t2
4694         dtx_11=cm1*(t1*dr+2.*dt2)+cm2*t2*dr
4695       end if
4696       if (nder.eq.3) then
4697         dr=2.*(cnpa*omega*xi)**2
4698         dtx_11=(cm1*t1+cm2*t2)*dr
4699       end if
4700       if (nder.eq.4) then
4701         dr=2.*ci*xi*omega*(1./vt**2-ci*xi*omega)-
4702      %     4.*(cnpe*omega)**2*(1.-dcos(xim))-
4703      %     2.*(cnpa*omega*xi)**2
4704         dt2=-2.*t2
4705         dtx_11=cm1*(t1*dr+2.*dt2)+cm2*t2*dr
4706       end if
4707       if (nder.eq.5) then
4708         dr=-(2./vt)**2*(1./vt**2-ci*xi*omega)
4709         dtx_11=(cm1*t1+cm2*t2)*dr
4710       end if
4711 C
4712 C
4713       return
4714       end
4715 C
4716 C---------------------------------------------------------------------
4717 C
4718       double complex function dtx_12(xi)
4719 C
4720       implicit double precision (a-h, o-z)
4721 C
4722       double complex dcmplx,cdexp,cdsqrt
4723       double complex cnpe
4724       double complex ci,cz,r,sr,t2
4725       double complex ck3,ck4,dr,dt2
4726       double complex cm1,cm2
4727       double complex ce,cfac,cbesk_a
4728 C
4729       double complex ckbs(2)
4730 C
4731       common /para9/ alpha,omega,vt,besk2
4732       common /wave0/ cnpa
4733       common /wave1/ cnpe
4734       common /trubd/ nder
4735 C
4736       external dcbks,umach,erset
4737 C
4738       ci=dcmplx(0.d0,1.d0)
4739       cz=dcmplx(0.d0,0.d0)
4740 C
4741       twopi=6.283185307179586476925287
4742 C
4743 C..... xim is xi modulo two*pi (needed in the sines and cosines) .....
4744 C
4745       xim=dmod(xi,twopi)
4746 C
4747 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
4748 C.... this takes care of the normalization factors in the tensor .....
4749 C... sr is the square-root of r (argument of the Bessel functions) ...
4750 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
4751 C
4752       dom=omega*vt**2
4753       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
4754      %  (cnpa*dom*xi)**2
4755       sr=cdsqrt(r)/vt**2
4756       r=r*besk2
4757 C
4758       if (dreal(sr).lt.50.d0) then
4759 C
4760         numb=2
4761         xnu=3.
4762         call umach(-3,9)
4763         call dcbks(xnu,sr,numb,ckbs)
4764         call erset(0,1,1)
4765 C
4766         ck3=ckbs(1)*dexp(1./vt**2)
4767         ck4=ckbs(2)*dexp(1./vt**2)
4768 C
4769       else
4770 C
4771         ce=1./vt**2-sr
4772         cfac=cdexp(ce)
4773         ck3=cbesk_a(3,sr)*cfac
4774         ck4=cbesk_a(4,sr)*cfac
4775 C
4776       end if
4777 C
4778       cm1=-ck3/(2.*sr*r)
4779       cm2= ck4/(2.*sr**2*r)
4780 C
4781       t1=-dsin(xim)
4782       t2=-(cnpe*omega)**2*dsin(xim)*(1.d0-dcos(xim))
4783 C
4784 C.... derivatives of r (un-normalized) w.r.t. various parameters .....
4785 C.. derivatives of t2 and tensor element w.r.t. various parameters ...
4786 C
4787       if (nder.eq.1) then
4788         dr=-2.*ci*xi*omega*(1./vt**2-ci*xi*omega)
4789         dtx_12=(cm1*t1+cm2*t2)*dr
4790       end if
4791       if (nder.eq.2) then
4792         dr=4.*(cnpe*omega)**2*(1.-dcos(xim))
4793         dt2=2.*t2
4794         dtx_12=cm1*(t1*dr+2.*dt2)+cm2*t2*dr
4795       end if
4796       if (nder.eq.3) then
4797         dr=2.*(cnpa*omega*xi)**2
4798         dtx_12=(cm1*t1+cm2*t2)*dr
4799       end if
4800       if (nder.eq.4) then
4801         dr=2.*ci*xi*omega*(1./vt**2-ci*xi*omega)-
4802      %     4.*(cnpe*omega)**2*(1.-dcos(xim))-
4803      %     2.*(cnpa*omega*xi)**2
4804         dt2=-2.*t2
4805         dtx_12=cm1*(t1*dr+2.*dt2)+cm2*t2*dr
4806       end if
4807       if (nder.eq.5) then
4808         dr=-(2./vt)**2*(1./vt**2-ci*xi*omega)
4809         dtx_12=(cm1*t1+cm2*t2)*dr
4810       end if
4811 C
4812 C
4813       return
4814       end
4815 C
4816 C---------------------------------------------------------------------
4817 C
4818       double complex function dtx_13(xi)
4819 C
4820       implicit double precision (a-h, o-z)
4821 C
4822       double complex dcmplx,cdexp,cdsqrt
4823       double complex cnpe
4824       double complex ci,cz,r,sr,t2
4825       double complex ck3,ck4,dr,dt2
4826       double complex cm1,cm2
4827       double complex ce,cfac,cbesk_a
4828 C
4829       double complex ckbs(2)
4830 C
4831       common /para9/ alpha,omega,vt,besk2
4832       common /wave0/ cnpa
4833       common /wave1/ cnpe
4834       common /trubd/ nder
4835 C
4836       external dcbks,umach,erset
4837 C
4838       ci=dcmplx(0.d0,1.d0)
4839       cz=dcmplx(0.d0,0.d0)
4840 C
4841       twopi=6.283185307179586476925287
4842 C
4843 C..... xim is xi modulo two*pi (needed in the sines and cosines) .....
4844 C
4845       xim=dmod(xi,twopi)
4846 C
4847 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
4848 C.... this takes care of the normalization factors in the tensor .....
4849 C... sr is the square-root of r (argument of the Bessel functions) ...
4850 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
4851 C
4852       dom=omega*vt**2
4853       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
4854      %  (cnpa*dom*xi)**2
4855       sr=cdsqrt(r)/vt**2
4856       r=r*besk2
4857 C
4858       if (dreal(sr).lt.50.d0) then
4859 C
4860         numb=2
4861         xnu=3.
4862         call umach(-3,9)
4863         call dcbks(xnu,sr,numb,ckbs)
4864         call erset(0,1,1)
4865 C
4866         ck3=ckbs(1)*dexp(1./vt**2)
4867         ck4=ckbs(2)*dexp(1./vt**2)
4868 C
4869       else
4870 C
4871         ce=1./vt**2-sr
4872         cfac=cdexp(ce)
4873         ck3=cbesk_a(3,sr)*cfac
4874         ck4=cbesk_a(4,sr)*cfac
4875 C
4876       end if
4877 C
4878       cm1=-ck3/(2.*sr*r)
4879       cm2= ck4/(2.*sr**2*r)
4880 C
4881       t2=cnpe*cnpa*omega**2*xi*dsin(xim)
4882 C
4883 C.... derivatives of r (un-normalized) w.r.t. various parameters .....
4884 C.. derivatives of t2 and tensor element w.r.t. various parameters ...
4885 C
4886       if (nder.eq.1) then
4887         dr=-2.*ci*xi*omega*(1./vt**2-ci*xi*omega)
4888         dtx_13=cm2*t2*dr
4889       end if
4890       if (nder.eq.2) then
4891         dr=4.*(cnpe*omega)**2*(1.-dcos(xim))
4892         dt2=t2
4893         dtx_13=2.*cm1*dt2+cm2*t2*dr
4894       end if
4895       if (nder.eq.3) then
4896         dr=2.*(cnpa*omega*xi)**2
4897         dt2=t2
4898         dtx_13=2.*cm1*dt2+cm2*t2*dr
4899       end if
4900       if (nder.eq.4) then
4901         dr=2.*ci*xi*omega*(1./vt**2-ci*xi*omega)-
4902      %     4.*(cnpe*omega)**2*(1.-dcos(xim))-
4903      %     2.*(cnpa*omega*xi)**2
4904         dt2=-2.*t2
4905         dtx_13=2.*cm1*dt2+cm2*t2*dr
4906       end if
4907       if (nder.eq.5) then
4908         dr=-(2./vt)**2*(1./vt**2-ci*xi*omega)
4909         dtx_13=cm2*t2*dr
4910       end if
4911 C
4912 C
4913       return
4914       end
4915 C
4916 C---------------------------------------------------------------------
4917 C
4918       double complex function dtx_22(xi)
4919 C
4920       implicit double precision (a-h, o-z)
4921 C
4922       double complex dcmplx,cdexp,cdsqrt
4923       double complex cnpe
4924       double complex ci,cz,r,sr,t2
4925       double complex ck3,ck4,dr,dt2
4926       double complex cm1,cm2
4927       double complex ce,cfac,cbesk_a
4928 C
4929       double complex ckbs(2)
4930 C
4931       common /para9/ alpha,omega,vt,besk2
4932       common /wave0/ cnpa
4933       common /wave1/ cnpe
4934       common /trubd/ nder
4935 C
4936       external dcbks,umach,erset
4937 C
4938       ci=dcmplx(0.d0,1.d0)
4939       cz=dcmplx(0.d0,0.d0)
4940 C
4941       twopi=6.283185307179586476925287
4942 C
4943 C..... xim is xi modulo two*pi (needed in the sines and cosines) .....
4944 C
4945       xim=dmod(xi,twopi)
4946 C
4947 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
4948 C.... this takes care of the normalization factors in the tensor .....
4949 C... sr is the square-root of r (argument of the Bessel functions) ...
4950 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
4951 C
4952       dom=omega*vt**2
4953       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
4954      %  (cnpa*dom*xi)**2
4955       sr=cdsqrt(r)/vt**2
4956       r=r*besk2
4957 C
4958       if (dreal(sr).lt.50.d0) then
4959 C
4960         numb=2
4961         xnu=3.
4962         call umach(-3,9)
4963         call dcbks(xnu,sr,numb,ckbs)
4964         call erset(0,1,1)
4965 C
4966         ck3=ckbs(1)*dexp(1./vt**2)
4967         ck4=ckbs(2)*dexp(1./vt**2)
4968 C
4969       else
4970 C
4971         ce=1./vt**2-sr
4972         cfac=cdexp(ce)
4973         ck3=cbesk_a(3,sr)*cfac
4974         ck4=cbesk_a(4,sr)*cfac
4975 C
4976       end if
4977 C
4978       cm1=-ck3/(2.*sr*r)
4979       cm2= ck4/(2.*sr**2*r)
4980 C
4981       t1=dcos(xim)
4982       t2=-(cnpe*omega*(1.d0-dcos(xim)))**2
4983 C
4984 C.... derivatives of r (un-normalized) w.r.t. various parameters .....
4985 C.. derivatives of t2 and tensor element w.r.t. various parameters ...
4986 C
4987       if (nder.eq.1) then
4988         dr=-2.*ci*xi*omega*(1./vt**2-ci*xi*omega)
4989         dtx_22=(cm1*t1+cm2*t2)*dr
4990       end if
4991       if (nder.eq.2) then
4992         dr=4.*(cnpe*omega)**2*(1.-dcos(xim))
4993         dt2=2.*t2
4994         dtx_22=cm1*(t1*dr+2.*dt2)+cm2*t2*dr
4995       end if
4996       if (nder.eq.3) then
4997         dr=2.*(cnpa*omega*xi)**2
4998         dtx_22=(cm1*t1+cm2*t2)*dr
4999       end if
5000       if (nder.eq.4) then
5001         dr=2.*ci*xi*omega*(1./vt**2-ci*xi*omega)-
5002      %     4.*(cnpe*omega)**2*(1.-dcos(xim))-
5003      %     2.*(cnpa*omega*xi)**2
5004         dt2=-2.*t2
5005         dtx_22=cm1*(t1*dr+2.*dt2)+cm2*t2*dr
5006       end if
5007       if (nder.eq.5) then
5008         dr=-(2./vt)**2*(1./vt**2-ci*xi*omega)
5009         dtx_22=(cm1*t1+cm2*t2)*dr
5010       end if
5011 C
5012 C
5013       return
5014       end
5015 C
5016 C---------------------------------------------------------------------
5017 C
5018       double complex function dtx_23(xi)
5019 C
5020       implicit double precision (a-h, o-z)
5021 C
5022       double complex dcmplx,cdexp,cdsqrt
5023       double complex cnpe
5024       double complex ci,cz,r,sr,t2
5025       double complex ck3,ck4,dr,dt2
5026       double complex cm1,cm2
5027       double complex ce,cfac,cbesk_a
5028 C
5029       double complex ckbs(2)
5030 C
5031       common /para9/ alpha,omega,vt,besk2
5032       common /wave0/ cnpa
5033       common /wave1/ cnpe
5034       common /trubd/ nder
5035 C
5036       external dcbks,umach,erset
5037 C
5038       ci=dcmplx(0.d0,1.d0)
5039       cz=dcmplx(0.d0,0.d0)
5040 C
5041       twopi=6.283185307179586476925287
5042 C
5043 C..... xim is xi modulo two*pi (needed in the sines and cosines) .....
5044 C
5045       xim=dmod(xi,twopi)
5046 C
5047 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
5048 C.... this takes care of the normalization factors in the tensor .....
5049 C... sr is the square-root of r (argument of the Bessel functions) ...
5050 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
5051 C
5052       dom=omega*vt**2
5053       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
5054      %  (cnpa*dom*xi)**2
5055       sr=cdsqrt(r)/vt**2
5056       r=r*besk2
5057 C
5058       if (dreal(sr).lt.50.d0) then
5059 C
5060         numb=2
5061         xnu=3.
5062         call umach(-3,9)
5063         call dcbks(xnu,sr,numb,ckbs)
5064         call erset(0,1,1)
5065 C
5066         ck3=ckbs(1)*dexp(1./vt**2)
5067         ck4=ckbs(2)*dexp(1./vt**2)
5068 C
5069       else
5070 C
5071         ce=1./vt**2-sr
5072         cfac=cdexp(ce)
5073         ck3=cbesk_a(3,sr)*cfac
5074         ck4=cbesk_a(4,sr)*cfac
5075 C
5076       end if
5077 C
5078       cm1=-ck3/(2.*sr*r)
5079       cm2= ck4/(2.*sr**2*r)
5080 C
5081       t2=cnpe*cnpa*omega**2*xi*(1.d0-dcos(xim))
5082 C
5083 C.... derivatives of r (un-normalized) w.r.t. various parameters .....
5084 C.. derivatives of t2 and tensor element w.r.t. various parameters ...
5085 C
5086       if (nder.eq.1) then
5087         dr=-2.*ci*xi*omega*(1./vt**2-ci*xi*omega)
5088         dtx_23=cm2*t2*dr
5089       end if
5090       if (nder.eq.2) then
5091         dr=4.*(cnpe*omega)**2*(1.-dcos(xim))
5092         dt2=t2
5093         dtx_23=2.*cm1*dt2+cm2*t2*dr
5094       end if
5095       if (nder.eq.3) then
5096         dr=2.*(cnpa*omega*xi)**2
5097         dt2=t2
5098         dtx_23=2.*cm1*dt2+cm2*t2*dr
5099       end if
5100       if (nder.eq.4) then
5101         dr=2.*ci*xi*omega*(1./vt**2-ci*xi*omega)-
5102      %     4.*(cnpe*omega)**2*(1.-dcos(xim))-
5103      %     2.*(cnpa*omega*xi)**2
5104         dt2=-2.*t2
5105         dtx_23=2.*cm1*dt2+cm2*t2*dr
5106       end if
5107       if (nder.eq.5) then
5108         dr=-(2./vt)**2*(1./vt**2-ci*xi*omega)
5109         dtx_23=cm2*t2*dr
5110       end if
5111 C
5112 C
5113       return
5114       end
5115 C
5116 C---------------------------------------------------------------------
5117 C
5118       double complex function dtx_33(xi)
5119 C
5120       implicit double precision (a-h, o-z)
5121 C
5122       double complex dcmplx,cdexp,cdsqrt
5123       double complex cnpe
5124       double complex ci,cz,r,sr,t2
5125       double complex ck3,ck4,dr,dt2
5126       double complex cm1,cm2
5127       double complex ce,cfac,cbesk_a
5128 C
5129       double complex ckbs(2)
5130 C
5131       common /para9/ alpha,omega,vt,besk2
5132       common /wave0/ cnpa
5133       common /wave1/ cnpe
5134       common /trubd/ nder
5135 C
5136       external dcbks,umach,erset
5137 C
5138       ci=dcmplx(0.d0,1.d0)
5139       cz=dcmplx(0.d0,0.d0)
5140 C
5141       twopi=6.283185307179586476925287
5142 C
5143 C..... xim is xi modulo two*pi (needed in the sines and cosines) .....
5144 C
5145       xim=dmod(xi,twopi)
5146 C
5147 C... r is the actual r multiplied by (vte/c)**4 and K_2({c/vt}**2) ...
5148 C.... this takes care of the normalization factors in the tensor .....
5149 C... sr is the square-root of r (argument of the Bessel functions) ...
5150 C............ dom = omega*(vte/c)**2 = (w/wce)*(vte/c)**2 ............
5151 C
5152       dom=omega*vt**2
5153       r=(1.-ci*xi*dom)**2+2.*(cnpe*dom)**2*(1.-dcos(xim))+
5154      %  (cnpa*dom*xi)**2
5155       sr=cdsqrt(r)/vt**2
5156       r=r*besk2
5157 C
5158       if (dreal(sr).lt.50.d0) then
5159 C
5160         numb=2
5161         xnu=3.
5162         call umach(-3,9)
5163         call dcbks(xnu,sr,numb,ckbs)
5164         call erset(0,1,1)
5165 C
5166         ck3=ckbs(1)*dexp(1./vt**2)
5167         ck4=ckbs(2)*dexp(1./vt**2)
5168 C
5169       else
5170 C
5171         ce=1./vt**2-sr
5172         cfac=cdexp(ce)
5173         ck3=cbesk_a(3,sr)*cfac
5174         ck4=cbesk_a(4,sr)*cfac
5175 C
5176       end if
5177 C
5178       cm1=-ck3/(2.*sr*r)
5179       cm2= ck4/(2.*sr**2*r)
5180 C
5181       t1=1.d0
5182       t2=(cnpa*omega*xi)**2
5183 C
5184 C.... derivatives of r (un-normalized) w.r.t. various parameters .....
5185 C.. derivatives of t2 and tensor element w.r.t. various parameters ...
5186 C
5187       if (nder.eq.1) then
5188         dr=-2.*ci*xi*omega*(1./vt**2-ci*xi*omega)
5189         dtx_33=(cm1+cm2*t2)*dr
5190       end if
5191       if (nder.eq.2) then
5192         dr=4.*(cnpe*omega)**2*(1.-dcos(xim))
5193         dtx_33=(cm1+cm2*t2)*dr
5194       end if
5195       if (nder.eq.3) then
5196         dr=2.*(cnpa*omega*xi)**2
5197         dt2=2.*t2
5198         dtx_33=cm1*(dr+2.*dt2)+cm2*t2*dr
5199       end if
5200       if (nder.eq.4) then
5201         dr=2.*ci*xi*omega*(1./vt**2-ci*xi*omega)-
5202      %     4.*(cnpe*omega)**2*(1.-dcos(xim))-
5203      %     2.*(cnpa*omega*xi)**2
5204         dt2=-2.*t2
5205         dtx_33=cm1*(dr+2.*dt2)+cm2*t2*dr
5206       end if
5207       if (nder.eq.5) then
5208         dr=-(2./vt)**2*(1./vt**2-ci*xi*omega)
5209         dtx_33=(t1+cm2*t2)*dr
5210       end if
5211 C
5212 C
5213       return
5214       end
5215 C
5216 C=====================================================================
5217 C
5218 C.......... solving the power flow and the power dissipated ..........
5219 C
5220       subroutine pwr_rel(roots,data)
5221 C
5222       implicit double precision (a-h, o-z)
5223 C
5224       double complex dconjg,dcmplx,ci
5225       double complex roots,cnpe
5226       double complex x11h,x12h,x13h,x22h,x23h,x33h
5227       double complex x11a,x12a,x13a,x22a,x23a,x33a
5228       double complex d11h,d12h,d13h,d22h,d23h,d33h
5229       double complex dx11h,dx12h,dx13h,dx22h,dx23h,dx33h
5230       double complex ex,ey,ez
5231       double complex sx,sy,sz,pdis
5232 C
5233       double complex ep(6),dep(6)
5234 C
5235       common /para1/ w0,rp,r0
5236       common /wave0/ cnpa
5237       common /wave1/ cnpe
5238       common /pwrdr/ d11h,d12h,d13h,d22h,d23h,d33h
5239 C
5240       ci=dcmplx(0.d0,1.d0)
5241 C
5242       cnpe=roots
5243 C
5244 C........... getting the tensor elements (this is 4*pi*xi) ...........
5245 C............. ( 4*pi*xi = epsilon_{i,j} - delta_{i,j} ) .............
5246 C.... (kherm = 1 is for getting the hermitian part of the tensor) ....
5247 C
5248       kherm=1
5249       call trub_herm(kherm,ep)
5250 C
5251       x11h=ep(1)
5252       x12h=ep(2)
5253       x13h=ep(3)
5254       x22h=ep(4)
5255       x23h=ep(5)
5256       x33h=ep(6)
5257 C
5258       d11h=1.-cnpa**2+x11h
5259       d12h=x12h
5260       d13h=cnpe*cnpa+x13h
5261       d22h=1.-cnpe**2-cnpa**2+x22h
5262       d23h=x23h
5263       d33h=1.-cnpe**2+x33h
5264 CCCCCCCCCCCCCCCCCCC
5265 CC      write (10,423)
5266 CC  423 format("inside pwr-rel")
5267 CC      write (10,424) d11h,d12h,d13h
5268 CC  424 format("d11h,d12h,d13h = ",6g18.10)
5269 CC      write (10,425) d22h,d23h,d33h
5270 CC  425 format("d22h,d23h,d33h = ",6g18.10)
5271 CCCCCCCCCCCCCCCCCCC
5272 C
5273 C............. evaluating the electric field components ..............
5274 C
5275       call wave_polar_pwr_rel(roots,ex,ey,ez)
5276 C
5277       kherm=1
5278       kdiff=2
5279       call d_drel(kherm,kdiff,ep,dep)
5280 C
5281       dx11h=dep(1)
5282       dx12h=dep(2)
5283       dx13h=dep(3)
5284       dx22h=dep(4)
5285       dx23h=dep(5)
5286       dx33h=dep(6)
5287 C
5288 C.......................... speed of light ...........................
5289       c=2.99792458d10
5290 C
5291 C........... power flow in the x-direction (in cgs units) ............
5292 C
5293       sx=0.5*c*dreal(cnpe*(cdabs(ey)**2+cdabs(ez)**2)-
5294      %               cnpa*dconjg(ex)*ez)-
5295      %   0.25*c/cnpe*(cdabs(ex)**2*dx11h+cdabs(ey)**2*dx22h+
5296      %                cdabs(ez)**2*dx33h+
5297      %                2.*(ci*dimag(dconjg(ex)*ey)*dx12h+
5298      %                       dreal(dconjg(ex)*ez)*dx13h+
5299      %                    ci*dimag(dconjg(ey)*ez)*dx23h))
5300 C
5301 C........... power flow in the y-direction (in cgs units) ............
5302 C
5303       sy=-0.5*c*dreal(cnpa*dconjg(ey)*ez+cnpe*ex*dconjg(ey))
5304 C
5305 C........... power flow in the y-direction (in cgs units) ............
5306 C
5307       ep(1)=x11h
5308       ep(2)=x12h
5309       ep(3)=x13h
5310       ep(4)=x22h
5311       ep(5)=x23h
5312       ep(6)=x33h
5313 C
5314       kherm=1
5315       kdiff=3
5316       call d_drel(kherm,kdiff,ep,dep)
5317 C
5318       dx11h=dep(1)
5319       dx12h=dep(2)
5320       dx13h=dep(3)
5321       dx22h=dep(4)
5322       dx23h=dep(5)
5323       dx33h=dep(6)
5324 C
5325       sz=0.5*c*dreal(cnpa*(cdabs(ex)**2+cdabs(ey)**2)-
5326      %               cnpe*ex*dconjg(ez))-
5327      %   0.25*c/cnpa*(cdabs(ex)**2*dx11h+cdabs(ey)**2*dx22h+
5328      %                cdabs(ez)**2*dx33h+
5329      %                2.*(ci*dimag(dconjg(ex)*ey)*dx12h+
5330      %                       dreal(dconjg(ex)*ez)*dx13h+
5331      %                    ci*dimag(dconjg(ey)*ez)*dx23h))
5332 C
5333 C.............. converting the power flow to mks units ...............
5334 C..... this gives the power flow density (in watts/meter**2) per .....
5335 C.... unit electric field amplitude squared (in {volts/meter}**2) ....
5336 C...... the conversion factor is = 1.e-11/(4*pi*2.99792458**2) .......
5337 C
5338       fact=8.854187817620391d-14
5339 C
5340       sxr=dreal(sx)
5341       sxi=dimag(sx)
5342       if (dabs(sxr).lt.1.d-20) sxr=0.d0
5343       if (dabs(sxi).lt.1.d-20) sxi=0.d0
5344       syr=dreal(sy)
5345       syi=dimag(sy)
5346       if (dabs(syr).lt.1.d-20) syr=0.d0
5347       if (dabs(syi).lt.1.d-20) syi=0.d0
5348       szr=dreal(sz)
5349       szi=dimag(sz)
5350       if (dabs(szr).lt.1.d-20) szr=0.d0
5351       if (dabs(szi).lt.1.d-20) szi=0.d0
5352 C
5353       sx=dcmplx(sxr,sxi)
5354       sy=dcmplx(syr,syi)
5355       sz=dcmplx(szr,szi)
5356       sx=sx*fact
5357       sy=sy*fact
5358       sz=sz*fact
5359 C
5360 C......................... power disspiated ..........................
5361 C
5362       kherm=0
5363       call trub_herm(kherm,ep)
5364 C
5365       x11a=ep(1)
5366       x12a=ep(2)
5367       x13a=ep(3)
5368       x22a=ep(4)
5369       x23a=ep(5)
5370       x33a=ep(6)
5371 C
5372       pdis=-0.5*ci*w0*(x11a*cdabs(ex)**2+x22a*cdabs(ey)**2+
5373      %                 x33a*cdabs(ez)**2+
5374      %                 2.*ci*x12a*dimag(dconjg(ex)*ey)+
5375      %                 2.*x13a*dreal(dconjg(ex)*ez)+
5376      %                 2.*ci*x23a*dimag(dconjg(ey)*ez))
5377 C
5378 C...................... converting to mks units ......................
5379 C.... this gives the power dissipated density (in watts/meter**3) ....
5380 C.... per electric field amplitude squared (in {volts/meter}**2) .....
5381 C........ the conversion factor = 1.e-9/(4.*pi*2.99792458**2) ........
5382 C
5383       fact=8.854187817620391d-12
5384       pdisr=dreal(pdis)
5385       pdisi=dimag(pdis)
5386       if (dabs(pdisr).lt.1.d-20) pdisr=0.d0
5387       if (dabs(pdisi).lt.1.d-20) pdisi=0.d0
5388       pdis=dcmplx(pdisr,pdisi)
5389       pdis=fact*pdis
5390 C
5391       write (13,101) data,cnpe,sx,sy,sz,pdis,ex,ey,ez
5392   101 format(17g18.10)
5393 C
5394 C
5395       return
5396       end
5397 C
5398 C---------------------------------------------------------------------
5399 C
5400 C........... calculating polarizations with respect to E_x ...........
5401 C
5402       subroutine wave_polar_pwr_rel(cnpe,ex,ey,ez)
5403 C
5404       implicit double precision (a-h, o-z)
5405 C
5406       double complex dcmplx,ci,cnpe,cdummy
5407       double complex d11,d12,d13,d22,d23,d33
5408       double complex ex,ey,ez
5409       double complex pol_p,pol_m,pol_z
5410 C
5411       double complex det(3),pdum(3)
5412 C
5413       common /pwrdr/ d11,d12,d13,d22,d23,d33
5414 C
5415       ci=dcmplx(0.d0,1.d0)
5416 C
5417       det(1)=-d12*d23-d13*d22
5418       det(2)= d22*d33-d23*d23
5419       det(3)= d12*d33+d13*d23
5420 C
5421       if (cdabs(det(1)).gt.cdabs(det(2))) then
5422         if (cdabs(det(1)).gt.cdabs(det(3))) then
5423           nn=1
5424         else
5425           nn=3
5426         end if
5427       else
5428         if (cdabs(det(2)).gt.cdabs(det(3))) then
5429           nn=2
5430         else
5431           nn=3
5432         end if
5433       end if
5434 C
5435       if (nn.eq.1) then
5436         ey=(d11*d23-d13*d12)/det(1)
5437         ez=(d12*d12+d11*d22)/det(1)
5438       end if
5439       if (nn.eq.2) then
5440         ey=(d12*d33-d13*d23)/det(2) 
5441         ez=(d12*d23-d13*d22)/det(2)
5442       end if
5443       if (nn.eq.3) then
5444         ey= (d13*d13-d11*d33)/det(3)
5445         ez=-(d12*d13+d11*d23)/det(3)
5446       end if
5447 C
5448 C....... zeroing out the small components of the polarizations .......
5449 C
5450       pdum(1)=(1.d0+ci*ey)/dsqrt(2.d0)
5451       pdum(2)=(1.d0-ci*ey)/dsqrt(2.d0)
5452       pdum(3)=ez
5453 C
5454       do 1 i=1,3
5455         if (dabs(dreal(pdum(i))).lt.1.0d-15) 
5456      %      pdum(i)=dcmplx(0.d0,dimag(pdum(i)))
5457         if (dabs(dimag(pdum(i))).lt.1.0d-15) 
5458      %      pdum(i)=dcmplx(dreal(pdum(i)),0.d0)
5459     1 continue
5460 C
5461       ex=dcmplx(1.0d0,0.d0)
5462       ey=ci*(dsqrt(2.d0)*pdum(2)-1.d0)
5463       ez=pdum(3)
5464 C
5465 C
5466       return
5467       end
5468 C
5469 C---------------------------------------------------------------------
5470 C
5471     subroutine length_jd(string,l_string)
5472 
5473     character string*(*), c
5474     integer l_string
5475 
5476     l_string = 1
5477 1003    continue
5478         c = string(l_string + 1:l_string + 1)
5479         if(c.eq.' ') then
5480             goto 1004
5481         end if
5482     l_string = l_string + 1
5483     goto 1003
5484 
5485 1004    continue
5486 
5487     return
5488     end
5489