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

          Line data    Source code
       1             :       MODULE m_kernel1
       2             : C..............................................................kernel1
       3             : c solution of 2 coupled equations (logariphmic mesh) for (l,\mu:|mu|>=l)
       4             : c 1->k1
       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]
       7             : c notations:
       8             : c Pi=exp(x)*gi
       9             : c Qi=exp(x)*cfi
      10             : c \gamma=sqrt(k^2-(Z/c)^2)
      11             : c v=0.5[V(1/2)+V(-1/2)]
      12             : c b=0.5[V(1/2)-V(-1/2)]
      13             : c Rydberg units: in charge
      14             : c hartree units: com.
      15             : c NSOL= 1 -  2 equations
      16             : c ------------                                     a. shick KFA 1996
      17             :       CONTAINS
      18           0 :       SUBROUTINE kernel1(mrad,xmj,kap1,xx1,e,v,b,ri,dx,nmatch,nstart,
      19           0 :      +                   dp,dq,wp,wq)
      20             : 
      21             : c
      22             :       USE m_constants, ONLY : c_light
      23             :       USE m_diff
      24             :       IMPLICIT NONE
      25             : C     ..
      26             : C     .. Scalar Arguments ..
      27             :       INTEGER, INTENT (IN) :: mrad
      28             :       REAL dx,e,xmj
      29             :       INTEGER kap1,nmatch,nstart
      30             : C     ..
      31             : C     .. Array Arguments ..
      32             :       REAL b(mrad),dp(2,2,mrad),dq(2,2,mrad),ri(mrad),v(mrad),
      33             :      +     wp(2,2,mrad),wq(2,2,mrad),xx1(4)
      34             : C     ..
      35             : C     .. Local Scalars ..
      36             :       REAL abz,bh,bmn,bmnp1,bn,bnp1,c1,c2,c3,cc,csq,dp11,dp12,
      37             :      +     dp13,dp14,dq11,dq12,dq13,dq14,dxd8,expdxh,gk,gmk,p1c,
      38             :      +     q1c,r,vh,vmn,vmnp1,vn,vnp1,xk,xm
      39             :       INTEGER i,ir,j,jri,n,nrk
      40             : C     ..
      41             : C     .. Local Arrays ..
      42           0 :       REAL bm(mrad),p1(mrad),p1s(mrad),q1(mrad),q1s(mrad),vm(mrad)
      43             : C     ..
      44             : C     .. Intrinsic Functions ..
      45             :       INTRINSIC exp,real
      46             : C     ..
      47             : C     .. Data statements ..
      48             :       DATA abz/0.5/
      49             : C     ..
      50           0 :       cc = c_light(2.0)
      51             : c
      52           0 :       DO ir = 1,mrad
      53           0 :          p1(ir) = 0.0
      54           0 :          q1(ir) = 0.0
      55           0 :          p1s(ir) = 0.0
      56           0 :          q1s(ir) = 0.0
      57             : c
      58           0 :          DO i = 1,2
      59           0 :             DO j = 1,2
      60             : c
      61           0 :                wp(i,j,ir) = 0.0
      62           0 :                wq(i,j,ir) = 0.0
      63             : c
      64           0 :                dp(i,j,ir) = 0.0
      65           0 :                dq(i,j,ir) = 0.0
      66             :             END DO
      67             :          END DO
      68             :       END DO
      69             : C
      70           0 :       csq = cc*cc
      71           0 :       expdxh = exp(dx/2.0)
      72           0 :       dxd8 = dx/8.0
      73           0 :       jri = nmatch
      74           0 :       nrk = jri
      75             :       CALL diff(
      76             :      >          mrad,v,dx,jri,
      77           0 :      <          vm)
      78             :       CALL diff(
      79             :      >          mrad,b,dx,jri,
      80           0 :      <          bm)
      81           0 :       xk = real(kap1)
      82           0 :       xm = xmj
      83           0 :       gk = -xm/ (xk+abz)
      84           0 :       gmk = xm/ (xk-abz)
      85             : c
      86           0 :       r = ri(nstart)
      87           0 :       vn = v(nstart)
      88           0 :       vmn = vm(nstart)
      89           0 :       bn = b(nstart)
      90           0 :       bmn = bm(nstart)
      91           0 :       p1(nstart) = xx1(1)
      92           0 :       q1(nstart) = xx1(2)
      93             : c---->f.po solution of Dirac eq.
      94           0 :       c1 = vn/r**2 - e
      95           0 :       c2 = 1.0 - c1/csq
      96             : c      C2=2.0-C1/CSQ !H
      97           0 :       c3 = bn/r/r
      98           0 :       p1s(nstart) = -xk*p1(nstart) + r* (c2+c3/csq*gmk)*q1(nstart)
      99           0 :       q1s(nstart) = xk*q1(nstart) + r* ((c1+c3*gk)*p1(nstart))
     100             : C   *     THE RUNGE-KUTTA METHOD IS USED IN THE <NKR> FIRST STEPS    *  WAB02520
     101           0 :       n = nstart
     102             : c-> y1
     103           0 :    10 p1c = p1(n)
     104           0 :       q1c = q1(n)
     105             : c-> k1=f(x1,y1)
     106           0 :       dp11 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
     107           0 :       dq11 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
     108             : c-> y11=y1+0.5*k1
     109           0 :       p1c = p1c + 0.5*dp11
     110           0 :       q1c = q1c + 0.5*dq11
     111             : c-> x1'=x1+0.5h
     112           0 :       r = r*expdxh
     113             : c-> interpolation of V and B for x1'
     114           0 :       vnp1 = v(n+1)
     115           0 :       vmnp1 = vm(n+1)
     116           0 :       bnp1 = b(n+1)
     117           0 :       bmnp1 = bm(n+1)
     118           0 :       vh = (vn+vnp1)*0.5 + (vmn-vmnp1)*dxd8
     119           0 :       bh = (bn+bnp1)*0.5 + (bmn-bmnp1)*dxd8
     120             : c-> k2=f(x1',y11)
     121           0 :       c1 = vh/r/r - e
     122           0 :       c2 = 1.0 - c1/csq
     123             : c      C2=2.0-C1/CSQ !H
     124           0 :       c3 = bh/r/r
     125             : C                                                                       WAB02850
     126           0 :       dp12 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
     127           0 :       dq12 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
     128             : c-> y12=y1+0.5*k2 (y12=y11-0.5*k1+0.5*k2)
     129           0 :       p1c = p1c + 0.5* (dp12-dp11)
     130           0 :       q1c = q1c + 0.5* (dq12-dq11)
     131             : c-> k3=f(x1',y12)
     132           0 :       dp13 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
     133           0 :       dq13 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
     134             : c-> y13=y1+k3 (y13=y12-0.5*k2+k3)
     135           0 :       p1c = p1c + dp13 - 0.5*dp12
     136           0 :       q1c = q1c - 0.5*dq12 + dq13
     137             : c-> x2=x1+h
     138           0 :       n = n + 1
     139           0 :       r = ri(n)
     140             : c-> k4=f(x2,y13)
     141           0 :       c1 = vnp1/r/r - e
     142           0 :       c2 = 1.0 - c1/csq
     143             : c      C2=2.0-C1/CSQ !H
     144           0 :       c3 = bnp1/r/r
     145             : c
     146           0 :       dp14 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
     147           0 :       dq14 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
     148             : c-> y2=y1+1/6*(k1+2*k2+2*k3+k4)
     149           0 :       p1(n) = p1(n-1) + (dp11+2.0* (dp12+dp13)+dp14)/6.0
     150           0 :       q1(n) = q1(n-1) + (dq11+2.0* (dq12+dq13)+dq14)/6.0
     151             : c-> f(x2,y2)
     152           0 :       p1s(n) = -xk*p1(n) + r* (c2+c3/csq*gmk)*q1(n)
     153           0 :       q1s(n) = xk*q1(n) + r* ((c1+c3*gk)*p1(n))
     154             : c-> redefinition of V ind B (2->1)
     155           0 :       vn = vnp1
     156           0 :       bn = bnp1
     157           0 :       vmn = vmnp1
     158           0 :       bmn = bmnp1
     159             : C                                                                       WAB03340
     160           0 :       IF (n-nrk.LT.0) GOTO 10
     161             : c---->Milne's method off now!
     162           0 :       DO ir = 1,mrad
     163           0 :          wp(1,1,ir) = p1(ir)
     164           0 :          wq(1,1,ir) = q1(ir)
     165             : c
     166           0 :          dp(1,1,ir) = p1s(ir)
     167           0 :          dq(1,1,ir) = q1s(ir)
     168             :       END DO
     169             : 
     170           0 :       END SUBROUTINE kernel1
     171             :       END MODULE m_kernel1

Generated by: LCOV version 1.13