c=====================================================================c

      subroutine plotd(lpr)

c=====================================================================c
C
C     prepares plot of densities in coordinate space
C
c---------------------------------------------------------------------c
      implicit real*8 (a-h,o-z)
      include 'dirhb.par'
c
      
      logical lpr
c
      dimension pn(nox), ppn(nox),pnn(nox),pnp(nox),pntot(nox)
      dimension dtot(-100:100,-100:100)
      dimension roSLy4(mg)
c
      common /baspar/ hom,hb0,b0
      common /defbas/ beta0,q,bp,bz
      common /mathco/ zero,one,two,half,third,pi
      common /optopt/ itx,icm,icou,ipc,inl,idd
      common /dens  / ro(MG,4),dro(MG,4)
      common /densjp  / rojp(MG,3)
      common /rocan/ roshca(NHHX,NB2X),rocan(MG,8)    
      common /tapes / l6,lin,lou,lwin,lwou,lplo,laka,lvpp
      common /physco/ hbc,alphi,r0
      common /masses/ amu,amsig,amome,amrho   
      common /gaussh/ xh(NGH),wh(NGH),zb(NGH)
      common /gaussl/ xl(NGL),wl(NGL),sxl(NGL),rb(NGL)      
      common /gaucor/ ww(MG)    
      common /dloc  / dl(MG,3)
c     common /dloc/ rvltot(MG),clocpup(MG) 
c      common /roloc/ roloc(MG,2,2)
c      common /dcanjp/ dcanjp(mg),dloc(mg,4),ddcanjp(mg)
c      common /roljpt / roljpt(mg),droljpt(mg)
c      common /dlocloc/ dlocloc(mg)
c      common /dlljp /dlljp(mg),dlljpud(mg),ddlljp(mg),ddlljpud(mg)
c      common /dlljp2/ dllpu(mg),dllpd(mg),dllnu(mg),dllnd(mg)
c
c
      if (lpr)
     &write(l6,*) ' ****** BEGIN PLOTD ********************************'
c
c     number of points for the plot
      mxr  = 100
      mxz  = 100
      mxy  = 100
c     plot step in (fm)
      dz=0.10
      dr=0.10
      dy = 0.10
      am = amu*hbc
      do ir = -mxr,mxr
         do iz = -mxz,mxz
            dtot(iz,ir) = 0.0d0
         enddo
      enddo
      ih=1
      il=1
      i=ih+(il-1)*ngh
      z=zb(ih)
      r=rb(il)
c---- plot for densities:
      open(lplo,file='dis.plo',status='unknown')
      open(20,file='r.res',status='unknown')
      open(21,file='z.res',status='unknown')
      open(22,file='rho.res',status='unknown')
      open(23,file='rhoIFRIT.txt',status='unknown')
      open(24,file='rhoPARADDME2.csv',status='unknown')
      open(25,file='rhoPARASLy4.csv',status='unknown')
      open(26,file='rhoPARAD1S.csv',status='unknown')
c$$$      open(30,file='ropup.res',status='unknown')
c$$$      open(31,file='ropdown.res',status='unknown')
c$$$      open(32,file='ronup.res',status='unknown')
c$$$      open(33,file='rondown.res',status='unknown')
c$$$      open(34,file='rototup.res',status='unknown')
c$$$      open(35,file='rototdown.res',status='unknown')
c$$$      open(36,file='rotot.res',status='unknown')
c$$$      open(37,file='rolocpup.res',status='unknown')
c$$$      open(38,file='rolocpdown.res',status='unknown')
c$$$      open(39,file='rolocnup.res',status='unknown')
c$$$      open(40,file='rolocndown.res',status='unknown')
c$$$      open(41,file='rholocn.res',status='unknown')
c$$$      open(42,file='rholocp.res',status='unknown')
c$$$      open(43,file='rholocalpha.res',status='unknown')
c$$$      open(44,file='rhon.res',status='unknown')
c$$$      open(45,file='rhop.res',status='unknown')
c$$$      open(46,file='densrad_n.res',status='unknown')
c$$$      open(47,file='densrad_p.res',status='unknown')
      z=zero
      do ir=-mxr,mxr
         write(20,*) ir*dr
      enddo
      do iz=-mxz,mxz
         z=iz*dz
         write(21,*) z
      enddo
