      program PlotDistributions
c The purpose is to save the charge distributions to disk, once normalized
      implicit none
      integer i,N
      parameter(N=1000)
      double precision a0,R10,R20,R21,R30,R31,R32,R40,h,Rcut,r
c Bohr radius in nanometers
      a0=0.053d0  
c Maximum extent of the integration grid
      Rcut=2d0/a0
c integration step in r
      h=Rcut/dble(N+1)
****************************************************************
c files where charge distributions will be saved
      open(unit=7,file='DATA.DIR/1s.dat')
      open(unit=8,file='DATA.DIR/2s.dat')
      open(unit=9,file='DATA.DIR/2p.dat')
      open(unit=10,file='DATA.DIR/3s.dat')
      open(unit=11,file='DATA.DIR/3p.dat')
      open(unit=12,file='DATA.DIR/3d.dat')
      open(unit=13,file='DATA.DIR/4s.dat')
      open(unit=15,file='DATA.DIR/ChargeDist.dat')
      do i=1,N 
       r=dble(i)*h
       call HydLikeFunctions(r,R10,R20,R21,R30,R31,R32,R40)
       write(7,*)  r,r**2*R10**2
       write(8,*)  r,r**2*R20**2
       write(9,*)  r,r**2*R21**2
       write(10,*) r,r**2*R30**2
       write(11,*) r,r**2*R31**2
       write(12,*) r,r**2*R32**2 
       write(13,*) r,r**2*R40**2
c The inner charge distribution weighted by the number of electrons
       write(15,*) r,
     & r**2*2d0*(R10**2+R20**2+3d0*R21**2+R30**2+3d0*R31**2)
      enddo
*************************************************************
      close(7)
      close(8)
      close(9)
      close(10)   
      close(11)
      close(12) 
      close(13)
      close(15)
      end

