-- read the radial wave functions
-- they are stored in the file "NiO_Radial"
-- order of functions
-- r	1S	2S	2P	3S	3P	3D
file = io.open("NiO_Radial", "r")
Rnl = {}
for line in file:lines() do
  RnlLine={}
  for i in string.gmatch(line, "%S+") do
    table.insert(RnlLine,i)
  end
  table.insert(Rnl,RnlLine)
end

-- some constants
a0      =  0.52917721092
Rydberg = 13.60569253
Hartree = 2*Rydberg

-- dd transitions from 3d (index 7 in Rnl) to 
-- 3d (index 7 in Rnl)
-- define a function that calculates
-- <R(r) | j_k(q r) | R(r)>
function RjRdd (q)
  Rj0R = 0
  Rj2R = 0
  Rj4R = 0
  dr = Rnl[3][1]-Rnl[2][1]
  r0 = Rnl[2][1]-2*dr
  for ir = 2, #Rnl, 1 do
    r = r0 + ir * dr
    Rj0R = Rj0R + Rnl[ir][7] * 
      math.SphericalBesselJ(0,q*r) * Rnl[ir][7] * dr
    Rj2R = Rj2R + Rnl[ir][7] * 
      math.SphericalBesselJ(2,q*r) * Rnl[ir][7] * dr
    Rj4R = Rj4R + Rnl[ir][7] * 
      math.SphericalBesselJ(4,q*r) * Rnl[ir][7] * dr
  end
  return Rj0R, Rj2R, Rj4R
end

-- the angular part is given as 
-- C(theta_q, phi_q)^* C(theta_r, phi_r)
-- which is a potential expanded on 
-- spherical harmonics
function ExpandOnClm(k,theta,phi,scale)
  ret={}
  for m=-k, k, 1 do
    table.insert(ret,{k,m,scale * math.SphericalHarmonicC(k,m,theta,phi)})
  end
  return ret
end

-- define nIXS transition operators
function TnIXS_dd(q, theta, phi)
  Rj0R, Rj2R, Rj4R = RjRdd(q)
  k=0
  A0 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj0R)
  T0 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A0)
  k=2
  A2 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj2R)
  T2 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A2)
  k=4
  A4 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj4R)
  T4 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A4)
  T = T0+T2+T4
  T.Chop()
  return T
end

-- q in units per a0 (if you want in units per A 
-- take 5*a0 to have a q of 5 per A)
q=4.5

print("for q=",q," per a0 (",q / a0," per A) The ratio of k=0, k=2 and k=4 transition strength is:",RjRdd(q))

-- define some transition operators
qtheta=0
qphi=0
T = TnIXS_dd(q,qtheta,qphi)

-- calculate the spectra
nIXSSpectra = CreateSpectra(H, T, psi)