kernel1.f 4.64 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171
      MODULE m_kernel1
C..............................................................kernel1
c solution of 2 coupled equations (logariphmic mesh) for (l,\mu:|mu|>=l)
c 1->k1
c d/dxP1=[-k1]P1+exp(x)[e-v+2mc^2]/c^2Q1+exp(x)B/c^2\sigma_z[-1,-1]Q1
c d/dxQ1=[k1]Q1-exp(x)[e-v]P1+exp(x)B[\sigma_z[1,1]P1]
c notations:
c Pi=exp(x)*gi
c Qi=exp(x)*cfi
c \gamma=sqrt(k^2-(Z/c)^2)
c v=0.5[V(1/2)+V(-1/2)]
c b=0.5[V(1/2)-V(-1/2)]
c Rydberg units: in charge
c hartree units: com.
c NSOL= 1 -  2 equations
c ------------                                     a. shick KFA 1996
      CONTAINS
      SUBROUTINE kernel1(mrad,xmj,kap1,xx1,e,v,b,ri,dx,nmatch,nstart,
     +                   dp,dq,wp,wq)

c
      USE m_constants, ONLY : c_light
      USE m_diff
      IMPLICIT NONE
C     ..
C     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: mrad
      REAL dx,e,xmj
      INTEGER kap1,nmatch,nstart
C     ..
C     .. Array Arguments ..
      REAL b(mrad),dp(2,2,mrad),dq(2,2,mrad),ri(mrad),v(mrad),
     +     wp(2,2,mrad),wq(2,2,mrad),xx1(4)
C     ..
C     .. Local Scalars ..
      REAL abz,bh,bmn,bmnp1,bn,bnp1,c1,c2,c3,cc,csq,dp11,dp12,
     +     dp13,dp14,dq11,dq12,dq13,dq14,dxd8,expdxh,gk,gmk,p1c,
     +     q1c,r,vh,vmn,vmnp1,vn,vnp1,xk,xm
      INTEGER i,ir,j,jri,n,nrk
C     ..
C     .. Local Arrays ..
      REAL bm(mrad),p1(mrad),p1s(mrad),q1(mrad),q1s(mrad),vm(mrad)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC exp,real
C     ..
C     .. Data statements ..
      DATA abz/0.5/
C     ..
      cc = c_light(2.0)
c
      DO ir = 1,mrad
         p1(ir) = 0.0
         q1(ir) = 0.0
         p1s(ir) = 0.0
         q1s(ir) = 0.0
c
         DO i = 1,2
            DO j = 1,2
c
               wp(i,j,ir) = 0.0
               wq(i,j,ir) = 0.0
c
               dp(i,j,ir) = 0.0
               dq(i,j,ir) = 0.0
            END DO
         END DO
      END DO
C
      csq = cc*cc
      expdxh = exp(dx/2.0)
      dxd8 = dx/8.0
      jri = nmatch
      nrk = jri
      CALL diff(
     >          mrad,v,dx,jri,
     <          vm)
      CALL diff(
     >          mrad,b,dx,jri,
     <          bm)
      xk = real(kap1)
      xm = xmj
      gk = -xm/ (xk+abz)
      gmk = xm/ (xk-abz)
c
      r = ri(nstart)
      vn = v(nstart)
      vmn = vm(nstart)
      bn = b(nstart)
      bmn = bm(nstart)
      p1(nstart) = xx1(1)
      q1(nstart) = xx1(2)
c---->f.po solution of Dirac eq.
      c1 = vn/r**2 - e
      c2 = 1.0 - c1/csq
c      C2=2.0-C1/CSQ !H
      c3 = bn/r/r
      p1s(nstart) = -xk*p1(nstart) + r* (c2+c3/csq*gmk)*q1(nstart)
      q1s(nstart) = xk*q1(nstart) + r* ((c1+c3*gk)*p1(nstart))
C   *     THE RUNGE-KUTTA METHOD IS USED IN THE <NKR> FIRST STEPS    *  WAB02520
      n = nstart
c-> y1
   10 p1c = p1(n)
      q1c = q1(n)
c-> k1=f(x1,y1)
      dp11 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
      dq11 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
c-> y11=y1+0.5*k1
      p1c = p1c + 0.5*dp11
      q1c = q1c + 0.5*dq11
c-> x1'=x1+0.5h
      r = r*expdxh
c-> interpolation of V and B for x1'
      vnp1 = v(n+1)
      vmnp1 = vm(n+1)
      bnp1 = b(n+1)
      bmnp1 = bm(n+1)
      vh = (vn+vnp1)*0.5 + (vmn-vmnp1)*dxd8
      bh = (bn+bnp1)*0.5 + (bmn-bmnp1)*dxd8
c-> k2=f(x1',y11)
      c1 = vh/r/r - e
      c2 = 1.0 - c1/csq
c      C2=2.0-C1/CSQ !H
      c3 = bh/r/r
C                                                                       WAB02850
      dp12 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
      dq12 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
c-> y12=y1+0.5*k2 (y12=y11-0.5*k1+0.5*k2)
      p1c = p1c + 0.5* (dp12-dp11)
      q1c = q1c + 0.5* (dq12-dq11)
c-> k3=f(x1',y12)
      dp13 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
      dq13 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
c-> y13=y1+k3 (y13=y12-0.5*k2+k3)
      p1c = p1c + dp13 - 0.5*dp12
      q1c = q1c - 0.5*dq12 + dq13
c-> x2=x1+h
      n = n + 1
      r = ri(n)
c-> k4=f(x2,y13)
      c1 = vnp1/r/r - e
      c2 = 1.0 - c1/csq
c      C2=2.0-C1/CSQ !H
      c3 = bnp1/r/r
c
      dp14 = dx* (-xk*p1c+r* (c2+c3/csq*gmk)*q1c)
      dq14 = dx* (xk*q1c+r* ((c1+c3*gk)*p1c))
c-> y2=y1+1/6*(k1+2*k2+2*k3+k4)
      p1(n) = p1(n-1) + (dp11+2.0* (dp12+dp13)+dp14)/6.0
      q1(n) = q1(n-1) + (dq11+2.0* (dq12+dq13)+dq14)/6.0
c-> f(x2,y2)
      p1s(n) = -xk*p1(n) + r* (c2+c3/csq*gmk)*q1(n)
      q1s(n) = xk*q1(n) + r* ((c1+c3*gk)*p1(n))
c-> redefinition of V ind B (2->1)
      vn = vnp1
      bn = bnp1
      vmn = vmnp1
      bmn = bmnp1
C                                                                       WAB03340
      IF (n-nrk.LT.0) GOTO 10
c---->Milne's method off now!
      DO ir = 1,mrad
         wp(1,1,ir) = p1(ir)
         wq(1,1,ir) = q1(ir)
c
         dp(1,1,ir) = p1s(ir)
         dq(1,1,ir) = q1s(ir)
      END DO

      END SUBROUTINE kernel1
      END MODULE m_kernel1