LCOV - code coverage report
Current view: top level - core - kernel2.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 130 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             :       MODULE m_kernel2
       2             : C.............................................................kernel2
       3             : c solution of 4 coupled equations (logariphmic mesh) for (l,\mu)
       4             : c 1->k1, 2->k2
       5             : c d/dxP1=[-k1]P1+exp(x)[e-v+2mc^2]/c^2Q1+exp(x)B/c^2\sigma_z[-1,-1]Q1
       6             : c d/dxQ1=[k1]Q1-exp(x)[e-v]P1+exp(x)B[\sigma_z[1,1]P1+\sigma_z[1,2]P2]
       7             : c d/dxP2=[-k2]P2+exp(x)[e-v+2mc^2]/c^2Q2+exp(x)B/c^2\sigma_z[-1,-1]Q2
       8             : c d/dxQ2=[k2]Q2-exp(x)[e-v]P2+exp(x)B[\sigma_z[2,1]P1+\sigma_z[2,2]P2]
       9             : c notations:
      10             : c Pi=exp(x)*gi
      11             : c Qi=exp(x)*cfi
      12             : c \gamma=sqrt(k^2-(Z/c)^2)
      13             : c v=0.5[V(1/2)+V(-1/2)]
      14             : c b=0.5[V(1/2)-V(-1/2)]
      15             : c Rydberg units: in charge
      16             : c Hartree units: com.
      17             : c NSOL= 2 -  4 equations
      18             : c ------------                                     a. shick KFA 1996
      19             :       CONTAINS
      20           0 :       SUBROUTINE kernel2(mrad,nsol,xmj,k1,k2,xx1,xx2,e,v,b,ri,dx,
      21           0 :      +                   nmatch,nstart,dp,dq,wp,wq)
      22             : 
      23             : 
      24             :       USE m_constants, ONLY : c_light
      25             :       USE m_diff
      26             :       IMPLICIT NONE
      27             : C     ..
      28             : C     .. Scalar Arguments ..
      29             :       INTEGER, INTENT (IN) :: mrad
      30             :       REAL dx,e,xmj
      31             :       INTEGER k1,k2,nmatch,nsol,nstart
      32             : C     ..
      33             : C     .. Array Arguments ..
      34             :       REAL b(mrad),dp(2,2,mrad),dq(2,2,mrad),ri(mrad),v(mrad),
      35             :      +     wp(2,2,mrad),wq(2,2,mrad),xx1(4),xx2(4)
      36             : C     ..
      37             : C     .. Local Scalars ..
      38             :       REAL abz,bh,bmn,bmnp1,bn,bnp1,c1,c2,c3,cc,csq,dp11,dp12,
      39             :      +     dp13,dp14,dp21,dp22,dp23,dp24,dq11,dq12,dq13,dq14,dq21,dq22,
      40             :      +     dq23,dq24,dxd8,expdxh,gk,gk1,gk1k,gkk1,gmk,gmk1,p1c,p2c,
      41             :      +     q1c,q2c,r,vh,vmn,vmnp1,vn,vnp1,xk,xk1,xm
      42             :       INTEGER i,iny,ir,j,jri,n,nrk
      43             : C     ..
      44             : C     .. Local Arrays ..
      45           0 :       REAL bm(mrad),p1(mrad),p1s(mrad),p2(mrad),p2s(mrad),q1(mrad),
      46           0 :      +     q1s(mrad),q2(mrad),q2s(mrad),vm(mrad)
      47             : C     ..
      48             : C     .. Intrinsic Functions ..
      49             :       INTRINSIC exp,real,sqrt
      50             : C     ..
      51             : C     .. Data statements ..
      52             :       DATA abz/0.5/
      53             : C     ..
      54           0 :       cc = c_light(2.0)
      55             : c
      56           0 :       DO ir = 1,mrad
      57           0 :          p1(ir) = 0.0
      58           0 :          q1(ir) = 0.0
      59           0 :          p1s(ir) = 0.0
      60           0 :          q1s(ir) = 0.0
      61             : c
      62           0 :          p2(ir) = 0.0
      63           0 :          q2(ir) = 0.0
      64           0 :          p2s(ir) = 0.0
      65           0 :          q2s(ir) = 0.0
      66             : c
      67           0 :          DO i = 1,2
      68           0 :             DO j = 1,2
      69             : c
      70           0 :                wp(i,j,ir) = 0.0
      71           0 :                wq(i,j,ir) = 0.0
      72             : c
      73           0 :                dp(i,j,ir) = 0.0
      74           0 :                dq(i,j,ir) = 0.0
      75             :             END DO
      76             :          END DO
      77             :       END DO
      78             : C
      79           0 :       csq = cc*cc
      80           0 :       expdxh = exp(dx/2.0)
      81           0 :       dxd8 = dx/8.0
      82           0 :       jri = nmatch
      83           0 :       nrk = jri
      84             : c-derivatives of potential
      85             :       CALL diff(
      86             :      >          mrad,v,dx,jri,
      87           0 :      <          vm)
      88             :       CALL diff(
      89             :      >          mrad,b,dx,jri,
      90           0 :      <          bm)
      91             : c- begin to solve coupled Dirac equation in "magnetic" field
      92           0 :       xk = real(k1)
      93           0 :       xk1 = real(k2)
      94           0 :       xm = xmj
      95           0 :       gk = -xm/ (xk+abz)
      96           0 :       gk1 = -xm/ (xk1+abz)
      97           0 :       gkk1 = -sqrt(1.0-gk*gk)
      98           0 :       gk1k = gkk1
      99           0 :       gmk = xm/ (xk-abz)
     100           0 :       gmk1 = xm/ (xk1-abz)
     101             : c
     102             : c     NB: V=R**2*POTENTIAL
     103             : c     NB: B=R**2*FIELD
     104             : c
     105           0 :       DO 30 iny = 1,nsol
     106             : C                                                                       WAB02140
     107           0 :          r = ri(nstart)
     108           0 :          vn = v(nstart)
     109           0 :          vmn = vm(nstart)
     110           0 :          bn = b(nstart)
     111           0 :          bmn = bm(nstart)
     112           0 :          IF (iny.EQ.1) THEN
     113           0 :             p1(nstart) = xx1(1)
     114           0 :             q1(nstart) = xx1(2)
     115           0 :             p2(nstart) = xx1(3)
     116           0 :             q2(nstart) = xx1(4)
     117             :          ELSE
     118           0 :             p1(nstart) = xx2(1)
     119           0 :             q1(nstart) = xx2(2)
     120           0 :             p2(nstart) = xx2(3)
     121           0 :             q2(nstart) = xx2(4)
     122             :          END IF
     123             : c---->f.po solution of Dirac eq.
     124           0 :          c1 = vn/r**2 - e
     125           0 :          c2 = 1.0 - c1/csq
     126             : c      C2=2.0-C1/CSQ !H
     127           0 :          c3 = bn/r/r
     128           0 :          p1s(nstart) = -xk*p1(nstart) + r* (c2+c3/csq*gmk)*q1(nstart)
     129             :          q1s(nstart) = xk*q1(nstart) + r*
     130           0 :      +                 ((c1+c3*gk)*p1(nstart)+c3*gkk1*p2(nstart))
     131             : C
     132           0 :          p2s(nstart) = -xk1*p2(nstart) + r* (c2+c3/csq*gmk1)*q2(nstart)
     133             :          q2s(nstart) = xk1*q2(nstart) +
     134           0 :      +                 r* ((c1+c3*gk1)*p2(nstart)+c3*gk1k*p1(nstart))
     135             : C   ******************************************************************
     136             : C   *                                                                *
     137             : C   *     THE RUNGE-KUTTA METHOD IS USED IN THE <NKR> FIRST STEPS    *
     138             : C   *                                                                *
     139             : C   ******************************************************************
     140             : C
     141           0 :          n = nstart
     142             : c-> y1
     143           0 :    10    p1c = p1(n)
     144           0 :          q1c = q1(n)
     145           0 :          p2c = p2(n)
     146           0 :          q2c = q2(n)
     147             : c-> k1=f(x1,y1)
     148           0 :          dp11 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
     149           0 :          dq11 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
     150           0 :          dp21 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
     151           0 :          dq21 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
     152             : c-> y11=y1+0.5*k1
     153           0 :          p1c = p1c + 0.5*dp11
     154           0 :          q1c = q1c + 0.5*dq11
     155           0 :          p2c = p2c + 0.5*dp21
     156           0 :          q2c = q2c + 0.5*dq21
     157             : c-> x1'=x1+0.5h
     158           0 :          r = r*expdxh
     159             : c-> interpolation of V and B for x1'
     160           0 :          vnp1 = v(n+1)
     161           0 :          vmnp1 = vm(n+1)
     162           0 :          bnp1 = b(n+1)
     163           0 :          bmnp1 = bm(n+1)
     164           0 :          vh = (vn+vnp1)*0.5 + (vmn-vmnp1)*dxd8
     165           0 :          bh = (bn+bnp1)*0.5 + (bmn-bmnp1)*dxd8
     166             : c-> k2=f(x1',y11)
     167           0 :          c1 = vh/r/r - e
     168           0 :          c2 = 1.0 - c1/csq
     169             : c      C2=2.0-C1/CSQ !H
     170           0 :          c3 = bh/r/r
     171             : C                                                                       WAB02850
     172           0 :          dp12 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
     173           0 :          dq12 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
     174           0 :          dp22 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
     175           0 :          dq22 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
     176             : c-> y12=y1+0.5*k2 (y12=y11-0.5*k1+0.5*k2)
     177           0 :          p1c = p1c + 0.5* (dp12-dp11)
     178           0 :          q1c = q1c + 0.5* (dq12-dq11)
     179           0 :          p2c = p2c + 0.5* (dp22-dp21)
     180           0 :          q2c = q2c + 0.5* (dq22-dq21)
     181             : c-> k3=f(x1',y12)
     182           0 :          dp13 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
     183           0 :          dq13 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
     184           0 :          dp23 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
     185           0 :          dq23 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
     186             : c-> y13=y1+k3 (y13=y12-0.5*k2+k3)
     187           0 :          p1c = p1c + dp13 - 0.5*dp12
     188           0 :          q1c = q1c - 0.5*dq12 + dq13
     189           0 :          p2c = p2c + dp23 - 0.5*dp22
     190           0 :          q2c = q2c - 0.5*dq22 + dq23
     191             : c-> x2=x1+h
     192           0 :          n = n + 1
     193           0 :          r = ri(n)
     194             : c-> k4=f(x2,y13)
     195           0 :          c1 = vnp1/r/r - e
     196           0 :          c2 = 1.0 - c1/csq
     197             : c      C2=2.0-C1/CSQ !H
     198           0 :          c3 = bnp1/r/r
     199             : c
     200           0 :          dp14 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
     201           0 :          dq14 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c+c3*gkk1*p2c))
     202           0 :          dp24 = dx* (-xk1*p2c+r* (c2+c3/csq*gmk1)*q2c)
     203           0 :          dq24 = dx* (xk1*q2c+r* ((c1+c3*gk1)*p2c+c3*gk1k*p1c))
     204             : c-> y2=y1+1/6*(k1+2*k2+2*k3+k4)
     205           0 :          p1(n) = p1(n-1) + (dp11+2.0* (dp12+dp13)+dp14)/6.0
     206           0 :          q1(n) = q1(n-1) + (dq11+2.0* (dq12+dq13)+dq14)/6.0
     207           0 :          p2(n) = p2(n-1) + (dp21+2.0* (dp22+dp23)+dp24)/6.0
     208           0 :          q2(n) = q2(n-1) + (dq21+2.0* (dq22+dq23)+dq24)/6.0
     209             : c-> f(x2,y2)
     210           0 :          p1s(n) = -xk*p1(n) + r* (c2+c3/csq*gmk)*q1(n)
     211           0 :          q1s(n) = xk*q1(n) + r* ((c1+c3*gk)*p1(n)+c3*gkk1*p2(n))
     212           0 :          p2s(n) = -xk1*p2(n) + r* (c2+c3/csq*gmk1)*q2(n)
     213           0 :          q2s(n) = xk1*q2(n) + r* ((c1+c3*gk1)*p2(n)+c3*gk1k*p1(n))
     214             : c-> redefinition of V ind B (2->1)
     215           0 :          vn = vnp1
     216           0 :          bn = bnp1
     217           0 :          vmn = vmnp1
     218           0 :          bmn = bmnp1
     219             : C                                                                       WAB03340
     220           0 :          IF (n-nrk.LT.0) GOTO 10
     221             : c---->Milne's method off now!
     222             : c
     223             : c storage of wave functions and their derivatives.
     224           0 :          DO ir = 1,jri
     225           0 :             wp(1,iny,ir) = p1(ir)
     226           0 :             wp(2,iny,ir) = p2(ir)
     227           0 :             wq(1,iny,ir) = q1(ir)
     228           0 :             wq(2,iny,ir) = q2(ir)
     229             : c derivatives
     230           0 :             dp(1,iny,ir) = p1s(ir)
     231           0 :             dp(2,iny,ir) = p2s(ir)
     232           0 :             dq(1,iny,ir) = q1s(ir)
     233           0 :             dq(2,iny,ir) = q2s(ir)
     234             :          END DO
     235             : C***********************************************************            WAB05210
     236           0 :    30 CONTINUE
     237             : c      write(6,*) wp
     238             :       
     239           0 :       END SUBROUTINE kernel2
     240             :       END MODULE m_kernel2

Generated by: LCOV version 1.13