-- The calculations are done using a set of basis spin-orbitals
-- These are very often given by some radial function times a spherical harmonic
-- (or given by sum of these)
-- Here I show for the very simple example of an H atom how one can calculate the 
-- appropriate integrals

-- a gnuplot script to make the plots
gnuplotInput = [[
set autoscale 
set xtic auto 
set ytic auto  
set style line  1 lt 1 lw 1 lc rgb "#FF0000"
set style line  2 lt 1 lw 1 lc rgb "#0000FF"
set style line  3 lt 1 lw 1 lc rgb "#00C000"
set style line  4 lt 1 lw 1 lc rgb "#FF00FF"
set style line  5 lt 1 lw 1 lc rgb "#00FFFF"
set style line  6 lt 1 lw 1 lc rgb "#000000"

set ylabel "Rnl" font "Times,12"
set xlabel "r(A)" font "Times,12"

set out 'HWaveFunctions.ps'
set size 1.0, 1.0
set terminal postscript portrait enhanced color  "Times" 8

plot "Radial" using 1:2  title "R_{1s}(r)" with lines ls 1
plot "Radial" using 1:3  title "R_{2s}(r)" with lines ls 1,\
     "Radial" using 1:4  title "R_{2p}(r)" with lines ls 2
plot "Radial" using 1:5  title "R_{3s}(r)" with lines ls 1,\
     "Radial" using 1:6  title "R_{3p}(r)" with lines ls 2,\
     "Radial" using 1:7  title "R_{3d}(r)" with lines ls 3
plot "Radial" using 1:8  title "R_{4s}(r)" with lines ls 1,\
     "Radial" using 1:9  title "R_{4p}(r)" with lines ls 2,\
     "Radial" using 1:10 title "R_{4d}(r)" with lines ls 3,\
     "Radial" using 1:11 title "R_{4f}(r)" with lines ls 4
plot "Radial" using 1:12 title "R_{5s}(r)" with lines ls 1,\
     "Radial" using 1:13 title "R_{5p}(r)" with lines ls 2,\
     "Radial" using 1:14 title "R_{5d}(r)" with lines ls 3,\
     "Radial" using 1:15 title "R_{5f}(r)" with lines ls 4,\
     "Radial" using 1:16 title "R_{5g}(r)" with lines ls 5

set ylabel "<R|j|R>" font "Times,12"
set xlabel "q(A^{-1})" font "Times,12"

i=1
do for [n1=1:3] {
  do for [l1=0:n1-1] {
    do for [n2=1:3] {
      do for [l2=0:n2-1] {
        do for [k=abs(l1-l2):l1+l2:2] {
          i=i+1
          plot "Bessel" using 1:(column(i)) title "<R_{".n1." ".l1."}(r)|j[qr]|R_{".n2." ".l2."}(r)" with lines ls 6
        }
      }
    }
  }
}

]]

-- bohr radius in units of Angstrom
a0 = 0.52917721092;

-- Hydrogen wave functions with effective Z
function Rnl (n, l, r, Zeff)
  return 2^(l+1) * e^(-r*Zeff/(a0 * n)) * n^(-2-l) * r^l * (Zeff/a0)^(l+3/2) * math.sqrt(math.fact(n-l-1)/math.fact(n+l)) * math.LaguerreL(n-l-1, 2*l+1, 2*r*Zeff/(a0*n))
end

-- write the radial functions to a file
file = io.open("Radial", "w")
dr=0.01
Nr=5000
for ir = 0, Nr, 1 do
  r = ir*dr
  file:write(string.format("%15.8E ",r))
  for n=1, 5, 1 do
    for l=0, n-1, 1 do
      file:write(string.format("%15.8E ",Rnl(n,l,r,1)))
    end
  end
  file:write("\n")
end
file:close()

-- normalization
io.write("\n<Rnl|Rnl>\n");
dr=0.01
Nr=5000
for n=1, 5, 1 do
  for l=0, n-1, 1 do
    sum=0
    for ir = 0, Nr, 1 do
      r = ir*dr
      sum = sum + dr * r^2 * Rnl(n,l,r,1)^2
    end
    io.write(string.format("%15.8E ",sum))
  end
  io.write("\n")
end

-- average distance
io.write("\n<Rnl|r|Rnl> [units of A]\n");
dr=0.01
Nr=5000
for n=1, 5, 1 do
  for l=0, n-1, 1 do
    sum=0
    for ir = 0, Nr, 1 do
      r = ir*dr
      sum = sum + dr * r^2 * r * Rnl(n,l,r,1)^2
    end
    io.write(string.format("%15.8E ",sum))
  end
io.write("\n")
end

