program densxz
c
c reads 3-dim mesh density data and prints the density in the x-z plane
c
c 17/apr/96
c 9,17/nov/94
c 8/nov/94 created by modifying readdens.for
c
      parameter(ndi=11,ndo=6)
      parameter(mx=13,my=13,mz=14)
      dimension irho(mx,my,mz,2),nx(2),ny(2),nz(2),unit(2)
      dimension rho(mx,mz,2)
      character buff*132,fmt*6,neupro(2)*7
      data neupro / 'neutron' , ' proton ' /
c
c-------------------- (1) reading the input data ---------------------
c
c Change the input-file name in the following open statement
c to the one you want to see.
c
      open(ndi,file='/dat/ftp/hfs3/output/n126082g.res',status='OLD'
     & ,action='READ')
c
      do 10 it=1,2
      do 10 k=1,mz
      do 10 j=1,my
      do 10 i=1,mx
   10 irho(i,j,k,it)=0.0d0
c
   20 read(ndi,'(A)') buff
      if(buff(1:22).ne.' NEUTRON DENSITY:MESH=') go to 20
      do 30 it=1,2
        read(buff(23:31),'(3I3)') nx(it),ny(it),nz(it)
        read(buff(38:52),'(F15.8)') unit(it)
        read(buff(58:63),'(A6)') fmt
        read(ndi,'('//fmt//')') 
     &  (((irho(i,j,k,it),i=1,nx(it)),j=1,ny(it)),k=1,nz(it))
        read(ndi,'(A)') buff
   30 continue   
      close(ndi)
c
c--------- (2) interpolation of the density to the xz-plane ----------
c
c Note that unit(it)*irho(i,j,k,it) is the nucleon density [fm**(-3)]
c at (|x|,|y|,|z|)=(i-0.5, j-0.5, k-0.5) [fm]. 
c 'it' means the type of nucleon: it=1 for neutron, it=2 for proton.
c
c The following interpolation procedure gives rho(i,k,it)
c which is the nucleon density at (|x|,|y|,|z|)=(i-0.5, 0, k-0.5) [fm].
c
      do 40 it=1,2
      do 40 k=1,mz
      do 40 i=1,mx
        dens1=irho(abs(i),1,abs(k),it)*unit(it)
        dens2=irho(abs(i),2,abs(k),it)*unit(it)
        rho(i,k,it)=(9.0d0*dens1-dens2)/8.0d0
   40 continue
c
c---------------- (3) checking the normalization ---------------------
c
c If the precision of the noramization is unsatisfactory, renormalize 
c the density by yourself. We do not do it here because the
c renormalization factor depends on the method of numerical integration
c you use.
c
c (1) integration with the mid-point formula 
c     in the Cartesian coordinates
c
      do 44 it=1,2
        sum=0.0d0
        do 42 k=1,mz
        do 42 j=1,my
        do 42 i=1,mx
          sum=sum+irho(i,j,k,it)
   42   continue
c The factor 8 of the following line is to account for the parts
c where at least one of {x,y,z} are negative.
        sum=sum*unit(it)*8
        write(6,900) neupro(it),sum,' (Cartesian)'
   44 continue
c
c (2) integration with the mid-point formula 
c     in the cylindrical coordinates
c
      pi=4.0d0*atan(1.0d0)
      do 80 it=1,2
        sum=0
        do 70 i=1,mx
          x=i-0.5d0
          do 60 k=1,mz
            sum=sum+rho(i,k,it)*(2*pi*x)
   60     continue
   70   continue
c The factor 2 of the following line is to account for the z<0 part.
        sum=sum*2
        write(6,900) neupro(it),sum,' (cylindrical)'
  900   format(' ',a,' number =',f10.5,a
     &  ,' :see the comments in the program')
   80 continue
c
c---------------------- (4) printing the result ----------------------
c
c*    open(ndo,file='dens_xz.dat',status='UNKNOWN',form='FORMATTED')
c
      do 100 it=1,2
        write(ndo,910) neupro(it)
  910   format(' ',A,' density in unit of 1E-4 fm**(-3) in y=0 plane')
        write(ndo,920) (i-0.5d0,i=1,mx) 
  920   format(' z/x ',13f5.1)
        do 90 k=1,mz
          z=k-0.5d0
          write(ndo,930) z,(int(rho(i,k,it)/1.0d-4),i=1,mx)
  930     format(f5.1,13i5)
   90   continue
  100 continue
c
c*    close(ndo)
c
        end