c         do ir=-mxr,mxr
c            r=ir*dr
c            rho = rdens(1,ro(1,2),pn,abs(z)/bz/b0,abs(r)/bp/b0)
c            rho1 = rdens(1,rocan(1,1),pn,abs(z)/bz/b0,abs(r)/bp/b0)
c            rho2 = rdens(1,rocan(1,2),pn,abs(z)/bz/b0,abs(r)/bp/b0)
c            rhotot1 = rdens(1,ro(1,1),pn,abs(z)/bz/b0,abs(r)/bp/b0)
c            rhotot2 = rdens(1,ro(1,2),pn,abs(z)/bz/b0,abs(r)/bp/b0)
c            write(lplo,100) z,r,rho1+rho2
c            
c  100       format(2f10.3,f15.6)             
c         enddo
c      enddo
      ijp = 0
      ijpp = 1
      if(ijpp.eq.1) then
         do ir = -mxr,mxr
            r=ir*dr
            do iz = -mxz,mxz
               z=iz*dz
               ijp = ijp + 1
!            rhotot1 = rdens(1,ro(1,1),pn,abs(z)/bz/b0,abs(r)/bp/b0)
               rhotot2 = rdens2(1,ro(1,2),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
               rhotest = rdens2(1,rojp(1,3),pntot,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
               rhon    = rdens2(1,rojp(1,1),pnn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
               rhop    = rdens2(1,rojp(1,2),pnp,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
               ijp = ijp+1
c            rhocan = rdens2(1,rocan(1,1),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
!           rhocan = rdens2(1,dl(1,1),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
!           rhocan2 = rdens2(1,dl(1,2),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
!           rhocan3 = rdens2(1,dl(1,3),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
c            rhon = rdens2(1,rojp(1,1),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
c            rhop = rdens2(1,rojp(1,2),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
c            rhotot2 = rdens(1,rvltot,pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
c            rholpup  = rdens2(1,dllpu,pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
c            rholpdown= rdens2(1,dllpd,pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
c            rholnup  = rdens2(1,dllnu,pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
c            rholndown= rdens2(1,dllnd,pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)
c            rhoup  = rdens2(1,rocan(1,4),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)            
c            rhodown  = rdens2(1,rocan(1,3),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp)            
!            rhotot2jp = rdensjp(1,ro(1,2),pn,abs(z)/bz/b0,abs(r)/bp/b0)
c            dtot(iz,ir) = rhotot2
!              write(22,*) ir,iz,rhotot2
               write(22,'(2i5,4f18.9)') ir,iz,rhon,rhop,rhotest,abs(rhon-rhop)
c            write(37,*) ir,iz,rholpup
c            write(38,*) ir,iz,rholpdown
c            write(39,*) ir,iz,rholnup
c            write(40,*) ir,iz,rholndown
c            write(34,*) ir,iz,rhoup
c            write(35,*) ir,iz,rhodown
c            write(41,*) ir,iz,rhocan
c            write(42,*) ir,iz,rhocan2
c            write(43,*) ir,iz,rhocan2
c            write(44,*) ir,iz,rhon               
c            write(45,*) ir,iz,rhop   
c            if(ir.ge.0.and.iz.eq.0) then
c               write(46,*) r,rhon 
c               write(47,*) r,rhop 
c            endif
            enddo
         enddo
      endif                     !ijpp
      ijpok=0
      if(ijpok.eq.1) then
c         do i = 1,mg
c            read(95,*) roSLy4(i)
!      read(94,*) roD1S(i)
!      roSLy4(i) = roSLy4(i)*ww(i)
c         enddo
!        bzSLy4 = 1.50206407294644 !1.15673359384027!1.82562524170824 !1.77098448898468
!        bpSLy4 = 1.65112633190644 !1.15673359384027!1.44105985404351 !1.3979290989760
!      bzD1S = 1.77098448898468
!      bpD1S = 1.3979290989760
         write(23,*) 2*mxz+1, 2*mxz+1, 2*mxz+1
!     y0 = 0.5
!      do iy = -mxy,mxy
!         y = iy*dy
!         if(y.ge.-y0.and.y.le.y0) then
!            mxz = 25
!            dz = 0.20
!         else
!            mxz = 10
!            dz = 0.50
!         endif
         ijp2 = 0
         do ix = -mxz,mxz
            x = ix*dz
            do iy = -mxy,mxy
               y = iy*dy
               r = sqrt(x**2+y**2)
               do iz=-mxz,mxz
                  z=iz*dz
                  ijp2 = ijp2 + 1
                  rhototjpp1 = rdens(1,ro(1,2),pn,abs(z)/bz/b0,abs(r)/bp/b0,ijp2)
c                  rad = 2.0d0*(2.0d0+x**2+y**2+z**2+0.5*(x**2-y**2)/(x**2+y**2+z**2))
!           rhototjpp2 = rdens(1,roSLy4,pn,abs(z)/bzSLy4,abs(r)/bpSLy4)
!             rhototjpp3 = rdens(1,roD1S,pn,abs(z)/bzD1S,abs(r)/bpD1S)
                  write(23,*) rhototjpp1
c                  write(23,*) rad
!           if(rhototjpp1.ge.0.005) write(24,110) x,y,z,rhototjpp1
!             if(rhototjpp2.ge.0.005)  write(25,110) x,y,z,rhototjpp2
!             write(26,110) x,y,z,rhototjpp3
               enddo
            enddo
         enddo
      endif
 110  format(f6.3,';',f6.3,';',f6.3,';',f20.18)
      close(lplo)
      close(20)
      close(21)
      close(22)
      close(23)
      close(24)
      close(25)
      close(26)
c----- Integrals
      dint = 0.0d0
      anrjp = (40.0d0+40.0d0) +1.0d0 
      arjp = (40.0d0+40.0d0)*dr/anrjp
      anzjp = (74.0d0-11.0d0) +1.0d0 
      azjp = (74.0d0-11.0d0)*dz/anzjp
      do ir = -40, 40
         r=abs(ir*dr)
         do iz = 11,74
            z=iz*dz
            dint = dint + dtot(iz,ir)*r
         enddo
      enddo
      dint = arjp*azjp*pi*dint
      write(*,*)'=========> DINT =',dint
      write(*,*)'b0,bp,bz = ',b0,bp,bz
c
c$$$      close(30)
c$$$      close(31)
c$$$      close(32)
c$$$      close(33)
c$$$      close(34)
c$$$      close(35)
c$$$      close(36)
c$$$      close(37)
c$$$      close(38)
c$$$      close(39)
c$$$      close(40)
c$$$      close(41)
c$$$      close(42)
c$$$      close(43)
c$$$      close(44)
c$$$      close(45)
c$$$      close(46)
c$$$      close(47)
      if (lpr)
     &write(l6,*) ' ****** END PLOTD **********************************'
      return
C-end-PLOT
      end







c======================================================================c

      double precision function rdens2(is,ro,pn,z,r,ijp)

c======================================================================c
c
c     calculation of the density ro at arbitrary meshpoint x
c     x is given in units of the oscillator lenght: x = r/b0f 
c
c     is = 1: density given at Gauss-Meshpoints ro(ih)
c          2: density given through oscillator expansion pn(n)
c
c----------------------------------------------------------------------c
      implicit real*8 (a-h,o-z)
      include 'dirhb.par'
c

c
      dimension ro(MG),pn(nox),qhb(0:nbx),qlb(0:nbx),phi(nox)
C
      common /baspar/ hom,hb0,b0
      common /basnnn/ n0f,n0b
      common /propag/ gg(NOX,NOX,3),psi(NOX,MG)
      common /bosqua/ nzb(NOX),nrb(NOX),NO
      common /gaucor/ ww(MG) 
      common /gfvsq / sq(0:igfv)
      common /gfvsqi/ sqi(0:igfv)
      common /gfvsqh/ sqh(0:igfv)
      common /gfvshi/ shi(0:igfv)
      common /gfvwgi/ wgi(0:igfv)
      common /gfvwfi/ wfi(0:IGFV)
      common /mathco/ zero,one,two,half,third,pi
      common /tapes / l6,lin,lou,lwin,lwou,lplo,laka,lvpp

c
c---- calculation of the oscillator expansion for ro
      if (is.eq.1) then
         do n = 1,no
               s = zero
               do ih = 1,MG
                  s = s + ro(ih)*psi(n,ih)*ww(ih)
               enddo
               pn(n) = s
         enddo
      endif
c
      zz=z**2
      w4pii = pi**(-0.25d0)
      w0  = w4pii*exp(-half*zz)
      wb  = w0
      qhb(0) = wb
      qhb(1) = sq(2)*wb*z
      do n = 2,N0B
         qhb(n) = sqi(n)*(sq(2)*z*qhb(n-1)-sq(n-1)*qhb(n-2))
      enddo   ! n      
 
      x=r**2
      w0 = sq(2)*exp(-half*x)
      wb = sq(2)*exp(-half*x) / sqrt(2*pi) 
      qlb(0) = wb*wfi(0)
      qlb(1) = wb*(1-x)*wfi(1)
      do n = 2,N0B
         qlb(n) = ((2*n-1-x)*qlb(n-1)-(n-1)*qlb(n-2))/n 
      enddo       
      do n = 1,NO
         nz = nzb(n)
         nr = nrb(n)
	     phi(n) = qhb(nz)*qlb(nr) 
      enddo   ! no      
c
      s = zero
      do n = 1,no
         s = s + pn(n)*phi(n)
      enddo
      rdens2 = s/b0**3
      if (rdens2 .lt. 0.d0) rdens2 =0.d0
c
      return
c-end-RDENS
      end










c======================================================================c

      double precision function rdens(is,ro,pn,z,r,ijp)

c======================================================================c
c
c     calculation of the density ro at arbitrary meshpoint x
c     x is given in units of the oscillator lenght: x = r/b0f 
c
c     is = 1: density given at Gauss-Meshpoints ro(ih)
c          2: density given through oscillator expansion pn(n)
c
c----------------------------------------------------------------------c
      implicit real*8 (a-h,o-z)
      include 'dirhb.par'
c

c
      dimension ro(MG),pn(nox),qhb(0:nbx),qlb(0:nbx),phi(nox)
C
      common /baspar/ hom,hb0,b0
      common /basnnn/ n0f,n0b
      common /propag/ gg(NOX,NOX,3),psi(NOX,MG)
      common /bosqua/ nzb(NOX),nrb(NOX),NO
      common /gaucor/ ww(MG) 
      common /gfvsq / sq(0:igfv)
      common /gfvsqi/ sqi(0:igfv)
      common /gfvsqh/ sqh(0:igfv)
      common /gfvshi/ shi(0:igfv)
      common /gfvwgi/ wgi(0:igfv)
      common /gfvwfi/ wfi(0:IGFV)
      common /mathco/ zero,one,two,half,third,pi
      common /tapes / l6,lin,lou,lwin,lwou,lplo,laka,lvpp
      common /denstool/ pnjp(nox)

c
c---- calculation of the oscillator expansion for ro
      if (is.eq.1) then
         do n = 1,no
            if(ijp.eq.1) then
               s = zero
               do ih = 1,MG
                  s = s + ro(ih)*psi(n,ih)*ww(ih)
               enddo
               pn(n) = s
               pnjp(n) = s
            else
               pn(n) = pnjp(n)
            endif
         enddo
      endif
c
      zz=z**2
      w4pii = pi**(-0.25d0)
      w0  = w4pii*exp(-half*zz)
      wb  = w0
      qhb(0) = wb
      qhb(1) = sq(2)*wb*z
      do n = 2,N0B
         qhb(n) = sqi(n)*(sq(2)*z*qhb(n-1)-sq(n-1)*qhb(n-2))
      enddo   ! n      
 
      x=r**2
      w0 = sq(2)*exp(-half*x)
      wb = sq(2)*exp(-half*x) / sqrt(2*pi) 
      qlb(0) = wb*wfi(0)
      qlb(1) = wb*(1-x)*wfi(1)
      do n = 2,N0B
         qlb(n) = ((2*n-1-x)*qlb(n-1)-(n-1)*qlb(n-2))/n 
      enddo       
      do n = 1,NO
         nz = nzb(n)
         nr = nrb(n)
	     phi(n) = qhb(nz)*qlb(nr) 
      enddo   ! no      
c
      s = zero
      do n = 1,no
         s = s + pn(n)*phi(n)
      enddo
      rdens = s/b0**3
      if (rdens .lt. 0.d0) rdens =0.d0
c
      return
c-end-RDENS
      end


c======================================================================c

      double precision function rdensjp(is,ro,pn,z,r)

c======================================================================c
c
c     calculation of the density ro at arbitrary meshpoint x
c     x is given in units of the oscillator lenght: x = r/b0f 
c
c     is = 1: density given at Gauss-Meshpoints ro(ih)
c          2: density given through oscillator expansion pn(n)
c
c----------------------------------------------------------------------c
      implicit real*8 (a-h,o-z)
      include 'dirhb.par'
c

c
      dimension ro(MG),pn(nox),qhb(0:nbx),qlb(0:nbx),phi(nox)
C
      common /baspar/ hom,hb0,b0
      common /basnnn/ n0f,n0b
      common /propag/ gg(NOX,NOX,3),psi(NOX,MG)
      common /bosqua/ nzb(NOX),nrb(NOX),NO
      common /gaucor/ ww(MG) 
      common /gfvsq / sq(0:igfv)
      common /gfvsqi/ sqi(0:igfv)
      common /gfvsqh/ sqh(0:igfv)
      common /gfvshi/ shi(0:igfv)
      common /gfvwgi/ wgi(0:igfv)
      common /gfvwfi/ wfi(0:IGFV)
      common /mathco/ zero,one,two,half,third,pi
      common /tapes / l6,lin,lou,lwin,lwou,lplo,laka,lvpp

c
c---- calculation of the oscillator expansion for ro
      if (is.eq.1) then
         do n = 1,no
            s = zero
            do ih = 1,MG
               s = s + ro(ih)*psi(n,ih)*ww(ih)
            enddo
            pn(n) = s
         enddo
      endif
c
      zz=z**2
      w4pii = pi**(-0.25d0)
      w0  = w4pii*exp(-half*zz)
      wb  = w0
      qhb(0) = wb
      qhb(1) = sq(2)*wb*z
      do n = 2,N0B
         qhb(n) = sqi(n)*(sq(2)*z*qhb(n-1)-sq(n-1)*qhb(n-2))
      enddo   ! n      
 
      x=r**2
      w0 = sq(2)*exp(-half*x)
      wb = sq(2)*exp(-half*x) / sqrt(2*pi) 
      qlb(0) = wb*wfi(0)
      qlb(1) = wb*(1-x)*wfi(1)
      do n = 2,N0B
         qlb(n) = ((2*n-1-x)*qlb(n-1)-(n-1)*qlb(n-2))/n 
      enddo       
      do n = 1,NO
         nz = nzb(n)
         nr = nrb(n)
	     phi(n) = qhb(nz)*qlb(nr) 
      enddo   ! no      
c
      s = zero
      do n = 1,no
         s = s + pn(n)*phi(n)
      enddo
      rdensjp = s/b0**3
      if (rdens .lt. 0.d0) rdens =0.d0
c
      return
c-end-RDENS
      end