-- average distance
io.write("\n<Rnl|r|Rnl> [units of a0]\n");
dr=0.01
Nr=5000
for n=1, 5, 1 do
  for l=0, n-1, 1 do
  sum=0
    for ir = 0, Nr, 1 do
    r = ir*dr
    sum = sum + dr * r^2 * r * Rnl(n,l,r,1)^2 / a0
  end
  io.write(string.format("%15.8E ",sum))
end
io.write("\n")
end

-- dipole moments
io.write("\n<Rn1l1|r|Rn2l2> [units of a0]\n");
dr=0.01
Nr=5000
io.write("        R(1s)           R(2s)           R(2p)           R(3s)           R(3p)            R(3d)\n")
for n1=1, 3, 1 do
  for l1=0, n1-1, 1 do
    io.write(string.format("R(%i %i) ",n1,l1))
    for n2=1, 3, 1 do
      for l2=0, n2-1, 1 do
        sum=0
        for ir = 0, Nr, 1 do
          r = ir*dr
          sum = sum + dr * r^2 * Rnl(n1,l1,r,1) * r * Rnl(n2,l2,r,1) / a0
        end
        io.write(string.format("%15.8E ",sum))
      end
    end
    io.write("\n")
  end
end

-- Slater Integrals
io.write("\n<Rn1l1(r1) Rn2(r2)l2(r)| e^2/(|r1-r2|) |Rn1l1(r1) Rn2(r2)l2(r)> and <Rn1l1(r1) Rn2(r)l2(r2)| e^2/(|r1-r2|) |Rn1l1(r2) Rn2(r)l2(r1)> [units of Rydberg]\n");
dr=0.05
Nr=1000
for n1=1, 3, 1 do
  for l1=0, n1-1, 1 do
    for n2=1, 3, 1 do
      for l2=0, n2-1, 1 do
        io.write(string.format("R1(%i %i) R2(%i %i) ",n1,l1,n2,l2))
        for k=0, 2*math.min(l1,l2), 2 do
          io.write(string.format("F[%i]= ",k))
          sum=0
          for ir1 = 1, Nr, 1 do
            r1 = ir1*dr
            for ir2 = 1, Nr, 1 do
              r2 = ir2*dr
              sum = sum + a0 * dr^2* r1^2 * r2^2 * Rnl(n1,l1,r1,1)^2 * Rnl(n2,l2,r2,1)^2 * math.min(r1,r2)^k / math.max(r1,r2)^(k+1)
            end
          end
          io.write(string.format("%15.8E ",sum))
        end
        io.write("\n")
        io.write(string.format("                "))
        for k=math.abs(l1-l2), l1+l2, 2 do
          io.write(string.format("G[%i]= ",k))
          sum=0
          for ir1 = 1, Nr, 1 do
            r1 = ir1*dr
            for ir2 = 1, Nr, 1 do
              r2 = ir2*dr
                 sum = sum + a0 * dr^2* r1^2 * r2^2 * Rnl(n1,l1,r1,1) * Rnl(n1,l1,r2,1) * Rnl(n2,l2,r1,1) * Rnl(n2,l2,r2,1) * math.min(r1,r2)^k / math.max(r1,r2)^(k+1)
              end
            end
            io.write(string.format("%15.8E ",sum))
          end
          io.write("\n")
      end
    end
  end
end

-- Expectation value of bessel functions
-- q in units of inverse bohr radii
file = io.open("Bessel", "w")
file:write("<Rn1l1|j_k(qr)|Rn2l2> ");
for n1=1, 3, 1 do
  for l1=0, n1-1, 1 do
    for n2=1, 3, 1 do
      for l2=0, n2-1, 1 do
        for k=math.abs(l1-l2), l1+l2, 2 do
          file:write(string.format("%2i %2i %2i %2i %2i         ",n1,l1,n2,l2,k));
        end
      end
    end
  end
end
file:write("\n")
dr=0.05
Nr=1000
dq=0.1
Nq=100
for iq=0, Nq, 1 do
  q = iq*dq
  file:write(string.format("%15.8E ",q))
  for n1=1, 3, 1 do
    for l1=0, n1-1, 1 do
      for n2=1, 3, 1 do
        for l2=0, n2-1, 1 do
          for k=math.abs(l1-l2), l1+l2, 2 do
            sum=0
            for ir = 0, Nr, 1 do
              r = ir*dr
              sum = sum + dr * r^2 * Rnl(n1,l1,r,1) * math.SphericalBesselJ(k,q*r) * Rnl(n2,l2,r,1)
            end
            file:write(string.format("%15.8E ",sum))
          end
        end
      end
    end
  end
file:write("\n")
end

-- write the gnuplot script to a file
file = io.open("HWaveFunctions.gnuplot", "w")
file:write(gnuplotInput)
file:close()

-- call gnuplot to execute the script
os.execute("gnuplot HWaveFunctions.gnuplot ; ps2pdf HWaveFunctions.ps")