(E) Spin.f program for XMCD and SPXAFS

From FEFF

'spin.f’:
      program spin_sum
c    written in grandmotherly fortran-77
      implicit double precision (a-h,o-z)
c     This program read two xmu.dat files for spin -up and -down,
c     calculated with Feff8.20 for the SAME paths list.
c     spin-up file is fort.1, spin-down file is fort.2
c     Both have to be edited: All lines should be deleted except
c       1) line: xsedge+100, used to normalize mu           1.3953E-04
c          leave only on this line:  1.3953E-04
c       2) 6-column data lines
c     The output will be written in fort.3 in 6 columns
c     E+shift1  E(edge)+shift2  xk cmd_total cmd_background  cmd_fs
c     where total = atomic background + fine structure
c     There are 3 possibilities
c     case 1) you want XMCD signal and used SPIN  1 or -1
c     case 2) you want XMCD signal and used SPIN 2 or -2, in order
c       to use non-relativistic formula for XMCD
c       factor li/2j+1 which was not convenient to do in a program
c     case 3) you want SPXAFS and used  SPIN \pm 2
c     ENTER your case here (icase is positive integer only)
icase = 2
c     if icase=2 ENTER factor=(-1)**(L+1/2-J) * L/(2*J+1)
c     where L,J are for your edge (ex. for L3 L=1 J=3/2, for L2 L=1 J=1/2)
c     for L3
c 
c
c
c
c factor = 0.25
c for L2
c factor = -0.5
c ENTER the energy shift you want for columns 1 and 2 in xmu.dat
shift1 = 0
shift2 = 0
c everything below is automated further
read (1,*,end=10) ap
read (2,*,end=10) am
xnorm = 0.5 *(ap+am)
c read the data
3    read(1,*,end=10)   x1, x2, ek, y1, y2, y3
   read(2,*,end=10)  x1, x2, ek, z1, z2, z3
   if (icase.eq.1) then
c      no xafs in this case:xfs - atomic part of XMCD
      t1 = (y1*ap + z1*am)/xnorm
      t2 = (y2*ap + z2*am)/xnorm
      t3 = (y3*ap + z3*am) /xnorm
   elseif (icase.eq.2) then
      t1 = (y1*ap - z1*am)*factor /xnorm
      t2 = (y2*ap - z2*am)*factor /xnorm
      t3 = (y3*ap - z3*am)*factor /xnorm
   elseif (icase.eq.3) then
 c     factor=0.5 always for SPXAFS
      t1 = (y1*ap - z1*am)/2.0/xnorm
      t2 = (y2*ap - z2*am)/2.0/xnorm
      t3 = (y3*ap - z3*am)/2.0/xnorm
c      you may want average total XAS as output in last column
c      t3 = (y1*ap + z1*am)/2.0/xnorm
   endif
   x1 =x1 + shift1
   x2 =x2 + shift2
   write(3,5)    x1, x2, ek, t1, t2, t3
5   format (6e13.5)
c
    goto 3
10  continue
stop 
end


developer's resources