      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 i,N
      parameter(N=)
      double precision a0,R10,R20,R21,R30,R31,R32,R40,h,Rcut,
     & int10,int20,int21,int30,int31,int32,int40,auxJac,r

c Bohr radius in nanometers
      a0=0.053d0  
c Maximum extent of the integration grid
      Rcut=  /a0
c integration step in r
      h=Rcut/dble(N+1)
****************************************************************
c This block to check the normalization of all the orbitals
****************************************************************
c initialize all integrals to zero
      int10=0d0
      int20=0d0
      int21=0d0
      int30=0d0
      int31=0d0
      int32=0d0
      int40=0d0
      do i=1,N 
       r=dble(i)*h
c The radial wavefunctions are given in HydLikeFunctions in another file
       call HydLikeFunctions(r,R10,R20,R21,R30,R31,R32,R40)
c      print*, R10,R20,R21,R30,R31,R32,R40
       if(i.eq.1.or.i.eq.N) then
c with the trapezoidal integral rule the first and last point have
c to be treated differently
        auxJac=r**2*h*1.5d0
       else
        auxJac=r**2*h
       endif
       int10=int10+R10**2*auxJac
       int20=int20+R20**2*auxJac
       int21=int21+R21**2*auxJac 
       int30=int30+R30**2*auxJac
       int31=int31+R31**2*auxJac
       int32=int32+R32**2*auxJac
       int40=int40+R40**2*auxJac
      enddo
      print*, int10, int20, int21, int30, int31, int32, int40
*************************************************************


      end


