      program checkNormalization
c The purpose is to decide which orbital has the smallest energy, 3d or 4s,
c in a simple model of the atom where filled orbitals are taken from hydrogen
      implicit none
      integer i1,i2,j,Nr,Nx,lpass
      parameter(Nr=50)
      parameter(Nx=16)
      double precision a0,R10,R20,R21,R30,R31,R32,R40,Rcut,Z,
     & auxJacr1,auxJacr2,r1,r2,x,intr1,intr2,intx,chargeDens,pi,
     & Q10,Q20,Q21,Q30,Q31,Q32,Q40,SquaredSphHarmonic,auxJacx,hx
      double precision Vscreened,hr1,hr2
      external SquaredSphHarmonic
      pi=3.14159265d0
c Bohr radius in nanometers
      a0=0.053d0  
c Maximum extent of the integration grid
      Rcut=  /a0
c integration steps in r and x
      hr1=Rcut/dble(Nr+1)
      hx=2d0/dble(Nx+1)
c Z is 19 but it can be set to 0 to test the code without the nucleus
      Z=19d0
c a file to save the screened potential
      open(unit=15,file='DATA.DIR/Vscreened.dat')
c here is the choice between 3d or 4s (pass 2 or 0 argument)
      lpass=0
****************************************************************
c This block to compute the triple integral
****************************************************************
c initialize all integrals to zero
      intr1=0d0
c     Start integral in r1
      do i1=1,Nr 
       r1=dble(i1)*hr1
       call HydLikeFunctions(r1,R10,R20,R21,R30,R31,R32,R40)
       if(i1.eq.1.or.i1.eq.Nr) then
        auxJacr1=r1**2*hr1*1.5d0
       else
        auxJacr1=r1**2*hr1
       endif
c      Start integral in r2, here choose interval 0 to r1 or infty.  
       hr2=r1/dble(Nr+1)
c       hr2=hr1
       intr2=0d0
       do i2=1,Nr
        r2=dble(i2)*hr2
        call HydLikeFunctions(r2,Q10,Q20,Q21,Q30,Q31,Q32,Q40)
        chargeDens=2d0*(Q10**2+Q20**2+3d0*Q21**2+Q30**2+3d0*Q31**2)
        if(i2.eq.1.or.i2.eq.Nr) then
         auxJacr2=r2**2*hr2*1.5d0
        else
         auxJacr2=r2**2*hr2
        endif
c       Start integral in cos(theta)
        intx=0d0
c We separate the first and last point of the loop due to the 1.5
c to avoid having to evaluate an "if" condition in the inner loop
        auxJacx=hx*chargeDens
        do j=2,Nx-1
         x=-1d0+dble(j)*hx
         intx=intx+auxJacx/dsqrt(r1**2+r2**2-2d0*r1*r2*x)
     &             *SquaredSphHarmonic(lpass,x)
        enddo
c here are the border points to complete the trapezoidal rule
        auxJacx=auxJacx*1.5d0
        x=-1d0+hx
        intx=intx+auxJacx/dsqrt(r1**2+r2**2-2d0*r1*r2*x)
     &            *SquaredSphHarmonic(lpass,x)
        x=-1d0+dble(Nx)*hx
        intx=intx+auxJacx/dsqrt(r1**2+r2**2-2d0*r1*r2*x)
     &            *SquaredSphHarmonic(lpass,x)
c Can use this line to set the integral to zero and leave only the nuclear 
c charge for studying the pieces separately
c        intx=0d0
c       Ended integral in cos(theta), continue with r' integral
        intr2=intr2+2d0*pi*intx*auxJacr2
       enddo
       print*, i1

c      Ended integral in r', continue with r integral
c      Choose between the two lines depending on 3d or 4s
       Vscreened=(intr2-Z/r1)*auxJacr1
       if(lpass.eq.0) then
        intr1=intr1+Vscreened*R40**2
       elseif(lpass.eq.2) then
        intr1=intr1+Vscreened*R32**2
       else
        print*, 'only set for 3d and 4s'
        stop
       endif
c this line saves to disk the effective potential
       write(15,*) r1, Vscreened/auxJacr1
c can also use the next line to save the effective charge
c       write(15,*) r1, Z-intr2*r1
      enddo
c     Ended integral in r1. Print energy in Hartree and in eV to screen
      print*, intr1, intr1*27.2d0
      close(15)
*************************************************************
      end



      double precision function SquaredSphHarmonic(l,costheta)
      implicit none
      integer l
      double precision costheta,pi
      pi=3.14159265d0
      if(l.eq.0) then
       SquaredSphHarmonic=1d0/(4d0*pi)
       return
      endif
      if(l.eq.2) then
       SquaredSphHarmonic=5d0/(16d0*pi)*(3d0*costheta**2-1d0)**2
       return
      endif
      print*, 'only set for 0,2'
      stop
      end

